MathematicS In Action Optimal strokes for the 4-sphere swimmer at low Reynolds number in the regime of small deformations

The paper deals with the optimal control problem that arises when one studies the 4 sphere artiﬁcial swimmer at low Reynolds number. Composed of four spheres at the end of extensible arms, the swimmer is known to be able to swim in all directions and orientations in the 3D space. In this paper, optimal strokes, in terms of the energy expended by the swimmer to reach a prescribed net displacement, are fully described in the regime of small strokes. In particular, we introduce a bivector formalism to model the displacements that turns out to be elegant and practical. Numerical simulations are also provided that conﬁrm the theoretical predictions. Two trajectories are the same displacement and rotation δp such that Λ − , one in solid and the other in dotted line. The three picture correspond respectively to the front, side and top views of those two trajectories.


Introduction
Since the seminal paper by E. M. Purcell [19], the subject of understanding the swimming at low Reynolds number has known a growing interest in the Physics community. A comprehensive bibliography was given and reviewed in [14], and a more recent overview in [11]. The interested reader will find there a wide view of the different directions that have been considered so far together with future topics that need to be more thoroughly studied, to better understand the peculiarities of this topic.
Among others is the subject of proposing and studying mechanisms that overcome the celebrated "Scallop Theorem" stated by Purcell and generalizing the approach of his "three-link swimmer". Several devices have been proposed (see for instance [7,10,18]).
In the meantime, the problem of "being able to swim", i.e., does there exist a deformation strategy that produces, when the interaction with the fluid is taken into account, a motion of the swimmer, has been rephrased in terms of Control Theory, and tools have been introduced to rigorously prove the controllability of several systems (see e.g. [3,8,16]).
Finding the best swimming strategies, or, in other words (using Lighthill's efficiency definition [15]) the strokes that produce a given displacement for the least amount of expended energy, is a problem of optimal control that has also been studied. A clear link is done with geometry since optimal strokes can be seen as geodesics in a suitable sub-Riemannian space [2,3].
with arms whose length can be controlled and it has been shown that it is possible to control both the position and the orientation of the swimmer in 2D or 3D. For instance, a "plane swimmer", called SPr3, composed of three rigid spheres at the extremities of three arms making 120°o ne to another can be controlled to move freely, both in translation and rotation, in the 2D plane. Finding optimal periodic strokes for this plane swimmer, i.e. the ones that induce a given displacement for the smallest amount of expanded energy, was realized in [4,5] in the regime of small deformations around a symmetric configuration. It is shown in particular that they are planar ellipses in the space of shapes.
In this paper, we pursue the approach by considering the 4-sphere swimmer SPr4, still proposed in [2], moving in the 3D space. This artificial swimmer possesses four rigid spheres at the extremities of four arms that are free to elongate and retract. We show that, as a consequence of the symmetric initial shape where all the arms have an equal length, the system can be completely understood, still in the regime of small deformations. In particular, we explicitly find the optimal strokes, i.e. the ones that consume the least energy and that achieve a given displacement (both translations and rotations). In general, the strokes are no longer planar ellipses contrarily to what happens with control systems in lower dimensions. This is nevertheless the case for very special constrained displacements that are also identified. Eventually, we give an example of a displacement that can be achieved with non-unique optimal strokes, and provide the reader with explicit formulas for the remaining unknown constants appearing in the system, in the limit of large arms.
The paper is organized as follows. In Section 2 we describe the swimmer and write the control problem making use of the symmetries of the artificial swimmer. We show in particular how the formalism of bivectors can be used to express the optimal control problem in a concise and elegant way. In Section 3, we give the main result of the paper, namely a complete understanding of optimal periodic strokes for a given displacement. Strikingly enough, there is a subtle distinction to make depending on whether the prescribed displacement is a simple bivector or not. Eventually, Section 4 provides the reader with a series of numerical simulations and plots of optimal strokes in different situations. These numerical results are obtained with parameters that are given via an expansion of the system when the arms of the artificial swimmer are large when compared to the diameter of the balls.

The 4-sphere swimmer
This section is devoted to the description of the 4-sphere swimmer SPr4, together with the writing of the control problem to leading order in the regime of small deformations. Symmetries of the swimmer are used to unveil the canonical structure of the control problem that, as we shall see, only depends on a small number of parameters. In all this paper, we use the notation N 3 = {1, 2, 3} and N 4 = {1, 2, 3, 4}.

The model
Let us consider the swimmer SPr4 proposed in [2]. To that end, let (S 1 , S 2 , S 3 , S 4 ) be a regular reference tetrahedron centered at c ∈ R 3 such that dist(c, S i ) = 1 for all i ∈ N 4 . The swimmer consists of four balls (B i ) i∈N 4 in R 3 all of radius a > 0. Each ball B i for i ∈ N 4 is centered at b i ∈ R 3 , which can freely move along the ray starting at c and passing through S i , see This reflects the situation where the balls are linked to the center c by thin jacks that are able to elongate and retract. However, the viscous resistance of these jacks is neglected and therefore the fluid is assumed to permeate the entire open set R 3 \ 4 i=1 B i . The balls do not rotate around their arms which implies that the shape of the swimmer is completely determined by the four lengths ζ 1 , ζ 2 , ζ 3 , and ζ 4 of its arms measured from c to the center b i of each ball. Meanwhile, there are no restrictions on the global rotation of the swimmer around the center c, i.e. for fixed arm lengths, the swimmer is considered to be a rigid body in a Stokesian fluid.
The geometrical configuration of the swimmer can be described by two sets of variables: (1) The vector of shape variables ζ := (ζ 1 , ζ 2 , ζ 3 , ζ 4 ) ∈ S := ( 3/2a, +∞) 4 ⊆ R 4 + , each value ζ i being the length of the i−th arm, that we denote by i in what follows. The lower bound in the open interval is chosen such that the balls cannot overlap.
(2) The vector of position variables p = (c, R) ∈ P := R 3 × SO(3), which encodes the global position and orientation of the swimmer in space.
To be more precise, we consider the reference tetrahedron convexly spanned by the four unit vectors The position and orientation of the swimmer in R 3 are then respectively described by the coordinates of the center c ∈ R 3 and the rotation R ∈ SO(3) of the swimmer with respect to the reference orientation induced by the reference tetrahedron. Hence, if the arms are aligned with the (z i ) i∈N 4 , then this rotation matrix equals the identity matrix Id 3 ∈ SO(3). We therefore set b i := c + ζ i Rz i for the center of the ball B i . In [2] it is shown that the system SPr4 is fully controllable, i.e. both in shape ζ and position p, using only the rate of shape changeζ. Using the assumptions of self-propulsion and negligible inertia of the swimmer, the total viscous force and torque exerted by the surrounding fluid on the swimmer must vanish. This permits to write the system as (for details see [2]) In the preceding, denoting by L(V, W ) the set of linear maps between the vector spaces V and W , we respectively have 3) denotes the tangent space 1 to SO(3) at R. (We have denoted by Skew 3 (R) the space of real skew-symmetric 3 × 3 matrices.) Notice that the first component F c of F may be conveniently represented as a 3 × 4 matrix, while an analogous writing for F R leads to a 3 × 3 × 4 tensor.
In the spirit of [2] and [4], regarding the energy consumption of a swimming stroke we follow the notion of swimming efficiency suggested by Lighthill [15] and we adopt the following notion of optimality: energy minimizing strokes are the ones that minimize the kinematic energy dissipated while trying to reach a given net displacement δp. Mathematically speaking, strokes are periodic changes of shapes ζ ∈ H 1 (I, S), where the index stands for periodic functions, and we can take, without further restriction I = [0, 2π]. The total energy dissipation due to the stroke ζ given above can be evaluated through an adequate quadratic energy functional, c.f. [2], where the energy density G ∈ C 1 (S, M 4×4 (R)) is a function with values in the space of symmetric and positive definite matrices. Finally, the optimal stroke problem can be written as: a displacement δp = (δc, δR) T being given, find inf

