1 Introduction

The general gravitational three-body problem has well known sets of equilibrium solutions, classified through the Euler (collinear) and Lagrange (triangular) solutions, representing periodic orbits in a central configuration (Musielak and Quarles 2014). A specific straight line solution also arises as a special case of the Euler problem where the masses are collinear and fixed in a rotating frame of reference (Battin 1999). In the limit of the restricted three-body problem, five well-known equilibrium solutions are found where an infinitesimal mass is fixed in a frame of reference rotating with two primary masses. The addition of external forces has been considered in the context of the restricted three-body problem. For example, adding continuous low thrust, either continuous reaction propulsion or a solar sail, or adding discrete impulses to the infinitesimal mass will displace the location of the equilibria (Dusek 1966; McInnes et al. 1994; Morimoto et al. 2007; McInnes 2011).

Moreover, the acceleration of the centre-of-mass of the gravitational two-body problem due to the application of an external force has been noted by a number of authors. Shkadov (1987) pointed out that the radiation pressure exerted on a reflector in static equilibrium relative to the Sun can be used to accelerate the centre-of-mass of the Sun-reflector system. McInnes (1991) demonstrated from first principles that the centre-of-mass of a Sun-solar sail system will accelerate since the two masses are gravitationally coupled. Similarly, McInnes (2002) demonstrated that a reflector in static equilibrium relative to the Earth will accelerate the Earth-reflector system due to the same underlying dynamics. Lu and Love (2005) proposed the ‘gravity tractor’ to use reaction propulsion to accelerate an asteroid-spacecraft system, while McInnes (2007) noted that displaced orbits can be in some cases be used more effectively for such applications. Fahnestock and Scheeres (2008) investigated the stability and control of gravity-tractors, including both position and attitude control. Orbit manipulation of a binary asteroid pair using an external force applied to one of the asteroids was investigated by Liu (2019), including the acceleration of the centre-of-mass and control of the relative motion of the two asteroids.

Other related work includes the use of Coulomb interactions to generate regular geometric configurations of masses. For example, Coulomb interactions can be used to form a triangular configuration of three charged spacecraft in Earth orbit (Hussein and Schaub 2006; Felicetti and Palerini 2015). Similarly, the restricted three-body problem can be formulated with Coulomb interactions between two of the masses, and so the location of the equilibrium points can be manipulated via the product of their charge to mass ratio (Yamakawa and Bando 2010). An ‘electrostatic tractor’ can also be envisaged which augments the gravitational interaction between two masses with a Coulomb interaction (Murdock et al. 2008; Bengston et al. 2018).

An accelerated two-body problem has also been considered by Namouni and Guzzo (2007) to model problems such as the influence of asymmetric astrophysical jets on the dynamics of planetary discs. Here, the momentum loss of the parent star through an asymmetric jet from an accretion disc acts as a perturbing acceleration to the two-body problem. It can be shown that the centre-of-mass of the two-body system accelerates due to such momentum loss. The dynamics of jet-induced excitation has been used to investigate the origin of the eccentricity of exoplanet systems (Namouni 2005). Hovering orbits, see also McInnes and Simmons (1989), Dankowicz (1994) and McInnes (1997), were also found where the geometric centre of the orbit is displaced away from the parent star due to the acceleration of the star itself via momentum loss.

In this paper, three gravitationally interacting masses are considered with an external force applied to one of the masses. As will be seen, the centre-of mass of the system accelerates and it is found that a straight line (but accelerating) equilibrium solution exists. By analogy, the general gravitational three-body problem (Battin 1999) possesses straight line equilibrium solutions, where the masses are arranged in a collinear fashion fixed in a rotating frame of reference. In this configuration, the centripetal acceleration induced by the rotating frame counterbalances the mutual gravitational interaction between the masses. However, it will be shown that an accelerating straight line equilibrium solution can also exist if there is a counterbalancing linear acceleration generated by an external force applied to one of the masses. It will demonstrated that only a single, unstable equilibrium solution exists.

The paper is organised as follows. Section 2 derives the relative equations of motion for the problem and notes that the standard integrals of motions of the three-body problem are in general not available. Then, in Sect. 3 it is demonstrated that there exists a single, unique equilibrium solution to the accelerated three-body problem, whereby the three masses form a straight line in the presence of an external force. Section 4 demonstrates that this single equilibrium solution is always unstable, while Sect. 5 determines the conditions for the instability to be fully controllable. Finally, Sect. 6 provides an application of the problem to manoeuvring chains of small, gravitationally coupled asteroids for space resource utilisation.

2 Accelerated three-body problem

The dynamics of three gravitationally interacting masses can be described in a fixed inertial frame of reference I, as shown in Fig. 1. The masses \({m}_{i}\) (i = 1, 2, 3) are located at position \({{\varvec{r}}}_{i}\) (i = 1, 2, 3) in the inertial frame and have relative position \({{\varvec{r}}}_{ij}\)=\({{\varvec{r}}}_{j}-{{\varvec{r}}}_{i}\) (i, j = 1, 2, 3), scalar separation \({r}_{ij}=\left|{{\varvec{r}}}_{j}-{{\varvec{r}}}_{i}\right|\) (i, j = 1, 2, 3) with unit direction vectors \({\widehat{{\varvec{r}}}}_{ij}={{\varvec{r}}}_{ij}/{r}_{ij}\) (i, j = 1, 2, 3). The centre-of-mass C of the system of three masses is located at position \({{\varvec{r}}}_{c}\). The masses are acted on by their mutual gravitational interaction, while mass 3 is also acted on by an external force \({\varvec{f}}\). The equations of motion for each mass can then be written as:

$${m}_{1}{\ddot{{\varvec{r}}}}_{1}\left(t\right)=\frac{G{m}_{1}{m}_{2}}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}+\frac{G{m}_{1}{m}_{3}}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}$$
(1a)
$${m}_{2}{\ddot{{\varvec{r}}}}_{2}\left(t\right)=\frac{G{m}_{2}{m}_{1}}{{{r}_{21}}^{2}}{\widehat{{\varvec{r}}}}_{21}+\frac{G{m}_{2}{m}_{3}}{{{r}_{23}}^{2}}{\widehat{{\varvec{r}}}}_{23}$$
(1b)
$${m}_{3}{\ddot{{\varvec{r}}}}_{3}\left(t\right)=\frac{G{m}_{3}{m}_{2}}{{{r}_{32}}^{2}}{\widehat{{\varvec{r}}}}_{32}+\frac{G{m}_{3}{m}_{1}}{{{r}_{31}}^{2}}{\widehat{{\varvec{r}}}}_{31}+{\varvec{f}}$$
(1c)

where \(G\) is the universal gravitational constant. Since there is an external force \({\varvec{f}}\ne 0\) acting on the system of isolated masses, in general neither the total energy, linear momentum or angular momentum of the system will be conserved. Indeed, adapting Musielak and Quarles (2014), it can be shown that:

Fig. 1
figure 1

Accelerated three-body problem in an inertial frame \(I\) with masses \(\left({m}_{1}, {m}_{2}{,m}_{3}\right)\) at positions \(\left({{\varvec{r}}}_{1}, {{\varvec{r}}}_{2}{, {\varvec{r}}}_{3}\right)\), centre-of-mass \(C\) at position \({{\varvec{r}}}_{C}\) and external force \({\varvec{f}}\)

$$\frac{dE}{dt}={\dot{{\varvec{r}}}}_{3}\cdot {\varvec{f}}$$
(2)

