Solar wind dynamics around a comet: The paradigmatic inverse-square-law model

Aims. Observations of solar protons near comet 67P/Churyumov-Gerasimenko (67P) by the Rosetta spacecraft can be modelled by the planar motion in an effective magnetic field proportional to 1/r^2. We aim to provide a thorough study of such dynamics, with a clear description of the behaviour of an incoming flux of particles. We will be able, then, to calibrate the free parameters of the model to Rosetta observations. Methods. Basic tools of dynamical analysis are used. They lead to a definition of the relevant parameters for the system and a classification of the possible types of trajectories. Using the so-obtained formalism, the structures formed by a flux of particles coming from infinity can be studied. Results. All the trajectories are parametrised by two characteristic radii, rE and rC, derived from first integrals. There are three different types of motion possible divided by a separatrix corresponding to rE = rC. An analytical expression of the trajectories, defined by an integral, is developed. Using this formalism, the application to a flux of particles coming from infinity (modelling the incident solar wind) gives one free parameter only, the radius rE, which scales the problem. A circular cavity of radius 0.28 rE is created, as well as an overdensity curve (analogous to a caustic in optics). At each observation time, rE can be calibrated to Rosetta plasma measurements, giving a qualitative understanding of the solar particle dynamics (incoming direction, cavity and density map). We also deduce that, in order to properly capture the essence of the dynamics, numerical simulations of the solar wind around a comet must use simulation boxes much larger than rE and grids much finer than rE.


