# The reflex instability: exponential growth of a large-scale $m = 1$ mode in astrophysical discs

Aurélien Crida, Clément Baruteau, Jean-François Gonzalez, Frédéric Masset, Paul Segretain, Philippine Griveaud, Héloïse Méheut, Elena Lega

► To cite this version:

Aurélien Crida, Clément Baruteau, Jean-François Gonzalez, Frédéric Masset, Paul Segretain, et al. The reflex instability: exponential growth of a large-scale $m = 1$ mode in astrophysical discs. Open Journal of Astrophysics, 2025, 8, 10.33232/001c.143230. hal-05222976

HAL Id: hal-05222976

https://hal.science/hal-05222976v1

Submitted on 25 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
THE REFLEX INSTABILITY:
EXPONENTIAL GROWTH OF A LARGE-SCALE *M* = 1 MODE IN ASTROPHYSICAL DISCS

AURÉLIEN CRIDA

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, laboratoire Lagrange UMR7293, Nice, France

CLÉMENT BARUTEAU

IRAP, Université de Toulouse, CNRS, Université Paul Sabatier, CNES, Toulouse, France

JEAN-FRANÇOIS GONZALEZ

Université Claude Bernard Lyon 1, CRAL UMR5574, ENS de Lyon, CNRS, Villeurbanne, F-69622, France

FRÉDÉRIC MASSET

Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México, Av. Universidad s/n, 62210 Cuernavaca, Mor., Mexico

PAUL SEGRETAIN

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, laboratoire Lagrange UMR7293, Nice, France

PHILIPPINE GRIVEAUD

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, laboratoire Lagrange UMR7293, Nice, France and
Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

HÉLOÏSE MÉHEUT

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, laboratoire Lagrange UMR7293, Nice, France

ELENA LEGA

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, laboratoire Lagrange UMR7293, Nice, France

ABSTRACT

We report the finding of a linear, non-axisymmetric, global instability in gas discs around stars, which may be relevant to other astrophysical discs. It takes the form of an *m* = 1 mode that grows in the disc density distribution while the star-barycentre distance rises exponentially with a characteristic timescale that is orders of magnitude longer than the orbital period. We present results of hydrodynamical simulations with various codes and numerical methods, using either barycentric or stlocentric reference frames, with or without the disc's self gravity: all simulations consistently show an unstable mode growing exponentially.

The instability disappears if, and only if, the reflex motion of the star due to the disc's asymmetry is not
taken into account in the simulations. For this reason we refer to this instability as the *reflex instability*. We
identify a feedback loop as a possible origin, whereby the acceleration of the star excites the eccentricity of the
disc, yielding an *m* = 1 mode in the density distribution which, in turn, pulls the star. The growth timescale
of the instability decreases with increasing disc mass and is a few hundred orbits for disc-to-star mass ratios of
a few percent. If truly physical, and not due to a numerical artifact that would be common to all the codes
we have employed, the reflex instability could have a dramatic impact on protoplanetary discs evolution and
planetary formation.

*Subject headings:* accretion, accretion discs — hydrodynamics — methods: numerical — planetary systems: formation — planetary systems: protoplanetary discs

1. INTRODUCTION

Gravitational systems in astrophysics often comprise a body
– the primary – that far outweighs the others. A typical exam-
ple is protoplanetary discs (with or without planets) around
single stars. In this Letter, we show that such a situation
leads to an instability which had so far not been reported
(to our knowledge), probably because its growth time was
too long for the simulations times and disc masses considered
usually, or because the strange behaviour of the gas disc was
attributed to some other phenomenon. However, when run-
ning simulations on thousands of orbits with disc-to-star mass
ratios larger than 10<sup>-3</sup>, its appearance is inevitable. We ar-
gue that this instability is of physical, not numerical, origin,
and occurs as soon as the disc’s axisymmetry is broken. The

key ingredient is the reflex motion of the star in reaction to
an asymmetry in the disc, often called the Indirect Term in
stlocentric codes (and sometimes discarded). A complete
discussion of the meaning, the role, and the relevance of the
Indirect Terms in general is proposed in our companion paper
(Crida et al. 2025, hereafter Paper I).

