Issue 
A&A
Volume 617, September 2018



Article Number  A99  
Number of page(s)  15  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201832742  
Published online  28 September 2018 
Solar wind dynamics around a comet
The paradigmatic inversesquarelaw model
^{1}
IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, UPMC Univ. Paris 06, LAL, Université de Lille, 75014 Paris, France
email: melaine.saillenfest@obspm.fr
^{2}
LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, UPMC Univ. Paris 06, 75014 Paris, France
^{3}
Swedish Institute of Space Physics, Kiruna, Sweden
^{4}
Luleå University of Technology, Department of Computer Science, Electrical and Space Engineering, Kiruna, Sweden
Received:
31
January
2018
Accepted:
22
May
2018
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 soobtained 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, r_{E} and r_{C}, derived from first integrals. There are three different types of motion possible divided by a separatrix corresponding to r_{E} = r_{C}. 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 r_{E}, which scales the problem. A circular cavity of radius 0.28 r_{E} is created, as well as an overdensity curve (analogous to a caustic in optics). At each observation time, r_{E} 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 r_{E} and grids much finer than r_{E}.
Key words: solar wind / comets: general / magnetic fields
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. 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, hereafter 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 GiacobiniZinner: 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 selfconsistent numerical models have tackled the interaction between the solar 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 insitu 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 inversesquarelaw “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, 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 selfconsistent physical models of comet 67P are presented in companion papers (Behar et al. 2018a,b).
The planar dynamics in an inversecubelaw 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 1907, 1930; Graef & 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 inversesquare 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: Sect. 2 presents a general study of the inversesquarelaw 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 orderofmagnitude 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 inversesquare law for solar protons around comets
The analytical model introduced by Behar et al. (2018b) is built from three submodels. 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 photoionisation and electronimpact 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, newborn cometary ions, which are steadily created and lost. The second submodel is a description of the magnetic field piling up due to the local decrease in the average velocity of the electrons (as slow newborn 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 submodel is a description of the electric field, which is reduced to its main component, the socalled motional electric field,
where 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 inversesquare 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 cometSunelectric 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.
2.2. 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 , 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
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 zaxis.
From dimensionality arguments, the conservation of υ 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
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}).
2.3. Geometry of the trajectories
Fig. 1. Panel a: 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)). Panel b: angular velocity as a function of ρ. Panel 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 panel a. 
Introducing ρ_{C} (Eq. (17)) into the energy expression (Eq. (16)), we get
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 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. Figure B.2 provides details of the phase portrait of the system, where the different types of trajectories can be easily identified (panel n = 2). The extreme values of ρ reachable by the particle (i.e. ) can be obtained from (Eq. (18)) by solving the equation U(ρ) = 1. This equation can be rewritten as
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 υ = 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 get
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 . A possible parametrisation (θ(s), ρ(s)) of the three types of nonsingular trajectories is given by
where Δ = ρ_{2} − ρ_{1}. With this parametrisation, s(τ) increases with τ and ρ(s) is minimum at s = 0.
2.4. Time information
Fig. 2. Left: ratio of periods T_{C}/T_{E} as a function of the parameter ρ_{C}, computed from (Eq. (29)). Some examples of rational values, corresponding to periodic orbits, are plotted as horizontal lines. Right: same periodic orbits plotted in the physical plane. The corresponding parameter ρ_{C} was obtained by a Newton method applied to (Eq. (29)). 
By expressing the term 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 quasiperiodic, 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 tend to be equal at ρ_{C}, that is, for the circular unstable trajectory (for which 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 can be used to avoid this double definition:
with
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.
3. Application to an incoming flux of particles
Our first motivation for studying the inversesquarelaw magnetic field is the deflection of solar wind protons as a result of their interactions with a cometarytype 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 zcomponent 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 xaxis ( with υ > 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.
3.1. The cavity
Since the particles have all the same velocity υ, 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 thus ρ_{0}(Y_{i}) ≡ ρ_{C}(Y_{i})exp(−1). 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
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}). Figure 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 colourcoded 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 , 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 (23). 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 . 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 soformed 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.
Fig. 3. Type of orbit followed by a particle as a function of its initial position (D, Y_{i}). The regions are coloured according to the value of ρ_{C}(D,Y_{i}), and the level ρ_{C} = 1 is represented by the black line. The types of orbits are labelled as in (Eq. (20)). The magenta dashed line shows the initial position of the trajectory reaching the minimum radius over the whole vertical line. For D < D_{crit}, it is a trajectory of type T_{2}. For D > D_{crit}, it is a trajectory of type T_{1} but infinitely close to the ρ_{C} = 1 curve (see text). 
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 welldefined 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
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.
3.2. 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 welldefined 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 Dexp(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 welldefined 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.
Fig. 4. Top: simulated density map of the particles around the origin for D → ∞ (the particles come from the right). The inner cavity of radius ρ_{cav} = W_{0}(exp[−1]) is visible, as well as a caustic (line of overdensity). Bottom: some trajectories evenly sampled along ln(ρ_{C}) are shown. The cavity is represented by the white disc. 
Fig. 5. Value of ∂θ/∂ρ_{0} in the (ρ_{C}, s) plane for infinite D. 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 T_{1} trajectories (ρ_{C} > 1). See Fig. 6 for its shape in the physical plane. 
Fig. 6. Form of the caustic obtained numerically by finding the root of ∂θ/∂ρ_{0} (Eq. (41)). Three different zoom level are used, which can be interpreted as three levels of cometary activity. The left panel presents the same scale as Fig. 4, in which the density structure is clearly visible. 
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 socalled variational equations. Let us consider a smooth function f of time t, depending on one parameter α ∈ ℝ (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
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 Newtontype algorithm. Its shape in the (x,y) plane is presented in Fig. 6 (it should be compared to the density 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 highactivity 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 finitesized 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}).
Fig. 7. Form of the caustic for different starting distances D ⩾ D_{crit}. The distortion caused by the use of a finite value of D is shown by the difference with the D = ∞ curve (black line). 
4. 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/υ, 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 υ, 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). Orderofmagnitude 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 υ = 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^{−1}.
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 cometSunelectric 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.
5. Conclusion
During most of their trajectory around the Sun, comets are in a lowactivity 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 magneticfieldlike 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 trajectories: two of them gather unbounded orbits (r_{C} > r_{E} and r_{C} < r_{E}), and the other one contains quasiperiodic 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 turnoff 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. (2018a,b).
Acknowledgments
We thank the anonymous referee for her or his thorough reading, which led to a much better version of the article.
The axes in the plane of motion are labelled (x, z) in Behar et al. (2017, 2018b), which is the traditional convention used in 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.
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.
References
 Avrett, E. H. 1962, J. Geophys. Res., 67, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Bagdonat, T., & Motschmann, U. 2002, Earth, Moon, and Planets, 90, 305 [NASA ADS] [CrossRef] [Google Scholar]
 Behar, E., Lindkvist, J., Nilsson, H., et al. 2016, A&A, 596, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Behar, E., Nilsson, H., Alho, M., Goetz, C., & Tsurutani, B. 2017, MNRAS, 469, S396 [CrossRef] [Google Scholar]
 Behar, E., Nilsson, H., Henri, P., et al. 2018a, A&A, 616, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Behar, E., Tabone, B., Saillenfest, M., et al. 2018b, A&A, in press, DOI: 10.1051/00046361/201832736 [Google Scholar]
 Cowley, S. W. H. 1987, Philos. Trans. Roy. Soc. London Ser. A, 323, 405 [NASA ADS] [CrossRef] [Google Scholar]
 de Vogelaere, R. 1950, Canadian J. Math., 2, 440 [CrossRef] [Google Scholar]
 Deca, J., Divin, A., Henri, P., et al. 2017, Phys. Rev. Lett., 118, 205101 [NASA ADS] [CrossRef] [Google Scholar]
 Graef, C., & Kusaka, S. 1938, J. Math. Phys., 17, 43 [CrossRef] [Google Scholar]
 Grewing, M., Praderie, F., & Reinhard, R., eds. 1988, Exploration of Halley’s Comet (Berlin: Springer) [CrossRef] [Google Scholar]
 Hamlin, D. A., Karplus, R., Vik, R. C., & Watson, K. M. 1961, J. Geophys. Res., 66, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, K. C., Bagdonat, T., Motschmann, U., et al. 2007, Space Sci. Rev., 128, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, Z., Tóth, G., Gombosi, T. I., et al. 2018, MNRAS, 475, 2835 [NASA ADS] [CrossRef] [Google Scholar]
 Koenders, C., Glassmeier, K.H., Richter, I., Motschmann, U., & Rubin, M. 2013, Planet. Space Sci., 87, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Koenders, C., Goetz, C., Richter, I., Motschmann, U., & Glassmeier, K.H. 2016a, MNRAS, 462, S235 [CrossRef] [Google Scholar]
 Koenders, C., Perschke, C., Goetz, C., et al. 2016b, A&A, 594, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lifshitz, J. 1942, J. Math. Phys., 21, 94 [CrossRef] [Google Scholar]
 Milani, A., & Gronchi, G. F. 2010, Theory of Orbit Determination (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Rubin, M., Combi, M. R., Daldorff, L. K. S., et al. 2014a, ApJ, 781, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Rubin, M., Koenders, C., Altwegg, K., et al. 2014b, Icarus, 242, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Shaikhislamov, I. F., Posukh, V. G., Melekhov, A. V., et al. 2015, Plasma Phys. Control. Fusion, 57, 075007 [NASA ADS] [CrossRef] [Google Scholar]
 Störmer, C. 1907, Le Radium (Paris), 4, 2 [CrossRef] [Google Scholar]
 Störmer, C. 1930, ZAp, 1, 237 [Google Scholar]
 Williams, D. J. 1971, Adv. Geophys., 15, 137 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Complementary figures
Fig. A.1. Characteristic lengths along the line of initial position Y_{i} for different distances D. Each curve is labelled above the panels: the parameter ρ_{C} (Eq. (33)) is drawn in yellow; the initial starting distance 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 Y_{i}, the parameter ρ_{C} is monotonous if D ⩾ D_{lim} (panels a, b), it crosses 1 in one additional point if D = D_{crit} (panel c), and in two additional points if D < D_{crit} (panel d–f). The type of trajectory of the particle with initial position Y_{i} 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 T_{2} when ρ_{1} < ρ_{i} < ρ_{2}; T_{3} when ρ_{i} > ρ_{3}; and T_{1} when ρ_{i} > ρ_{1} (with no ρ_{2}, ρ_{3}). We note that trajectories of type T_{2} are only possible for D < D_{crit} (panel d–f). 
Fig. A.2 Flux of particles coming from the line X = D, with D = 1.25 > D_{crit}. The two types of possible orbits are represented separately, with the same colour code as in Fig. 3. Top: unbounded trajectories of type T_{3} 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 T_{1} are represented in red. They cross the characteristic radius ρ_{0}, at which their angular velocity is inverted (see the small loops). For initial positions Y_{i} tending to the black point, the minimal distance reached by the particle tends to ρ_{cav} = W_{0}(exp[−1]) (inner dashed circle, see (Eq. 37)). If we consider an infinite number of particles, this forms a cavity with radius W_{0}(exp[−1]) (grey disc). 
Fig. A.3 Same as Fig. A.2 but for D = 0.25 < D_{crit}. An interval of initial conditions produces bounded orbits (in green), which loop 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.2 are represented. We note that the red trajectories are still limited by the W_{0}(exp[−1]) radius (inner dashed circle). 
Appendix B: Comparison to other powers of 1/r
Starting from the pioneering work by Störmer (1907) applied to the geomagnetic field, there is a vast literature about the trajectories of particles in the equatorial plane of a magnetic dipole, for which the magnetic field is perpendicular to the plane and proportional to 1/r^{3}. The present paper reveals that the 1/r^{2} field has numerous similarities, so we propose here to compare the dynamics driven by the different powers of 1/r.
Let us introduce a positive integer n ∈ ℕ^{*}, and a physical constant k_{n} (in unit length to the power n per unit time). The equations of motion in polar coordinates for a magnetic field proportional to 1/r^{n} are
As before, these equations imply the conservation of the velocity norm and a generalised angular momentum c_{n} obtained by direct integration of (Eq. (B.2)).
Fig. B.1 Effective potential and angular velocity as a function of ρ in the three possible cases occurring for n > 2. The unit level on the vertical axis gives the intervals of ρ allowed for the particle, such that U_{n}(ρ) < 1. 
As for any planar problem with rotational symmetry, there exists an expression of the solutions (θ, t) as a function of r defined by an integral. Indeed, the conservation of c_{n} allows us to express as a function or r, which can be injected in the velocity norm. The solution is finally obtained by quadrature (see Formulas 26 and 31 obtained for n = 2).
Despite this general way of resolving the equations, the case n = 2 is special. Indeed, it is the only one for which (Eq. (B.1)) is also directly integrable, by expressing the term from the energy. This allowed us to express explicitly the time as a function of θ and r (Eq. (30)). In practice, this means the case n = 2 is the only one for which the “drift” proper frequency of all the trajectories is ω_{E} (see Sect. 2.4). For any other power of n, the quantity ω_{E} is only the frequency of the unstable circular orbit (see below); the drift frequency of the other trajectories is a function of c_{n} and thus different for all of them (Hamlin et al. 1961; Avrett 1962). This somewhat complicates the search for periodic trajectories (Graef & Kusaka 1938).
Fig. B.2 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 T_{1}, T_{2}, and T_{3} are plotted respectively in red, green, and blue. The thick black curve represents the unit level (separatrix). Along it lie the homoclinic orbit and the two branches of the 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. 
The case n = 1 must also be taken separately, because since k_{1} has the dimension of a velocity, the constant υ cannot be turned into a characteristic length analogous to r_{E} (Eq. (11)). The dynamics can though be studied by using the “effective potential” method. When considering an incoming flux of particles as in Sect. 3, we show that the launch distance d becomes the scaling parameter of the system. There is thus no limit when d → ∞. This means that the case n = 1 has no physical meaning for this setting.
For n ⩾ 2, the parameter k_{n} naturally defines a characteristic length and a characteristic frequency:
Examples for n = 2 and n = 3 can be found in (Eq. (11)) and Störmer (1930). As shown in Sect. 2, they can be used to define dimensionless coordinates ρ = r/r_{E} and dτ = ω_{E} dt. The case n = 2 is studied in detail above, so we will now suppose that n > 2. Using the dimensionless coordinates, the equations of motion and conserved quantities rewrite as
where, as before, the dot now means derivative with respect to the normalised time τ. We note that C_{n} is now the angular momentum at infinity, or equivalently, the impact parameter times the constant velocity norm. Using the “effective potential” method as in (Eq. (18)), we get
The ρ_{0} and ρ_{C}like characteristic lengths associated to this potential would be
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. (B3)) 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 inversesquarelaw 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
From Descartes‘ rule of sign, we obtain that has exactly one positive root whatever the value of γ_{n} (which defines ρ_{1}). On the other hand, has zero or two positive roots, and exactly zero if γ_{n} < 0. For γ_{n} > 0, has only one local maximum for ρ > 0, equal to . Given that and , we deduce as expected that 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 welldefined from (Eq. (B8)) 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 and can be expressed in terms of ψ, leading to the following expressions:
The corresponding phase portraits are shown in Fig. B2, 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
so all the possible values of γ_{n} are spanned by the initial positions Y_{i}, including the critical one γ_{n} = 1. The study of γ_{n} as a function of Y_{i} and D shows that the behaviour of the trajectories is qualitatively similar to what we obtained for n = 2 (Sect. 3). First of all, there is a limiting distance D_{lim} above which γ_{n} is monotonous with respect to Y_{i}. It can be written in a very general way as
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. (B8)) 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.
Fig. B.3 Comparison of the caustics obtained for the first few values of n. We recognise the sizes of the central cavities given in Table B.2. 
Critical starting distance below which bounded trajectories appear, given for the first few powers of 1/r.
Eventually, one can use the implicit solution obtained by quadrature in order to compute the shape of the caustic, as we did in Sect. 3. For n > 2, we get
These functions can be used directly as in Eqs. (26) and (31), with the same parametrisation s ∈ ℝ. 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 infinite: 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.
Radius of the cavity formed by a flux of particles coming from a distance , for the first few powers of 1/r.
All Tables
Critical starting distance below which bounded trajectories appear, given for the first few powers of 1/r.
Radius of the cavity formed by a flux of particles coming from a distance , for the first few powers of 1/r.
All Figures
Fig. 1. Panel a: 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)). Panel b: angular velocity as a function of ρ. Panel 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 panel a. 