Introduction
Plasma instruments on board the Rosetta mission have provided invaluable information about the dynamics of solar and cometary ions in a comet neighbourhood (comet 67P/Churyumov-Gerasimenko, or 67P). For the first time, these dynamics were followed as the comet nucleus activity was evolving, from more than 3.8 astronomical units (au) to its perihelion at 1.2 au. Because of the lower production rate of 67P's nucleus together with these large heliocentric distances, the interaction between the solar wind and the cometary atmosphere (coma), was fundamentally different from what was previously observed at more active comets closer to the Sun (comet Halley or comet Giacobini-Zinner: Grewing et al. 1988;Cowley 1987). Because the gyration scale of the cometary ions is larger than the interaction region, the classical fluid treatment of the plasmas does not apply at 67P, and a kinetic description of the interaction is necessary. In Behar et al. (2018b) the terms "fluid comet" and "kinetic comet" were introduced to separate these two different regimes. As of now, only self-consistent numerical models have tackled the interaction between the so-lar wind and the coma of a "kinetic" comet (Bagdonat & Motschmann 2002;Hansen et al. 2007;Rubin et al. 2014b;Koenders et al. 2016a,b;Behar et al. 2016;Deca et al. 2017;Huang et al. 2018). All these models result in a highly asymmetric plasma environment, in contrast with the classical symmetric picture obtained for more active comets closer to the Sun (Rubin et al. 2014a).
In the context of comet 67P and based on in-situ data, Behar et al. (2017) recently showed that a cavity completely free of solar particles is created around the comet's nucleus, surrounded by a region where they are focused in a specific direction. Moreover, the measured velocity of the solar protons is almost constant in norm throughout the mission, indicating that they are deflected without significant loss of energy. Remarkably, the main plasma observations are very well reproduced by a simple inverse-square-law "effective magnetic field" applied to the incoming flux of solar protons. Using this model, the density and velocity profiles become natural geometrical effects, also in qualitative agreement with numerical simulations (Behar et al. 2016(Behar et al. , 2017. Due to the striking success of this empirical approach, it became necessary to outline the intrinsic properties of this force field: this will allow us to state clearly what it would imply for the dynamics of solar wind protons, and hopefully to link the observables to physical quantities. The aim of this paper is to provide a full characterisation of the planar dynamics of charged particles in a magnetic field proportional to 1/r 2 . This way, the appropriate formalism will be available for further applications to the Rosetta mission or any analogous physical modelling. In particular, the behaviour of an incoming flux of charged particles in a 1/r 2 magnetic field has only been explored by numerical means so far, and its precise characteristics are still missing. Consequently, this paper is mainly devoted to dynamical aspects, and we will only hint at the physical considerations regarding its application to comets. Crucial discussions about the nature of this force and comparisons to self-consistent physical models of comet 67P are presented in companion papers (Behar et al. 2018b,a).
The planar dynamics in an inverse-cube-law magnetic field has been thoroughly studied by physicists because it describes the motion of charged particles in the equatorial plane of a magnetic dipole (e.g. Störmer 1907Störmer , 1930Graef & Kusaka 1938;Lifshitz 1942;de Vogelaere 1950;Avrett 1962). The interest for such dynamics was greatly enhanced by its direct applications to the geomagnetic field. Unbounded and bounded solutions exist and trapped particles are indeed observed around the Earth (Williams 1971). The deflection of an incoming flux of particles seems to produce similar structures as for an inverse-square law (compare Fig. 1 by Shaikhislamov et al. 2015 with Fig. 3 by Behar et al. 2017). Besides, complex plasma interactions in other physical contexts could possibly be modelled as well by such simple laws. A comparative study of the different powers of 1/r would thus be also valuable.
This paper is organised as follows: Section 2 presents a general study of the inverse-square-law magnetic field. After having summarised the model developed by Behar et al. (2018b), we define the different types of possible orbits, outline their properties, and give for them an analytical expression defined by an integral. In Sect. 3, we apply this formalism to an incoming flux of particles, similar to the solar protons. The properties of the cavity and of the overdensity region reported by Behar et al. (2017) are fully characterised. Then, Sect. 4 presents order-of-magnitude estimates of the characteristic quantities of the model calibrated on the plasma observations realised by the Rosetta spacecraft.
Additionally, the comparison of dynamics produced by magnetic fields proportional to an arbitrary power of 1/r is given in Appendix B: it could serve as reference when dealing with analogous problems.
2. General study of the dynamics 2.1. The inverse-square law for solar protons around comets The analytical model introduced by Behar et al. (2018b) is built from three sub-models. We outline here their main characteristics (readers mainly interested in dynamical aspects can safely go to Sect. 2.2).
Steady state is always assumed, implying that the change of heliocentric distance of the comet is slow enough to be considered as an adiabatic process. The first submodel is a description of the ionised coma and its density distribution. The cometary atmosphere is assumed to have a spherical symmetry: the neutral elements are produced at a rate Q and expand radially in all directions with constant velocity u 0 . The cometary ions are created from these neutral elements with a rate ν i (mainly by photo-ionisation and electron-impact ionisation). They initially have the radial velocity u 0 , but they are accelerated by the local electric and magnetic fields and lost from the system. This is taken into account by a "destruction" rate ν ml (where ml stands for mass loading). Therefore, in the regime of the system under study, the local density of cometary ions can be written as where R is the radial distance from the nucleus (see Behar et al. 2018b, for details). In this description, the ionised component of the coma is essentially made of the slow, new-born cometary ions, which are steadily created and lost. The second sub-model is a description of the magnetic field piling up due to the local decrease in the average velocity of the electrons (as slow new-born ions are added to the flow). The magnetic field B is considered frozen in the electron fluid, the latter coming from infinity on parallel trajectories. The third sub-model is a description of the electric field, which is reduced to its main component, the so-called motional electric field, where u i is the average velocity of all charges carried by solar and cometary ions. Considering only the Lorentz force, this results in a generalised gyromotion for both populations, where the two gyroradii depend strongly on the density ratio. This generalised gyromotion is the core of the model, giving a mechanism through which energy and momentum are transferred from one population to the other. For simplicity, we finally consider that the cometary particles are mainly composed of water, resulting in the same charge +e for the cometary (H 2 O + ) and solar wind (H + ) ions. Putting all things together, the force applied to the solar protons is In this expression, x is the position vector of the proton, m is its mass, and the dot means time derivative. The magnetic field B ∞ is the one carried by the solar wind before its encounter with the comet, and n sw is the average density of solar wind protons. Injecting the expression of n com from Eq. (1), we finally get an inverse-square law like the one studied in the rest of this article. It should be noted that the force applied here to solar wind protons is not a magnetic field as such, but it behaves like one. Hence, even if we speak generically of "magnetic field" throughout this article, the reader should understand "effective magnetic field" to mean a vector field behaving as a magnetic field but possibly produced as a result of more complex interactions. In our case, the relevant dynamics takes place in a plane, written (x, y) in the following. In the comet-Sun-electric frame (CSE) used by Behar et al. (2017), this plane contains the comet, the Sun, and the electric field vector produced by the incoming solar wind. 1 During its operating phase around comet 67P, the Rosetta spacecraft was not far from this plane (since it was not far from the comet itself). In other contexts, the (x, y) plane used here could be the equatorial plane of some source of magnetic field.