where the total energy of the three-body system \(E=T+U\) is defined by the kinetic energy \(T=1/2{\sum }_{i=1}^{3}{m}_{i}{\dot{{\varvec{r}}}}_{i}.{\dot{{\varvec{r}}}}_{i}\) and potential energy \(U=G/2{\sum }_{i=1}^{3}{\sum }_{j=1}^{3}{m}_{i}{m}_{j}/{r}_{ij}\) for \(i\ne j\). Similarly, it can be shown that:

$$\frac{d{\varvec{H}}}{dt}={{\varvec{r}}}_{3}\times {\varvec{f}}$$
(3)

where the total angular momentum \({\varvec{H}}\) of the three-body system is defined by \({\varvec{H}}={\sum }_{i=1}^{3}{m}_{i}{{\varvec{r}}}_{i}\times {\dot{{\varvec{r}}}}_{i}\).

A further property of key importance is the centre-of-mass \({{\varvec{r}}}_{c}\) of the three-body system defined as:

$${{\varvec{r}}}_{c}=\frac{{\sum }_{i=1}^{3}{m}_{i}{{\varvec{r}}}_{i}}{M}$$
(4)

where \(M={\sum }_{i=1}^{3}{m}_{i}\) is the total mass. Then, differentiating \({{\varvec{r}}}_{c}\) and substituting from Eqs. (1), while noting that that \({{\varvec{r}}}_{ji}={-{\varvec{r}}}_{ij}\), it can be shown that:

$${\ddot{{\varvec{r}}}}_{c}=\frac{1}{M}\boldsymbol{ }{\varvec{f}}$$
(5)

It can be seen that the centre-of-mass of the three-body system accelerates; hence, if the masses remain bound, all three masses will experience an acceleration relative to the fixed inertial frame of reference I. The external force \({\varvec{f}}\) is coupled to each mass via their mutual gravitational interaction.