In the text 
Fig. 2. Left: ratio of periods T_{C}/T_{E} as a function of the parameter ρ_{C}, computed from (Eq. (29)). Some examples of rational values, corresponding to periodic orbits, are plotted as horizontal lines. Right: same periodic orbits plotted in the physical plane. The corresponding parameter ρ_{C} was obtained by a Newton method applied to (Eq. (29)). 

In the text 
Fig. 3. Type of orbit followed by a particle as a function of its initial position (D, Y_{i}). The regions are coloured according to the value of ρ_{C}(D,Y_{i}), and the level ρ_{C} = 1 is represented by the black line. The types of orbits are labelled as in (Eq. (20)). The magenta dashed line shows the initial position of the trajectory reaching the minimum radius over the whole vertical line. For D < D_{crit}, it is a trajectory of type T_{2}. For D > D_{crit}, it is a trajectory of type T_{1} but infinitely close to the ρ_{C} = 1 curve (see text). 

In the text 
Fig. 4. Top: simulated density map of the particles around the origin for D → ∞ (the particles come from the right). The inner cavity of radius ρ_{cav} = W_{0}(exp[−1]) is visible, as well as a caustic (line of overdensity). Bottom: some trajectories evenly sampled along ln(ρ_{C}) are shown. The cavity is represented by the white disc. 