Equations of motion
A particle of mass m, charge q, and position x = (x, y, z) T is subject to a magnetic field of the form From the classical Lorentz force, the equations of motion are mẍ = qẋ × B(x), that is, in which the constant k has the dimension of length time velocity and the dot means derivative with respect to the time t. From Eq. (5), the vertical velocity is constant and imposed by the initial conditions. We are interested here in the dynamics in the (x, y) plane. Let us introduce the polar coordinates (r, θ). The equations of motion rewrite 2 as Since the force is always perpendicular to the velocity vector, its norm is constant (conservation of the total energy E). This leads to the first integral v = ṙ 2 + r 2θ2 , equal to the norm of the velocity projected in the (x, y) plane. Moreover, Eq. (7) is directly integrable: where the arbitrary constant r is added for dimensionality reasons. This leads to a second first integral, It can be thought as the conservation of a generalised angular momentum, coming from the symmetry of rotation around the z-axis.
solar wind studies (the effective magnetic field is thus oriented along the y axis). We think that the notation (x, y) is more appropriate for the present paper, focussed on dynamics only. This should not be too confusing for the reader. 2 We get here the same equations as Graef & Kusaka (1938). This comes from a mistake in their paper: they begin with the equations of a 1/r 2 field; they introduce the conserved quantities of a 1/r 3 field; they write down equations mixing both types of fields, and eventually, they study the 1/r 3 one for the rest of the paper. Since they deal with the motion in the equatorial plane of a magnetic dipole, 1/r 3 is the correct field to use.
From dimensionality arguments, the conservation of v makes a characteristic length and a characteristic frequency of the system appear: In the same way, the generalised angular momentum c can be turned into the characteristic length r C = r exp(c/k + 1) = r exp(r 2θ /k + 1) .
As we will see, the dynamics of the particle is entirely contained inside the independent constants r C and r E . Their physical meaning will appear later. Similarly to Störmer (1930), it is convenient at this point to use the normalised quantities We note that if k > 0, the direction of time is reversed. In the new coordinates, the equations of motion (6-7) become where this time the dot means derivative with respect to the normalised time τ (this double use of the dot should not be confusing for the reader, since t is only used with the dimensional r coordinate, while τ is only used with the dimensionless ρ coordinate). Equations (14-15) are equivalent to:ρ coming respectively from the energy and the generalised angular momentum (we have in particular ρ C = r C /r E ).

Geometry of the trajectories
Introducing ρ C (Eq. 17) into the energy expression (Eq. 16), we geṫ where we define We will see below that both ρ 0 and ρ C have a precise dynamical meaning. Since they are directly proportional, the problem can be indifferently parametrised by one or the other. For simplicity, we will consider either ρ 0 or ρ C in the following, according to the dynamical feature under discussion.
The function U can be interpreted as an effective potential, which counterbalances the kinetic term at all times. Its general form gives directly the values of ρ allowed as a A&A proofs: manuscript no. EBM2 Effective potential as a function of ρ. Changing the value of parameter ρC is equivalent to rescaling the axes. The unit level on the vertical axis gives the intervals of ρ allowed for the particle, such that U (ρ) < 1. These intervals are delimited by ρ1, ρ2, and ρ3, given at Eq. (22). b) Angular velocity as a function of ρ. c) Examples of trajectories for three orbit types, obtained by using the expression from Eq. (26) with parameters ρC = (1.01, 0.5, 0.9) for (red, green, blue). The axes are rescaled such that ρC appears the same, as in graph a.
function of the parameter ρ C (Fig. 1a). Noting {T i } the different types of trajectories, we have 3 We note that ρ 2 and ρ 3 are only defined if ρ C < 1.
The extreme values of ρ reachable by the particle in the different cases are thus where W 0 and W −1 are the Lambert functions. In accordance with Eq. (20), ρ 2 and ρ 3 are only defined if ρ 0 < exp(−1), that is, ρ C < 1.
On the other hand, the conservation of ρ C (Eq. 17) allows us to write the angular velocity as a function of ρ only (Fig. 1b). The stable equilibrium point at ρ = ρ 0 corresponds to v = 0 (motionless particle), and the unstable equilibrium point at ρ = ρ C corresponds to a circular orbit with constant angular velocity. Whatever the trajectory, the angular velocity vanishes at ρ = ρ 0 and changes sign. The inner part of T 1 trajectories shows thus a unique loop away from the origin, whereas T 2 trajectories continuously rotate around the radius ρ 0 . On the contrary, T 3 trajectories always rotate in the same direction around the origin (see Fig. 1c for some examples).
A parametric expression of the trajectories can be easily obtained from the first integrals. Indeed, from Eqs. (17) and (18) we geṫ which gives with in which the initial conditions are written (ρ i , θ i ). One can note that the integrand is singular in the extrema of ρ (Eq. 22), but the integral itself is always convergent (except for ρ C = 1, since in this case the particle makes an infinite number of loops before eventually reaching ρ = ρ C ). The ± sign in Eq. (24) stands for the branches approaching (−) and leaving (+) the minimum radius. This definition by parts can be avoided by parametrising the trajectories by a parameter s ∈ R. A possible parametrisation θ(s), ρ(s) of the three types of non-singular trajectories is given by where ∆ = ρ 2 − ρ 1 . With this parametrisation, s(τ ) increases with τ and ρ(s) is minimum at s = 0.