Since in general neither the total energy, linear momentum or angular momentum of the system will be conserved, a reduction in order of Eq. (1) is not available using the standard methods of the three-body problem. However, Eq. (1) can be reduced using relative coordinates \({{\varvec{r}}}_{12}\) and \({{\varvec{r}}}_{13}\) since the identity \({{\varvec{r}}}_{23}={{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\) exists (Cruz 2020; Broucke and Lass 1973). Therefore, again noting that \({{\varvec{r}}}_{ji}={-{\varvec{r}}}_{ij}\), from Eq. (1) it can be shown by subtraction that:

$${\ddot{{\varvec{r}}}}_{12}\left(t\right)=-\frac{G\left({m}_{1}+{m}_{2}\right)}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}+\frac{G{m}_{3}}{{\left|{{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right|}^{3}}\left({{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right)-\frac{G{m}_{3}}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}$$
(6a)
$${\ddot{{\varvec{r}}}}_{13}\left(t\right)=-\frac{G\left({m}_{1}+{m}_{3}\right)}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}+\frac{G{m}_{2}}{{\left|{{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right|}^{3}}\left({{\varvec{r}}}_{12}-{{\varvec{r}}}_{13}\right)-\frac{G{m}_{2}}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}+\frac{1}{{m}_{3}}{\varvec{f}}$$
(6b)

where \({r}_{23}=\left|{{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right|\). The initial set of 3 vector equations of motion in the inertial frame of reference I have therefore been reduced to relative equations of motion in coordinates \({{\varvec{r}}}_{12}\) and \({{\varvec{r}}}_{13}\). These reduced equations of motion can now be used to determine the conditions for an equilibrium configuration of the three masses.

3 Conditions for equilibrium

Now that the relative equations of motion have been defined, the conditions for equilibrium solutions can be investigated. In order to proceed, an equilibrium solution to Eqs. (6) can be sought with \({\ddot{{\varvec{r}}}}_{12}={\ddot{{\varvec{r}}}}_{13}=0\). Intuitively, it can be expected that this condition can be satisfied if \({\widehat{{\varvec{r}}}}_{12}={\widehat{{\varvec{r}}}}_{13}=\widehat{{\varvec{f}}}\), where \(\widehat{{\varvec{f}}}={\varvec{f}}/\left|{\varvec{f}}\right|\), so that the three masses are collinear with the external force \({\varvec{f}}\) directed along the unit vectors connecting the masses. To demonstrate that this is the only possible equilibrium configuration, Eqs. (6) can in principle be written as:

$${\ddot{{\varvec{r}}}}_{12}\left(t\right)={f}_{12}{\widehat{{\varvec{r}}}}_{12}+{f}_{13}{\widehat{{\varvec{r}}}}_{13}$$
(7a)
$${\ddot{{\varvec{r}}}}_{13}\left(t\right)={g}_{12}{\widehat{{\varvec{r}}}}_{12}+{g}_{13}{\widehat{{\varvec{r}}}}_{13}+\nu \widehat{{\varvec{f}}}$$
(7b)

for scalar functions \({f}_{12}\), \({f}_{13}\), \({g}_{12}\) and \({g}_{13}\) defined from Eq. (6) and \(\nu =\left|{\varvec{f}}\right|/{m}_{3}\). Equations (7) represent simultaneous linear equations which can be solved to verify the conditions for equilibrium. First, the equilibrium condition \({\ddot{{\varvec{r}}}}_{12}={\ddot{{\varvec{r}}}}_{13}=0\) is applied. Collecting the coefficients, Eqs. (7) can then be written as:

$${\varvec{R}}\left[\begin{array}{c}{\widehat{{\varvec{r}}}}_{12}\\ {\widehat{{\varvec{r}}}}_{13}\end{array}\right]+\left[\begin{array}{c}0\\ \nu \widehat{{\varvec{f}}}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right],\boldsymbol{ }{\varvec{R}}=\left[\begin{array}{cc}{f}_{12}& {f}_{13}\\ {g}_{12}& {g}_{13}\end{array}\right]$$
(8)

whose solution is given by:

$$\left[\begin{array}{c}{\widehat{{\varvec{r}}}}_{12}\\ {\widehat{{\varvec{r}}}}_{13}\end{array}\right]=-{\Vert {\varvec{R}}\Vert }^{-1}\left[\begin{array}{cc}{g}_{13}& {-f}_{13}\\ {-g}_{12}& {f}_{12}\end{array}\right]\left[\begin{array}{c}0\\ \nu \widehat{{\varvec{f}}}\end{array}\right],\Vert {\varvec{R}}\Vert \ne 0$$
(9)

It can be seen from Eq. (9) that \({\widehat{{\varvec{r}}}}_{12}={{\Vert {\varvec{R}}\Vert }^{-1}f}_{13}\nu \widehat{{\varvec{f}}}\) and \({\widehat{{\varvec{r}}}}_{13}=-{\Vert {\varvec{R}}\Vert }^{-1}{f}_{12}\nu \widehat{{\varvec{f}}}\) so that \({\widehat{{\varvec{r}}}}_{12}={\widehat{{\varvec{r}}}}_{13}=\widehat{{\varvec{f}}}\), demonstrating that the three masses must be collinear with the external force \({\varvec{f}}\), as shown in Fig. 2. As noted in Sect. 1, straight line equilibria exist in the general three-body problem where the three masses are collinear in a frame of reference rotating about their common centre-of-mass. Again, the centripetal acceleration induced by the rotating frame counterbalances the mutual gravitational attraction between the masses. The accelerated equilibrium configuration investigated here requires a counterbalancing linear acceleration generated by the external force \({\varvec{f}}\) acting on \({m}_{3}\). While that mass will accelerate, the other two masses \({m}_{1}\) and \({m}_{2}\) will also accelerate due to the mutual gravitational interaction between \({m}_{1}\), \({m}_{2}\) and \({m}_{3}\).

Fig. 2
figure 2

Collinear equilibrium configuration of the three masses \(\left({m}_{1}, {m}_{2}{,m}_{3}\right)\) with \({\widehat{{\varvec{r}}}}_{12}={\widehat{{\varvec{r}}}}_{13}=\widehat{{\varvec{f}}}\)

Given that it has been demonstrated that only a collinear equilibrium configuration is possible, without loss of generality, the 3 masses can now considered to be translating along the \({{\varvec{e}}}_{1}\) axis in the inertial frame \(I\), as shown in Fig. 3. Again, the centre-of mass of the system will accelerate such that \({\ddot{{\varvec{r}}}}_{c}=\left(1/M\right){\varvec{f}}\). It can be seen from Eq. (2) that \(\mathrm{d}E/\mathrm{d}t\ne 0\), since the external force \({\varvec{f}}\) still does work on the 3 gravitationally coupled masses, with \(E\) scaling as \({\Vert {\varvec{f}}\Vert }^{2}{t}^{2}/2M\). However, from Eq. (3) it can be seen that \(\mathrm{d}{\varvec{H}}/\mathrm{d}t=0\), so that the angular momentum of 3 masses is now fixed, with \({\varvec{H}}=0\).

Fig. 3
figure 3

Collinear equilibrium configuration of the three masses \(\left({m}_{1}, {m}_{2}{,m}_{3}\right)\) with spacing parameter \(\lambda \), total separation \(d\) and centre-of-mass acceleration \({\ddot{{\varvec{r}}}}_{c}=\boldsymbol{ }{\varvec{f}}/M\)

In order to proceed, the spacing between\({m}_{1}\), \({m}_{2}\) and \({m}_{3}\) required for equilibrium will now be determined, again with the masses translating along the \({{\varvec{e}}}_{1}\) axis. If the length-scale of the problem is defined as \(d\) such that \({r}_{13}=d\), with \({r}_{12}=\lambda d\) for some spacing parameter\(\lambda \), such that \(0<\lambda \le 1\), then\({r}_{23}=\left(1-\lambda \right)d\), as shown in Fig. 3. Then, it can be shown from Eqs. (6) that the conditions for equilibrium can be written as:

$$-\frac{G\left({m}_{1}+{m}_{2}\right)}{{\lambda }^{2}{d}^{2}}+\frac{G{m}_{3}}{{\left(1-\lambda \right)}^{2}{d}^{2}}-\frac{G{m}_{3}}{{d}^{2}}=0$$
(10a)
$$-\frac{G\left({m}_{1}+{m}_{3}\right)}{{d}^{2}}-\frac{G{m}_{2}}{{\left(1-\lambda \right)}^{2}{d}^{2}}-\frac{G{m}_{2}}{{\lambda }^{2}{d}^{2}}+\frac{1}{{m}_{3}}f=0$$
(10b)

These conditions can be simplified by defining mass ratios \({m}_{i}={\mu }_{i}M\) (i = 1 − 3), where \({\sum }_{i=1}^{3}{\mu }_{i}=1\), while a gravitational scaling parameter \(\Lambda =G{M}^{2}/{d}^{2}\) can also be defined. In non-dimensional form, Eqs. (10) then become:

$$1-\frac{1}{{\left(1-\lambda \right)}^{2}}+\frac{\kappa }{{\lambda }^{2}}=0$$
(11a)
$$\gamma ={\mu }_{3}\left(1-{\mu }_{2}\left(1-\frac{1}{{\lambda }^{2}}-\frac{1}{{\left(1-\lambda \right)}^{2}}\right)\right)$$
(11b)

where \(\kappa =\left({\mu }_{1}+{\mu }_{2}\right)/{\mu }_{3}\) and \(\gamma =f/\Lambda \) is the non-dimensional external force required to generate a collinear equilibrium solution with spacing parameter \(\lambda \). The spacing \(\lambda \) is therefore a function of the ratio of \({m}_{3}\), which is acted on by the external force, and the sum of the remaining masses \({m}_{1}+{m}_{2}\). It can be shown that Eq. (11a) can be rewritten as a quartic such that:

$${\lambda }^{4}-2{\lambda }^{3}+{\kappa \lambda }^{2}-2\kappa \lambda +\kappa =0$$
(12)

Since \(\kappa \ge 0\), from Descartes’ rule of signs it can be seen that Eq. (11) possess at most 4 positive roots (or 2 or zero), no negative roots and so at most 4 imaginary roots (or 2 or zero). In general, it is found that there are 2 positive real roots and 2 imaginary roots, with a single real root in the range \(0\le \lambda \le 1\), corresponding to a single equilibrium configuration of the masses. The existence of a single real root in the range \(0\le \lambda \le 1\) can be now demonstrated using Sturm’s theorem.

The quartic represented by Eq. (12) can be written as \({f}_{0}\left(\lambda ;\kappa \right)={\lambda }^{4}-2{\lambda }^{3}+{\kappa \lambda }^{2}-2\kappa \lambda +\kappa \). Then, a Sturm chain of polynomials (Thomas 1941) can be defined firstly by calculating the derivatives of the quartic such that \({f}_{i}\left(\lambda ;\kappa \right)={{d}^{i}f}_{i}\left(\lambda ;\kappa \right)/d{\lambda }^{i}\) (i = 1 − 4), and so with first derivative \({f}_{1}\left(\lambda ;\kappa \right)={4\lambda }^{3}-6{\lambda }^{2}+2\kappa \lambda -2\kappa \). Then, a sequence of polynomials can be defined using the relationship \({f}_{i+1}\left(\lambda ;\kappa \right)=-{f}_{i-1}\left(\lambda ;\kappa \right)\boldsymbol{ }\mathrm{rem }{f}_{i}\left(\lambda ;\kappa \right)\), where \(\mathrm{rem}\) indicates the remainder of the polynomial division of \({f}_{i-1}\left(\lambda ;\kappa \right)\) and \({f}_{i}\left(\lambda ;\kappa \right)\). The full Strum chain is found to be:

$${f}_{0}\left(\lambda ;\kappa \right)={\lambda }^{4}-2{\lambda }^{3}+{\kappa \lambda }^{2}-2\kappa \lambda +\kappa $$
(13a)
$${f}_{1}\left(\lambda ;\kappa \right)={4\lambda }^{3}-6{\lambda }^{2}+2\kappa \lambda -2\kappa $$
(13b)
$${f}_{2}\left(\lambda ;\kappa \right)=\frac{1}{4}\left({\lambda }^{2}\left(3-2\kappa \right)+5\kappa \lambda -3\kappa \right)$$
(13c)
$${f}_{3}\left(\lambda ;\kappa \right)=\frac{8\kappa ({\kappa }^{2}+9-({\kappa }^{2}-\kappa +18)\lambda )}{{(3-2\kappa )}^{2}}$$
(13d)
$${f}_{4}\left(\lambda ;\kappa \right)=-\frac{{\left(2\kappa -3\right)}^{2}({\kappa }^{2}+27)}{4{({\kappa }^{2}-\kappa +18)}^{2}}$$
(13e)

The number of real roots in the interval \(0\le \lambda \le 1\) can then be determined from the difference in the number of sign changes in the Sturm sequence at the limits of the internal. First, the Sturm chain can be evaluated at \(\lambda =0\) such that:

$${f}_{0}\left(0;\kappa \right)=\kappa $$
(14a)
$${f}_{1}\left(0;\kappa \right)=-2\kappa $$
(14b)
$${f}_{2}\left(0;\kappa \right)=-\frac{3}{4}\kappa $$
(14c)
$${f}_{3}\left(0;\kappa \right)=\frac{8\kappa ({\kappa }^{2}+9)}{{(3-2\kappa )}^{2}}$$
(14d)
$${f}_{4}\left(0;\kappa \right)=-\frac{{\left(3-2\kappa \right)}^{2}({\kappa }^{2}+27)}{4{({\kappa }^{2}-\kappa +18)}^{2}}$$
(14e)

It can be seen from Eqs. (14) that since \(\kappa >0\), there are 3 sign changes along the chain with the sequence \(\left[+, -,-,+,-\right]\). Similarly, the Sturm chain can be evaluated at \(\lambda =1\) such that:

$${f}_{0}\left(1;\kappa \right)=-1$$
(15a)
$${f}_{1}\left(1;\kappa \right)=-2$$
(15b)
$${f}_{2}\left(1;\kappa \right)=\frac{3}{4}$$
(15c)
$${f}_{3}\left(1;\kappa \right)=\frac{8\kappa (\kappa -9)}{{(3-2\kappa )}^{2}}$$
(15d)
$${f}_{4}\left(1;\kappa \right)=-\frac{{\left(3-2\kappa \right)}^{2}({\kappa }^{2}+27)}{4{({\kappa }^{2}-\kappa +18)}^{2}}$$
(15e)

Again, that since \(\kappa >0\), it can be seen from Eqs. (15) that there are 2 sign changes along the chain with the sequence \(\left[-, -,+,\pm ,-\right]\). This is independent of the sign of \({f}_{3}\left(1;\kappa \right)\), whose sign changes at \(\kappa =9\), since the total number of sign changes is unchanged. The number of real roots in the interval \(0\le \lambda \le 1\) is then given by the difference in the number of sign changes at the limits of the interval, which is one. Therefore, there is a single, unique equilibrium configuration of the 3 masses. In principle, a chain of N masses with an external force applied to the final mass in the chain can also form an equilibria configuration, although the problem is not considered here.

Now that only a single equilibrium configuration has been found, the limiting cases for \(\kappa \) can be investigated. First, it can be seen that in the limit \(\kappa \to 0\) that Eq. (12) reduces to \({\lambda }^{4}-2{\lambda }^{3}=0\) and so \(\lambda \to 0\) since \(0\le \lambda \le 1\). This limiting case corresponds to \({m}_{1}+{m}_{2}\ll {m}_{3}\) with \({r}_{12}\to 0\) and \({r}_{23}\to d\). Here, \({m}_{1}\) and \({m}_{2}\) become closely spaced and so are strongly coupled to ensure equilibrium with their smaller relative masses. Similarly in the limit \(\kappa \to \infty \) Eq. (12) reduces to \({\lambda }^{2}-2\lambda +1=0\) and so \(\lambda \to 1\). This limiting case corresponds to \({m}_{1}+{m}_{2}\gg {m}_{3}\) with \({r}_{12}\to d\) and \({r}_{23}\to 0\). Here, \({m}_{1}\) and \({m}_{2}\) become separated with \({m}_{2}\) and \({m}_{3}\) closely spaced, and so strongly coupled, to ensure equilibrium with the relatively smaller mass \({m}_{3}\). Furthermore, for the equal mass problem \(\kappa =2\), with \({\mu }_{1}={\mu }_{2}={\mu }_{3}=1/3\), although since \(\kappa =\left({\mu }_{1}+{\mu }_{2}\right)/{\mu }_{3}\) in principle \(\kappa =2\) only requires \({\mu }_{1}+{\mu }_{2}=2/3\). Solving Eq. (12) with \(\kappa =2\) yields a single real equilibrium solution in the range \(0\le \lambda \le 1\) given by \(\lambda \cong 0.606\). The equal mass equilibrium configuration therefore corresponds to the central mass \({m}_{2}\) being located closer to \({m}_{3}\). Other cases are considered later in Sect. 6.

For a given parameter \(\kappa \), a pair of parameters representing a collinear equilibrium solution to the accelerated three-body system can then be found, denoted by \(\left(\lambda ,\gamma \right)\), representing the equilibrium spacing \(\lambda \) and the required non-dimensional force \(\gamma \) to maintain the configuration. The equilibrium spacing parameter \(\lambda \) is shown in Fig. 4 as a function of \(\kappa \). The asymptotic limiting cases noted above with \(\kappa \to 0\) and \(\kappa \to \infty \) can be seen, along with the equal mass solution at \(\kappa =2\) with \(\lambda \cong \) 0.606 from Eq. (12) and \(\gamma \cong \) 1.241 from Eq. (11b). Example configurations in dimensional units again be considered later in Sect. 6.

Fig. 4
figure 4

Equilibrium separation \(\lambda \) as a function of \(\kappa \), with a separation of \(\lambda \cong 0.606\) for \(\kappa =2\) (\({\mu }_{3}=1/3\)) and non-dimensional force \(\gamma =1.241\)

4 Stability properties

In order to determine the stability properties of the single equilibrium configuration found in Sect. 3, the relative equations of motion defined by Eqs. (6) can be linearised. For ease of illustration, and without loss of generality, it will be assumed that the three masses are again translating along the inertial e1 axis, shown in Fig. 3. Then, an equilibrium solution \({\widetilde{\mathbf{r}}}_{12}=\left(\mathrm{\lambda d},\mathrm{0,0}\right)\) and \({\widetilde{\mathbf{r}}}_{13}=\left(\mathrm{d},\mathrm{0,0}\right)\) can be defined and perturbations \(\delta {{\varvec{r}}}_{12}\) and \(\delta {{\varvec{r}}}_{13}\) added such that \({{\varvec{r}}}_{12}\to {\widetilde{{\varvec{r}}}}_{12}+\delta {{\varvec{r}}}_{12}\) and \({{\varvec{r}}}_{13}\to {\widetilde{{\varvec{r}}}}_{13}+\delta {{\varvec{r}}}_{13}\). From Eqs. (6), the resulting coupled linear equations can be written as:

$$\delta {\ddot{{\varvec{r}}}}_{12}\left(t\right)={{\varvec{\Gamma}}}_{11}\delta {{\varvec{r}}}_{12}+{{\varvec{\Gamma}}}_{12}\delta {{\varvec{r}}}_{13}$$
(16a)
$$\delta {\ddot{{\varvec{r}}}}_{13}\left(t\right)={{\varvec{\Gamma}}}_{21}\delta {{\varvec{r}}}_{12}+{{\varvec{\Gamma}}}_{22}\delta {{\varvec{r}}}_{13}$$
(16b)

where, for a constant external force \({\varvec{f}}\), it can seen that:

$${{\varvec{\Gamma}}}_{11}=\frac{\partial }{\partial {{\varvec{r}}}_{12}}{\left[-\frac{G\left({m}_{1}+{m}_{2}\right)}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}+\frac{G{m}_{3}}{{\left|{{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right|}^{3}}\left({{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right)-\frac{G{m}_{3}}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}\right]}_{{{\varvec{r}}}_{12}={\widetilde{{\varvec{r}}}}_{12},{\boldsymbol{ }{\varvec{r}}}_{13}={\widetilde{{\varvec{r}}}}_{13}}$$
(17a)
$${{\varvec{\Gamma}}}_{12}=\frac{\partial }{\partial {{\varvec{r}}}_{13}}{\left[-\frac{G\left({m}_{1}+{m}_{2}\right)}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}+\frac{G{m}_{3}}{{\left|{{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right|}^{3}}\left({{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right)-\frac{G{m}_{3}}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}\right]}_{{{\varvec{r}}}_{12}={\widetilde{{\varvec{r}}}}_{12},{\boldsymbol{ }{\varvec{r}}}_{13}={\widetilde{{\varvec{r}}}}_{13}}$$
(17b)
$${{\varvec{\Gamma}}}_{21}=\frac{\partial }{\partial {{\varvec{r}}}_{12}}{\left[-\frac{G\left({m}_{1}+{m}_{3}\right)}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}+\frac{G{m}_{2}}{{\left|{{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right|}^{3}}\left({{\varvec{r}}}_{12}-{{\varvec{r}}}_{13}\right)-\frac{G{m}_{2}}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}\right]}_{{{\varvec{r}}}_{12}={\widetilde{{\varvec{r}}}}_{12},{\boldsymbol{ }{\varvec{r}}}_{13}={\widetilde{{\varvec{r}}}}_{13}}$$
(17c)
$${{\varvec{\Gamma}}}_{22}=\frac{\partial }{\partial {{\varvec{r}}}_{13}}{\left[-\frac{G\left({m}_{1}+{m}_{3}\right)}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}+\frac{G{m}_{2}}{{\left|{{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right|}^{3}}\left({{\varvec{r}}}_{12}-{{\varvec{r}}}_{13}\right)-\frac{G{m}_{2}}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}\right]}_{{{\varvec{r}}}_{12}={\widetilde{{\varvec{r}}}}_{12},{\boldsymbol{ }{\varvec{r}}}_{13}={\widetilde{{\varvec{r}}}}_{13}}$$
(17d)

Evaluating the derivatives, substituting for the equilibrium configuration \({\widetilde{\mathbf{r}}}_{12}=\left(\mathrm{\lambda d},\mathrm{0,0}\right)\) and \({\widetilde{\mathbf{r}}}_{13}=\left(\mathrm{d},\mathrm{0,0}\right)\), and defining \(\widetilde{\omega }=\sqrt{GM/{d}^{3}}\), it can then be shown that the matrices are of the form:

$${{\varvec{\Gamma}}}_{11}=\left[\begin{array}{ccc}-2a& 0& 0\\ 0& a& 0\\ 0& 0& a\end{array}\right]{,{\varvec{\Gamma}}}_{12}=\left[\begin{array}{ccc}-2b& 0& 0\\ 0& b& 0\\ 0& 0& b\end{array}\right]$$
(18a)
$${{\varvec{\Gamma}}}_{21}=\left[\begin{array}{ccc}-2c& 0& 0\\ 0& c& 0\\ 0& 0& c\end{array}\right]{,{\varvec{\Gamma}}}_{22}=\left[\begin{array}{ccc}-2d& 0& 0\\ 0& d& 0\\ 0& 0& d\end{array}\right]$$
(18b)

where the coefficients of the matrices are given by:

$$a={-\widetilde{\omega }}^{2}\left(\frac{1-{\mu }_{3}}{{\lambda }^{3}}+\frac{{\mu }_{3}}{{\left(1-\lambda \right)}^{3}}\right)$$
(19a)
$$b={-\widetilde{\omega }}^{2}\left(1-\frac{1}{{\left(1-\lambda \right)}^{3}}\right){\mu }_{3}$$
(19b)
$$c={-\widetilde{\omega }}^{2}\left(\frac{1}{{\lambda }^{3}}+\frac{1}{{\left(1-\lambda \right)}^{3}}\right){\mu }_{2}$$
(19c)
$$d={-\widetilde{\omega }}^{2}\left(1-{\mu }_{2}\left(1+\frac{1}{{\left(1-\lambda \right)}^{3}}\right)\right)$$
(19d)

In order to determine the stability properties of the equilibrium configuration, a trial solution can be substituted into Eqs. (16) of the form:

$$\left[\begin{array}{c}\delta {{\varvec{r}}}_{12}\\ \delta {{\varvec{r}}}_{13}\end{array}\right]=\left[\begin{array}{c}\delta {\overline{{\varvec{r}}} }_{12}\\ \delta {\overline{{\varvec{r}}} }_{13}\end{array}\right]{e}^{\omega t}$$
(20)

for some constant vectors \(\delta {\overline{{\varvec{r}}} }_{12}\) and \(\delta {\overline{{\varvec{r}}} }_{13}\) and exponent \(\omega \). This leads to the characteristic polynomial of the problem given, which is determined from:

$$\mathrm{det}\Vert {\varvec{\Gamma}}-{\upomega }^{2}{\varvec{I}}\Vert =0$$
(21)

and where the matrix \({\varvec{\Gamma}}\) is assembled from Eqs. (18) as:

$${\varvec{\Gamma}}=\left[\begin{array}{cc}{{\varvec{\Gamma}}}_{11}& {{\varvec{\Gamma}}}_{12}\\ {{\varvec{\Gamma}}}_{21}& {{\varvec{\Gamma}}}_{22}\end{array}\right]$$
(22)

Then, evaluating Eq. (21) it can be shown that the characteristic polynomial can be written as:

$$\left({\upxi }^{2}+2\left(a+d\right)\upxi +4\left(ad-bc\right)\right){\left({\upxi }^{2}-\left(a+d\right)\upxi +\left(ad-bc\right)\right)}^{2}=0$$
(23)

where \({\upxi =\upomega }^{2}\).

In order to determine stability properties of the single equilibrium configuration of the 3 masses, Eq. (23) can be partitioned into two quadratics such that:

$${\psi }_{1}\left(\upxi \right)={\upxi }^{2}+2\left(a+d\right)\upxi +4\left(ad-bc\right)$$
(24a)
$${{\psi }_{2}\left(\upxi \right)=\upxi }^{2}-\left(a+d\right)\upxi +\left(ad-bc\right)$$
(24b)

First, while the signs of \(a+d\) and \(ad-bc\) are determined by the mass ratios \({\mu }_{2}\) and \({\mu }_{3}\) and equilibrium spacing \(\lambda \), Descartes’ rule of signs can be used to investigate the number of positive roots of \({\psi }_{1}\) and \({\psi }_{2}\) by enumerating all 4 combinations of the signs of the terms \(a+d\) and \(ad-bc\):

\({\varvec{a}}+{\varvec{d}}>0\) and \({\varvec{a}}{\varvec{d}}-{\varvec{b}}{\varvec{c}}>0\): It can be seen that \({\psi }_{1}\left(\upxi \right)\) has no changes of sign and \({\psi }_{2}\left(\upxi \right)\) has two changes of sign, while \({\psi }_{1}\left(-\upxi \right)\) has 2 changes of sign and \({\psi }_{2}\left(-\upxi \right)\) has no changes of sign, so the exact number of real positive roots is undetermined. However, if it assumed that \(a+d>0\) and \(ad-bc>0\), and the two roots of \({f}_{1}\) are defined as \({\xi }_{1}\) and \({\xi }_{2}\), then \({\xi }_{1}{+\xi }_{2}=a+d>0\) and \({\xi }_{1}{\xi }_{2}=ad-bc>0\) which implies that \({\xi }_{1}>0\) and \({\xi }_{2}>0\). Hence, there are at least two real positive roots, with associated eigenvalues \(\pm \sqrt{{\xi }_{1}}\) and \(\pm \sqrt{{\xi }_{2}}\), and so the equilibrium configuration is unstable.

\({\varvec{a}}+{\varvec{d}}<0\) and \({\varvec{a}}{\varvec{d}}-{\varvec{b}}{\varvec{c}}<0\): It can be seen that \({\psi }_{1}\left(\upxi \right)\) has one change of sign and \({\psi }_{2}\left(\upxi \right)\) has one change of sign, while \({\psi }_{1}\left(-\upxi \right)\) has one change of sign and \({\psi }_{2}\left(-\upxi \right)\) has one change of sign. Therefore, both \({\psi }_{1}\left(\upxi \right)\) and \({\psi }_{2}\left(\upxi \right)\) each have a positive real root, and so the equilibrium configuration is unstable.

\({\varvec{a}}+{\varvec{d}}>0\) and \({\varvec{a}}{\varvec{d}}-{\varvec{b}}{\varvec{c}}<0\): It can be seen that \({\psi }_{1}\left(\upxi \right)\) has one change of sign and \({\psi }_{2}\left(\upxi \right)\) has one change of sign, while \({\psi }_{1}\left(-\upxi \right)\) has one change of sign and \({\psi }_{2}\left(-\upxi \right)\) has one change of sign. Therefore, both \({\psi }_{1}\left(\upxi \right)\) and \({\psi }_{2}\left(\upxi \right)\) each have a positive real root, and so the equilibrium configuration is unstable.

\({\varvec{a}}+{\varvec{d}}<0\) and \({\varvec{a}}{\varvec{d}}-{\varvec{b}}{\varvec{c}}>0\): It can be shown that \({\psi }_{1}\left(\upxi \right)\) has two changes of sign, and \({\psi }_{2}\left(\upxi \right)\) has no changes of sign, while \({\psi }_{1}\left(-\upxi \right)\) has no changes of sign and \({\psi }_{2}\left(-\upxi \right)\) has two changes of sign, so the exact number of real positive roots is undetermined. However, if it assumed that \(a+d<0\) and \(ad-bc>0\), and the two roots of \({\psi }_{1}\) are defined as \({\xi }_{1}\) and \({\xi }_{2}\), then \({\xi }_{1}{+\xi }_{2}=-2\left(a+d\right)>0\) and \({\xi }_{1}{\xi }_{2}=4\left(ad-bc\right)>0\) which implies that \({\xi }_{1}>0\) and \({\xi }_{2}>0\). Hence, there are at least two real positive roots, with associated eigenvalues \(\pm \sqrt{{\xi }_{1}}\) and \(\pm \sqrt{{\xi }_{2}}\), and so the equilibrium configuration is unstable.

As an example, consider now the equal mass problem with \({\mu }_{1}=1/3\), \({\mu }_{2}=1/3\) and \(\kappa =2\), and so from Eq. (9) the equilibrium spacing is again given by \(\lambda \cong 0.606\). The matrix coefficients are then found to be \(a=-8.448\widetilde{\omega }\), \(b=5.120\widetilde{\omega }\), \(c=-6.951\widetilde{\omega }\) and \(d=4.787\widetilde{\omega }\). The eigenvalues are therefore given by \({\upomega }_{\mathrm{1,2}}=\pm 1.016\widetilde{\omega }\), \({\upomega }_{3.4}=\pm 3.064\widetilde{\omega }\), \({\upomega }_{\mathrm{5,6}}=\pm 2.167\mathrm{i}\widetilde{\omega }\) and \({\upomega }_{\mathrm{7,8}}=\pm 1.437\mathrm{i}\widetilde{\omega }\) where the repeated eigenvalues from the second quadratic in Eq. (23) have not been listed. Clearly, the equal mass equilibrium configuration is unstable as expected, since there are two positive real eigenvalues.

In order to further investigate the stability properties of the single equilibrium configuration, the full eigenvalue spectrum is shown in Fig. 5 as a function of \(\kappa \). This includes overlapping repeated roots. For each value of \(\kappa \), Eq. (11a) is solved to determine the equilibrium spacing \(\lambda \) and the coefficients defined by Eqs. (19) are determined, assuming \({\mu }_{1}=1/3\) and \({\mu }_{2}=1/3\) and so \({\mu }_{3}=\left({\mu }_{1}+{\mu }_{1}\right)/\kappa \). Again, as expected from the analysis above, a positive real eigenvalue is always present.

Fig. 5
figure 5

Eigenvalue spectrum of the single equilibrium configuration with \(\widetilde{\omega }=1\) as a function of \(\kappa \) (● \(\kappa =1 (\lambda =0.531)\), ■\(\kappa =2 (\lambda =0.606)\), × \(\kappa =30 (\lambda =0.847)\), ▼\(\kappa =100 (\lambda =0.909)\))

5 Controllability properties

In order to investigate active control of the single, unstable equilibrium configuration of the three masses discussed in Sects. 3 and 4, it will be assumed that a small control force \(\delta {{\varvec{f}}}_{i}\) (i = 1,2,3) is applied to each mass, in addition to the main external force \({\varvec{f}}\) applied to \({m}_{3}\). Then, from Eqs. (1), the equation of motion for each mass can now be written as:

$${m}_{1}{\ddot{{\varvec{r}}}}_{1}\left(t\right)=\frac{G{m}_{1}{m}_{2}}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}+\frac{G{m}_{1}{m}_{3}}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}+\delta {{\varvec{f}}}_{1}$$
(25a)
$${m}_{2}{\ddot{{\varvec{r}}}}_{2}\left(t\right)=\frac{G{m}_{2}{m}_{1}}{{{r}_{21}}^{2}}{\widehat{{\varvec{r}}}}_{21}+\frac{G{m}_{2}{m}_{3}}{{{r}_{23}}^{2}}{\widehat{{\varvec{r}}}}_{23}+\delta {{\varvec{f}}}_{2}$$
(25b)
$${m}_{3}{\ddot{{\varvec{r}}}}_{3}\left(t\right)=\frac{G{m}_{3}{m}_{2}}{{{r}_{32}}^{2}}{\widehat{{\varvec{r}}}}_{32}+\frac{G{m}_{3}{m}_{1}}{{{r}_{31}}^{2}}{\widehat{{\varvec{r}}}}_{31}+{\varvec{f}}+\delta {{\varvec{f}}}_{3}$$
(25c)

By reducing the equations of motion using the same procedure as detailed in Sect. 2, it can be shown using relative coordinates \({{\varvec{r}}}_{12}\) and \({{\varvec{r}}}_{13}\) that:

$${\ddot{{\varvec{r}}}}_{12}\left(t\right)=-\frac{G\left({m}_{1}+{m}_{2}\right)}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}+\frac{G{m}_{3}}{{\left|{{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right|}^{3}}\left({{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right)-\frac{G{m}_{3}}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}+\delta {{\varvec{u}}}_{12}$$
(26a)
$${\ddot{{\varvec{r}}}}_{13}\left(t\right)=-\frac{G\left({m}_{1}+{m}_{3}\right)}{{{r}_{13}}^{2}}{\widehat{{\varvec{r}}}}_{13}+\frac{G{m}_{2}}{{\left|{{\varvec{r}}}_{13}-{{\varvec{r}}}_{12}\right|}^{3}}\left({{\varvec{r}}}_{12}-{{\varvec{r}}}_{13}\right)-\frac{G{m}_{2}}{{{r}_{12}}^{2}}{\widehat{{\varvec{r}}}}_{12}+\frac{1}{{m}_{3}}{\varvec{f}}+\delta {{\varvec{u}}}_{13}$$
(26b)

where the two control terms \(\delta {{\varvec{u}}}_{12}\) and \(\delta {{\varvec{u}}}_{13}\) are defined by:

$$\delta {{\varvec{u}}}_{12}=\frac{1}{{m}_{2}}\delta {{\varvec{f}}}_{2}-\frac{1}{{m}_{1}}\delta {{\varvec{f}}}_{1}$$
(27a)
$$\delta {{\varvec{u}}}_{13}=\frac{1}{{m}_{3}}\delta {{\varvec{f}}}_{3}-\frac{1}{{m}_{1}}\delta {{\varvec{f}}}_{1}$$
(27b)

It can be seen that \(\delta {{\varvec{u}}}_{12}\ne 0\) and \(\delta {{\varvec{u}}}_{13}\ne 0\) if \(\delta {{\varvec{f}}}_{1}=0\), so in principle control accelerations can be applied to \({m}_{2}\) and \({m}_{3}\) only, or the control accelerations can be distributed between all 3 masses.

In order to determine the controllability properties of the single unstable equilibrium configuration, the rank of the controllability matrix can be evaluated. In order to proceed, the linearised dynamics of the problem defined by Eqs. (16) can be written in state space form as:

$$\frac{d}{dt}\left[\begin{array}{c}\begin{array}{c}\delta {{\varvec{r}}}_{12}\\ \delta {{\varvec{r}}}_{13}\\ \delta {\dot{{\varvec{r}}}}_{12}\end{array}\\ \delta {\dot{{\varvec{r}}}}_{13}\end{array}\right]=\left[\begin{array}{cccc}0& 0& {\varvec{I}}& 0\\ 0& 0& 0& {\varvec{I}}\\ {{\varvec{\Gamma}}}_{11}& {{\varvec{\Gamma}}}_{12}& 0& 0\\ {{\varvec{\Gamma}}}_{21}& {{\varvec{\Gamma}}}_{22}& 0& 0\end{array}\right]\left[\begin{array}{c}\begin{array}{c}\delta {{\varvec{r}}}_{12}\\ \delta {{\varvec{r}}}_{13}\\ \delta {\dot{{\varvec{r}}}}_{12}\end{array}\\ \delta {\dot{{\varvec{r}}}}_{13}\end{array}\right]+\left[\begin{array}{cc}\begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}\\ {\mathbf{B}}_{12}& 0\\ 0& {\mathbf{B}}_{13}\end{array}\right]\delta {\varvec{u}}$$
(28)

where \({\varvec{I}}\) is the 3 × 3 identity matrix, \(0\) is the 3 × 3 null matrix and the 3 × 1 input matrices \({\mathbf{B}}_{12}\) and \({\mathbf{B}}_{13}\) define the connection between the control vector \(\delta {\varvec{u}}={\left(\delta {{\varvec{u}}}_{12},\delta {{\varvec{u}}}_{13}\right)}^{T}\) and the dynamics of the problem. The controllability matrix \(\mathbf{C}\) is defined in terms of the system matrix \(\mathbf{M}\) and input matrix \(\mathbf{B}\) where:

$$\mathbf{M}=\left[\begin{array}{cccc}0& 0& {\varvec{I}}& 0\\ 0& 0& 0& {\varvec{I}}\\ {{\varvec{\Gamma}}}_{11}& {{\varvec{\Gamma}}}_{12}& 0& 0\\ {{\varvec{\Gamma}}}_{21}& {{\varvec{\Gamma}}}_{22}& 0& 0\end{array}\right]$$
(29)
$$\mathbf{B}=\left[\begin{array}{cc}\begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}\\ {\mathbf{B}}_{12}& 0\\ 0& {\mathbf{B}}_{13}\end{array}\right]$$
(30)

such that:

$$\mathbf{C}=\left[\begin{array}{ccc}\mathbf{B}& \mathbf{M}\mathbf{B}& {\mathbf{M}}^{2}\mathbf{B}{...\mathbf{M}}^{N-1}\mathbf{B}\end{array}\right]$$
(31)

where N = 12 is the size of the NxN system matrix \(\mathbf{M}\). In order to investigate the controllability of the problem, combinations of control accelerations can be considered and the rank of the controllability matrix \(r\left({\varvec{C}}\right)\) determined. Using the input matrices \({\mathbf{B}}_{12}\) and \({\mathbf{B}}_{13}\), different example combinations of control acceleration can be considered, as shown in Table 1. It can be seen that in general the system is fully controllable if there exists a control acceleration for each axis of motion, either applied to each mass (case a), applied to a single mass only (cases b, c), or distributed between the masses (cases d-g). For distributed control (cases d-g), a control acceleration is required for each axis of motion.

Table 1 Rank of the controllability matrix \({\varvec{C}}\) with a range of input vector configurations \({\mathbf{B}}_{12}\) and \({\mathbf{B}}_{13}\)

In order to illustrate the use of active control for stabilisation of the single unstable equilibrium configuration, the control vector \(\delta {\varvec{u}}\) will now be defined with \(\delta {{\varvec{f}}}_{1}=0\) such that \(\delta {{\varvec{u}}}_{12}=\delta {{\varvec{f}}}_{2}/{m}_{2}\) and \(\delta {{\varvec{u}}}_{13}=\delta {{\varvec{f}}}_{3}/{m}_{3}\). Moreover, it will be assumed that \({\mathbf{B}}_{12}={\varvec{I}}\) and \({\mathbf{B}}_{13}={\varvec{I}}\), so that there are control accelerations along all 3 axes on masses \({m}_{2}\) and \({m}_{3}\), in addition to the main external force \({\varvec{f}}\) acting on mass \({m}_{3}\). Again for illustration, the individual controllers will be defined simply as:

$$\delta {{\varvec{u}}}_{12}={\mathbf{G}}_{r}\delta {{\varvec{r}}}_{12}+{\mathbf{G}}_{v}\delta {\dot{{\varvec{r}}}}_{12}$$
(32a)
$$\delta {{\varvec{u}}}_{13}= {\mathbf{G}}_{r}\delta {{\varvec{r}}}_{13}+{\mathbf{G}}_{v}\delta {\dot{{\varvec{r}}}}_{13}$$
(32b)

with gain vectors \({\mathbf{G}}_{r}\) and \({\mathbf{G}}_{v}\) defined by:

$${\mathbf{G}}_{r}=\left[\begin{array}{ccc}{G}_{r}& 0& 0\\ 0& {G}_{r}& 0\\ 0& 0& {G}_{r}\end{array}\right]{,\mathbf{G}}_{{\varvec{v}}}=\left[\begin{array}{ccc}{G}_{v}& 0& 0\\ 0& {G}_{v}& 0\\ 0& 0& {G}_{v}\end{array}\right]$$
(33)

for scalar gains \({G}_{r}\) and \({G}_{v}\). For illustration, an example of the unstable and controlled responses of the three mass equilibrium configuration is shown in Figs. 6 and 7, with non-dimensional units used such that \(G=M=d=1\). In this example, an equal mass problem is considered such that \(\kappa =2\) and so \(\lambda \cong 0.606\) with \(\gamma =1.241\). Figure 6 shows the instability of the equal mass configuration due the positive real eigenvalue(s) identified in Sect. 4. The controlled response is shown in Fig. 7, where It can be seen that the instability is supressed, as expected from the controllability analysis. The control actions are shown in Fig. 8.

Fig. 6
figure 6

Evolution of the instability of the equal mass problem with initial conditions \(\delta {{\varvec{r}}}_{12}=\left(\varepsilon ,2\varepsilon ,3\varepsilon \right)\) and \(\delta {{\varvec{r}}}_{13}=\left(\mathrm{0,2}\varepsilon ,3\varepsilon \right)\) for \(\varepsilon ={10}^{-3}\) (\(\delta {{\varvec{r}}}_{12}\left(1\right)\) and \(\delta {{\varvec{r}}}_{13}(1)\)–solid line, \(\delta {{\varvec{r}}}_{12}(2)\) and \(\delta {{\varvec{r}}}_{13}(2)\)–dashed line, \(\delta {{\varvec{r}}}_{12}(3)\) and \(\delta {{\varvec{r}}}_{13}(3)\)–dotted line)

Fig. 7
figure 7

Control of the instability of the equal mass problem with initial conditions \(\delta {{\varvec{r}}}_{12}=\left(\varepsilon ,2\varepsilon ,3\varepsilon \right)\) and \(\delta {{\varvec{r}}}_{13}=\left(\mathrm{0,2}\varepsilon ,3\varepsilon \right)\) for \(\varepsilon ={10}^{-3}\) and gains \({G}_{r}={50, G}_{v}=100\) (\(\delta {{\varvec{r}}}_{12}\left(1\right)\) and \(\delta {{\varvec{r}}}_{13}(1)\)–solid line, \(\delta {{\varvec{r}}}_{12}(2)\) and \(\delta {{\varvec{r}}}_{13}(2)\)–dashed line, \(\delta {{\varvec{r}}}_{12}(3)\) and \(\delta {{\varvec{r}}}_{13}(3)\)–dotted line)

Fig. 8
figure 8

Control effort required for the equal mass problem with initial conditions \(\delta {{\varvec{r}}}_{12}=\left(\varepsilon ,2\varepsilon ,3\varepsilon \right)\) and \(\delta {{\varvec{r}}}_{13}=\left(\mathrm{0,2}\varepsilon ,3\varepsilon \right)\) for \(\varepsilon ={10}^{-3}\) and gains \({G}_{r}={50, G}_{v}=100\) (\(\delta {{\varvec{u}}}_{12}\left(1\right)\) and \(\delta {{\varvec{u}}}_{13}(1)\)–solid line, \(\delta {{\varvec{u}}}_{12}(2)\) and \(\delta {{\varvec{u}}}_{13}(2)\)–dashed line, \(\delta {{\varvec{u}}}_{12}(3)\) and \(\delta {{\varvec{u}}}_{13}(3)\)–dotted line)

6 Applications

The use of two-body equilibria with an external force has been considered for a range of applications, as discussed in Sect. 1 (Shkadov 1987; McInnes 2002; Liu and Love 2005), including manoeuvring small near Earth asteroids. Similarly, the accelerated three-body equilibria investigated in this paper can be considered as a means of manoeuvring chains of small asteroids. Although the equilibria are unstable, as demonstrated in Sect. 4, they are controllable, as demonstrated in Sect. 5. An illustrative example is shown in Fig. 9, with three assumed uniform spherical masses \(\left({m}_{1}, {m}_{2}{,m}_{3}\right)\) with density 2 g cm−3 and a total separation \(d\) of 10 km, where the absolute size of each mass is illustrated. At 2.5 astronomical units (asteroid main belt), the Hill radius of a 50 m asteroid is approximately 20 km, and so is much larger than the maximum separation of the asteroids considered here. While the Hill radius indicates that an isolated three-body problem can be considered, a full multi-body simulation with solar perturbations (and indeed the harmonics of the asteroid gravitational potentials) could investigate the details of the dynamics and control of the asteroid chain problem.

Fig. 9
figure 9

Equilibrium configurations of three small asteroids, density 2 g cm−3 and a total separation \(d\) of 10,000 m a equal mass case with asteroid radii of 50 m b \({m}_{1}\) and \({m}_{3}\) with a radius of 50 m and \({m}_{2}\) with a radius of 150 m c \({m}_{1}\) and \({m}_{2}\) with a radius of 50 m and \({m}_{3}\) with a radius of 150 m

Firstly, the equal mass case is considered with each asteroid having a radius of 50 m. The equilibrium spacing is again found from Eq. (12) to be \(\lambda \cong 0.606\), so that the spacing of \({m}_{1}\) and \({m}_{2}\) is approximately 600 m and the spacing of \({m}_{2}\) and \({m}_{3}\) is approximately 400 m, as shown in Fig. 9. The required external force \(f\) acting on \({m}_{3}\) is found from Eq. (11b) to be 8.2 N. Moreover, from Sect. 4, the largest positive, real eigenvalue is found to be \(3.064\widetilde{\omega }\), where \(\widetilde{\omega }=\sqrt{GM/{d}^{3}}\). The instability timescale is therefore of order \(1/3.064\widetilde{\omega }\), which for the equal mass case defined above is approximately 8.25 days.

Furthermore, if \({m}_{1}\) and \({m}_{3}\) are now assumed to have a radius of 50 m, while \({m}_{2}\) has a radius of 150 m, the equilibrium spacing found to be \(\lambda \cong 0.843\), shown in Fig. 9. The required external force \(f\) acting on \({m}_{3}\) is again found to be 827.7 N. Finally, if \({m}_{1}\) and \({m}_{2}\) are assumed to have a radius of 50 m, while \({m}_{3}\) has a radius of 150 m, the equilibrium spacing found to be \(\lambda \cong 0.281\), again shown in Fig. 9. The required external force \(f\) acting on \({m}_{3}\) is found to be 840.8 N.

It should be noted that while a relatively small acceleration is required to establish an equilibrium configuration, the acceleration of the centre-of-mass of the asteroid three-body system will be small. For example, for the case above where \({m}_{1}\) and \({m}_{2}\) are assumed to have a radius of 50 m, while \({m}_{3}\) has a radius of 150 m, the acceleration of the centre-of-mass is 2.8 × 10–8 ms−2. This is equivalent to a change in speed, \(\Delta v\), of 0.87 ms−1 per year. It is noted that while the \(\Delta v\) is small, it is comparable to that considered in the literature for asteroid orbit modification (Lu and Love 2005).

7 Conclusions

An accelerated three-body problem has been presented whereby three masses interact gravitationally, while one mass is also acted on by an external force. It has been shown that the centre-of-mass of this three-body system accelerates. A single equilibrium solution has been identified where the masses form a straight line configuration. The linear stability properties of this single equilibrium configuration have been investigated, and it has been demonstrated that the configuration is always unstable. However, it has been demonstrated that the equilibrium configuration, although unstable, is in principle controllable. An application of such accelerated three-body equilibria for manoeuvring chains small asteroids has also been presented.