arXiv:1603.05244v1 [math.OC] 16 Mar 2016 Three-Dimensional Multi-Tethered Satellite Formation with the Elements Moving Along Lissajous Curves D. Yarotsky∗ V. Sidorenko† D. Pritykin‡ Abstract This note presents a novel approach to maintain three-dimensional multi-tethered satellite formation in space. For a formation consisting of a main body connected by tethers with several deputy satellites (the so-called “hub-and-spoke” configuration) we demonstrate that under proper choice of the system’s parameters the deputy satellites can move along Lissajous curves in the plane normal to the local vertical with all tethers stretched; the total force due to the tension forces acting on the main satellite is balanced in a way allowing it to be in relative equilibrium strictly below or strictly above the system’s center of mass. We analyze relations between the system’s essential parameters and obtain conditions under which the proposed motion does take place. We also study analytically the motion stability for different configurations and whether the deputy satellites can collide or the tethers can entangle. Our theoretical findings are corroborated and validated by numerical experiments. Keywords: Tethered Satellite System, Satellite Formation, Dynamics, Control, Stability 1 Introduction Three-dimensional satellite formations are often discussed in connection with multi-point measurements needed for atmospheric, geodetic or plasma physics studies. To simplify control strategies and to minimize fuel consumption, tethers can be used to maintain desired relative positions of satellites in the formation flying. For the first time three-dimensional multi-tethered formations were discussed probably by [6], who proposed double-pyramid configurations. It seems that the most straightforward way to obtain a multi-tethered formation is to deploy from the main satellite several tethers with deputy satellites at their ends. To ∗ Institute for Information Transmission Problems, Russian Academy of Sciences, Bolshoy Karetny per. 19, 127051 Moscow, Russian Federation † Keldysh Institute of Applied Mathematics, Russian Academy of Sciences, Miusskaya Sq., 4, 125047 Moscow, Russian Federation ‡ Moscow Institute of Physics and Technology, Institutskiy per., 9, 141700, Dolgoprudny, Moscow Region, Russian Federation 1 specify such formations [13] introduced the term “hub-and-spoke”. Behavior of “hub-andspoke” multi-tethered formations has been studied for different dynamical environments: in circular orbit [1, 4], in elliptic orbit [5], in halo-orbit [7, 18] and near collinear Lagrangian points [17]. To keep the tethers taut the combination of rotation with gravity-gradient forces is usually proposed. Among other opportunities the relatively new concept of the Tethered Coulomb Structure (TCS) is worth mentioning [12,14,15]. In this case the satellites are electrostatically charged to produce repulsive forces between them. Nevertheless, it looks as if Coulomb repulsive forces can be used to prevent the slack of short enough tether: in [12, 14, 15] the discussed length is 10 m by order of magnitude. To give an idea of our approach we begin with the system of two bodies connected by a single tether aligned along the local vertical; the mass center of the system moves in a circular orbit. Small oscillations of this system around local vertical are a combination of in-plane and out-of-plane natural oscillations with incommensurable frequencies [8, 16]. The motion of end bodies in these oscillations can be roughly described as a motion along curves densely filling certain areas on planes normal to the local vertical. Then let us consider the degenerate “hub-and-spoke” configuration in the relative equilibrium with all tethers aligned along the vertical (Fig. 1, left) The displacement of a single deputy satellite from the relative equilibrium position causes oscillations whose frequencies differ from those inherent in the preceding case. With the proper tuning of the system’s parameters these frequencies can be made commensurable resulting in the motion of the deputy satellite along a Lissajous curve in the plane, normal to the local vertical (Fig. 1, right). Naturally all other deputy satellites can also be put in motion along similar curves. The choice of initial conditions allows to avoid collisions among them and to ensure the balance of the tension forces applied to the main satellite so as to preserve its relative equilibrium. We suppose that the described structure can be useful for some applications or at least become the starting point for the development of new approaches to maintain 3D multitethered satellite formations in space. In Section 2 we start our study with a simplified dynamical model of multi-tethered formation (point masses + weightless tethers). In Section 3 we consider the linearized dynamics near vertical equilibrium and describe the oscillations in the system. In Section 4 we consider deputy satellite formations moving along Lissajous curves. In Section 5 we present the results obtained by numerical simulation of the system’s dynamics. 2 Deputy satellite dynamics at small deviations from the relative equilibrium As mentioned in the Introduction, we will consider the multi-tethered satellite formation comprising N +1 bodies: the main satellite C of mass mC and N deputy satellites D1 , ..., DN (each of mass mD ) linked to the main satellite (but not to each other) by identical extensible tethers; the tethers’ masses are ignored. 2 D2 D2 D3 D1 D3 D1 O y x z C C Figure 1: Multi-tethered formation consisting of the main body C and three deputy satellites D1 , D2 , D3 . On the left the system is in relative equilibrium with all tethers aligned along the local vertical. On the right the deputy satellites move along a Lissajous curve (see also the animation provided as Electronic Supplementary Material). To write down the equations of motion, we shall introduce a Local Vertical Local Horizontal (LVLH) reference frame Oxyz, centered on the position of the system’s center of mass (CoM) in its nominal orbital motion: Oz axis is aligned with the local vertical and oriented towards Earth’s center, Ox runs tangentially to the orbit in the direction of the CoM motion, and Oy axis is directed along the normal line to the orbit plane (Fig. 1). It is assumed that nominally the system’s center of mass moves in circular orbit with the mean motion ω0 . The unit vectors in the directions of the axes x, y, z will be denoted as ex , ey , ez respectively. In the LVLH frame the motion of any satellite in the considered formation can be described by the Hill-Clohessy-Wiltshire (HCW) equations [9]: ẍ = 2ω0 ż + Tx /m, ÿ = ω02 y + Ty /m, z̈ = −2ω0 ẋ + 3ω02 z + Tz /m, (1) where Tx , Ty , Tz denote the components of the (total) tether tension T applied to a given body of mass m. We shall adopt the usual visco-elastic model of massless tether, hence the tension force applied to the ith deputy satellite is   d (rC − ri ) Ti = 1(|rC −ri |>l0 ) · k(|rC − ri | − l0 ) + b |rC − ri | , (2) dt |rC − ri | 3 where l0 is the slack tether length, k the elastic coefficient, b the damping coefficient, and ( 1, if |rC − ri | > l0 1(|rC −ri |>l0 ) = 0, otherwise The vertical equilibrium of the system takes place if all tethers rest (ṙC = ṙi = 0) aligned along the local vertical: yi∗ = yC∗ = 0,  −1 mC ∗ N m D + mC mD ω02 ∗ zi = − . z = l0 −3 N mD C mC k x∗i = x∗C , (3) Here and below equilibrium quantities are marked by the asterisk *. In contrast to the y and z coordinates, the x coordinate in this configuration, though shared by all satellites, is arbitrary, since the equations are invariant with respect to translations along the orbit (i.e., along the x axis in the approximation of the orbital dynamics, provided by the HCW equations). For the above configuration to be valid the denominator in the last formula must be positive, i.e. the tethers must be sufficiently rigid to counteract the microgravity: mC mD . (4) N m D + mC Denoting the tether’s length in the equilibrium configuration by l∗ = |r∗C − r∗i |, we obtain a useful relation k > 3ω02 l∗ − l0 3ω 2 mC mD = 0 . l∗ k N m D + mC (5) Denote 3ω 2 mC mD l∗ − l0 = 0 m r , mr = . l∗ k N m D + mC The tension forces in all the tethers in the equilibrium configuration (3) have the same value T∗ given by the obvious formula λ∗ = T∗ = k(l∗ − l0 )ez . (6) It is not difficult to derive an approximate expression for the tension forces in case of small displacements of the main satellite and the ith deputy satellite with respect to equilibrium (3): Ti =T∗ + kλ∗ (∆xC − ∆xi )ex + kλ∗ (∆yC − ∆yi )ey + + [k(∆zC − ∆zi ) + b(∆żC − ∆żi )] ez . Here ∆xC , ∆yC , ∆zC are the components of the main satellite’s displacements, and ∆xi , ∆yi , ∆zi are the components of the ith deputy satellite displacements. 4 With the linearized tension the HCW equations read k (∆xC − ∆xi ), mD k (∆yC − ∆yi ), ∆ÿi = −ω02 ∆yi + λ∗ mD b k (∆zC − ∆zi ) + (∆żC − ∆żi ) ∆z¨i = −2ω0 ∆ẋi + 3ω02 ∆zi + mD mD ∆ẍi = 2ω0 ∆z˙i + λ∗ (7) for the deputy satellites, and N k X ∆x¨C = 2ω0 ∆z˙C − λ∗ (∆xC − ∆xi ), mC i=1 ∆y¨C = −ω02 ∆yC N k X − λ∗ (∆yC − ∆yi ), mC i=1 ∆z¨C = −2ω0 ∆x˙C + 3ω02 ∆zC (8) N N b X k X (∆zC − ∆zi ) − (∆żC − ∆żi ) − mC i=1 mC i=1 for the main satellite. 3 Decoupling the equations of motion The system of equations derived above can be split up into three independent groups by taking appropriate linear combinations, separately for x, y and z components: 1. By taking the sum of the equation for the main satellite with weight mC and all the respective equations for deputy satellites with weights mD we obtain a triple of scalar equations for ! N X 1 mC ∆rC + mD ∆ri , N m D + mC i=1 i.e., for the motion of the whole system’s center of mass. 2. By taking the sum of the equations for deputy satellites with coefficients 1/N and subtracting from it the equation for the main satellite we obtain a triple of scalar equations for N 1 X ∆ri − ∆rC , N i=1 i.e. describing the relative motion between the main satellite and the auxiliary satellites’ center of mass. 5 3. Finally, if we take any linear P combination of equations for deputy satellites with N some coefficients σi such that i=1 σi = 0, we obtain a triple of scalar equations PN for i=1 σi ∆ri , which can be interpreted as a partial description of the relative motion between deputy satellites. To obtain the full description, we need to consider all N − 1 linearly independent such assignments of coefficients, thus giving the total of 3(N − 1) scalar equations. The simplest example of the set with suitable combinations of coefficients is {−1, 1, 0, . . . , 0}, {0, −1, 1, 0, . . . , 0}, . . . , {0, . . . , −1, 1}. Physically it means that we use relative displacements ∆ri+1 − ∆ri (i = 1, N − 1) of consecutive objects to describe the dynamics of the subsystem composed of the deputy satellites. One more opportunity is provided by the set     1 1 1 1 1 1 ,..., − ,...,1 − ,− , 1 − ,− ,...,− N N N N N N characterizing relative displacements of the deputy satellites D1 , . . . , D(N −1) with respect to the center of mass of all deputy satellites. Clearly, the collection of these three groups of equations is equivalent to the original system of 3(N + 1) scalar equations for the main and deputy satellites. Let us deal with the three groups one by one. 3.1 System CoM motion equations The motion of the system’s center of mass is described by free HCW equations [9]. It is well-known that these equations describe oscillations with frequency ω0 in the y component, and oscillations with frequency ω0 combined with a linear drift in the orbital xz plane. 3.2 Relative motion between the main satellite and the deputy satellites’ center of mass Let ∆x, ∆y, ∆z denote the scalar components of 1 N PN i=1 ∆ri − ∆rC . Then using relation (5) k ∆x = 2ω0 ∆ż − 3ω02 ∆x, mr k ∆ÿ = −ω02 ∆y − λ∗ ∆y = −4ω02 ∆y, mr b k ∆z − ∆ż. ∆z̈ = −2ω0 ∆ẋ + 3ω02 ∆z − mr mr ∆ẍ = 2ω0 ∆ż − λ∗ (9) The equation for ∆y is independent from the equations for ∆x, ∆z and describes harmonic oscillations with the frequency 2ω0 . To analyze the remaining equations for ∆x, ∆z we consider the corresponding first-order system 6   0 1 0 ∆x 2    d  ∆ẋ   −3ω0 0 0 =   0 0 0 ∆z dt 2 0 −2ω0 3ω0 − ∆ż  k mr  0 ∆x   2ω0   ∆ẋ 1   ∆z − mbr ∆ż   .  (10) Let us denote by A the matrix on the right-hand side of the system (10). To examine the stability property of this system, we write down the characteristic equation     b 4 3bω02 k k 4 2 2 2 2 det(A − ρI4 ) = ρ + ρ + + 4ω0 ρ + ρ + 3ω0 − 3ω0 = 0. mr mr mr mr The symbol I4 is used here for the identity matrix of the fourth order. Applying the Routh-Hurwitz stability criterion [10], it is not difficult to establish that all eigenvalues ρ of A belong to the left half-plane Reρ < 0 iff the condition (4) is satisfied and b > 0. Summarizing, under the stability assumption (4) the center of mass of deputy satellites oscillates with respect to the main satellite; in the linear approximation the energy dissipation in tethers does not affect the oscillations along the y axis. 3.3 Relative motion between deputy satellites Let ∆x, ∆y, ∆z denote the scalar components of P ficients σi subject to N i=1 σi = 0. Then PN i=1 σi ∆ri with some assignment of coef- k ∆x, mD k ∆ÿ = −ω02 ∆y − λ∗ ∆y, mD ∆ẍ = 2ω0 ∆ż − λ∗ ∆z̈ = −2ω0 ∆ẋ + 3ω02 ∆z − (11) k b ∆z − ∆ż. mD mD The equation for the y component describes oscillations with frequency r r 4mC + N mD k 2 = ω0 . ωy = ω0 + λ∗ mD mC + N m D The remaining equations for ∆x, ∆z can be written in the first-order form as      0 1 0 0 ∆x ∆x k     0 0 2ω0  d    ∆ẋ  .  ∆ẋ  =  −λ∗ mD 0 0 0 1   ∆z  dt  ∆z   k 2 ∆ż ∆ż 0 −2ω0 3ω0 − mD − mbD The matrix A1 on the right-hand side of (13) has the characteristic polynomial     b 3 k λ∗ kb λ∗ k (λ∗ + 1)k 4 2 2 2 det(A1 − ρI4 ) = ρ + ρ + + ω0 ρ + 2 ρ + − 3ω0 . mD mD mD mD mD 7 (12) (13) Applying again Routh-Hurwitz stability criterion, we establish that the system (13) is asymptotically stable iff b > 0 and k ≥ 3ω02 mD . (14) If the condition (14) is fulfilled, but the dissipation is absent (b = 0), then the spectrum of A1 is purely imaginary. If the condition is violated, matrix A1 has a positive eigenvalue. Note that the stability condition (14) is stronger than the earlier condition (4). Assuming b = 0, in the limit of large rigidity k the linearized dynamics of ∆x, ∆z approximately decouples into independent oscillations of ∆x and ∆z with frequencies r  
mD ω02 λ∗ k −2 1−2 = ωx = +O k mD k r  
3mC mD ω02 −2 (15) 1−2 +O k = ω0 , N m D + mC k r  
k mD ω02 −2 . 1+ +O k ωz = mD 2k 3.4 Lyapunov function Our solution of the linearized dynamics implies, in particular, that in the linear approximation of tether tension the vertical equilibrium is stable. This conclusion can also be shown in a stronger sense – without linearization of tether tension – by directly providing a Lyapunov function. Specifically, let r′ = r − rCoM , ṙ′ = ṙ − ṙCoM denote position and velocity relative to the system’s center of mass. Consider the energy E = T + V of relative motion, where N mD X ′ 2 mC ′ 2 ṙ ṙ + T = 2 i=1 i 2 C and N  2 mD ω02 X ′ 2 mC ω02  ′ 2 2 V = (y i − 3z ′ i ) + y C − 3z ′ C + 2 i=1 2 N + kX 2 1(|r′ C −r′ i |>l0 ) · (|r′ C − r′ i | − l0 ) . 2 i=1 If b = 0, then Ė = 0, otherwise Ė ≤ 0. Taylor expansion of V near the equilibrium yields 8 N   m ω2 mD ω0 2 X  ′ 2 C 0 ′ 2 ′ 2 ′2 ∆yC − 3∆zC + ∆yi − 3∆zi + V = Vequilibrium + 2 2 i=1 N io h k Xn 2 2 2 + (∆zi′ − ∆zC′ ) + λ∗ (∆x′ i − ∆x′ C ) + (∆x′ i − ∆x′ C ) 2 i=1 (16) 3 + O(|∆r′ | ) Taking into account the identity ′ ∆r C N mD X ′ ∆r i , =− mC i=1 one can check that this form is positive definite exactly if stability conditions (4) and (14) hold. 4 Motion of deputy satellites along Lissajous curves Having described oscillations of the “hub-and-spoke” system near the equilibrium, we consider now the possibility of the satellites moving in such a way that the system elements (satellites and tethers) never collide. Given our assumption of small deviations from the vertical equilibrium, we formulate this as the requirement that the projections (xi , yi ) of deputy satellites to the xy plane never come close to each other. Results of Section 3.3 imply that the position of a deputy satellite relative to the center of mass of all deputy satellites, N 1 X rk , r′ i = ri − N k=1 oscillates with frequency r 4mC + N mD mC + N m D in the y direction and, for sufficiently rigid tethers, with frequency r 3mC ωx ≈ ω0 N m D + mC ωy = ω0 (17) (18) in the x direction (with a more accurate value P given by (15)). Since, as shown in Section 3.2, the center of mass of deputy satellites N1 N k=1 rk performs independent oscillations, we may ignore these latter and assume without loss of generality that N N 1 X 1 X xk = yk ≡ 0 N k=1 N k=1 9 (19) at all times. By linearity, for each pair i, j of deputy satellites their relative position ri − rj also oscillates with the same frequencies ωx , ωy . If the frequencies are incommensurate, i.e., ωx /ωy is irrational, then the trajectory of the oscillation is aperiodic and it comes arbitrarily close to the origin, i.e. the two satellites come arbitrarily close to each other. We are thus naturally led to consider commensurate oscillations: ωx p = , ωy q (20) where p, q are co-prime natural numbers. In this case the two oscillations have a common period of 2πp 2πq TL = = . ωx ωy It is convenient to introduce the non-dimensional time τ= t . TL The xy-motion of a single deputy satellite can then be written as x = x0 sin (2πpτ + ϕx ) , y = y0 sin(2πqτ + ϕy ), (21) with some initial phases ϕx , ϕy , and has a period 1. The trajectory of this motion is known as a Lissajous curve [3, Sec. 25] Substituting expressions (17), (18) for ωx , ωy into formula (20), we obtain the following relation between the frequency ratio and the satellite mass ratio: 3q 2 N mD = 2 − 4. mC p In particular, positivity of the left-hand side entails √ p 3 < . q 2 (22) We remark in passing that this condition excludes the usual elliptic (or circular) oscillations corresponding to p = q = 1. The frequencies ωx , ωy can then be expressed in terms of p and q: pω0 ωx = p , q 2 − p2 qω0 ωy = p . q 2 − p2 10 4.1 Balanced formations avoiding collisions We seek now formations of deputy satellites moving according to (21) subject to the following conditions: A. The arrangement of deputy satellites must satisfy the balance condition (19); B. The satellites must never collide, i.e. (xi (t), yi (t)) 6= (xj (t), yj (t)) for all t and i 6= j; C. Optionally, we may wish to ensure that (xi (t), yi (t)) 6= (0, 0) for all i and t – in this case an additional satellite can be added at the center of the “hub-and-spoke” system without collisions with this system. We consider two types of uniform arrangement of deputy satellites. Type I is a uniform arrangement of N satellites along a single Lissajous curve: the position of the i’th deputy satellite is given by         i i xi (τ ) = x0 sin 2πp τ + yi (τ ) = y0 sin 2πq τ + + ϕx , + ϕy , N N where the phases ϕx , ϕy are the same for all satellites. Type II is a uniform arrangement of N satellites along several Lissajous curves: the position of the i’th deputy satellite is given by         i i xi (τ ) = x0 sin 2π pτ + yi (τ ) = y0 sin 2π qτ + + ϕx , + ϕy , N N where the phases ϕx , ϕy are the same for all satellites. The following proposition summarizes properties of such formations with respect to the above three conditions. Proposition 1 Let N = 2, 3, . . . Denote ϕ0 = qϕx −pϕy . π a) For a Type I formation, the balance condition A is fulfilled if and only if neither p nor q is divisible by N . For a Type II formation, the balance condition A is fulfilled for all N. b) For a Type I formation, the no-collision condition B is fulfilled iff the number ϕ0 + (p − q)/2 is not an integer and N is co-prime with p and q. 11 c) For a Type II formation, in case N ≥ 3 the no-collision condition B is fulfilled if and only if (ϕ0 + (p − q)/2) N is not divisible by the greatest common divisor of N and q − p. In case N = 2 condition B is fulfilled iff ϕ0 is not an integer. d) Lissajous curve (21) goes through the origin (0, 0) iff ϕ0 is an integer. It follows that a Type I formation satisfies condition C iff ϕ0 is not an integer. A Type II formation satisfies condition C iff ϕ0 is not of the form a + 2b(q − p)/N with integer a, b. The proof of Proposition 1 is given in Appendix A. Note that, given p, q, admissible formations (satisfying all three conditions A-C) exist for all N = 2, 3, . . . in case of Type II, but not for all N in case of Type I. In Table 1 we list parameters of admissible formations for all p, q ≤ 4 subject to condition (22). p/q(= ωx /ωy ) N mD /mC 1/2 1/3 2/3 1/4 3/4 8 23 11/4 44 4/3 Admissible N for Type I 3, 2, 5, 3, 5, 5, 4, 7, 5, 7, 7, . . . 5, . . . 11, . . . 7, . . . 11, . . . Table 1: Admissible formations for small values of p, q. In Fig. 2 we show several examples of formations of Types I and II satisfying conditions A,B,C with p, q subject to condition (22). 4.2 Entanglement Practical implementation of the introduced formations requires to resolve another issue. The tethers’ ends are not attached to the main satellite at exactly the same point. Even if the tethers are connected very close to each other, they still have nonzero thickness. Relative motion of satellites in formations of Type I and Type II may cause not only contacts between the tethers, but also tethers entanglement. We will distinguish two kinds of entanglement: one that can be canceled out by rotating the main satellite about the z axis as shown in Fig. 3a, and one that can not be eliminated by such rotations (Fig. 3b). These two kinds will be referred to as weak and strong entanglement, respectively. We can define entanglement rigorously by applying simple topological concepts [2] to the dynamics of the deputy satellites. Consider the motion of a deputy satellites formation as a 12 y Type I, p=1, q=2, N=3, phi0=1/4 Type I, p=1, q=2, N=5, phi0=1/4 1.0 1.0 0.5 0.5 0.5 0.0 0.0 0.0 −0.5 −0.5 −0.5 −1.0 −1.0 −1.0 −1.0 −0.5 0.0 0.5 1.0 −1.0 Type I, p=1, q=3, N=2, phi0=1/2 y Type II, p=1, q=2, N=2, phi0=1/2 1.0 −0.5 0.0 0.5 1.0 −1.0 Type I, p=3, q=4, N=5, phi0=1/4 1.0 1.0 0.5 0.5 0.5 0.0 0.0 0.0 −0.5 −0.5 −0.5 −1.0 −1.0 −0.5 0.0 0.5 1.0 x 0.0 0.5 1.0 Type II, p=2, q=3, N=2, phi0=1/2 1.0 −1.0 −0.5 −1.0 −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 x 0.0 0.5 1.0 x Figure 2: Examples of formations of Types I and II. continuous map M0 (τ ) =  x1 (τ ) y1 (τ )      xN (τ ) x2 (τ ) , τ ∈ [0, 1]. ,... , yN (τ ) y2 (τ ) By periodicity, M0 (0) = M0 (1). Entanglement means the impossibility to “straighten” this map subject to topological constraints. Precisely, let us say that M0 is strongly homotopic to the constant map       xN (0) x2 (0) x1 (0) , τ ∈ [0, 1] , ,... , M1 (τ ) = yN (0) y2 (0) y1 (0) if there is a continuous family of maps Mα (τ ) , τ ∈ [0, 1] , α ∈ [0, 1] , that deforms M0 into M1 without collisions and subject to the periodicity constraint Mα (0) = Mα (1) , α ∈ [0, 1]. (23) If M0 is not strongly homotopic to M1 , then we call M0 weakly entangled. In order to define strong entanglement, we say that M0 is weakly homotopic to the constant map M1 if the constraint (23) is replaced by the relaxed constraint Mα (0) = Rϕ(α) Mα (1) , α ∈ [0, 1] , 13 D2 D3 D2 D1 D2 D1 D3 C D1 o o 360 turn of the D1D2D3 subsystem D3 360 turn of the main satellite C C D2 C (a) D2 D1 D3 D1 D3 C C (b) Figure 3: Examples of weak and strong entanglement. In the case of weak entanglement (a) tethers can be disentangled by rotation of the main satellite. Tethers in strong entanglement (b) can not be disentangled by rotation of the main satellite. where Rϕ(α) is the rotation of the satellite positions in the xy plane by continuously varying angles ϕ(α) common to all deputy satellites, for some choice of the continuous function ϕ(α). We then say that M0 is strongly entangled if M0 is not weakly homotopic to M1 . The above definitions do not involve the main satellite. We assume that the initial state is non-entangled, and the geometry of the tether attachment points to the main satellite is reflected in the initial positions of the deputy satellites as shown in Fig. 3. An obvious example of a weak, but not strong entanglement results from a circular motion as shown in Fig. 3a. Also note that a motion of N = 2 deputy satellites is never strongly entangled but may be weakly entangled. By considering a few simple examples, it is easy to see that entanglement does take place in some formations of Type I and Type II and does not in others. For instance, there is a strong braid-like entanglement for the Type I formation with p = 1, q = 2, and N = 3 (top left plot in Fig. 2), weak entanglement for the Type I formation with p = 1, q = 3, and N = 2 (bottom left plot in Fig. 2), and no entanglement for the Type II formations with N = 2 and p = 1, q = 2 or p = 2, q = 3, and N = 2 (top and bottom right plots in Fig. 2). Complete analysis of entanglement in our formations appears to be relatively complex mathematically, hence we will not attempt it here. We will, however, state a simple proposition involving pairwise relations between satellites. Namely, for any pair (i, j) of deputy satellites consider the winding number wi,j defined as the number of turns that the j’th satellite makes about the i’th satellite in the xy-plane over the period: wi,j 1 = arg 2π  xj (τ ) − xi (τ ) yj (τ ) − yi (τ )  1 τ =0 Clearly, wi,j = wj,i .The winding numbers are obviously invariant under a strong homotopy. 14 Under a weak homotopy, they are incremented by an amount common to all pairs and equal 1 to the number of full turns 2π ϕ(α)|1α=0 . Since winding is absent in the straightened map M1 , the numbers wi,j can serve to establish sufficient (but not necessary) conditions of entanglement. Proposition 2 a) Consider a formation of Type I or Type II moving without collisions. If at least one of p, q is even, then wi,j = 0 for all pairs. Otherwise, wi,j = ±1 for all pairs.Consequently, if both p, q are odd, then there is at least a weak entanglement. b) Moreover, consider a formation of Type I and suppose that both p, q are odd so that wi,j = ±1 by a). Then both values +1 and −1 are encountered among the winding numbers wi,j iff neither of the numbers q−p, q+p is divisible by 2N . As a consequence, if neither of q − p, q + p is divisible by 2N , then the entanglement is not only weak, but also strong. Proposition 2 is proved in Appendix B. It leads to an expected result if applied to a circular or elliptic motion (p = q = 1): part a) confirms a weak entanglement, and part b) ensures that all winding numbers are equal since p − q = 0, so the sufficient condition for a strong entanglement is not fulfilled. Note also that for odd p, q one of the numbers q − p, q + p is always divisible by 4, so statement b) agrees with our earlier remark that there can be no strong entanglement for N = 2. 4.3 Second order perturbation theory In contrast to the general stability of the “hub-and-spoke” system established in Section 3.4 in the sense of smallness of deviation from the vertical equilibrium, we do not expect the motion of deputy satellites along Lissajous curves to be stable. This motion is a subtle phenomenon which, in particular, is strongly affected by nonlinearities and can be maintained without additional control action only at relatively small oscillation amplitudes. In this section we examine the nonlinearity effects by deriving second order corrections to the evolution equations in case k/3ω02 mD ≫ 1 when one can neglect the tethers’ extensibility and put l0 = l∗ . If we assume that the center of mass of the system rests at the origin, then the configuration can be parametrized by the 2N coordinates xi , yi , i = 1, .., N, of the 15 auxiliary satellites. Specifically, the remaining coordinates are given by N N mD X mD X xi , yC = − yi , mC i=1 mC i=1 N q mr X 2 zC = l − (xi − xC )2 − (yi − yC )2 ≈ mC i=1 ) ( N 1 X mr 2 2  , N l0 − (xi − xC ) + (yi − yC ) ≈ mC 2l0 i=1 q mr zi = zC − l2 − (xi − xC )2 − (yi − yC )2 ≈ − l0 + mD ) ( N X    1  m r . + (xi − xC )2 + (yi − yC )2 − (xk − xC )2 + (yk − yC )2 2l0 mC k=1 xC = − Here and below we consistently keep terms only up to second order in xi , yi , ẋi , ẏi , ẍi , ÿi . The x, y components Ti,x , Ti,y of the tension acting on the i’th deputy satellite can be expressed through the z component by Ti,x = xi − xC xi − xC yi − yC yi − yC Ti,z ≈ − Ti,z , Ti,y = Ti,z ≈ − Ti,z . zi − zC l0 zi − zC l0 The z component Ti,z , in turn, is found from the corresponding HCW equation: Ti,z mr ≈ 2ω0 ẋi − 3ω02 zi ≈ 2ω0 ẋi + 3ω02 l0 . m mD The remaining two HCW equations then yield   xi − xC 2 mr ẍi ≈ 2ω0 żi − 2ω0 ẋi + 3ω0 l0 , l0 mD   yi − yC 2 2 mr ÿi ≈ −ω0 yi − 2ω0 ẋi + 3ω0 l0 . l0 mD Differentiating zi and retaining our notation for ωx , ωy from (17)-(18), we finally obtain  2ω0 2 ẍi + ωx (xi − xC ) ≈ − ẋC (xi − xC ) + (ẏi − ẏC ) (yi − yC ) − l0  N  mr X  (ẋj − ẋC ) (xj − xC ) + (ẏj − ẏC ) (yj − yC ) , − mC j=1 ÿi + ωy2 yi − 3ω02 2ω0 mr yC ≈ − ẋi (yi − yC ) , mD l0 where we have placed first order terms on the left and second order terms on the right. 16 The obtained equations can be used to find anharmonic corrections to a particular small harmonic oscillation. The usual procedure is to substitute the harmonic oscillation in the right-hand side and find the correction, to leading order, by solving the resulting nonhomogeneous linear equation [11]. If the right-hand side contains secular terms, i.e., those whose frequencies match some of the eigenfrequencies of the linear equation, then, additionally, the solution’s frequencies need to be adjusted to prevent its non-physical growth. In the case at hand we take the motion of a Type I or II formation with a small amplitude as a base harmonic oscillation that we denote xi,0 (t) , yi,0 (t). From the balance condition (19) we have xC (t) = yC (t) ≡ 0, so that the evolution equations simplify to # " N X 2ω m 0 r ẍi + ωx2 (xi − xC ) ≈ (ẋj,0 xj,0 + ẏj,0 yj,0 ) , ẏi,0 yi,0 − l0 mC j=1 ÿi + ωy2 yi − 3ω02 2ω0 mr yC ≈ − ẋi,0 yi,0 . mD l0 The linear terms describe oscillations with frequencies √ ωx , ωy , ωCx = 3ω0 , ωCy = 2ω0 , where the latter two correspond to the motion of the center of mass of auxiliary satellites or, equivalently, to the motion of the main satellite. The terms ẋk,0 xk,0 , ẏk,0 yk,0 , and ẋi,0 yi,0 on the right-hand side result in oscillations with frequencies 2ωx , 2ωy , and ωy ± ωx , respectively. Though the second order correction of motion is obviously present for each deputy satellite for any choice of system parameters, it is possible to choose parameters so as to make the second order correction completely vanish for the main satellite. The evolution equations for the main satellite are derived by adding up the equations for the deputy satellites: ! N N X X
2ω0 N mD mC mC 2 ẏj,0 yj,0 − ẋj,0 xj,0 , − ẍC + ωCx xC ≈ mD l0 mC + N mD j=1 mC + N mD j=1 − N
2ω0 X mC 2 ÿC + ωCy yC ≈ − ẋj,0 yj,0 . mD l0 j=1 Note that the N mD tends to be much larger than mC for a system satisfying our assumptions (see Table 1), so the motion of the main satellite tends to be generally affected by second order corrections much stronger than the motion of the deputy satellites. However, contributions from different j’s here may cancel out. Proposition 3 Consider a Type I formation such that neither of 2p, 2q, q ± p is divisible by N . Then N X k=1 ẋk,0 (t)xk,0 (t) = N X ẏk,0 (t)yk,0 (t) = k=1 N X k=1 ẋk,0 (t)yk,0 (t) ≡ 0 for all t, so that the right-hand sides of the above equations vanish identically. 17 The proof is elementary, and we omit it. Examples of parameter sets fulfilling the proposition’s hypothesis are N = 5, p = 1, q = 2 and N = 5, p = 3, q = 4. Our numerical experiments below confirm that in these cases the main satellite is indeed stable in contrast to the generic settings at the same oscillation amplitude of deputy satellites. 5 5.1 Examples and numerical simulations Setup and general observations In our numerical simulations we consider oscillations of the system with rigid tether at small angles and with equal amplitudes in x and y: x0 = y0 = a = κrad l0 Here a is the linear amplitude and κrad is the corresponding angular amplitude expressed in radians; the same angular amplitude expressed in degrees is denoted κdeg . In all our experiments the dimensionless coefficient k/3ω02 mD characterizing tether rigidity falls in the range [3 · 102 , 103 ] and κdeg ≤ 6◦ (to justify the above linear relation between a and κrad ). In order to numerically examine the stability of the Lissajous motion we introduce quantities characterizing relative deviations of the main and deputy satellites from their theoretical positions obtained in the linear approximation. Specifically, we consider the relative deviation of the main satellite’s numerically computed position from the z axis 1q 2 2 δC (t) = (t), xC,num (t) + yC,num a and the mean relative deviation of the deputy satellites’ numerically found trajectories from the theoretical Lissajous curves: N 1 X δD (t) = N a i=1 q (xi,num (t) − xi,Liss (t))2 + (yi,num (t) − yi,Liss (t))2 These quantities can be compared with the minimum distance between different deputy satellites on their theoretical trajectories: q 1 δmin = min (xi,Liss (t) − xj,Liss (t))2 + (yi,Liss (t) − yj,Liss (t))2 a t,i6=j Here we neglect the difference in z coordinates of the deputy satellites since for small amplitude oscillations this difference has higher order of smallness. If each satellite remains at all times within a distance of δmin /2 from its theoretical position on the Lissajous curve, then all satellites are guaranteed to avoide collisions. We will consider a slightly relaxed condition of stability δD (t) < 18 δmin 2 (24) that constrains only the mean deviation δD (t) of the deputy satellites. The time interval during which this condition holds can be roughly considered as a “system stability interval”. We restrict ourselves to the simplest ratio ωx : ωy = 1 : 2 achieved at p = 1, q = 2. In this case, the theoretical mass ratio is N mD = 8, mC and the full period TL of system oscillations is related to the orbital period T0 = 1/(2πω0 ) by p √ TL = q 2 − p2 T0 = 3T0 , so that the non-dimensional time t . 3T0 We consider the three different configurations of the system shown in the first row of Fig. 2. For each configuration, we perform simulations for κdeg = 1◦ and κdeg = 3◦ for 10 orbital periods. In our simulations, the central satellite moves along a geostationary orbit, and tether length l0 = 10000 m. δ τ=√ 0.30 0.25 0.20 0.15 0.10 0.05 0.00 0 0.5 δ 0.4 ϰdeg =1 ◦ δC δD δmin/2 2 4 6 8 ϰdeg =3 ◦ δC δD δmin/2 0.3 0.2 0.1 0.0 0 10 2 4 orbit number 6 8 10 Figure 4: Numerically observed relative deviations of the satellites from the theoretical trajectories for the Type I formation with N = 3 deputy satellites. In this case the minimum relative distance between deputy satellites on the theoretical trajectories is δmin = 0.60. At κdeg = 1◦ stability condition (24) holds with a large margin for the whole simulation interval of 10 orbital periods, while at κdeg = 3◦ it breaks down after six orbital periods. 19 The results for Type I formation at N = 3 (Fig. 4) show that the deviations of the main satellite are initially much larger than those of the deputy satellites. This is not surprising, since the main satellite is eight times lighter. However, deviations of the deputy satellites approximately linearly accumulate with time, and eventually catch up with those of the main satellite. In case κdeg = 3◦ relative deviations are much larger than in case κdeg = 1◦ : approximately three times larger for δC and six times larger for δD . The second order perturbation theory in Section 4.3 suggests a linear dependence of relative deviations on κdeg , but for κdeg = 3◦ deviations of the main satellite from the equilibrium position are already comparable to the oscillation amplitude, so this perturbation theory is not truly applicable here. In case κdeg = 1◦ the stability condition (24) holds with a large margin for the whole simulation interval, while at κdeg = 3◦ it breaks down after six orbital periods. 0.20 δ 0.15 ϰdeg =1 ◦ δC δD δmin/2 0.10 0.05 0.00 0 0.20 δ 0.15 2 4 6 8 ϰdeg =3 ◦ δC δD δmin/2 0.10 0.05 0.00 0 10 2 4 orbit number 6 8 10 Figure 5: Deviations for Type I and N = 5. In this case δmin = 0.43. Stability condition (24) holds with a large margin for the whole simulation interval of 10 orbital periods for both κdeg = 1◦ and κdeg = 3◦ . The results for Type I formation at N = 5 (Fig. 5) are drastically different due to the cancellation of second order corrections of the main satellite’s motion pointed out in Proposition 3. Not only is the main satellite almost immobile, but also the deviations of the deputy satellites are from two to four times smaller than in the previous case. In particular, the no-collision condition (24) holds at κdeg = 3◦ throughout the whole simulation interval. Note also that in the κdeg = 3◦ case the deviations of deputy satellites are approximately three times as large as in the κdeg = 1◦ case, in good agreement with the second order perturbation theory which is now applicable since the deviations are small. Results in the case of Type II formation at N = 2 (Fig. 6) are on the whole similar to 20 δ 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0 1.0 δ 0.8 ϰdeg =1 ◦ δC δD δmin/2 2 4 6 8 ϰdeg =3 ◦ 0.6 δC δD δmin/2 0.4 0.2 0.0 0 10 2 4 orbit number 6 8 10 Figure 6: Deviations for Type II and N = 2. In this case δmin = 1.32. At κdeg = 1◦ stability condition (24) holds with a large margin for the whole simulation interval of 10 orbital periods, while at κdeg = 3◦ it breaks down after six orbital periods. those obtained in the first case (Type I, N = 3). The deviations are in fact now higher than in that case, but this is somewhat compensated by the larger δmin , so that the no-collision condition at κdeg = 3◦ again holds about up to the half of the simulation interval. 5.2 Frequency adjustment by tuning the mass ratio The specific pattern of deviation growth observed in the above examples (linear growth in time with superimposed periodicity) strongly suggests that this growth is largely due to the gradual shift of the trajectories occurring because the ratio of the true x and y system’s frequencies does not exactly match the approximate value α0 = N mD /mC = 8. We can expect to negate these shifts by adjusting the mass ratio: N mD = α0 − ∆α. mC We find the appropriate ∆α numerically, by minimizing the maximum deviation of deputy satellites for 10 orbital periods: max δD (t) → min ∆α t∈[0,10T0 ] We perform these optimizations at different angles κdeg for the second, most stable configuration from the previous section (Type I, N = 5). The results are shown in Table 2. The 21 κdeg 1◦ 2◦ 3◦ 4◦ 5◦ 6◦ ∆αopt 0.021 0.039 0.071 0.116 0.169 0.230 max δD (t) t∈[0,10T0 ] max δD (t) t∈[0,10T0 ] max δD (t) t∈[0,30T0 ] w/o adjustment with adjustment with adjustment 0.0408 0.0766 0.1230 0.1780 0.2600 0.3700 0.0210 0.0368 0.0541 0.0702 0.0874 0.1110 0.0375 0.0530 0.0814 0.1280 0.1650 0.2020 Table 2: Results of the mass ratio adjustment experiments. At each angle κdeg we numerically find the optimal adjustment ∆αopt . The maximum deviations of the main satellite from the vertical without adjustment or with the optimized adjustment are shown in the third and fourth columns, respectively. The last column shows results of simulations with adjustment spanning 30 orbital periods. results clearly show a big improvement over the earlier results obtained without mass ratio adjustment. In Fig. 7 the deviations of the adjusted systems with different values of κdeg are plotted for 30 orbital periods. 6 Conclusions We found that subject to appropriate choice of “hub-and-spoke” system parameters deputy satellites can move along Lissajous curves so that the system stays in free motion (i.e. no fuel is consumed in the nominal operation mode). The main satellite is in a state of relative equilibrium on the local vertical, passing through the system CoM. Our analysis shows the existence of rather nontrivial relations between the system’s parameters that, when satisfied, produce a well-balanced system without collisions between the deputy satellites or tethers. Certain configurations determined by the proper choice of parameters allow placing an additional satellite at the center of the “hub-and-spoke” system without collisions with the other parts of the system. Of course, the intricate way in which the deputy satellites change their positions is challenging for technical implementation. In particular, the tethers must be attached to the main satellite so as not to intertwine. Analysis of the tether entanglement in terms of homotopy provides another set of constraints that the system parameters should satisfy. One of the most curious results in this paper is the study of nonlinear effects, which allows to formulate yet another set of conditions for the systems parameters to cancel out the second order corrections in the systems’ equations of motion. Our numerical experiments with two parameter sets (for 5 deputy satellites) satisfying these conditions corroborate the theoretically predicted system’s stability The application of the proposed motion pattern is limited to small angular deviations of 22 0.20 δ 0.15 ϰdeg =1 ◦ δC δD δmin/2 0.10 0.05 0.00 0 0.20 δ 0.15 5 10 15 20 25 ϰdeg =3 ◦ δC δD δmin/2 0.10 0.05 0.00 0 0.20 δ 0.15 5 10 15 20 25 30 ϰdeg =6 ◦ δC δD δmin/2 0.10 0.05 0.00 0 30 5 10 15 orbit number 20 25 30 Figure 7: Deviations in the Type I, N = 5 formation with the optimally adjusted satellite mass ratio. Stability condition (24) holds for the whole simulation interval of 30 orbital periods for all three considered values of κdeg . the tethers from the local vertical. The possibility to extend it to large deviations from the vertical could be a subject of further investigation. Acknowledgment The authors first conceived the idea of the motion described in this paper during the dynamical analysis of the rotating multi-tethered satellite system [1], and the authors sincerely thank their collaborators Didier Alary, Kirill Andreev, Pavel Boyko, Elena Ivanova, and Cyrille Tourneur for the warm and stimulating atmosphere of that study. The work of one of the authors (DY) on the present paper was supported by Russian Science Foundation (project 14-50-00150). 23 Appendix A Proof of Proposition 1 Statements a) and d) are very simple, so we only provide proofs for b) and c). b) Suppose that the satellite i collides with the satellite j at the time moment τ . The condition xi (τ ) = xj (τ ) admits two series of solutions: 