Time information
By expressing the term ρ 2θ2 from the energy constant (Eq. 16) and by injecting it in the first equation of motion (Eq. 14), we get which can be directly integrated to give The polar angle is thus composed of a linear part plus a term proportional toρ. The physical meaning of the frequency ω E (Eq. 11) is now clear: it is the drift angular velocity of every particle. This is of particular interest for bounded trajectories. Indeed, they are quasi-periodic, with two proper frequencies: the "drift" frequency (rotation around the origin) and the "loop" frequency (small loop around the ρ 0 radius). Sinceρ vanishes at the extreme values of ρ (Eq. 22), the period T C of the loops is while the period of the overall rotation around the origin is simply T E = 2π (that is 2π/ω E in dimensional coordinates). Bounded trajectories represented in a frame rotating with τ consist thus only in the small loop around ρ 0 . Periodic orbits are produced when the fraction T C /T E is a rational number. Figure 2 shows the behaviour of T C /T E as a function of the parameter ρ C , along with some examples of periodic trajectories. We note that the two frequencies tends to be equal at ρ C , that is, for the circular unstable trajectory (for whichρ = 0 at all time). More generally, Eq. (28) can be used to express the time as a function of ρ just like θ in Eq. (24). Expressing ρρ from Eq. (18), we get where ± means (−) when the particle gets closer to the origin, and (+) when it goes back. As before, a parameter s ∈ R can be used to avoid this double definition: This equation could have been obtained also directly from the energy constant (see Appendix B). It can be added among the parametrisation given by Eq. (26) in order to compute the time at every position. From Eq. (30), one can note that in order to compute θ and τ at a given value of s, it is enough to compute only one integral.

Application to an incoming flux of particles
Our first motivation for studying the inverse-square-law magnetic field is the deflection of solar wind protons as a result of their interactions with a cometary-type atmosphere. At very large distances from the comet, they can be considered as following parallel trajectories. In this section, we thus consider a permanent flux of particles initially evolving on parallel trajectories. As before, the z-component of the dynamics is trivial. We choose the orientation of the reference frame such that the initial velocity of the particles projected in the (x, y) plane is along the x-axis (ẋ i = −v with v > 0). At the position x i = d, the magnetic field from Eq. (4) is activated. We aim to determine how the particles are distributed in the plane (x, y) in the permanent regime, and in particular when d → ∞. For finite d distances, as we will consider in a first step, one can think of a continuous source of particles with the shape of a vertical infinite "wall". A similar setup was studied numerically by Shaikhislamov et al. (2015) in a dipole 1/r 3 field, but with the addition of a magnetopause (the particles were launched from a curved line instead of a fixed horizontal distance). As we will see, the two situations create similar features.

