© ESO 2018

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 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 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 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, 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². 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² 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. 2018a,b).

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 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 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 sub-model 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₀. 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₀, 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

(1)

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,

(2)

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₂O⁺) and solar wind (H⁺) ions. Putting all things together, the force applied to the solar protons is

(3)

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¹. 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

(4)

From the classical Lorentz force, the equations of motion are , that is,

(5)

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² as

(6)

(7)

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

(8)

equal to the norm of the velocity projected in the (x, y) plane. Moreover, (Eq. (7)) is directly integrable:

(9)

where the arbitrary constant r_⋆ is added for dimensionality reasons. This leads to a second first integral,

(10)

It can be thought as the conservation of a generalised angular momentum, coming from the symmetry of rotation around the z-axis.

From dimensionality arguments, the conservation of υ makes a characteristic length and a characteristic frequency of the system appear:

(11)

In the same way, the generalised angular momentum c can be turned into the characteristic length

(12)

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

(13)

We note that if k > 0, the direction of time is reversed. In the new coordinates, the equations of motion (6–7) become

(14)

(15)

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:

(16)

(17)

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

(18)

where we define

(19)

We will see below that both ρ₀ 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 ρ₀ 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³

(20)

We note that ρ₂ and ρ₃ 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

(21)

The extreme values of ρ reachable by the particle in the different cases are thus

(22)

where W₀ and W_-1 are the Lambert functions. In accordance with (Eq. (20)), ρ₂ and ρ₃ are only defined if ρ₀ 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 ρ = ρ₀ 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 ρ = ρ₀ and changes sign. The inner part of T₁ trajectories shows thus a unique loop away from the origin, whereas T₂ trajectories continuously rotate around the radius ρ₀. On the contrary, T₃ 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

(23)

which gives

(24)

with

(25)

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 non-singular trajectories is given by

(26)

where Δ = ρ₂ − ρ₁. With this parametrisation, s(τ) increases with τ and ρ(s) is minimum at s = 0.

2.4. Time information

Fig. 2. Left: ratio of periods TC/TE 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

(27)

which can be directly integrated to give

(28)

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 ρ₀ radius). Since vanishes at the extreme values of ρ (Eq. (22)), the period T_C of the loops is

(29)

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 ρ₀. 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

(30)

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:

(31)

with

(32)

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 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 ( 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³ 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))

(33)

and thus ρ₀(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⁴. First of all, we note that

(34)

The particles have thus all the possible values of ρ_C, including the critical one ρ_C = 1 (Eq. (20)). Let us write

(35)

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,

(36)

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 colour-coded in Fig. 3 as a function of their initial position. Trajectories of types T₂ and T₃ can be distinguished by considering their initial radius , which should be smaller or larger than 1, respectively (since ρ₂ < 1 and ρ₃ > 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 ρ₁(Y_i) over the T₁ and T₂ 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 ρ₁ is reached in the T₂ zone. The value of the minimum is

(37)

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₂. 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₁ 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.

Fig. 3. Type of orbit followed by a particle as a function of its initial position (D, Yi). The regions are coloured according to the value of ρC(D,Yi), 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 < Dcrit, it is a trajectory of type T2. For D > Dcrit, it is a trajectory of type T1 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 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

(38)

where the function φ(ρ) is defined in (Eq. (25)). This improper integral being convergent, the trajectories parametrised by their ρ_C constant have indeed a well-defined 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 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 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 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³ 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 = W0(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 T1 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 so-called 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

(39)

(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 ρ₀, itself bijectively linked to the initial condition Y_i, plays the part of α. The variational equation can thus be written as

(40)

Using the chain rule, this partial derivative can be computed from (Eq. (38), considering s as a function of ρ, itself a function of ρ₀ via ρ₁ or ρ₃. We obtain

(41)

with

(42)

and

(43)

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 ∂θ/∂ρ₀ 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 ∂θ/∂ρ₀ becomes zero and changes sign. This means that neighbouring trajectories cross in this point, creating an overdensity. The curve along which ∂θ/∂ρ₀ 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 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 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).

Fig. 7. Form of the caustic for different starting distances D ⩾ Dcrit. 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). 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 υ = 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⁵ km² 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-Sun-electric 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 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² (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 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ⁿ, 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², 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.

References

Appendix A: Complementary figures

Fig. A.1. Characteristic lengths along the line of initial position Yi 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 Yi, the parameter ρC is monotonous if D ⩾ Dlim (panels a, b), it crosses 1 in one additional point if D = Dcrit (panel c), and in two additional points if D < Dcrit (panel d–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 (panel d–f).

Fig. A.2 Flux of particles coming from the line X = D, with D = 1.25 > Dcrit. The two types of possible orbits are represented separately, with the same colour code as in Fig. 3. Top: 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).

Fig. A.3 Same as Fig. A.2 but for D = 0.25 < Dcrit. 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 W0(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³. The present paper reveals that the 1/r² 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ⁿ are

(B.1)

(B.2)

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 Un(ρ) < 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 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 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₁ 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:

(B.3)

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

(B.4)

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

(B.5)

The ρ₀ and ρ_C-like characteristic lengths associated to this potential would be

(B.6)

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:

(B.7)

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⁵.

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₁, but for which the radius ρ₀ would be sent to infinity. Hence, their angular velocity is always negative.

The extreme reachable radii are the positive roots of the two polynomials

(B.8)

From Descartes‘ rule of sign, we obtain that has exactly one positive root whatever the value of γ_n (which defines ρ₁). 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 ρ₂ and ρ₃). 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. (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:

(B.9)

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

(B.10)

and it is enough to study the case k_n > 0. We note that

(B.11)

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

(B.12)

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 ρ₁ 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.

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

(B.13)

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.

All Tables

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 TC/TE 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, Yi). The regions are coloured according to the value of ρC(D,Yi), 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 < Dcrit, it is a trajectory of type T2. For D > Dcrit, it is a trajectory of type T1 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 = W0(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 T1 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 ⩾ Dcrit. 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 Yi 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 Yi, the parameter ρC is monotonous if D ⩾ Dlim (panels a, b), it crosses 1 in one additional point if D = Dcrit (panel c), and in two additional points if D < Dcrit (panel d–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 (panel d–f).

In the text

Fig. A.2 Flux of particles coming from the line X = D, with D = 1.25 > Dcrit. The two types of possible orbits are represented separately, with the same colour code as in Fig. 3. Top: 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).

In the text

Fig. A.3 Same as Fig. A.2 but for D = 0.25 < Dcrit. 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 W0(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 Un(ρ) < 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 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 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