Symmetries
In analogy to [4], we investigate the invariance of the system due to symmetries. For any initial condition p 0 := (c 0 , R 0 ) ∈ P and any shape curve ζ ∈ C 1 (J, S) with J an open interval containing zero, we denote by γ(c 0 , R 0 , ζ) : J → P the solution associated with the dynamical systeṁ as well as by γ c (c 0 , R 0 , ζ) and γ R (c 0 , R 0 , ζ) its projections on R 3 and SO(3), respectively. Therefore, for any t ∈ J, we havė (2.7)

Translational and rotational invariance
The Stokes equations are invariant to translations. Moreover, changing the initial orientation, i.e. applying a rotation to p 0 , in the dynamical system (2.6) leads to the same rotation applied to the resulting trajectory p(t). Therefore, we may state the following symmetry property of the control system (2.6) with respect to translations and rotations:

8)
and for the angular part of the solution

9)
at any point in time t ∈ J. 1 For an introduction to manifolds, we refer to [12].

Corollary 2.2.
The translational and rotational invariances given above immediately translate into symmetry properties of the dynamical system (2.1). Namely, we have the formulas In what follows, we therefore simplify the notation by defining (2.11)

Permutation of two arms
Previous symmetries were linked with the fact that the Stokes equations are invariant with respect to translations and rotations. They have nothing to do with the special parametrization of the swimmer itself. Here, we instead pay attention to the geometry of the swimmer and its symmetries to get further relations. To that aim, let us consider the effect of swapping two arms on the generic solution of the dynamical system (2.6). Take i, j ∈ {1, 2, 3, 4} with i = j, and let P ij ∈ L(R 4 , R 4 ) denote the map that interchanges the i-th and j-th coordinates. This physically corresponds to the swap of the arms i and j, if applied to the 4-dimensional shape space S. In addition, let S ij denote the reflection of R 3 that maps the arm i onto the arm j, in the reference orientation Id 3 . Geometrical inspection of the reference tetrahedron shows that S ij is always a reflection at a plane containing the remaining arms k and l. The condition that follows is easily understood when one realizes that for a given shape history ζ(t), γ(0, Id 3 , ζ)(t) and γ(0, Id 3 , P ij ζ)(t) can be deduced one from another by applying the symmetry S ij in space, as in Fig. 2.2. In other words, an observer watching the dynamics of γ(0, Id 3 , ζ)(t) of SPr4 in a mirror in the reflection plane of S ij sees the dynamics γ(0, Id 3 , P ij ζ)(t) of a micro-swimmer obtained from SPr4 by swapping arms i and j. Therefore, the translational part γ c (0, Id 3 , P ij ζ)(t) is obtained by applying the symmetry S ij to γ c (0, Id 3 , ζ)(t), while the rotational part γ R (0, Id 3 , P ij ζ)(t) is conjugated to γ R (0, Id 3 , P ij ζ)(t) under S ij .