The cavity
Since the particles have all the same velocity v, they have the same characteristic radius r E and drift frequency ω E given by Eq. (11). We are thus able to use the normalised variables ρ = r/r E and τ = ω E t (same as previous section) in order to describe their motions in a common way. However, the characteristic radius ρ C of each particle is a function of its initial position along the Oy axis. Using the normalised coordinates Y i = y i /r E and D = d/r E , we get from Eq. (17) and The problem is about determining the different types of orbits followed by the particles as a function of D and of their initial position Y i . In the following, we suppose that k is positive. 4 First of all, we note that The particles have thus all the possible values of ρ C , including the critical one ρ C = 1 (Eq. 20). Let us write A&A proofs: manuscript no. EBM2 the limiting distance above which ρ C (Y i ) is monotonous. For D > D lim , there is thus only one trajectory with ρ C = 1 among the initial positions Y i . For D < D lim , on the contrary, ρ C (Y i ) has a local maximum (larger than 1) and a local minimum. It is straightforward to show that there is a critical distance, such that if D < D crit , the local minimum of ρ C (Y i ) is smaller than 1. This produces two other critical trajectories with ρ C = 1 (or only one in the limiting case D = D crit ). In Appendix A, Fig. A.1 shows the behaviour of all the characteristic lengths as a function of D, where the meaning of D lim and D crit is obvious. As a summary, the types of orbits followed by the particles are colour-coded in Fig. 3 as a function of their initial position. Trajectories of types T 2 and T 3 can be distinguished by considering their initial radius ρ i = D 2 + Y 2 i , which should be smaller or larger than 1, respectively (since ρ 2 < 1 and ρ 3 > 1).
In order to determine the distribution of the particles in the plane, useful information is given by the extreme radii reached by the particles. Each of them can be expressed in terms of D and Y i by using Eqs. (22) and (33). For a fixed distance D, the minimum radius reached by the whole flux of particles is given by the minimum of ρ 1 (Y i ) over the T 1 and T 2 sets of orbits (red and green regions of Fig. 3). For D > D crit , it corresponds to the inner loop of the asymptotic trajectory, that is, at the very limit between the red and blue zones of Fig. 3. On the contrary, for D < D crit , the minimum value of ρ 1 is reached in the T 2 zone. The value of the minimum is where γ = 1 2 1 − √ 1 − 4D 2 . Interestingly, for D D crit , the value of ρ cav is independent of the starting distance D. Moreover, we note that whatever the value of D > 0, the minimum distance ρ cav is never zero. This implies that among the whole flux of particles, none reaches the origin. In other words, the magnetic field naturally creates a cavity around the origin, devoid of any particle. The shape of this cavity can be inferred as follows: • For D < D crit , the minimal radius is reached by a bounded orbit of type T 2 . Hence, particles following this orbit come back periodically in ρ cav and spread at all θ values (if we exclude periodic orbits as in Fig. 2). The cavity in the permanent regime is thus circular.
• For D D crit , the minimal radius is reached by an unbounded orbit of type T 1 for which ρ C → 1. We note that a particle following the exact critical trajectory (ρ C = 1) never reaches the minimum radius, because it would have to pass through the asymptotic circular orbit (around which it circles infinitely, see Fig. 1a). However, particles starting from a position Y i slightly larger than the critical one do reach their minimal radii (though slightly larger than ρ cav ) in a finite time. Moreover, these "neighbour" trajectories reach the latter with a different phase θ: the so-formed cavity is thus (asymptotically) circular with radius ρ cav .
Some examples of trajectories are presented in Appendix A (Figs. A.2 and A.3), showing the formation of the cavity in the two regimes. We insist on the fact that it is circular in both cases, contrary to what was primarily suggested by Behar et al. (2017). As we will see in the next section, the fact that it could seem elongated in the case D > D crit comes from important contrasts in particle densities.
In the case of solar wind protons deflected around an active comet, the starting distance can be considered as infinite (D D crit ). Hence, we need to verify that the trajectories produced by this simple model of the solar wind have a well-defined limit when D → ∞. The function ρ C (Y i ) presented in Eq. (33) being monotonous whenever D > D lim (Eq. 35) and spanning all the possible values (as shown by Eq. 34), the particles can be indifferently parametrised by their initial condition Y i or by their characteristic radius ρ C . For a fixed value of ρ C , the ratio Y i /D tends to 0 when D → ∞. This means that the initial angle θ of particles coming from infinity is 0. From Eq. (26), the expression of the trajectories is thus θ(s) = s +∞ ϕ ρ(s ) ds and T 1 : ρ(s) = ρ 1 + |s| where the function ϕ(ρ) is defined in Eq. (25). This improper integral being convergent, the trajectories parametrised by their ρ C constant have indeed a welldefined limit when D → ∞. It should be noted, though, that their initial position Y i tends to −∞ (even if the ratio D/Y i tends to zero). The notion of "impact parameter" has thus no physical meaning in this problem. This information is crucial when dealing with simulations based on more realistic models of the solar wind because they are necessarily performed in a limited region of space (that is, for a finite value of D and a finite range of Y i ). For now, we already know that the size of the simulation cells (considering a regular grid) should not exceed the radius of the cavity, which is the smallest scale of the system. As we will see in the next section, our simplistic model can also be used to infer the size of the simulation box required to obtain relevant results.

