EDP Sciences
Open Access
Issue
A&A
Volume 617, September 2018
Article Number A99
Number of page(s) 15
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201832742
Published online 28 September 2018

© ESO 2018

Licence Creative Commons
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 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/r2. 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/r2 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 u0. 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 u0, 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 (H2O+) 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 nsw is the average density of solar wind protons. Injecting the expression of ncom 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 wind1. 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 rewrite2 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 rC and rE. 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 (1415) are equivalent to:

(16)

(17)

coming respectively from the energy and the generalised angular momentum (we have in particular ρC = rC/rE).

2.3. Geometry of the trajectories

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

Open with DEXTER

Introducing ρC (Eq. (17)) into the energy expression (Eq. (16)), we get

(18)

where we define

(19)

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 (Ti) the different types of trajectories, we have3

(20)

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

(21)

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

(22)

where W0 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 T1 trajectories shows thus a unique loop away from the origin, whereas T2 trajectories continuously rotate around the radius ρ0. On the contrary, T3 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 Δ = ρ2ρ1. With this parametrisation, s(τ) increases with τ and ρ(s) is minimum at s = 0.

2.4. Time information

thumbnail 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)).

Open with DEXTER

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

(29)

while the period of the overall rotation around the origin is simply TE = 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 TC/TE is a rational number. Figure 2 shows the behaviour of TC/TE 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 xi = 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/r3 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 rE and drift frequency ωE given by (Eq. (11)). We are thus able to use the normalised variables ρ = r/rE 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 Yi = yi/rE and D = d/rE, we get from (Eq. (17))

(33)

and thus ρ0(Yi) ≡ ρC(Yi)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 Yi. In the following, we suppose that k is positive4. 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(Yi) is monotonous. For D > Dlim, there is thus only one trajectory with ρC = 1 among the initial positions Yi. For D < Dlim, on the contrary, ρC(Yi) 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 < Dcrit, the local minimum of ρC(Yi) is smaller than 1. This produces two other critical trajectories with ρC = 1 (or only one in the limiting case D = Dcrit). Figure A.1 shows the behaviour of all the characteristic lengths as a function of D, where the meaning of Dlim and Dcrit 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 T2 and T3 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 Yi 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(Yi) over the T1 and T2 sets of orbits (red and green regions of Fig. 3). For D > Dcrit, 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 < Dcrit, the minimum value of ρ1 is reached in the T2 zone. The value of the minimum is

(37)

where . Interestingly, for DDcrit, 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 < Dcrit, the minimal radius is reached by a bounded orbit of type T2. 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 DDcrit, the minimal radius is reached by an unbounded orbit of type T1 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 Yi 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 > Dcrit comes from important contrasts in particle densities.

thumbnail 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).

Open with DEXTER

In the case of solar wind protons deflected around an active comet, the starting distance can be considered as infinite (DDcrit). 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(Yi) presented in (Eq. (33)) being monotonous whenever D > Dlim (Eq. (35)) and spanning all the possible values (as shown by (Eq. 34)), the particles can be indifferently parametrised by their initial condition Yi or by their characteristic radius ρC. For a fixed value of ρC, the ratio Yi/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 Yi tends to −∞ (even if the ratio D/Yi 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 Yi). 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 Yi/D → 0 when D → ∞, that is, Yi becomes negligible compared to D. The function ρC(Yi) from (Eq. (33)) behaves thus like Dexp(Yi+1), so we can replace the uniform distribution of the particles along the Yi 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 Yi (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/r3 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.

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

Open with DEXTER

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

Open with DEXTER

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

Open with DEXTER

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 ρ0, itself bijectively linked to the initial condition Yi, 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 ρ0 via ρ1 or ρ3. 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 ∂θ/∂ρ0 in the (ρC, s) plane. Particles with ρC < 1 do not produce any accumulation (they rather spread). Particles with ρC > 1, on the contrary, arrive at a point where ∂θ/∂ρ0 becomes zero and changes sign. This means that neighbouring trajectories cross in this point, creating an overdensity. The curve along which ∂θ/∂ρ0 is zero can be obtained numerically using a Newton-type algorithm. Its shape in the (x,y) plane is presented in Fig. 6 (it should be compared to the density map of Fig. 4). For particles coming from infinity, the shape of the caustic only depends on the characteristic radius rE, 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 DDcrit).

thumbnail Fig. 7.

Form of the caustic for different starting distances DDcrit. The distortion caused by the use of a finite value of D is shown by the difference with the D = ∞ curve (black line).