Our physical setup and the various codes we use are pre-
sented in Sect. 2, before the results of our simulations and a
description of the instability are reported in Sect. 3. A phys-
ical interpretation is proposed and the role of the different
ingredients is analysed in Sect. 4. Finally we discuss the im-
plications of our work and compare with other results in the
literature in Sect. 5, notably in the field of N-body problems
and galactic dynamics, where similar results have been found.
## 2. PHYSICAL MODEL AND NUMERICAL SETUP

### 2.1. Physical equations and reference values for the disc parameters

The physical system that we consider is comprised of a protoplanetary gaseous disc orbiting around a star. Our fiducial disc model is deliberately simple: razor-thin (two-dimensional), unmagnetized, inviscid/non-turbulent, isothermal. The disc is described by polar coordinates $(R, \phi)$ and its evolution is modelled by the continuity equation, Euler's equation and the ideal gas law, with a temperature both uniform and stationary. The choice of a globally isothermal disc allows to avoid the Vertical Shear Instability when examining the onset of our instability in supplementary 3D simulations (Nelson et al. 2013). The disc's self-gravity is included in some of our 2D simulations, which further solve Poisson's equation.

Our disc model assumes a uniform initial radial profile for the Toomre-Q parameter (Toomre 1964, see also footnote 7 in Paper I). Since the disc is assumed to be isothermal, the aspect ratio $h(R) \propto R^{1/2}$ so a uniform $Q$ implies that the initial disc's surface density profile scales as $R^{-3/2}$. The instability mechanism examined in this Letter being non-axisymmetric, white noise is added to the initial density profile at the 0.1% level in order to break the initial symmetry.

Our simulations adopt the following code units: a reference radius $R_0$ as the length unit, the stellar mass $M_*$ as the mass unit, and the inverse Keplerian frequency at $R_0$ as the time unit. In the results shown in Sect. 3, time is expressed in orbital periods at $R_0$: $T_0 = 2\pi(GM_*/R_0^3)^{-1/2}$. Our reference model has $\Sigma(R_0) = 4 \times 10^{-3}$ and $h(R_0) = 0.05$, thus $Q \sim 4$ initially. The disc extends from $0.32R_0$ to $3.125R_0$, so our fiducial value for $\Sigma(R_0)$ implies an initial disc-to-star mass ratio of $\sim 0.06$. We stress that such a physical set-up should be stable according to the present common knowledge on Keplerian discs.

### 2.2. Numerics

We have used several codes to study the linear growth stage of the instability:

* **FARGO-2D1D** (Crida et al. 2007), a sibling of the original FARGO code (Masset 2000) which solves the governing equations in the barycentric frame. We did not make use of the 1D grid in the simulations presented here.

* **FARGO-ADSG** (Baruteau & Masset 2008a,b), an extension of FARGO used to study the impact of self-gravity on the growth of the instability. Contrary to FARGO-2D1D, FARGO-ADSG solves the governing equations in the stellocentric frame.

* **FARGO3D** (Benítez-Llambay & Masset 2016), the 3D successor of FARGO, which is used to examine the instability in 3D (the code also employs a stellocentric frame).

* **IDEFIX** (Lesur et al. 2023), which, contrary to the FARGO family codes, is a code based on a higher order Godunov scheme that uses a second-order Runge-Kutta time integrator to solve the governing equations in the stellocentric frame.

* **PHANTOM** (Price et al. 2018), a 3D SPH (Smoothed Particle Hydrodynamics) code which, like FARGO-2D1D, employs a barycentric frame but, unlike all previous grid-based codes, has the obvious advantage of having no boundaries. However, SPH particles moving inside an accretion radius centred on the star are removed and