The caustic
In this section, we are interested in the relative density of particles in the (x, y) plane in the permanent regime. As for the geometry of the trajectories (see Eq. 38 and text above), we should first determine if the density of particles in the (X, Y ) plane has a well-defined limit for D → ∞. For a fixed value of ρ C , we saw that Y i /D → 0 when D → ∞, that is, Y i becomes negligible compared to D. The function ρ C (Y i ) from Eq. (33) behaves thus like D exp(Y i + 1), so we can replace the uniform distribution of the particles along the Y i axis by a uniform distribution of ln(ρ C ). Since, as shown above, the geometry of the trajectory for a given ρ C has itself a unique limit (Eq. 38), the density map has also a well-defined limit for D → ∞.
The relative density of particles can be simulated by distributing points randomly along trajectories evenly sampled along Y i (for finite D) or evenly sampled along ln(ρ C ) (for infinite D). Since the particles have all the same velocity, we must use an homogeneous distribution in time τ . An illustration for infinite starting distance is given in Fig. 4. As already reported by Shaikhislamov et al. (2015) for the 1/r 3 magnetic field, a line of overdensity appears. This is a purely geometrical effect since, in the limit of this physical model, the particles do not interact which each other. This line is constituted of the points where two neighbouring trajectories cross each other. We will call it a "caustic" by analogy to light rays. For small values of D, several types of caustic appear. We will not go into details here, though, because small values of D have no physical interest.
In a general way, an overdensity region appears whenever the flux of particles is contracted, that is, when two trajectories of neighbouring initial conditions get closer to each other. This is quantified by the so-called variational equations. Let us consider a smooth function f of time t, depending on one parameter α ∈ R (which can be the initial condition f (t = 0)). At a given time t, the distance df between two curves with neighbouring values of the parameter α is at first order (see Milani & Gronchi 2010, for thorough details in the context of error propagations). Of course, the distance between the two curves vanishes if they cross, implying that ∂f /∂α = 0. In our case, the θ angle (Eq. 38) plays the part of f , the radial variable ρ plays the part of t, and the parameter ρ 0 , itself bijectively linked to the initial condition Y i , plays the part of α. The variational equation can thus be written as Using the chain rule, this partial derivative can be computed from Eq. (38), considering s as a function of ρ, itself a function of ρ 0 via ρ 1 or ρ 3 . We obtain with Article number, page 7 of 15 A&A proofs: manuscript no. EBM2 and For finite values of D, an analogous formula can be obtained from Eq. (26), containing additional terms due to the initial conditions. Figure 5 shows the general form of ∂θ/∂ρ 0 in the (ρ C , s) plane. Particles with ρ C < 1 do not produce any accumulation (they rather spread). Particles with ρ C > 1, on the contrary, arrive at a point where ∂θ/∂ρ 0 becomes zero and changes sign. This means that neighbouring trajectories cross in this point, creating an overdensity. The curve along which ∂θ/∂ρ 0 is zero can be obtained numerically using a Newton-type algorithm. Its shape in the (X, Y ) plane is presented in Fig. 6 (it should be compared to the den-sity map of Fig. 4). For particles coming from infinity, the shape of the caustic only depends on the characteristic radius r E , which acts as a scaling parameter. For solar wind protons deflected around a comet, this means that whatever the cometary activity (expressed in the k parameter), the structure formed by the proton trajectories is always exactly the same, though it is seen at a different "zoom level". This is illustrated in Fig. 6. As mentioned earlier, complex numerical simulations of the interaction of solar protons with cometary ions are always limited to finite simulation boxes. In practice, this means that solar particles, supposed unaffected yet by the comet, are launched from a finite distance D. This necessarily distorts the dynamical structures, as already pointed out by Koenders et al. (2013) for high-activity comets. In our case, by comparing the shape of the caustic for different starting distances D, our simplistic model can give an estimate of the error introduced by the finite-sized simulation boxes. This is shown in Fig. 7: simulations with a small box tend to underestimate the opening angle of the caustic. The error is thus larger at larger distances from the nucleus (but the cavity radius is unaffected as long as D D crit ).

Parameter values for a realistic comet
From Eq. (37), we know that solar wind protons are in the regime for which the radius of the circular cavity is independent of D. Switching back to dimensional quantities, it writes r cav ≈ 0.28 r E . The radius of the cavity depends thus only on r E = |k|/v, that is, on the incident velocity of the particles and on the k constant of the effective magnetic field. In particular, the cavity boundary was crossed by the Rosetta spacecraft: knowing v, its distance from the comet at time of crossing allows us to measure the k parameter (assuming that the solar wind protons did follow this simple model). Order-of-magnitude estimates can be obtained from Fig. 1 by Behar et al. (2017): the spacecraft crossed the boundary from inside to outside the cavity in December 2015, when the comet was at about 1.7 au from the Sun. The data give v = 300 km/s and r cav = 130 km at the time of crossing, resulting in a characteristic length r E ≈ 470 km and a k parameter of the order of 10 5 km 2 /s. As shown in Sect. 2.1, the value of k is proportional to the outgassing rate Q of the comet (see the companion paper by Behar et al. 2018b, for details). Actually, considering the very high velocity of the solar protons, the change of k due to the varying cometary activity can safely be modelled as an adiabatic process. Each time of an observation by Rosetta corresponds thus to a different value of k (or equivalently r E ). Still assuming that the solar wind protons did follow the dynamics described in this paper, the parameter k can be estimated at any time from the observed deflection of solar particles. Indeed, knowing the position of the spacecraft during each observation in the comet-Sunelectric frame (CSE), we just have to rescale the picture (that is, to find the unit length r E ), such that the incoming flux of particles is indeed deflected by the observed amount at this specific position. This method will be presented in detail in a forthcoming article, in which the data points will be systematically compared to theoretical values. It leads to the cavity radius being larger than 5 km when the comet is closer than 2.6 au from the Sun, and growing beyond 1500 km at perihelion. Particles come from s = ∞, they reach their minimum radii at s = 0, at which the last term of Eq. (41) diverges, and they go on with negative s. The white line shows the level curve ∂θ/∂ρ0 = 0, corresponding to the caustic (overdensity of particles). It is formed by the set of T1 trajectories (ρC > 1). See Fig. 6 for its shape in the physical plane.