In the text 
Fig. 5. Value of ∂θ/∂ρ_{0} in the (ρ_{C}, s) plane for infinite D. 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 T_{1} trajectories (ρ_{C} > 1). See Fig. 6 for its shape in the physical plane. 

In the text 
Fig. 6. Form of the caustic obtained numerically by finding the root of ∂θ/∂ρ_{0} (Eq. (41)). Three different zoom level are used, which can be interpreted as three levels of cometary activity. The left panel presents the same scale as Fig. 4, in which the density structure is clearly visible. 

In the text 
Fig. 7. Form of the caustic for different starting distances D ⩾ D_{crit}. The distortion caused by the use of a finite value of D is shown by the difference with the D = ∞ curve (black line). 

In the text 
Fig. A.1. Characteristic lengths along the line of initial position Y_{i} for different distances D. Each curve is labelled above the panels: the parameter ρ_{C} (Eq. (33)) is drawn in yellow; the initial starting distance 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 Y_{i}, the parameter ρ_{C} is monotonous if D ⩾ D_{lim} (panels a, b), it crosses 1 in one additional point if D = D_{crit} (panel c), and in two additional points if D < D_{crit} (panel d–f). The type of trajectory of the particle with initial position Y_{i} 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 T_{2} when ρ_{1} < ρ_{i} < ρ_{2}; T_{3} when ρ_{i} > ρ_{3}; and T_{1} when ρ_{i} > ρ_{1} (with no ρ_{2}, ρ_{3}). We note that trajectories of type T_{2} are only possible for D < D_{crit} (panel d–f). 