<table><caption>TABLE 1
DIFFERENCES BETWEEN THE CODES AND NUMERICAL SET-UPS USED IN OUR SIMULATIONS. COLUMN 5 LISTS THE ASSOCIATED EXPONENTIAL GROWTH TIMESCALE $\tau_g$.</caption><thead><tr><th>Code</th><th>Frame</th><th>2D/3D</th><th>Self-gravity</th><th>τ<sub>g</sub> (T<sub>0</sub>)</th></tr></thead><tbody><tr><td>FARGO-2D1D</td><td>barycentric</td><td>2D</td><td>No</td><td>300</td></tr><tr><td>PHANTOM</td><td>barycentric</td><td>3D</td><td>No</td><td>62</td></tr><tr><td rowspan="2">FARGO3D</td><td rowspan="2">stellocentric</td><td>2D</td><td rowspan="2">No</td><td>83</td></tr><tr><td>3D</td><td>130</td></tr><tr><td>IDEFIX</td><td>stellocentric</td><td>2D</td><td>No</td><td>118</td></tr><tr><td rowspan="2">FARGO-ADSG</td><td rowspan="2">stellocentric</td><td rowspan="2">2D</td><td>No</td><td>91</td></tr><tr><td>Yes</td><td>27</td></tr></tbody></table>

their mass and momentum are transferred to the sink particle representing the star. The discrete initial position of the SPH particles plays the equivalent role of white noise in our grid-based simulations to break the disc's axisymmetry (though not at the same 0.1% level).

By default, in the grid-based simulations carried out in a stellocentric frame, the disc's indirect term $ITdd$, which is the opposite of the total acceleration exerted by the disc on the star (see Paper I), is applied to each grid cell. The 2D grid-based simulations use 460 grid cells with a logarithmic spacing along the radial direction and 1256 cells along the azimuthal direction, ensuring a resolution $dR/R = d\phi = 0.005$. The 3D run performed with FARGO3D has $230 \times 628 \times 48$ cells in the radial, azimuthal and latitudinal directions, respectively, with latitudes spanning $\pm 0.15$ rad about the disc midplane. So-called wave-damping zones are used for $R \in [0.32 - 0.4]R_0 \cup [2.5 - 3.125]R_0$, where the disc fields are damped towards their initial axisymmetric profile. The PHANTOM simulation uses $5 \times 10^5$ particles, with the sink accretion radius taken equal to $r_{\text{sink}} = 0.1R_0$. In FARGO-ADSG simulations with gas self-gravity, the self-gravitating acceleration uses a smoothing length<fnref n="1" /> of $0.3H(R)$, with $H(R) = h(R) \times R$ the disc's pressure scale height.

## 3. OVERVIEW OF THE INSTABILITY

### 3.1. Results of the reference setup

We present in Fig. 1 the results of a simulation for our reference setup with FARGO-2D1D, which employs a barycentric frame. The left panel shows that the star undergoes an outward spiraling motion around the centre-of-mass of the system. This outward spiraling motion is associated with the growth of a global mode with azimuthal wavenumber $m = 1$, as can be seen in the perturbed density in the middle and right panels. Between the middle and right panels, which are separated by 100 $T_0$, the global mode has rotated by about $180^\circ$, and it is clear from the left panel that this corresponds to about half the period of the star's spiraling motion.

We point out that the growing perturbed density varies in phase with orbital radius; for instance, in the right panel, the bulk of the disc has a negative (positive) density perturbation in the southern (northern) half of the disc, but it is not the case close to the grid's inner edge. The disc's outer regions set the position of the centre-of-mass, so that the star is located opposite to the large, smooth positive perturbed density at both times. However, the disc's inner regions set the acceleration exerted by the disc on the star, as shown by the blue

<footnote n="1">¹ Since the disc is isothermal, the smoothing length in the self-gravitating acceleration does not scale linearly with $R$, which is yet required for the self-gravitating acceleration to be formally expressed as a convolution product and be computed with fast Fourier transforms (Baruteau & Masset 2008b). Using a smaller smoothing length, or a different disc model with uniform Toomre-Q parameter and uniform aspect ratio, does not qualitatively change the results of self-gravitating simulations presented in §4.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1.— Results of our reference setup for the FARGO-2D1D simulation which uses a barycentric frame. The left panel shows the position of the star in the barycentric frame. The open square symbols mark the star's position in the middle and right panels, and the arrows display the acceleration exerted by the disc on the star (arbitrary magnitude). The middle and right panels display the perturbation of the surface density of the disc relative to its azimuthally-averaged initial density, at two times separated by about half a pattern period. The solid curve shows the inner edge of the computational grid, the dashed curves show the edges of the inner and outer wave-damping zones.</figcaption></figure>

arrows in the left panel of the figure, which both point toward the positive perturbed density near the disc's inner edge.