Conclusion
During most of their trajectory around the Sun, comets are in a low-activity regime. When studying the dynamics of cometary and solar wind ions, this results in a gyration scale larger than the interaction region. In this situation, solar wind protons can be efficiently modelled by test particles subject to a magnetic-field-like force proportional to 1/r 2 (in a cometocentric reference frame). In this article, we provided a full characterisation of their trajectories in the plane perpendicular to this field.
As for every autonomous vector field with rotational symmetry, the system admits two conserved quantities: the kinetic energy E and a generalised angular momentum C. In our case, both of them can be turned into characteristic radii r E and r C , which entirely define the dynamics (throughout the text, we rather use the adimensional quantity ρ C = r C /r E ). There are three families of tra-jectories: two of them gather unbounded orbits (r C > r E and r C < r E ), and the other one contains quasi-periodic bounded orbits (r C < r E ). A bifurcation occurs at r C = r E , with a homoclinic orbit asymptotic to a circle of radius r C (hyperbolic equilibrium point) and two branches coming from and going to infinity. Generic analytical expressions of the trajectories (r, θ, t) are obtained, of the form θ(r) = ω E t(r) + f (r), where ω E is a constant, f (r) is an explicit function, and the time t(r) is defined by an integral.
When considering an incoming flux of particles coming from infinity on parallel trajectories and at the same velocity, a cavity is naturally created around the origin. This cavity, entirely free of particle, is circular with radius r cav ≈ 0.28 r E . Extending away from it, a curve of overdensity of particles spreads similarly to an optical caustic. This overdensity curve has no explicit expression but its shape in the plane can be computed at an arbitrary precision. The whole setting depends only on r E , which acts as a scaling parameter.
If we model the motion of solar wind protons around comet 67P by this simple dynamics, the radius r E can be calibrated from Rosetta plasma observations. From the arrival of Rosetta in the vicinity of the comet until the signal turn-off when reaching the cavity boundary, r E grew from a few kilometres up to about 470 km. This gives not only a qualitative understanding of the observed deflection of solar particles and the formation of the cavity, but also the relevant scales for the problem. In particular, if refined simulations of solar wind are used, we stress the need for simulation boxes much larger than the characteristic radius r E (to account for the caustic shape) and a grid much finer than r E (if one wants to resolve the cavity structure).
According to the results by Behar et al. (2018b), it is important to note that the capacity of this simple dynamics to account for the motion of solar wind protons around a comet is increasing with the distance to the nucleus: the farther away from the nucleus, the better the model. In other words, physical assumptions on which the physical model is based may start to crumble at the origin (the nucleus) first, leaving the modelled deflection far from the nucleus unaffected.
Comparisons to a generic magnetic field proportional to 1/r n , added in Appendix B, reveal similar features whenever n > 1. Albeit the radius of the cavity and the precise shape of the caustic are different for each n, the specific choice of 1/r 2 , if only taken on empirical grounds, would be difficult to justify. However, this law can now be retrieved from the physical modelling of solar wind protons and cometary activity, as it is presented by Behar et al. (2018b,a).
but since they become negative or complex numbers when C n < 0 (or even undefined when C n = 0), they would not have such a general physical meaning as for n = 2. The system is thus better parametrised by C n itself, or by a wisely chosen parameter γ n : This parameter, as well as the characteristic length from Eq. (B.3) has been introduced by Störmer (1907) in the particular context of n = 3. As we will see, this allows us to describe all the possible trajectories in a unified way. 5 We must now consider the cases of negative, positive, and zero values of C n . The effective potential and angular velocity as functions of ρ are represented in Fig. B.1 in the three cases. For C n > 0, the dynamics is pretty similar to the inverse-square-law field and we have the same types of trajectories. For C n 0, the only possible trajectories are of a type analogous to T 1 , but for which the radius ρ 0 would be sent to infinity. Hence, their angular velocity is always negative.
The extreme reachable radii are the positive roots of the two polynomials P ± n (ρ) = ±(n − 2) ρ n−1 + (n − 1) γ n ρ n−2 − 1 . (B.8) From Descartes' rule of sign, we obtain that P + n has exactly one positive root whatever the value of γ n (which defines ρ 1 ). On the other hand, P − n has zero or two positive roots, and exactly zero if γ n < 0. For γ n > 0, P − n has only one local maximum for ρ > 0, equal to γ n−1 n − 1. Given that P − n (0) = −1 and P − n (∞) = −∞, we deduce as expected Each curve is labelled above the graphs: the parameter ρC (Eq. 33) is drawn in yellow; the initial starting distance ρi = D 2 + Y 2 i is drawn in magenta; the extreme reachable radii ρ1, ρ2, and ρ3 (Eq. 22) are drawn in red, green, and blue. As a function of Yi, the parameter ρC is monotonous if D D lim (a, b), it crosses 1 in one additional point if D = Dcrit (c), and in two additional points if D < Dcrit (d, e, f ). The type of trajectory of the particle with initial position Yi is determined by the location of its initial distance ρi with respect to the characteristic lengths ρ1, ρ2, and ρ3: the trajectory is of type T2 when ρ1 < ρi < ρ2; T3 when ρi > ρ3; and T1 when ρi > ρ1 (with no ρ2, ρ3). We note that trajectories of type T2 are only possible for D < Dcrit (d, e, f ). that P − n has zero positive root if γ n < 1 and two if γ n 1 (which define ρ 2 and ρ 3 ). Since the polynomials are of order n − 1, these roots have necessarily an explicit expression for n 5 (Abel's impossibility theorem). For n > 5, we have no guarantee that an explicit expression of the three radii exists, but they are still well-defined from Eq. (B.8) and they can be determined numerically.
Whatever the value of n, the different types of trajectories can be easily distinguished by plotting a phase portrait of the system. In our case, the best option is to use the level curves γ n in the plane (ρ, ψ), where ψ is the angle between the position and velocity vectors. Indeed, bothρ anḋ θ can be expressed in terms of ψ, leading to the following expressions: The corresponding phase portraits are shown in Fig. B.2, showing that the dynamics in the cases n 2 are qualitatively similar. If we consider an incoming flux of particles as in Sect. 3, the γ n constant of the particles for n > 2 is and it is enough to study the case k n > 0. We note that This formula is also valid for n = 2. Then, there is a critical distance D crit < D lim below which bounded trajectories appear (as in Fig. 3). Finally the flux of incoming particles naturally creates a circular cavity similar to the case n = 2.
For D > D crit , the radius ρ cav of this cavity is also independent of D: it is equal to the ρ 1 radius (Eq. B.8) at γ n = 1. We give in Tables B.1 and B.2 the values of D crit and ρ cav for the first few n. We give also their analytical expression when we found one. Unbounded trajectories of type T3 are represented in blue. They approach a minimum distance equal to 1 but never reach it exactly (outer dashed circle, around which they can perform an arbitrary number of turns before going back). In the limiting case where ρC = 1 (black initial condition), the particle makes a infinite number of turns as ρ → 1. Bottom: Unbounded trajectories of type T1 are represented in red. They cross the characteristic radius ρ0, at which their angular velocity is inverted (see the small loops). For initial positions Yi tending to the black point, the minimal distance reached by the particle tends to ρcav = W0(exp[−1]) (inner dashed circle, see Eq. 37). If we consider an infinite number of particles, this forms a cavity with radius W0(exp[−1]) (grey disc).