Condition 2.3 (Swap of i and j)
). With the preceding notation, we have
This enables us to deduce symmetry relationships that F defined in (2.11) has to satisfy. Proposition 2.4. Let (i, j) ∈ {1, 2, 3, 4} 2 , with i = j. Let χ ∈ S be a given shape and η ∈ R 4 be any shape derivative. (2.14) and Proof. Let us consider J to be an open interval containing 0, and ζ ∈ C 1 (J, S) be a shape function satisfying ζ(0) = χ andζ(0) = η. We consider the trajectories γ(0, Id 3 , ζ)(t) and γ(0, Id 3 , P ij ζ)(t) and use Condition 2.3. Looking first at the translational component and differentiating with respect to time t gives on the one hanḋ while on the other hand, using Condition 2.3, we havė Equating both results, using (2.13), and taking t = 0 leads to Since η can be arbitrarily chosen, this gives (2.14). The proof of (2.15) follows the same lines. Namely, looking at the R-component of γ and differentiating again in time leads tȯ and, using Condition 2.3, tȯ Equation (2.15) follows again from equating both right-hand sides, using (2.13) and taking t = 0.
Eventually, we investigate the symmetry properties of the energy density through an arm swap P ij in the shape space. Due to the invariance of the energy dissipation under reflection, we deduce, with the same notation as before, the following proposition, whose proof, which follows the same arguments as in Proposition 2.4 is left to the reader. (2.16) or, equivalently, P ij G(P ij χ)P ij = G(χ) .