Open with DEXTER

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 rcav ≈ 0.28 rE. The radius of the cavity depends thus only on rE = |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 rcav = 130 km at the time of crossing, resulting in a characteristic length rE ≈ 470 km and a k parameter of the order of 105 km2 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 rE). 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 rE), 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/r2 (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 rE and rC, which entirely define the dynamics (throughout the text, we rather use the adimensional quantity ρC = rC/rE). There are three families of trajectories: two of them gather unbounded orbits (rC > rE and rC < rE), and the other one contains quasi-periodic bounded orbits (rC < rE). A bifurcation occurs at rC = rE, with a homoclinic orbit asymptotic to a circle of radius rC (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 rcav ≈ 0.28 rE. 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 rE, which acts as a scaling parameter.

If we model the motion of solar wind protons around comet 67P by this simple dynamics, the radius rE 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, rE 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 rE (to account for the caustic shape) and a grid much finer than rE (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/rn, 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/r2, 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.


1

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.

2

We get here the same equations as Graef & Kusaka (1938). This comes from a mistake in their paper: they begin with the equations of a 1/r2 field; they introduce the conserved quantities of a 1/r3 field; they write down equations mixing both types of fields, and eventually, they study the 1/r3 one for the rest of the paper. Since they deal with the motion in the equatorial plane of a magnetic dipole, 1/r3 is the correct field to use.

3

In dimensional coordinates, the three cases correspond to rC > rE, rC < rE, and rC = rE.

4

Since ρC(Yi)|k = ρC(−Yi)|k, it is enough to study the case k > 0. The case k < 0 is obtained by mirror symmetry Yi → −Yi.

5

In the case n = 2, the analogous parameter is γ2 = −lnρ0.

References

Appendix A: Complementary figures

thumbnail 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 DDlim (panels a, b), it crosses 1 in one additional point if D = Dcrit (panel c), and in two additional points if D < Dcrit (panel df). 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 df).

Open with DEXTER

thumbnail 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).

Open with DEXTER

thumbnail 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).

Open with DEXTER

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/r3. The present paper reveals that the 1/r2 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 kn (in unit length to the power n per unit time). The equations of motion in polar coordinates for a magnetic field proportional to 1/rn are

(B.1)

(B.2)

As before, these equations imply the conservation of the velocity norm and a generalised angular momentum cn obtained by direct integration of (Eq. (B.2)).

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

Open with DEXTER

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 cn 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 cn and thus different for all of them (Hamlin et al. 1961; Avrett 1962). This somewhat complicates the search for periodic trajectories (Graef & Kusaka 1938).

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

Open with DEXTER

The case n = 1 must also be taken separately, because since k1 has the dimension of a velocity, the constant υ cannot be turned into a characteristic length analogous to rE (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 kn 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/rE 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 Cn 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 ρ0 and ρC-like characteristic lengths associated to this potential would be

(B.6)

but since they become negative or complex numbers when Cn < 0 (or even undefined when Cn = 0), they would not have such a general physical meaning as for n = 2. The system is thus better parametrised by Cn 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 way5.

We must now consider the cases of negative, positive, and zero values of Cn. The effective potential and angular velocity as functions of ρ are represented in Fig. B.1 in the three cases. For Cn > 0, the dynamics is pretty similar to the inverse-square-law field and we have the same types of trajectories. For Cn ⩽ 0, the only possible trajectories are of a type analogous to T1, 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

(B.8)

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 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 kn > 0. We note that

(B.11)

so all the possible values of γn are spanned by the initial positions Yi, including the critical one γn = 1. The study of γn as a function of Yi 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 Dlim above which γn is monotonous with respect to Yi. 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 Dcrit < Dlim 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 > Dcrit, 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 Dcrit and ρcav for the first few n. We give also their analytical expression when we found one.

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

Open with DEXTER

Table B.1.

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

(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 Yi. It is thus much simpler than in Sect. 3, since Yi 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.

Table B.2.

Radius of the cavity formed by a flux of particles coming from a distance , for the first few powers of 1/r.

All Tables

Table B.1.

Critical starting distance below which bounded trajectories appear, given for the first few powers of 1/r.

Table B.2.

Radius of the cavity formed by a flux of particles coming from a distance , for the first few powers of 1/r.

All Figures

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

Open with DEXTER
In the text
thumbnail 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)).

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 7.

Form of the caustic for different starting distances DDcrit. The distortion caused by the use of a finite value of D is shown by the difference with the D = ∞ curve (black line).

Open with DEXTER
In the text
thumbnail 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 DDlim (panels a, b), it crosses 1 in one additional point if D = Dcrit (panel c), and in two additional points if D < Dcrit (panel df). 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 df).

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.