To illustrate this spiraling motion, Fig. 2 shows the distance between the star and the centre-of-mass of the star and disc system as a function of time. The red dots on the right of the left panel show that it grows exponentially with a timescale $\tau_g \approx 300 T_0$.

## 3.2. Numerical robustness

To confirm this puzzling result, we carried out simulations using our reference set-up with other codes. Table 1 lists their differences. All simulations presented in this section neglect the disc self-gravity. The simulation run with PHANTOM in 3D yields qualitatively the same result as shown in the previous subsection with FARGO-2D1D, with the exponential growth of an $m = 1$ mode in the disc, though the growth timescale is reduced to $\tau_g \approx 60 T_0$ (see the left panel of Fig. 2, purple dots on the left). The other codes, FARGO-ADSG, IDEFIX, FARGO3D (in 2D and in 3D), which use a stlocentric frame, all find the exponential growth of this distance too. It is shown in the left panel of Fig. 2, and is equivalent to the amplitude of the $m = 1$ mode in the disc. This effect is therefore independent of the reference frame, and appears in simulations using codes of various types and solvers.

To assess whether this phenomenon could be due to a numerical effect common to all the codes used above (like a bad implementation of the Indirect Term, or the boundary conditions), we ran additional simulations with the same physical set-up, by changing a few numerical parameters. The right panel in Fig. 2 illustrates the results obtained when increasing resolution in PHANTOM and FARGO-ADSG, or removing the FARGO algorithm (hence shortening the hydrodynamical time-step) in FARGO-ADSG, or modifying the numerical scheme of integration in IDEFIX. The results are consistently independent of these numerical parameters.

We performed further tests to assess the impact of the boundary conditions for the grid codes. For instance, we ran simulations where the disc's axisymmetry was broken by the inclusion of a planet on a fixed orbit instead of white noise, and using a non-reflecting boundary condition that avoids reflections of the planet's wake off the grid's edges instead of wave-damping zones. We also ran simulations where the disc's initial density was smoothly dropped by a large factor ($\approx 100$) near the edges of the grid. Finally, we ran simulations with FARGO-2D1D that used an additional 1D grid surround-

ing the 2D grid along with the use of open 'outflow' boundary conditions at the edges of the 1D grid (David Fang, private communication). In all cases, we always found an exponential growth of the star-barycentre distance. We further stress that PHANTOM has no boundaries stricto-sensu, but the sink particle would be equivalent to an open boundary condition at a radius ($0.1 R_0$ in the default case) smaller than the inner edge of the grid in our simulations. Besides, we note that real protoplanetary disks are limited in space, and their physical boundaries may be regions of wave reflection similar to possible wave reflection at the edge of the grid in the simulations.

Finally, we varied several disc parameters and found that the growth timescale is a non-trivial combination of the disc's surface density and aspect ratio at the inner edge. The use of wave-damping zones and their detailed implementation or, in the specific case of FARGO-2D1D, the fact that the inner edge of the grid is centred on the centre of mass of the system and not on the star, may affect the physical conditions at the effective inner edge radius. A last PHANTOM simulation, using $r_{\text{sink}} = 0.32 R_0$ and plotted in the right panel of Fig. 2, shows a longer $\tau_g \sim 110 T_0$ more in line with the grid simulations, which would support this interpretation. We think that this explains the quantitative differences observed in our simulations (in the results shown in Fig. 2, the growth timescales are $\tau_g = 113 T_0 \pm 0.2$ dex). In particular, we note that the two codes with the most extreme values of the growth timescale are PHANTOM and FARGO-2D1D, that is the ones with a different geometry of the inner edge. In the next section, we present arguments in favor of a physical origin for the qualitatively consistent appearance of the exponential growth of the star-barycentre distance in all our simulations.

## 4. PHYSICAL INTERPRETATION

**A positive feedback loop**— Previous sections strongly suggests that the instability found using either a stlocentric frame with ITdd or a barycentric frame is of physical origin. However, when the same simulations are run in the stlocentric frame without ITdd, the disc remains stable for ever. It is expected indeed for this set-up that the disc is intrinsically stable. Similarly, a barycentric simulation in which the star does not feel the disc acceleration remains stable for ever. This shows that the instability originates in all the simulations from a positive feedback loop between an asymmetry in the disc and the reflex motion of the star.