1. 2πp τ + Ni + ϕx = 2π p τ + Nj + nx + ϕx with some integer nx , that is i−j nx = . N p 2. 2πp τ + i N
+ ϕx = 2π
 + nx − p τ + Nj − ϕx with some integer nx , that is     i+j ϕx 1 1 1 − . + nx − τ= 2 p 2 π N 1 2 The condition yi (τ ) = yj (τ ) has similar series but with q instead of p, ϕy instead of ϕx and ny instead of nx . Thus there are four possibilities for the satellites to collide: an element in series 1 or 2 for the coordinate x must occur simultaneously with an element in series 1 or 2 for the coordinate y. Let us deal with these cases one by one. 1. Series 1 for x and series 1 for y – impossible because p and q are co-prime and −N < i − j < N. 2. Series 2 for x and series 1 for y. Series 1 for y has non-trivial solutions iff N and q are not co-prime. If a non-trivial solution exists, the time moment τ is found from series 2 for x. Thus in this series collisions occur iff N and q are not co-prime. 3. Series 1 for x and series 2 for y. Similarly, collisions occur iff N and p are not co-prime. 4. Series 2 for x and series 2 for y. Equating τ from both series, we obtain     1 1 1 1 ϕx ϕy = + nx − + ny − p 2 π q 2 π whence nx q−ny p= qϕx − pϕy p − q + . π 2 Since p and q are co-prime, with nx , ny running over all possible integers, the left side also runs over all possible integers, i.e. collisions in this case happen iff the right side is an integer. c) Like in the proof of b), we obtain two series of relations from the condition xi (τ ) =xj (τ ). However, the first series takes the form i−j =nx N 24 and has no non-trivial solutions, because −N <i − j<N . This leaves the second series, which has the form   1 1 ϕx i + j τ= . + nx − − 2p 2 π N Like before, we equate τ from the series for x and y and obtain     1 1 ϕx i + j ϕy i + j 1 1 = + nx − − + ny − − p 2 π N q 2 π N whence N (qnx − pny ) − (q − p) (i + j) =  qϕx − pϕy p − q + π 2  N. As nx , ny run over all integer values, the expression qnx −pny also takes all integer values. If N ≥ 3, the sum i + j for all possible pairs of different numbers from 1 to N assumes all possible integer values modulo N . Thus the expression in the left side takes all possible integer values divisible by the greatest common divisor of N and q − p. On the other hand, if N = 2, then i + j = 3 and the above relation simplifies to the requirement that qϕx − pϕy π be integer. Appendix B Proof of Proposition 2 a) For a closed curve (x (τ ) , y(τ ))τ ∈[0,1] not containing the origin, the number w of its rotations about the origin can be written as w= 1 X sign (x (τ ) ẏ (τ )), 2 τ :y(τ )=0 assuming that the intersections of the curve with the x coordinate axis are non-degenerate. We will apply this formula to xij (τ ) = xj (τ ) − xi (τ ) , yij (τ ) = yj (τ ) − yi (τ ) . We consider separately the two types of formations. Type I. In this case     πp(j − i) j+i xij (τ ) = 2x0 sin + ϕx , cos 2πp τ + N 2N     j+i πq(j − i) + ϕy , cos 2πq τ + yij (τ ) = 2y0 sin N 2N 25 so that wi,j    πq(j − i) πp(j − i) sign sin w∗ = ±w∗ , = sign (x0 y0 ) sign sin N N  where w∗ is the winding number of the Lissajous curve x∗ (τ ) = cos (2πpτ + ϕx ) , y∗ (τ ) = cos (2πqτ + ϕy ) . We have 1 y∗ (τ ) = 0 for τ = 2q  ϕy 1 s− − π 2  , s = 1, ..2q. Moreover, with this choice of τ we have ẏ∗ (τ ) > 0 for even s and ẏ∗ (τ ) < 0 for odd s. Applying the formula for the winding number, we obtain    2q 1X πps s + ϕ∗ (−1) sign cos w∗ = 2 s=1 q with some (unimportant) phase constant ϕ∗ . Now consider separately the cases when both p, q are odd and when one of them is even. Let both p, q be odd. Since p, q are co-prime, the values (πps/q)mod 2π run   over the 2q−1 πps values in the set Z2q = {πr/q}r=0 as s runs over 1, 2, .., 2q. Then the quantity cos q + ϕ∗ takes equally many positive and negative values as s runs over 1, 2, .., 2q. Now, if s runs only over even values 2, .., 2q, then (πps/q)mod 2π runs over the values in the set Zq ={2πr/q}q−1 r=0 . Then, since q is odd, the number of positive and negative values taken by cos as s runs only even values 2, .., 2q, differs by 1. It follows that w∗ = πps q + ϕ∗ , 1 (±1 − (∓1)) = ±1. 2 Now let one of p, q be even; without loss of generality we can assume that it is   q. We can πps then repeat the above argument, but since q is now even, we conclude that cos q + ϕ∗ takes equally many negative and positive values as s runs over the even or odd subset of 1, 2.., 2q. It follows that w∗ = 21 (0 − 0) = 0. We have thus proved statement a) for Type I formation. Type II. In this case     π(j − i) j+i xij (τ ) = 2x0 sin + ϕx , cos 2π pτ + N 2N     j+i π(j − i) cos 2π qτ + + ϕy , yij (τ ) = 2y0 sin N 2N which leads to wi,j      2q sign (x0 y0 ) X ps i + j q − p s (−1) sign cos π = + + ϕ∗ 2 q N q s=1 26 with some phase constant ϕ∗ not depending on i, j. The statement a) then follows just like in case of Type I. b) From the proof of a), wi,j = w∗ sign (x0 y0 ) sign(fj−i ), where we set  πps   πqs  fs = sin sin . N N Note first that if q − p is divisible by 2N , then fs = sin2 (πps/N ) > 0 for all s not divisible by N (recall that N and p are co-prime by the assumption of no collisions). Accordingly, wi,j = w∗ sign (x0 y0 ) for all i, j = 1, .., N with i 6= j. Similarly, wi,j = −w∗ sign (x0 y0 ) if q + p is divisible by 2N . Now suppose that neither of q − p, q + p is divisible by 2N . We need to show that fs takes both positive and negative values as s runs over 1, 2, .., N − 1. Since fN = f0 = 0 and P2N f2N −s = fs , it suffices to show that s=0−1 fs = 0. But that immediately follows from the hypothesis and the identity   πs(p − q) πs(p + q) 1 cos . − cos fs = 2 N N The presence of strong entanglement follows from the presence of different winding numbers, since under a weak homotopy the winding number (possibly nonzero) must be the same for all pairs. References [1] Alary, D., Andreev, K., Boyko, P., Ivanova, E., Pritykin, D., Sidorenko, V., Tourneur, C., Yarotsky, D.: Dynamics of multi-tethered pyramydal satellite formation. Acta Astronautica. 117, 222-232 (2015) [2] Armstrong, M.A.: Basic topology. Springer, New York (1983) [3] Arnold, V.I.: Ordinary differential equations. MIT Press, Cambridge, Massachussets (1978) [4] Avanzini, G., Fedi, M.: Refined dynamical analysis of multi-tethered satellite formations. Acta Astronautica. 84, 36-48 (2013). [5] Avanzini, G., Fedi, M.: Effects of eccentricity of the reference orbit on multi-tethered satellite formations. Acta Astronautica. 94, 338-350 (2014). [6] Bekey, I.: Tethers open new space options. Astronautics and Aeronautics, 21, 23-40 (1983). [7] Cai, Z., Zhao, J., Peng, H., Qi, Z.: Nonlinear control of rotating multi-tethered formations in halo orbits. Int. J. Comput. Methods, 11, 1344008 [26 pages](2014). 27 [8] Celletti, A., Sidorenko, V.V.: Some properties of dumbbell satellite attitude dynamics. Cel. Mech. Dyn. Astron., 101, 105-126 (2008). [9] Clohessy, W.H., Wiltshire, R.S.: Terminal guidance system for satellite rendezvous. J. Aerosp. Sci., 27, 653-658 (1960). [10] Gantmacher, F.: Applications of the theory of matrices, Interscience, New York (1959) [11] Nayfeh, A.H.: Perturbations methods, Wiley-Interscience, New York (1973) [12] Panosian, S., Seubert, C.R., Schaub, H.: Tethered Coulomb Structure Applied to Close Proximity Situational Awareness AIAA Journal of Spacecraft and Rockets, 49, 11831193 (2012). [13] Pizarro-Chong, A., Misra, A.K.: Dynamics of multi-tethered satellite formations containing a parent body. Acta Astronautica, 63, 1188-1202 (2008). [14] Seubert, C.R., Schaub, H.: Tethered Coulomb Structures: Prospects and Challenges, Journal of Astronautical Sciences, 57, 347-348 (2009) [15] Seubert, C.R., Panosian, S., Schaub, H.: Attitude and Power Analysis of MultiTethered, Two-Node Tethered Coulomb Structures AIAA Journal of Spacecraft and Rockets, 48, 1033-1045 (2011) [16] Sidorenko, V.V., Celletti, A.: A “spring-mass” model of tethered satellite systems: properties of planar periodic motions. Cel. Mech. Dyn. Astron., 107, 209-231 (2010). [17] Wong, B., Misra, A.: Planar dynamics of variable length multi-tethered spacecraft near collinear Lagrangian points. Acta Astronautica, 63, 1178-1187 (2008). [18] Zhao, J., Cai, Z.: Nonlinear dynamics and simulation of multi-tethered satellite formations in Halo orbits. Acta Astronautica, 63, 673-681 (2008). 28