In the text 
Fig. A.2 Flux of particles coming from the line X = D, with D = 1.25 > D_{crit}. The two types of possible orbits are represented separately, with the same colour code as in Fig. 3. Top: unbounded trajectories of type T_{3} 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 T_{1} are represented in red. They cross the characteristic radius ρ_{0}, at which their angular velocity is inverted (see the small loops). For initial positions Y_{i} tending to the black point, the minimal distance reached by the particle tends to ρ_{cav} = W_{0}(exp[−1]) (inner dashed circle, see (Eq. 37)). If we consider an infinite number of particles, this forms a cavity with radius W_{0}(exp[−1]) (grey disc). 

In the text 
Fig. A.3 Same as Fig. A.2 but for D = 0.25 < D_{crit}. An interval of initial conditions produces bounded orbits (in green), which loop 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.2 are represented. We note that the red trajectories are still limited by the W_{0}(exp[−1]) radius (inner dashed circle). 

In the text 
Fig. B.1 Effective potential and angular velocity as a function of ρ in the three possible cases occurring for n > 2. The unit level on the vertical axis gives the intervals of ρ allowed for the particle, such that U_{n}(ρ) < 1. 

In the text 
Fig. B.2 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 T_{1}, T_{2}, and T_{3} are plotted respectively in red, green, and blue. The thick black curve represents the unit level (separatrix). Along it lie the homoclinic orbit and the two branches of the 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. 

In the text 
Fig. B.3 Comparison of the caustics obtained for the first few values of n. We recognise the sizes of the central cavities given in Table B.2. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.