Further insight into the feedback loop can be gained by
<figure><img src="image_3.png" /><figcaption>FIG. 2.— Time evolution of the radial distance between the star and the centre-of-mass of the star-disc system for simulations using our reference setup, obtained for different codes and by varying numerical parameters (see caption and text).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3.— Results of our reference setup for FARGO-ADSG simulations using a stellocentric frame. The left panel shows the time evolution of the radial distance between the star and the centre-of-mass of the star-disc system for four simulations with or without self-gravity ('sg'), and with or without the indirect term of the disc on itself ('ITdd'). The middle and right panels display the perturbation of the surface density of the disc relative to its azimuthally-averaged initial density for the simulations with ITdd, without self-gravity (middle panel) and with self-gravity (right panel). The solid curve shows the inner edge of the computational grid, the dashed curves show the edges of the inner and outer wave-damping zones.</figcaption></figure>

integrating numerically the motion of a test particle that is initially on a circular orbit around a star. When applying a constant acceleration on the star, the eccentricity of the test particles rises linearly, with the longitude of pericentre shifted by $-\pi/2$ with respect to the acceleration vector (see also Namouni 2005). Now, if the acceleration is set proportional to the particle's eccentricity, and phase-shifted by $+\pi/2$ from the direction of pericenter, we find an exponential increase in the particle's eccentricity (see Appendix A). Similarly, in hydrodynamical simulations where we apply a constant and uniform acceleration to the whole disc (like the indirect term), the disc eccentricity grows and an $m=1$ mode appears in the perturbed density. In turn, this $m=1$ mode exerts a force on the star<fnref n="2" /> which is proportional to the amplitude of the mode

(thus to the eccentricity) and, as pointed out in § 3.1, may not point towards the centre of mass of the system if this mode has a phase that varies with radius. If the resulting acceleration has a component pointing at $+\pi/2$ from the longitude of the pericentre at a given radius, the loop is closed.

**Role of self-gravity** — All the simulations presented in section 3 suffer however – like most disc simulations in the literature – from an inconsistency which we put forward in Paper I: while the star feels the gravitational acceleration from the disc, the disc does not, because we do not include its self-gravity. This imbalance can be problematic, and typically makes tidal effects impossible to account for. Therefore, we have run our reference set-up with FARGO-ADSG, with the gas self-gravity on. The results are shown in Fig. 3, right panel, compared to the middle panel where the self-gravity is off, at different

<footnote n="2">² $m > 1$ modes result in zero acceleration of the star.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4.— Growth timescale $\tau_g$ versus the initial surface density $\Sigma_0$ at the reference radius $R_0$, with and without self-gravity. The dashed lines show suggestive trends.</figcaption></figure>

times corresponding to a similar amplitude of the mode. In both cases, the star-barycentre distance grows exponentially, as shown in the left panel, and only the timescale and the shape of the mode change. The growth timescale $\tau_g \approx 27 T_0$ for the self-gravitating run, that is $\approx 3.3$ times shorter than for the same run without self-gravity. Conversely, in a simulation with only self-gravity but no ITdd, the disc remains stable. Of course with neither self-gravity nor ITdd, nothing happens either. Therefore, the feedback loop described above is *not* caused by the absence of self-gravity in the disc. The stellar acceleration due to the disc is sufficient and necessary to trigger this instability. Self-gravity then modifies the response of the disc quantitatively.

*Dependence with the disc's surface density*— While a detailed characterisation of the instability is beyond the scope of this Letter and is deferred to a future, on going study, we briefly present the results of simulations where we varied the initial surface density at the reference radius ($\Sigma_0$). Such simulations were carried out with FARGO-ADSG, with and without self-gravity, and the uniform Toomre-Q parameter varies from 2 to 16 across simulations. For each simulation, the growth timescale $\tau$ is inferred from the time evolution of the star-barycentre distance. Figure 4 shows that $\tau \propto 1/\Sigma_0$ for non-self-gravitating runs, as expected for a linear perturbation where the feedback is proportional to $\Sigma_0$. Self-gravitating simulations show a slightly steeper dependence with $\Sigma_0$ ($\tau \propto \Sigma_0^{-3/2}$).