Small stroke expansion
In this section, we approximate the energy (2.4) and the dynamics (2.6) in the range of small strokes around a symmetric configuration, where all four arms have an identical length. In order to proceed, following [2], we first consider the first-order expansion of the system in ζ and take advantage of the symmetries described in the previous section to reveal the structure of the different terms. Then, we write the corresponding linearized control equations. We work in the neighborhood of the symmetric shape ξ 0 = (l 0 , l 0 , l 0 , l 0 ) T by considering ζ := ξ 0 + ξ. We consider deformations ξ around the initial and symmetric configuration ξ 0 , that belong to the spaceḢ 1 (I, R 4 ), i.e. the Sobolev space of 2π-periodic vector-valued functions having first order weak derivative in L 2 (I, R 4 ) and a vanishing average. In other words, we consider the deformations to be small and zero-averaged deviations from ξ 0 . Here, I denotes the closed interval [0, 2π] as before.
We first consider the energy (2.4) and set G ξ 0 (ξ) := G(ξ 0 + ξ). In this small stroke regime, we can approximate the energy density by is symmetric and positive definite, and ξ Ḣ1 = ξ L 2 . More precisely, Looking at the extra diagonal terms reveals that G 0 has the form for two parameters h and ρ. The matrix G 0 can then be diagonalized as where (τ i ) i∈N 4 is the following orthonormal basis of R 4 , made of eigenvectors of G and the (g i ) i∈N 4 are the corresponding eigenvalues Note that, G 0 being positive definite, we also have ρ > max(h, −3h).
Let us now focus on the equations of the dynamics (2.6). First, let us recall that from rotational invariance, we obtained in (2.10), (2.11) that F c and F R can be factorized as Similarly as before, we set F c,ξ 0 (ξ) := F c (ξ 0 + ξ) and F R,ξ 0 (ξ) := F R (ξ 0 + ξ). It has been shown in [2] that F is an analytic function, and we may therefore consider the first-order expansions in ξ: where we have denoted by ( e p ) p∈N 3 the canonical basis vectors of R 3 and by L := (L 1 , L 2 , L 3 ) the following basis of Skew 3 (R): (0)). In order to reveal the structure of (A p ) 0≤p≤3 and (B p ) 0≤p≤3 , we transfer the symmetry properties (2.14), (2.15) The following sections aim at finding the structures of both (A p ) 0≤p≤3 and (B p ) 0≤p≤3 that are consequences of these symmetry relations.

Zeroth order terms
The structure of A 0 and B 0 is investigated in the following, based on identities (2.23) and (2.24).
Let (i, j, k, l) ∈ N 4 be any permutation of (1, 2, 3, 4). Writing (2.23) with η = e k , we obtain: so that A 0 e k is an eigenvector associated with the eigenvalue 1 of S ij . Recall that S ij is the reflection that maps the arm i onto arm j, leaving arms k and l invariant. This implies that A 0 e k ∈ span{z k , z l }, but, since l can be otherwise arbitrarily chosen, we deduce that A 0 e k = α 0 k z k for some α 0 k ∈ R. This is true for all k ∈ N 4 . Writing now (2.23) with η = e i leads to which can also be written as 4 is the orthonormal basis of eigenvectors of G 0 defined in (2.19), and α 0 ∈ R.
As far as B 0 is concerned, we again consider any permutation (i, j, k, l) of (1, 2, 3, 4), and write (2.24) with η = e k . We obtain: Since B 0 e k ∈ Skew 3 (R), it may be decomposed using the basis L given in (2.22) as B 0 e k = p∈N 3 β k p L p so that we have, for all u ∈ R 3 : Let us define β k = p∈N 3 β k p e p ∈ R 3 (vector in R 3 whose coordinates are the ones of B 0 e k in L). Since L p u = e p × u for any u in R 3 , we obtain which implies that β k ∈ span{z k , z l } ⊥ . The index l = k being arbitrary, we obtain that β k = 0, and consequently, B 0 e k = 0. Eventually, since k is arbitrary, we finally deduce