(B.13)
These functions can be used directly as in Eqs. (26) and (31), with the same parametrisation s ∈ R. If we consider particles coming from infinity (D → ∞), there is a notable difference with respect to the case n > 2. Indeed, the parameter γ n of the particles (Eq. B.10) becomes directly proportional to Y i . It is thus much simpler than in Sect. 3, since Y i now keeps a clear meaning even when D is infi-  An interval of initial conditions produces bounded orbits (in green), which loops forever inside the unit circle. They approach closer to the origin than the red trajectories, producing a smaller circular cavity (grey disc, see Eq. 37). For comparison, the same circles as in Fig a) C n > 0 unit level for: γ n > 1 γ n = 1 0 < γ n < 1 angular velocityθ radial distance ρ b) C n ≤ 0 unit level for: γ n ≤ 0 radial distance ρ Phase portraits of the system in the plane (ρ, ψ) for different powers n, obtained in terms of the level curves of γn (with, in particular, γ2 = − ln ρ0). Trajectories of type T1, T2, and T3 are plotted respectively in red, green, and blue. The thick black curve represents the unit level (separatrix). Along it lie the homoclinic orbit T 2 and the two branches of the T 3 orbit, whereas the circular trajectory T is plotted in orange. The white level curve represents the zero level. For n > 2, it always remains in the ρ sin ψ < 0 side. n analytical D  − 7. The expression for n = 2 is taken from Sect. 3. nite: it becomes the impact parameter of the particles. For completeness, Fig. B.3 compares the shape of the caustics obtained for the first few values of n.