## 5. DISCUSSION

We have found that non-axisymmetric Keplerian discs are linearly unstable against the growth of an $m=1$ mode arising from the reflex motion of the star around the centre-of-mass of the star and disc system. For this reason we shall dub it *reflex instability* for future reference.

We identified a feedback loop which could be the cause of this instability, and described how the disc eccentricity and the acceleration of the star can excite each other. However, due to the lack of a fundamental linear analysis so far (work in progress), there is still a risk that this instability is of numerical origin. Even if it were the case, we stress that it seems to be ubiquitous in numerical codes, and we warn the community that they may encounter this exponential growth of an $m=1$ mode in the disc in their long-term simulations that include ITdd.

At first glance it may seem surprising that the reflex instability has not been reported before for typical protoplanetary disc masses. The explanation could be two-fold. On the one hand, for typical disc-to-star mass ratios below a few percent, the growth timescale can exceed $10^3$ dynamical timescales, thereby requiring rather long-term simulations in order not to fly under the radar. Actually, some of us (and other colleagues, e.g. G. Pichierri who observed the growth of an oscillation of the torque felt by a planet around its average value, private communication) have been surprised in the past by the sudden and unexplained excitation of the disc and/or crash of the code in very long term simulations; in retrospect, we think this strange behaviour can be attributed to the reflex instability, whose amplitude was negligible for most of the time. On the other hand, most simulations of star-disc interactions use a stellocentric frame and thus require the disc's indirect term to be activated for the instability to manifest itself, and it turns out that the inclusion of indirect terms has long remained nebulous (see Paper I). In fact, in most non-self-gravitating simulations of discs without planets, ITdd is not included (e.g. Robert et al. 2020).

This reflex instability shares characteristics with the work of Heemskerk et al. (1992). They did a linear stability analysis of $m=1$ disturbances showing that the shift of the star away from the system's centre-of-mass leads to the growth of an instability in large ($\gtrsim 1$) disc-to-star mass ratios. This instability was further confirmed by non-linear hydrodynamical simulations for similarly large disc-to-star mass ratios. The authors concluded that the instability disappears when the disc mass becomes lower than the star mass. We note however that their linear analysis actually discards collective effects (namely self-gravity and pressure), so that it does not apply to our framework. In their section 3.5, the authors state that if self-gravity and pressure were included in the disc response, an instability may be allowed to persist "at some level" for any disc-to-star mass ratio. Furthermore, in their numerical simulations, they report finite growth rates even for their lowest disc-to-star mass ratio, 1/9, which is barely twice the value we use in our simulations. A dedicated linear analysis in the regime of low disc-to-star mass ratios will be presented in a future study.

Interestingly, Peng & Batygin (2020) have observed a growth of the eccentricity in pure N-body simulations containing only small particles, that is particles that have a small mass but do not feel each other's direct gravity. As stated by the authors, this unexpected dynamical excitation stems from the reflex motion of the star due to the global motion of the cloud of small particles. However, since the small particles do not have a collective behaviour, no $m=1$ mode can grow so that in their case, there is no feedback loop and no exponential growth.

Also, Saha et al. (2007) have studied the persistence of a $m=1$ mode in galactic discs with an exponentially decreasing surface density profile and no central mass, but with self-gravity. Although the physical set-up is different from ours and the rotation velocity is far from Keplerian, they find that the mode precesses very slowly and that "the disc is never unstable to the lopsided [$m=1$] modes without the indirect term arising due to the lopsided perturbation itself". The reflex instability therefore seems to have astrophysical applications well beyond the field of protoplanetary discs and to be a quite generic phenomenon.

In this paper we have only explored the linear growth stage of the reflex instability, intentionally leaving its non-linear development and its full characterisation for future work. Note that the grid-based simulations of the present study did not show any saturation of the instability, until the disc reaches
quite extreme eccentricities, causing most often a crash of the simulation. The growth of a global $m = 1$ mode could limit the lifetime of massive discs and excite the disc's eccentricity with interesting implications for scenarios of planet formation and orbital evolution. It may also imply the formation of a large-scale spiral throughout the disc, with potential connections to spiral structures observed in the gas emission of protoplanetary discs.

## ACKNOWLEDGMENTS