First order terms
In order to understand the structure of the matrices (A p ) p∈N 3 and (B p ) p∈N 3 , we split them into their symmetric and skew-symmetric parts, namely In view of (2.25) and (2.26) they satisfy the identities (for all i, j ∈ N 4 , i = j) We also introduce for all m, n ∈ N 4 that store the (m, n)-entry of the three matrices (A + p ) p∈N 3 and similarly for ( For the latter, we also observe that for any u ∈ R 3 and p ∈ N 3 , since L p u = e p ×u, one has 3 , being skew-symmetric, have vanishing diagonal coefficients. For the extradiagonal coefficients (i, j) with i = j, we write (2.28) with ξ = e i and η = e j . Using the fact that A − p is skew-symmetric leads to

Study of the skew-symmetric parts
This is true for any choice of couple (i, j). Notice that α − ij is unique up to the choice of a sign that depends on the order that we choose for the two vectors z k and z l in the cross product. We may therefore fix this sign by assuming furthermore that sgn(i, j, k, l) = 1. Now, we rewrite (2.28) with ξ = e i and η = e k and obtain which we deduce that the α − ij are all equal to α − ∈ R. (Notice that sgn(j, k, i, l) = sgn(i, k, l, j) = sgn(i, j, k, l) = 1.) By explicitly calculating the cross products z i × z j , we find that the corresponding matrices (A − p ) p∈N 3 can be expressed with respect to α − as For the matrices B − p , p ∈ N 3 , we again consider any fixed permutation (i, j, k, l) of (1, 2, 3, 4) and write (2.29) with ξ = e k and η = e l . We obtain, using (2.30) for any u ∈ R 3 Similarly, taking ξ = e j and η = e l leads to B − il = −S ij B − jl giving that the coefficients β − kl are all identical. We thus deduce that B − kl = β − z k × z l , for some β − ∈ R and any k = l ∈ N 4 .
From a matricial point of view, this now yields Study of the symmetric parts A + p and B + p , p ∈ N 3 . The study of the symmetric parts follows the same lines as for the skew-symmetric ones. First, we write (2.28) with ξ = e i and η = e j and use the fact that A + p is symmetric. This leads to To prove that all the constants are equal, we finally write , so that there exists α + ∈ R such that α + ij = α + for all i = j. Concerning the diagonal terms A + ii , we find in a similar way that A + ii ∈ span{z i , z j } ∩ span{z i , z k } ∩ span{z i , z l } so that there existsᾱ + i ∈ R such that A + ii =ᾱ + i z i and we prove as before that all the constants are equal: . It now remains to identify the matrices B + p . To do so, we use the same symmetries to check that the vectors B + ij being proportional to both z k × z l and z i × z j must vanish. Similarly, we can prove that B + ii = 0 for all i. Finally, we obtain the following expressions for the symmetric parts of the matrices:

34)
To finish this section, we conclude that the dynamics of SPr4 to first order on the amplitude of the strokes around ξ 0 only depends on the five parameters α 0 , α − , α + ,ᾱ + and β − that appear in the preceding formulas. We will compute those parameters explicitly in Section 4.1, in the asymptotic regime of large arms l 0 compared to the diameter a of the balls.

The linearized control equations
We know from the previous section, that the control system (2.1) governing the evolution of SPr4 around ζ = ξ 0 (i.e. around ξ = 0), up to higher order terms, simplifies to (2.36) In particular, fixing ξ ∈Ḣ 1 (I, R 4 ) and defining Γ(t) := p∈N 3 (B pξ (t)·ξ(t))L p : I → Skew 3 (R), the dynamics of R can be written as an ordinary differential equation on the Lie group SO (3): (2.37) To simplify equations (2.36) further, we are interested in the solution of (2.37) in the regime of a small stroke ξ ∈Ḣ 1 (I, R 4 ) or equivalently in the regime of uniformly small in time matrices Γ. Intuitively, the solution R of (2.37) should not deviate too much from the initial value R 0 if the vector field Γ driving the differential equation stays small. More precisely, from [1, p. 31] we have Hence, choosing R 0 = Id 3 yields in particular the following approximations for any t ∈ I : for ξ ↓ 0. We now integrate both previous relations over I. Using that A 0ξ vanishes due to the periodicity of the stroke ξ, together with the fact that the term involving Γ is of order 2 in ξ, namely Denoting by f := (2π) −1 I f (s)ds the average of f ∈ L 2 (I), we obtain the following result: Proposition 2.7. For any ξ ∈Ḣ 1 (I, R 4 ), in a neighborhood of 0 ∈Ḣ 1 (I, R 4 ), the following estimates hold: We also note that A pξ · ξ = A − pξ · ξ and B pξ · ξ = B − pξ · ξ for all p ∈ N 3 , so that the constraint can be written, to leading order:

The bivector formalism
We now reformulate the optimization problem (2.5), and more precisely the constraints, to leading order, in the framework of bivectors. For details about bivectors, we refer to [17], but in a nutshell, one can think of (simple) bivectors of R 4 as oriented plane segments, which are obtained from the exterior product u ∧ v of two vectors u, v ∈ R 4 . Linear combinations of such simple bivectors form a linear space of dimension six denoted by 2 R 4 , to which we will later identify the space of (sufficiently small) net displacements R 3 × Skew 3 (R). This approach allows us to handle both translations and rotations in a single formalism and will prove extremely useful for the understanding of the structure of the solutions. Recall that we have denoted by (e 1 , e 2 , e 3 , e 4 ) the canonical basis of R 4 , then (e 12 , e 13 , e 14 , e 23 , e 24 , e 34 ) forms a basis of 2 R 4 , where we have set e ij := e i ∧ e j to simplify notation. The space 2 R 4 can be equipped with a scalar product by linearly extending the Gramian determinant, i.e. for two simple bivectors u ∧ v and u ∧ v we set and then we extend to the entire space by linearity. Subsequently, we also have a norm on the space 2 R 4 that we denote by | · |. In particular, we have the identities and with equality if and only if u and v are orthogonal. We emphasize at this point that not every bivector of R 4 is simple, i.e. can be represented as a single wedge product u ∧ v of two vectors u and v of R 4 . For example, the bivector e 1 ∧ e 2 + e 3 ∧ e 4 ∈ 2 R 4 cannot be written as a single wedge product. It is however true, as the preceding example suggests, that every bivector of R 4 can be decomposed into the sum of at most two orthogonal simple bivectors. Moreover, this decomposition is unique, meaning that the pair of simple bivectors forming the sum is uniquely defined, if and only if the norms of both summands are different. On the other hand, in the non-unique case there are infinitely many such decompositions, which is connected to the relationship between bivectors and rotations [17]. This fact will play a crucial role in the characterization of the optimal control curves. However, note that this decomposition property is not restricted to the standard scalar product but extends to any choice of scalar product in R 4 , e.g. the one induced by G 0 .
Bivectors can also be identified with skew-symmetric matrices. Indeed, we may define the bijective map Ω : 2 R 4 → Skew 4 (R) by linearly extending on all 2 R 4 the map defined on simple bivectors This gives us a way to represent any bivector by a skew-symmetric matrix, and vice versa. We deduce from (2.42) that A short calculation shows that the skew-symmetric parts A − p and B − p of the matrices A p and B p may be represented as simple bivectors through the mapping Ω as The constraints (2.39) then become and similarly using the notation τ ij := τ i ∧ τ j . Therefore, the isomorphism sending the basis ( e 1 , e 2 , e 3 , L 1 , Finally, the optimal control problem (2.5) can be rewritten, to first order, using the bivector formalism as: the displacement δp ∈ 2 R 4 being given, find inf

The main result
In this section, we present the main result of this paper, which solves the optimization problem (2.51), and we lay out the proof thereof.

Statement of the theorem
Theorem 3.1. Let ω ∈ 2 R 4 . We consider the following optimization problem: Then, we distinguish the following two cases: The bivector ω is simple. Then, any minimizer ξ ∈Ḣ 1 (I, R 4 ) of (3.1) is of the form where ω = u ∧ v, u and v are G 0 -orthogonal and of the same G 0 -norm.
(B) The bivector ω is not simple. Then, any minimizer ξ ∈Ḣ 1 (I, R 4 ) of (3.1) is of the form where ω = u 1 ∧ v 1 + u 2 ∧ v 2 , the vectors u 1 , v 1 , u 2 and v 2 are pairwise G 0 -orthogonal, the two pairs of vectors (u 1 , v 1 ) and (u 2 , v 2 ) are respectively of the same G 0 -norm and they satisfy furthermore the condition The proof will be decomposed in the remaining of this section. We split it in three different steps, namely, the existence of optimal periodic strokes, and the proofs of parts (A) and (B), respectively.

Proof of Theorem 3.1 -Existence of optimal periodic strokes
We first show that the solution of the minimization problem (3.1) exists. This is a consequence of the direct method of the calculus of variations. Indeed, let us consider the set of strokes satisfying the constraints Assuming for the time being that H = ∅, the existence of a minimizer for the energy is guaranteed since, from the Rellich theorem, the constraint is continuous under weak H 1 convergence.
In order to prove that the set H is non-empty, we recall that any bivector ω can be written as where the four vectors u 1 , v 1 , u 2 , v 2 may be furthermore assumed to be orthogonal. Note that in the case of a simple bivector, one has u 2 = v 2 = 0. Consider the stroke We now compute from which we deduce that ξ ∈ H.

Proof of Theorem 3.1 -Part (A)
We now consider the case where the displacement ω is a simple bivector that can be written as ω = u ∧ v where u and v are two G 0 −orthogonal vectors, with the same G 0 −norm.
Considering the change of variable η = G  0 v), we deduce that the minimization problem can equivalently be rewritten as : Now, we set ω = (G 0 v are orthogonal and possess the same norm. Let η ∈Ḣ 1 (I, R 4 ) be a stroke that satisfies the constraint Iη (t) ∧ η(t)dt = ω . We emphasize that since η is a 2π periodic and null averaged curve. The energy I |η| 2 dt is therefore bounded from below by the norm of the constraint. We shall see that it is possible to get equality. For this to happen, we must have equality in the Poincaré inequality, which is only possible if the periodic curve η has only the fundamental mode in Fourier series, i.e.
We deduce that the problem has a unique solution up to any rotation in the plane defined by ω , one of which being given by We finish the proof by writing back ξ in terms of η as

Proof of Theorem 3.1 -Part (B)
Here, we assume that ω = 0 is not a simple bivector, or, in other words, that the decomposition (3.2) of ω involves four non-colinear vectors ω = u 1 ∧ v 1 + u 2 ∧ v 2 , and we may further assume without loss of generality that Using the same change of variable η = G  Writing the Euler-Lagrange equations associated with this minimization problem under constraints (3.5) leads for any variation δ ∈Ḣ 1 (I, R 4 ) to Iη ·δ dt = λ : where λ ∈ 2 R 4 is the Lagrange multiplier associated with the constraints. This equation in only valid, i.e. λ exists, whenever the constraints are qualified [13]. This is the case if the gradients of the constraints are linearly independent, or, in other words if µ : Assuming the left-hand side, and noticing that (see (2.44)) µ : we deduce that this vanishes for all δ ∈Ḣ 1 (I, R 4 ) if and only if Ω(µ)η vanishes identically for all t ∈ I. Now, the matrix Ω(µ) being a 4 × 4 skew-symmetric, we may block diagonalize as where the matrix P µ is orthogonal. We therefore only have the following possibilities: • Either Ω(µ) is invertible (φ µ = 0 and ψ µ = 0), and thenη vanishes identically. The stroke η being null average, we deduce that η = 0 which is not possible, since ω = 0.
Furthermore, noticing that η is a 2π-periodic function, we deduce that φ λ and ψ λ must be non-vanishing integers that, without loss of generality, and up to a change of sign of u 1 and u 2 , we may furthermore assume to be both non-negative. They are also distinct otherwise there would be only two terms for η in (3.8), and the associated constraint would be a simple bivector which is not the case by assumption.
Computing the energy and the constraint in terms of η using (3.8), we obtain respectively and since φ λ and ψ λ are different integers. Now, the decomposition of a bivector as the sum of two simple bivectors is unique whenever the norm of both simple bivectors are different. We therefore have the following alternative, 0 v 2 |, and we may assume for instance that |G 0 v 2 |, the other case being treated similarly. Up to a possible exchange of φ λ and ψ λ , we deduce that Remember that ( u 1 , v 1 ) must be orthogonal and possess the same norm, as well as ( u 2 , v 2 ). But, we know, from (3.4) that (G This choice is unique up to possible rotations in the plane (G In that case, the decomposition is not unique and there may be several optimal curves among them the ones given by But other curves constructed from a different decomposition of the constraint and the method above are also satisfactory.

Back to the original problem
The first order problem (2.51) can be solved using the previous theorem with ω = Λ −1 δp where (2.50)). The distinction between parts (A) and (B) of the theorem can be rephrased in terms of the simplicity of the displacement δp directly: due to the block diagonal structure of Λ, we have that ω is simple if and only if δp is simple. Indeed, consider ω = u ∧ v a simple bivector and express u and v in the basis (τ i ) i∈N 4 Then, identifying vectors and bivectors with their coordinates in their respective bases, the following short calculation gives so that δp is a simple bivector. The converse holds writing ω = Λ −1 δp and using the same argument.
Finally, we deduce that Theorem 3.1 with ω = Λ −1 δp enables us to solve the original problem (2.51). In case δp = u 1 ∧ v 1 is simple, part (A) of the theorem applies for

Numerical simulations in the long arm regime
From the above, the solution to the first order optimization problem (2.51) can be explicitly computed provided that Λ, and thus α − and β − , are known. The optimal trajectories can then be reconstructed by computing the solutions to the first order dynamical system (2.36), which now requires to know the matrices A p , B p for p = 0, . . . , 3. From (2.27), (2.31), (2.32), (2.33), recall that these matrices only depend on a set of five parameters: α 0 , α − , α + ,ᾱ + and β − . In a similar fashion to [5], we explain in this section how to determine the asymptotic expansion of these parameters in the regime where the arms of the swimmer are assumed to be very long compared to the radii of the balls. This will enable us afterwards to develop relevant numerical test cases.

The long arm regime
In order to proceed we come back to the original optimization problem (from which the different matrices were deduced as first order terms for small strokes) and compute explicitly the asymptotic expansion of F (R, ζ) that appears in (2.1), i.e.
Let us recall the classical steps that lead to this linear relationship betweenṗ andζ. The starting point consists in taking advantage of the linearity of the Stokes equations, together with the absence of inertia. Indeed, this provides us with a linear relationship between the velocities u = (u i ) i∈N 4 of the balls (B i ) i∈N 4 and the forces f = (f i ) i∈N 4 that the fluid applies to them f = Res(ζ)u . (4.1) Note that the matrix Res(ζ) is not explicit and depends on the shape ζ through the resolution of the Stokes equations in the fluid domain. However, once Res(ζ) is known, the linear relationship linkingṗ andζ that we expressed as (2.1) is a consequence of the following explicit calculations.

Orders of convergence
For the convergence experiment, we have rescaled the net displacements with a parameter ∈ (10 −1 , 1), which yields a corresponding optimal control curve ξ . The preceding comparison of the final displacement is plotted with respect to the magnitude of the stroke in  In the left picture that shows the error in translation, the predicted order of convergence is recovered except for the case δp 1 , where we observe that the computer accuracy is already reached from the beginning. In the right picture that gives the error in rotation, the right order of convergence is also recovered, together with the same behavior as before for the case δp 1 . We strikingly also observe a sudden increase of the accuracy for the cases corresponding to simple bivectors (δp 2 and δp 3 ) whereas the non-simple case does not show this kind of superconvergence. This would need further investigation to be explained.

Trajectories
Let us discuss the two trajectories associated with the simple displacements δp 1 and δp 2 rescaled by = 10 −3 in order to be in the small displacement regime. The results are given in In the first case, where we impose a vertical displacement without rotation, we observe that the swimmer only moves along the vertical axis and does not experience any rotation during the whole stroke. We have represented the z-coordinate of the swimmer: the final displacement is achieved, through a round trip trajectory, which is a classical behavior for this kind of systems. For the pure z-rotation case δp 2 , the picture shows a circular motion in the xy-plane of the swimmer that goes back to its initial position at the end of the stroke. As far as the rotation is concerned, we have represented the swimmer along the trajectory, magnifying the rotations to get a clearer picture. As expected, a pure z-rotation is clearly visible between the initial and final states. Both trajectories so far reflect the typical behavior of a Stokesian micro-swimmer  Trajectory of SPr4 associated with the optimal control curve realizing a net displacement proportional to δp 2 , i.e. a rotation around the z-axis. The initial and final positions are marked by neon green and purple arms, respectively.
(cf. [5]). Finally, the trajectory corresponding to a net displacement proportional to δp 3 is very similar and thus omitted.
More involved trajectories can also be obtained. As an example, we show in Figure 4.4 the optimal trajectory found for a prescribed displacement δp = (2, 2, 0, 0, 0, 1). We notice that this displacement is still simple, but nevertheless produces a non obvious optimal trajectory.
Let us now consider the non-simple net displacement δp 4 (screw motion along z) that we have again rescaled, this time by a factor 10 −4 . Figure 4.5 shows the corresponding trajectory, which is exactly the superposition of the oscillation along the z-axis and the circular motion in the xy-plane corresponding to δp 1 and δp 2 .
As shown, the trajectory becomes quite more complicated for such imposed (non-simple) displacements. Notice also that there is no specific reason for the axis of the imposed rotation to be aligned with the translation. A general translation, in any direction, can be imposed together with a general rotation with any axis. In such a context, when the periodic stroke  is applied several times, the artificial swimmer experiences an helical trajectory. This kind of helical trajectory has been heavily discussed in the literature and observed in real systems (see for instance [9,20] and references therein).

A non-uniqueness result
Let us finish this discussion with one of the most interesting features of the optimal control problem (2.51). Namely, the fact that for a certain class of prescribed net displacements, the corresponding optimal control curve is not unique. Those displacements are precisely the ones which are associated with a non-simple bivector whose orthogonal summands are of the same norm. Indeed, since there are infinitely many decompositions of such a bivector into the sum of two orthogonal simple bivectors, Theorem 3.1 yields infinitely many optimal control curves in this special case.
A general construction of such decompositions is beyond the scope of this paper. However, let us consider the following basic example: consider a net displacement δp such that i.e. ω is proportional to the right-hand side. Then, we easily find a second decomposition of ω by hand. Indeed, one has for instance This yields two distinct optimal control curves and thus two distinct trajectories for the swimmer SPr4, which are presented in Figure 4.6. As one can see, the two trajectories are distinct, while the net displacement remains the same.
Notice that yet another decomposition (there exist infinitely many) is also possible such as that would give another optimal control curve.

Final comments
The method described in this paper has been incorporated in a Python code from which the figures given before were computed. This code is available under the MIT license, in the GitLab of EPFL at the address https://gitlab.epfl.ch/weder/spr.