We thank very much several colleagues for enlightening discussions on this subject: William Béthune, Alessandro Morbidelli, Gabriele Pichierri, Nathan Magnan, Fabiola Gerosa, Dylan Kloster, Kate Minker, Léa Marques, Thomas Rometsch and Lucas Jordan. We also thank the anonymous referee whose comments helped to improve this manuscript.

We would like to warmly thank the "Programme National de Physique Stellaire" (PNPS) and "Programme National de Planétologie" (PNP) of CNRS/INSU co-funded by CEA and

CNES for their financial support of our BiPhasique research group. JFG thanks the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) for its financial support within the Plan France 2030 of the French government operated by the ANR. FM acknowledges support from UNAM's grant PAPIIT 107723, UNAM's DGAPA PASPA program and the Laboratoire Lagrange at Observatoire de la Côte d'Azur for hospitality during a one-year sabbatical stay. Part of the numerical simulations were performed on the CALMIP Supercomputing Centre of the University of Toulouse, and others on "Mesocentre SIGAMM" hosted by Observatoire de la Côte d'Azur. The authors are also grateful to the Université Côte d'Azur's Center for High-Performance Computing (OPAL infrastructure) for providing resources and support. We acknowledge support by DFG-ANR supported GEPARD project (ANR-18-CE92-0044 DFG: KL 650/31-1). This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence Université Côte d'Azur under reference number ANR-15-IDEX-01.

# APPENDIX

## LINK BETWEEN ACCELERATION AND ECCENTRICITY

### Laplace vector, eccentricity, and acceleration

For a particle in orbit around a central object of mass $M$, we remind that the Laplace vector is defined as $\vec{P} = \frac{1}{GM}(\vec{v} \times \vec{h}) - \vec{u}_r$, where $\vec{v}$ is the velocity of the test particle, $\vec{h} = \vec{r} \times \vec{v}$ is its specific orbital angular momentum, and $\vec{u}_r$ is the unit vector in the radial direction. One can show that $\vec{P}$ is constant: its norm is equal to the eccentricity of the orbit and its direction points towards the pericentre. Without loss of generality, we choose here the longitude of the pericentre $\varpi$ to be zero, so that $\vec{P} = e \vec{u}_x$, where $\vec{u}_x$ is the unit vector in the $x$-direction; $\vec{u}_y$ and $\vec{u}_z$ complete an orthonormal base with $\vec{u}_z$ perpendicular to the orbital plane, such that $\vec{h} = \sqrt{GMa(1 - e^2)}\vec{u}_z$ (where $a$ is the semi-major axis).

From the definition of $\vec{P}$ above, one can show that

$$ \vec{v} = \sqrt{\frac{GM}{a(1 - e^2)}} (\vec{u}_\phi + e\vec{u}_y) = \sqrt{\frac{GM}{a}} (\vec{u}_\phi + e\vec{u}_y + \mathcal{O}(e^2)). \quad (A1) $$

In other words, for small eccentricities, the velocity is decomposed into a uniform circular motion ($\vec{u}_\phi$ is the unit vector in the azimuthal direction) and a uniform linear motion parallel to the minor axis of the ellipse.

Hence, applying a constant acceleration $\vec{a}_c = \varepsilon\vec{u}_y$ (e.g. an indirect term) to a particle on an initially circular orbit should lead to an increase of the linear component of the velocity $v_y = \sqrt{\frac{GM}{a}}e$, that is an increase in the eccentricity as $\dot{e} = \varepsilon/\sqrt{GM/a}$.

This is confirmed by a simple numerical simulation where the motion of a test particle is integrated through a Runge-Kutta-4 algorithm with an acceleration given by $-r^{-2}\vec{u}_r + \vec{a}_c$ with $\vec{a}_c = 0.004\vec{u}_y$ (units are chosen so that $G = M = 1$ and the particle starts on a circular orbit of radius 1). The top left panel of Fig. 5 shows the resulting trajectory; it goes from circular to eccentric with a linear increase in the eccentricity (bottom right panel) and $\varpi = 0$.

### Creating a feedback loop

In a second simulation, we set the additional acceleration to be $\vec{a}_c = 0.02 \cdot R_{\pi/2}(\vec{P})$, where $R_\phi$ is a rotation of angle $\phi$ around the $z$-axis (in the $x-y$ plane). Since the acceleration is now proportional to the eccentricity, this leads to an exponential growth of the eccentricity, as shown in the top right and bottom right panels of Fig. 5. The initial orbit was again circular, but rounding errors are enough to initiate a positive feedback loop whereby the eccentricity triggers an acceleration which increases the eccentricity (as seen above). As a consequence, the Laplace vector grows exponentially in norm, keeping its direction constant (here $\vec{u}_x$).

If now the angle of the rotation is not exactly $\pi/2$ but between 0 and $\pi$, we observe a slower, but still exponential growth of the eccentricity, with a linear precession of the longitude of the pericentre. The result is shown in the bottom left and bottom right panels of Fig. 5 for an angle of $\pi/8$. This can be interpreted as $\vec{a}_c$ having a component perpendicular to $\vec{P}$, which yields the exponential growth of $e$, and a component parallel to $\vec{P}$, which yields the precession of the orbit.

### Link with the reflex instability

The above results illustrate how a positive feedback loop can easily ignite in our set-up. The eccentricity growth will build an $m = 1$ mode in the gas disc. In turn, this asymmetry in the density will accelerate the star, causing an indirect term proportional to the eccentricity. For a simple eccentric gas disk (uniform eccentricity and longitude of periastron), the surface
<figure><img src="image_6.png" /><figcaption>FIG. 5.— Orbital evolution of a particle feeling an acceleration $\vec{a} = -(1/r^2)\vec{u}_r + \vec{a}_c$. Top left panel: $\vec{a}_c = 0.004 \vec{u}_y$. Top right: $\vec{a}_c = 0.02 R_{\pi/2}(\vec{P})$. Bottom left panel: $\vec{a}_c = 0.02 R_{\pi/8}(\vec{P})$, where $R_{\phi}$ is the rotation by an angle $\phi$ in the x–y plane. Bottom right panel: evolution of the eccentricity in the three simulations (same color).</figcaption></figure>

density perturbation due to the eccentric motion is symmetric with respect to major axis. Hence, the indirect term would be strictly parallel to the Laplace vector in this case. However, in Fig. 1 the direction of the acceleration is not parallel to the star–barycentre line. Therefore, it may have a component exciting the eccentricity and the *m* = 1 mode in some part of the disc, yielding an exponential growth of the star–barycentre distance, and a perpendicular component which makes the system precess linearly. Hence, the star describes a logarithmic spiral, as observed.

Describing analytically the response of a Keplerian gas disc to a uniform acceleration, and deriving the acceleration this response exerts on the star, are the missing ingredients to fully characterise the reflex instability. Nonetheless, this simple analysis in the two-body problem illustrates how the feedback loop may work.

REFERENCES

Baruteau C., Masset F., 2008a, ApJ, 672, 1054

Baruteau C., Masset F., 2008b, ApJ, 678, 483

Benítez-Llambay P., Masset F. S., 2016, ApJS, 223, 11

Crida A., Morbidelli A., Masset F., 2007, A&A, 461, 1173

Crida A., et al., 2025, The Open Journal of Astrophysics, 8

Heemskerk M. H. M., Papaloizou J. C., Savonije G. J., 1992, A&A, 260, 161

Lesur G. R. J., Baghdadi S., Wafflard-Fernandez G., Mauxion J., Robert C. M. T., Van den Bossche M., 2023, A&A, 677, A9

Masset F., 2000, A&AS, 141, 165

Namouni F., 2005, AJ, 130, 280

Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610

Peng S., Batygin K., 2020, ApJ, 898, L46

Price D. J., et al., 2018, PASA, 35, e031

Robert C. M. T., Méheut H., Ménard F., 2020, A&A, 641, A128

Saha K., Combes F., Jog C. J., 2007, MNRAS, 382, 419

Toomre A., 1964, ApJ, 139, 1217

This paper was built using the Open Journal of Astro-
physics LATEX template, which was slightly modified by the

authors. The OJA is a journal which provides fast and easy
peer review for new papers in the astro-ph section of the
arXiv, making the reviewing process simpler for authors and
referees alike. Learn more at http://astro.theoj.org.