© ESO 2022

1 Introduction

Comets are a formidable laboratory for plasma experiments. As the surface of the nucleus is heated by solar radiation, the ices sublimate, turning into a gas that flows away from the nucleus at several hundred m s⁻¹. This gas, in turn, is ionised by the extreme ultraviolet solar radiation and accelerated solar wind electrons such that an ionosphere forms directly from the surface. This cloud of cometary ions and electrons will then interact with the ambient interplanetary plasma, mainly made of protons and electrons, carrying a convective electric field. However, as comets have very elliptical or hyperbolic trajectories, they may go through different stages: as they get closer to the Sun, the outgassing activity increases as well as photoionisation, and therefore the cometary ionosphere becomes denser. Noteworthily, within the cometary ionosphere, a particular region forms around the nucleus that is characterised by an extremely low ambient magnetic field, lower than 1 nT, named the diamagnetic cavity.

A cometary diamagnetic cavity was observed for the first time on 13 and 14 March 1986 (Neubauer et al. 1986) during the Giotto flyby (Reinhard 1986). As Giotto got closer to the nucleus, the magnetic field strength slowly increased and then abruptly dropped to almost zero for 2 min at around 4600 km from the nucleus (Neubauer 1987). The surface that encased the cavity was named the ‘contact surface’. In the literature, this boundary was also referred to as the ‘ionopause’ (Ma et al. 2008). Cravens (1986) ruled out the use of the latter term as in the case of Venus, this was defined as the boundary where the magnetic pressure balances the thermal plasma pressure which might not be the case at 1P/Halley. The contact surface is also commonly referred to as the cavity boundary or simply the boundary (e.g. Goetz et al. 2016b,a; Gunell et al. 2017). A more appropriate name in light of the results presented here and used by Israelevich et al. (2003) is the term diamagnetic cavity boundary layer (DCBL). Throughout this paper, we use the latter. The origin of the diamagnetic cavity and the balance at play at the DCBL are still debated. At the time, a possible explanation brought by Cravens (1987) and Ip & Axford (1987) was that the magnetic pressure was counterbalanced by the ion neutral drag, that is to say, the force applied by the neutral on the plasma (mainly ions) by means of ion-neutral collisions as they move at different velocities. However, the calculation relied on assumptions about the geometry of the cavity which is difficult to gauge with single spacecraft observations: (i) the ion speed is close to zero at the boundary and outwards, or (ii) the cometary plasma is at photochemical equilibrium. To be applicable, this balance at the boundary requires that (Cravens 1986): $$-\nabla \left(\frac{B^2}{2{\mu }_0}\right).\left(V-U_n\right)>0,$$(1)

where V stands for the mean plasma velocity and U_n for that of the neutrals, showing that the equilibrium may not hold if ions travel faster than the neutrals along the gradient direction of the magnetic pressure.

Almost 30 years later, from 2014 to 2016, the Rosetta mission (Glassmeier et al. 2007a) orbiting around Comet 67P/Churyumov-Gerasimenko (Churyumov & Gerasimenko 1972) provided a unique opportunity to explore the plasma environment of a comet over an extended period of time and follow its evolution. The findings were a heliocentric distance from 1.24 au (perihelion) to 3.8 au (end of mission) and an outgassing rate from 10²⁵ s⁻¹ to slightly below 10²⁹ s⁻¹ (Simon Wedlund et al. 2019, 2020). Thanks to its onboard magnetometer (RPC-MAG, Glassmeier et al. 2007b), a diamagnetic cavity was also observed but this time over an extended period, as Rosetta, unlike Giotto, was not a flyby mission. In the middle of its escort phase, Rosetta crossed the DCBL multiple times from April 2015 to February 2016 (Goetz et al. 2016a; Götz 2019), and went in and out of the cavity more than 600 times for short periods of time (from several tens of seconds to tens of minutes). However, Rosetta was an extremely slow spacecraft, not exceeding a few m s⁻¹ except during manoeuvres and excursions, meaning that the DCBL was a moving boundary likely behaving like ebb and flow passing over the immobile spacecraft: the origin of this behaviour is unknown and the behaviour itself is not yet clearly understood. Finally, crossings occurred at much closer distances from the nucleus than at 1P, namely at between 50 and 400 km, clearly linked to different outgassing rates and heliocentric distances. Nevertheless, these crossings occurred farther out than anticipated by models (Rubin et al. 2012; Koenders et al. 2015).

Plasma modelling and experiments are required in order to understand this boundary and the birth of the diamagnetic cavity. Three different modelling approaches exist in the case of comets: magnetohydrodynamics (MHD), hybrid, and fully kinetic particle-in-cell (PIC). MHD models have the ability to model large-scale structures – especially at large outgassing rates (≳10²⁷ s⁻¹) – and the interaction between the cometary ionosphere and the Solar wind, for instance. However, ions and electrons are both treated as fluids, and ion gyroradius scales cannot be resolved in MHD. Additional assumptions are made regarding the electron pressure, the generalised Ohm’s law, and the treatment of Maxwell’s equations (e.g. assuming plasma quasi-neutrality and neglecting the displacement current in the low-frequency limit). Although MHD models in general do exhibit a diamagnetic cavity around the comet nucleus and seem to agree with single-spacecraft observations (e.g. see Rubin et al. 2014; Huang et al. 2016, 2018) with similar sizes, the origin and equilibrium at play at the DCBL are still not understood. For example, Maxwell’s equations are not self-consistently solved and the inclusion of the Hall term in the Maxwell-Faraday (induction) equation may lead to very different geometries or sizes of the diamagnetic cavity (Huang et al. 2018). Indeed, this determines whether the magnetic field is frozen-in with the ions or the electrons. Moreover, some of these models might ignore ion-neutral chemistry whereas in situ observations have shown that the latter is significant in particular near perihelion (Beth et al. 2020) when diamagnetic cavity crossings are more likely to be observed. Finally, the spatial discretisation of the domain introduces a numerical resistivity and therefore a numerical diffusion of the magnetic field which prevents reliable studies in a quasi-collisionless case.

Hybrid models still treat electrons as a fluid while ions are treated as clouds of macro-particles with their own probabilistic weight. The ion velocity is updated and pushed in time using the fundamental relation of the dynamics, which may include a Langevin term to reproduce the ion-neutral drag (e.g. Puhl-Quinn & Cravens 1995). This corresponds to a continuous treatment of the collisions. As the number of neutral molecules is much larger than that of ions, during each time step, there are enough collisions with the neutrals that the total transfer of momentum can be modelled as an average friction force. This approximation can hold only if the cell size is much larger than the mean free path of the ions. There are also hybrid models that treat the collisions with a probabilistic approach (e.g. Koenders et al. 2015; Simon Wedlund et al. 2017; Alho et al. 2019).

It is only recently that full-PIC simulations have been carried out at comets in a very limited case for an outgassing rate of 10²⁵ s⁻¹ where collisions may be negligible. Indeed, in contrast to hybrid and MHD models, PIC models are designed to primarily simulate collisionless plasmas (driven by the Klimontovich equation) and solve electromagnetic fields self-consistently. Deca et al. (2017) performed a four-species (cometary ions and electrons, solar-wind protons and electrons) simulation of 67P for the conditions met at ~3 au with iPIC3D (Markidis et al. 2010), including an implicit numerical scheme (Mason 1981; Brackbill & Forslund 1982). Indeed, implicit schemes alleviate constraints present in explicit PIC models: observing the Courant-Friedrich-Levy (CFL, Courant et al. 1928) condition (cΔt ≤ Δx) and preventing the so-called grid instability (Δx ≲ λ_e,De, Hockney & Eastwood 1988; Birdsall & Langdon 2004). Implicit PIC models allow us to model larger spatial scales with larger time-steps (Deca et al. 2017). However, the fluctuations of the electric field might be damped as the fast motion of the electrons is not resolved as it is in iPIC3D (Markidis et al. 2010). Nevertheless, these models are valuable in that they provide details of the complex interaction between the Solar wind and the cometary ionosphere in the collisionless regime. In particular, they contribute to the understanding of the energisation of the solar wind electrons that dive towards the comet nucleus (see Galand et al. 2020). Although PIC simulations allow us to access the physics at smaller scales, the main drawback is that only a few include collisions by means of a Monte-Carlo approach (so-called PIC-MMC). This should be kept in mind as observations point out that plasma boundaries are correlated with the ion exobase (Mandt et al. 2016) or the electron exobase regarding the DCBL (Henri et al. 2017).

Prior to the Giotto flyby, active experiments in space were performed in order to simulate how a plasma cloud expands into a magnetised environment, providing insight into how the cometary plasma interacts with the Solar wind. These experiments are known as AMPTE artificial comet experiments (Active Magnetospheric Particle Tracer Explorers Valenzuela et al. 1986) and originate from an idea put forward by Biermann et al. (1961). They consist in the release of primarily a barium cloud into the Solar wind with the outcome being probed in situ during the release. Lühr et al. (1986b) and Haerendel et al. (1986) focused on respectively the magnetic field observations and the plasma dynamics. A similar experiment was set later with lithium (Lühr et al. 1986a). During both releases, a diamagnetic cavity was formed around the cloud. However, it could not be maintained as ions were not replenished as they are at comets through ionisation of the continuously outward-expanding neutrals.

In recent years, other experiments have been performed to simulate diamagnetic cavities in the laboratory (e.g. Bonde et al. 2015, 2018). In these experiments, a plasma was produced by a laser pulse hitting a target, and a diamagnetic cavity formed as the plasma expanded in the surrounding magnetic field.

Although there are many approaches and attempts to tackle the origin of the diamagnetic cavity, little is known and its origin is still debated. However, no fully kinetic simulation had ever been employed to investigate its formation until now. This paper is the first reported attempt to use a PIC simulation in such a case. This method allows us to look at physical phenomena down to the electron scale and the electron velocity distribution function, that is, the behaviour of electrons at the kinetic scales through this transition region. The paper is organised as follows. In Sect. 2, we describe the numerical model used for the simulations as well as the setup. In Sect. 3, we present the results from the simulations such as electromagnetic fields, thermodynamics variables, and distribution functions, followed by a discussion focusing on the properties of the electrons through the DCBL in Sect. 4. Finally, we summarise our findings and propose future investigations in Sect. 5.

2 Method

2.1 Formalism

In this work, we simulate the DCBL using the PIC method. In this section, we review the basics of the method and its relation to the kinetic and fluid theory of plasmas that we use in our analysis. PIC simulations represent the plasma by a finite number of macro particles (so-called Klimontovich-Dupree representation) and solve the Klimontovich equation (Klimontovich 1958; Dupree 1963). In this representation, the velocity distribution function f_s of the species s is discrete and given by: $$f_s\left(r,v,t\right)={\displaystyle \sum _{j=1}^N{W}_j{\delta }_r\left(r-r_j\left(t\right)\right){\delta }_0\left(v-v_j\left(t\right)\right),}$$(2)

where N is the number of macroparticles in the simulation, r_j(t) is the position of the macroparticle j at time t, v_j(t) is its velocity, and W_j is its associated weight. However, numerically, due to space discretisation, the Dirac function in space δ(r − r_i(t)) should be replaced by a shape function ‘S’ which deposits the charge and the current of the macro-particle onto spatial grid points.

For plasmas, in the limit when the plasma parameter $\Lambda =4\pi n_e{\lambda }_D^3$ is large enough (i.e. a weakly coupled and uncorrelated plasma, with no pair interactions between particles), the distribution function becomes the solution of the Vlasov equation as the sources and losses have been ignored, which gives a continuous description of the plasma in space and velocity such that: $${\partial }_t{f}_s=-v\cdot {\nabla }_r{f}_s-\frac{q_s}{m_s}\left(E+v\times B\right)\cdot {\nabla }_v{f}_s,$$(3)

with the electric E and magnetic B fields solutions of Maxwell’s equations: $$\begin{array}{lll}{\nabla }_r\cdot E\hfill & =\hfill & q\frac{n_i-n_e}{{\epsilon }_0}\hfill \\ {\nabla }_r\cdot B\hfill & =\hfill & 0\hfill \\ {\partial }_{t}B\hfill & =\hfill & -{\nabla }_r\times E\hfill \\ {\partial }_{t}E\hfill & =\hfill & c^2{\nabla }_r\times B-{\mu }_0{c}^2\underset{J_i}{\underbrace{(q{n}_i{V}_i}}\underset{J_i}{\underbrace{-q{n}_e{V}_e)}},\hfill \end{array}$$(4)

where n_i (n_e) stands for the total ion (electron) number density, V_i (V_e) stands for the mean ion (electron) bulk velocity, and J_i (J_e). stands for the mean ion (electron) current. Each quantity is a function of space and time.

Due to the finite number of particles, it is not possible to obtain a continuous description of f_s in phase space. Nevertheless, f_s has moments, solutions of the different equations in MHD, and hybrid simulations in the limit of a large plasma parameter (Λ ≫ 1) to limit the correlation between particles. Fluid equations are not solved in a PIC simulation, but fluid quantities derived from the moments of the velocity distribution function (VDF) can be of interest when analysing the results, as we show in Sect. 3. For instance, the continuity equation for species s in the absence of source and loss is as follows: $${\partial }_t{n}_s=-\nabla \cdot \left(n_s{V}_s\right),$$(5)

and the momentum equation is: $$\begin{array}{lll}{m}_s{\partial }_t\left(n{}_s{V}_s\right)\hfill & =\hfill & -\nabla \cdot \stackrel{\text{dynamic\hspace{0.17em}pressure\hspace{0.17em}tensor\hspace{0.17em}}{\text{P}}_{\text{s,dyn}}}{\overbrace{\left(m_s{n}_s{V}_s\otimes V_s\right)}}\hfill \\ \hfill & \hfill & -\nabla \cdot \underset{\text{thermal\hspace{0.17em}pressure\hspace{0.17em}tensor}}{\underbrace{{\text{P}}_{\text{s,th}}}}+\underset{\text{Lorentz\hspace{0.17em}force}}{\underbrace{q_s{n}_s\left(E+V_s\times B\right),}}\hfill \end{array}$$(6)

which can both be written in conservation form. By adding both momentum equations from the ions (index i) and the electrons (index e), we end up with: $$\begin{array}{l}{m}_i{\partial }_t\left(n_i{V}_i\right)+m_e{\partial }_t\left(n_e{V}_e\right)+{\epsilon }_0{\partial }_t\left(E\times B\right)=\hfill \\ -\nabla \cdot {\text{P}}_{\text{i,dyn}}-\nabla \cdot {\text{P}}_{\text{e,dyn}}-\nabla \cdot {\text{P}}_{\text{i,th}}-\nabla \cdot {\text{P}}_{\text{e,th}}+\nabla \cdot \sigma ,\hfill \end{array}$$(7)

where $$\sigma ={\epsilon }_{0}E\otimes E+\frac{B\otimes B}{{\mu }_0}-\left(\frac{{\epsilon }_0{E}^2}{2}+\frac{B^2}{2{\mu }_0}\right){\text{I}}_3,$$(8)

is the Maxwell stress tensor (a matrix), P_i,th (P_e,th) is the ion (electron) thermal pressure tensor, and I₃ is the identity matrix. In the MHD limit, the energy stored in the electric field is negligible compared to the magnetic energy (i.e. c²E² ≪ B²). Finally, when the electrons are assumed massless, the electric field can be derived explicitly through the electron momentum equation, the so-called generalised Ohm’s law given by: $$E\approx -\frac{\nabla \cdot {\text{P}}_{\text{e,th}}}{q{n}_e}-V_e\times B.$$(9)

Equation (9) will be scrutinised and verified in the context of our simulation in Sect. 3.3. As the simulation is performed in 1D along x, Maxwell’s equations can be simplified as follows: $$\begin{array}{ll}{\partial }_x{E}_x=\frac{q\left(n_i-n_e\right)}{{\epsilon }_0}\hfill & {\partial }_x{B}_x=0\hfill \\ {\partial }_t{B}_y=+{\partial }_x{E}_z\hfill & {\partial }_t{B}_z=-{\partial }_x{E}_y\hfill \\ {\partial }_t{E}_x=-{\mu }_0{c}^2{J}_x\hfill & \hfill \\ {\partial }_t{E}_y=-{\mu }_0{c}^2{J}_y-c^2{\partial }_x{B}_z\hfill & {\partial }_x{E}_z=-{\mu }_0{c}^2{J}_z+c^2{\partial }_x{B}_y.\hfill \end{array}$$(10)

This set of equations, alongside the generalised Ohm’s law, is important to keep in mind in order to interpret our results presented in Sect. 3. Finally, a word of caution: as sources and losses have been ignored in particular in the continuity equation of both species, it is impossible to reach equilibrium, as this would mean that n_sV_s · e_x = constant, as the simulation is performed in one spatial dimension. As the initial number density is imposed and decreasing as a function of x, steady state would require V_s · e_x ∝ 1/n_s. In addition, the simulation time cannot be too long because, in the real world, ions and electrons are continuously produced through photoionisation and electron-impact, and therefore the initial reservoir of ions and electrons would be replaced and refurnished with these newborn ions and electrons, which is not taken into account here.

2.2 Setup and initial conditions

The PIC simulations were carried out with SMILEI, an explicit and Cartesian high-performance open-source code designed to simulate various plasma physics situations, from astrophysics to relativistic laser-plasma interactions (Derouillat et al. 2018). The simulation is setup in 1D3V configuration: thermodynamic quantities only depend on the spatial direction x while particles in the velocity phase space may still evolve in the three directions in velocity space. Firstly, as an explicit scheme, the simulation time-step should be small enough to prevent light waves from propagating more than one cell at any given time (the so-called Courant-Friedrichs-Lewy condition), that is, cΔt ≤ Δx for a strict stability. However, in practice, it may be necessary to restrain them (cΔt < Δx). Secondly, in order to prevent numerical heating, the grid resolution Δx should resolve the electron Debye length ${\lambda }_{e,\text{De}}=\sqrt{\frac{{\epsilon }_0{k}_{\text{B}}{T}_e}{q^2{n}_e}}\left(\text{Δ}x<{\lambda }_{e,\text{De}}\right)$. The SMILEI unit of length is the electron skin depth L_e,sd = c/ω_pe so that ${\lambda }_{e,\text{De}}/L_{e,\text{sd}}=\sqrt{\frac{k_{\text{B}}{T}_e}{m_e{c}^2}}$. As our investigations are in the frame of classical physics and ‘cold’ plasmas, we need k_BT_e ≪ m_ec². Nevertheless, we cannot use realistic electron temperatures of the order of a few tens of eVs. Indeed, we have the hierarchical relation: $$\stackrel{\text{no\hspace{0.17em}instability}}{\overbrace{\underset{\text{CFL}}{\underbrace{c\text{Δ}t\le \text{Δ}x}}\le \underset{k_{\text{B}}{T}_e\ll m_e{c}^2}{\underbrace{{\lambda }_{e,\text{De}}\ll L_{e,sd}}}}}\ll N_{\text{cell}}\Delta x,$$(11)

and $${\omega }_{\text{pe}}\text{Δ}t\le \sqrt{\frac{k_{\text{B}}{T}_e}{m_e{c}^2}}.$$

The initial temperature of the plasma constrains the time-step and the total run time of the simulations. Therefore, trade-offs must be made. For our simulations, we initialise the electrons with a thermal energy of k_BT_e = 0.01m_ec². This value appears relatively large compared with the reality (around k_BT_e ≈ 10 eV) but according to Eq. (11), applying a realistic electron temperature is unfeasible because of limited computational resources. Indeed, reducing the temperature would require reducing both Δx and drastically increasing the runtime for simulating the same spatial and temporal domain. In addition, as discussed in Sect. 4, fluctuations of the electric field in the unmagnetised field would be theoretically larger for lower Debye lengths if the number of particles per cell were kept constant between simulations. Moreover, it should be noted that because the initial number density is not uniform along x, ω_pe and are not either. However, the inputs are defined with respect to the highest plasma number density n₀ = n_i(0, 0) meaning that ω_pe(x, 0) < ω_pe(0, 0) and λ_e,De(0, 0) < λ_e,De(x, 0) ensuring stability and no numerical heating. For our simulations, we chose n₀ = 10⁹ m⁻³ which is typical of ion densities observed around perihelion during diamagnetic cavity crossings (Henri et al. 2017; Hajra et al. 2018). Here, n₀ corresponds to the plasma number density on the left side of the simulation box (x = 0). Nevertheless, at the DCBL (x = L_box/2), the number density is close to n₀ exp(−2) ≈ 0.13n₀, which is one order of magnitude less than observations.

To perform a fully kinetic simulation, two additional ingredients are required: a spatial grid and a particle pusher. We used the default grid in SMILEI, the Yee grid (Yee 1966), which is the only one available for 1D geometry. Regarding the particle pusher, several options are available in SMILEI. All have their pros and cons. To be accurate, a pusher should be ultimately symplectic. At the moment, none of the pushers provided by Boris, Vay, or Higuera and Cary – that are available in SMILEI – (Boris & Shanny 1970; Vay 2008; Higuera & Cary 2017) have this property. Each scheme may introduce errors at different levels (Ripperda et al. 2018). For this simulation, we chose the pusher of Higuera and Cary.

Finally, the simulation needs appropriate initial and boundary conditions both for the species and electromagnetic fields. For the ions and electrons, we start from an initial profile decreasing exponentially along x. Although at comets, ion and electron number densities should decrease as 1/r where r is the cometocentric distance (Gombosi 2015; Beth et al. 2019), this holds only for a spherical symmetry. Here, the simulation is performed in Cartesian geometry along x and therefore it is not perfectly representative of a comet. The choice of an exponential profile has one purpose and benefit: as we initialise with a constant electron temperature T_e, in regions where the magnetic field is constant, the electric field is purely ambipolar such that E_x ≈ (k_BT_e/q)∂_xn_e/n_e. If the electron number density n_e follows an exponential law, E_x is almost constant. This is very helpful because it can be easily estimated from inputs and should be constant through the simulation box. At the boundaries, we do not inject ions. However, for the lower bound (x = 0, closer to the comet), ions and electrons are reflected and thermalised while for the upper bound they are purely lost. Regarding electromagnetic forces, the best choice seems to be Silver-Müller boundary conditions (Barucq & Hanouzet 1997) to prevent trapping of waves (especially light waves) within the box.

3 Results

We set up a 1D3V PIC simulation of the DCBL. The spatial dimension is denoted by x, and the x axis crosses the DCBL. The left-hand side (x = 0) is located inside the diamagnetic cavity, and the right-hand side in plasma surrounding the cavity. The simulated region is initialised with a plasma whose density decreases with increasing x and which is unmagnetised for low x values and magnetised for high ones. The time-step is 0.081${\omega }_{\text{e},0}^{-1}$ ≈ 4.54 × 10⁻⁸ s where ω_pe,0 corresponds to the highest plasma frequency met at the start of the simulation (i.e. at x = 0). We run the simulation over more than 750 000 timesteps which corresponds to ~0.035 s of real time. Values of the electromagnetic fields (E, B) and of the thermodynamic quantities (number density n, current J, thermal pressure tensor P of both species, ions and electrons) have been recorded every 100 time-steps as the PIC simulation uses a substantial amount of memory. The initialisation of the simulation is summarised in Table 1 and Fig. 1. One word should be said regarding the magnetic field profile. Although the shape was enforced, several tests were made on the most suitable H_B. Interestingly, steeper values (i.e. ≲400L_De,0) generated large amplitude waves in the electromagnetic fields towards the unmagnetised region. In addition, if too steep, the magnetic profile relaxes quickly towards a gentler profile with a scale height close to H_B. In that respect, we initialise the simulation directly to this stabler value. Moreover, around the DCBL as indicated in Fig. 1, H_B corresponds to 140 local L_De (the Debye length increases with x as n_e decreases and T_e is constant at the onset) and 14c/ω_p,e. This is of the order of the value suggested by Grad (1961). Although unknown at that time, Grad (1961) analytically explored a DCBL-like configuration of the plasma and magnetic field and estimated the ‘thinnest’ possible magnetic field profile. He found that the minimal width should be around ~8c/ω_p,e.

Fig. 1 Initial setup and profiles for ions, electrons, and magnetic field in our simulation box. Ion temperature is not indicated as it is set to 0.

3.1 Energy density of electromagnetic fields E and B

Figure 2 shows the evolution over time of the energy stored in the different components of the electric field during the simulation in SMILEI units (here n₀m_ec²). We also over-plotted the different speeds of interest in the case where waves are present: the speed of light (c), the thermal speed of the electrons $\left(\sqrt{k_{\text{B}}{T}_e/m_e}\approx 0.1c\right)$, and the ion acoustic speed $\left(\sqrt{k_{\text{B}}{T}_e/m_i}\approx 0.001c\right)$. The ion initial speed (10⁻⁴c) is not displayed as it would be almost vertical. Distinct features may be seen in the different components. In the top panel $\left({\epsilon }_0{E}_x^2\right)$, there is a layer of disturbance on the left side (x ≈ 0) associated with box boundary effects. For the left side, we chose a boundary that reflects ions and electrons with the initial speed and temperature set at the beginning of the simulation. Interestingly, the boundary moves at different speeds during the simulation. The boundary appears around 0.003 s, moving at around the electron thermal speed, and then at ~0.005 s the boundary slows down and moves roughly at the ion acoustic speed. Another perturbation is propagating inwards from the right side of the box (x ≈ L_box). Although we set a Silver-Mueller boundary condition, particles are removed there. As electrons are propagating and leaving the box faster than the ions, an electric field is set to remove the ions at the same speed as the electrons. The characteristic speed of the perturbation is not identified as it is between the ion acoustic speed and the electron thermal speed. As long as both perturbations remain far from the diamagnetic cavity boundary layer, initialised around x ~ 200−300 km, it should not affect our results.

Regarding ${\epsilon }_0{E}_x^2$, perturbations are only associated with the current along x such that: $${\partial }_t{\epsilon }_0{E}_x^2=-2J_x{E}_x,$$(12)

where J_x is the plasma current along x. Away from the DCBL where the magnetic field can be assumed constant, we see that the fluctuations of $E_x^2$. are on average smaller in the magnetised part (B_z(x) ≈ B₀) than in the unmagnetised one (B_z(x) ≈ 0). The mean electric field in both parts, magnetised and unmagnetised alike, is the ambipolar field, of the same value. However, as ions and electrons are magnetised, this prevents or limits any spurious current J_x that drives the fluctuations. The largest values and fluctuations of $E_x^2$ are observed near the DCBL and within the unmagnetised part. The DCBL grows larger over time but this becomes clearer in $E_y^2$. which is discussed below. The flanks move away from the initial position of the boundary faster than the ion acoustic speed $\left(v_{\text{th},i}\approx \sqrt{k_{\text{B}}{T}_e/m_i}={10}^{-3}c\right)$ but slower than the electron thermal speed $\left(v_{\text{th},i}\approx \sqrt{k_{\text{B}}{T}_e/m_i}=0.1c\right)$. Further investigations showed that the DCBL grows larger at a speed closer to the local Alfvén speed. Within the layer, E_x is equal to almost zero because of a diamagnetic current; this is discussed in detail in Sect. 4.

The variations seen in the middle panel of Fig. 2 are of particular interest. At the beginning of the simulation, we observe the propagation of light waves from the left side of the simulation box which are absorbed by the right side, more likely as boundary effects. However, from 0.001 s to 0.005 s, light waves start from the boundary. The presence of large fluctuations of $E_y^2$ highlights the fact that our initial magnetic profile is not in equilibrium. Indeed, the temporal variations of $E_y^2$ are entangled with those of $B_z^2$ such that: $${\partial }_t\left({\epsilon }_0{E}_y^2+B_z^2/{\mu }_0\right)=-2J_y{E}_y-2{\epsilon }_0{c}^2{\partial }_x{E}_y{B}_z.$$(13)

At the beginning of the simulation and until 0.005 s, the fluctuations are located at the DCBL, more towards the unmagnetised side. After 0.005 s, the DCBL starts to form and the fluctuations split into two components, one along each flank. These components appear to differ in nature as they have different characteristics, such as timescale. However, this is perhaps biased by the fact that the sampling rate of the field is lower than the typical frequency (e.g. the plasma frequency) and therefore the perceived frequencies are aliased (i.e. modified by the sampling rate; we do not comply with the Nyquist–Shannon sampling theorem as we are limited by data storage) and inaccurate. In order to identify or understand the nature of these waves, we perform a fast Fourier transform (FFT). Figure 3 shows a 2D FFT of E_y around the cavity between 22.7 and 25.7 km (left panel) and in the magnetised region (right panel), from 0 to 0.005 s in time. The diagram shows aliased typical dispersion for light waves and plasma waves, preferentially propagating to the right. As Maxwell’s equations are solved using a finite difference time domain, the dispersion relation for waves is modified, mainly at large wave numbers k. For instance, the wave dispersion relation for modified light waves (ω² = c²k²) becomes: $${\sin }^2\left(\frac{\omega \text{Δ}t}{2}\right)=\frac{c^2\text{Δ}{t}^2}{\text{Δ}{x}^2}{\sin }^2\left(\frac{k\text{Δ}x}{2}\right),$$(14)

and for plasma waves (${\omega }^2={\omega }_{pe}^2+c^2{k}^2$) becomes: $${\sin }^2\left(\frac{\omega \text{Δ}t}{2}\right)=\frac{{\omega }_{pe}^2\text{Δ}{t}^2}{4}+\frac{c^2\text{Δ}t}{\text{Δ}{x}^2}{\sin }^2\left(\frac{k\text{Δ}x}{2}\right).$$(15)

In addition, they are aliased. The sampling rate of the fields is 100Δt and thus the range of ωΔt covered here is [−π/100; π/100]. The numerically modified wave dispersion becomes: $$\omega \text{Δ}t\equiv \pm 2\arcsin \left(\sqrt{\frac{{\omega }_{pe}^2\text{Δ}{t}^2}{4}+\frac{c^2\text{Δ}{t}^2}{\text{Δ}{x}^2}\sin \left(\frac{k\text{Δ}x}{2}\right)}\right)\left(\text{mod}\frac{2\pi }{100}\right).$$(16)

Similarly, the same relation can be used for light waves by setting ω_pe = 0.

Plasma waves appear and seem likely trapped: they struggle to propagate through the magnetised region and are unlikely to propagate towards the left as the plasma density and are higher. Indeed, as there is a gradient in the plasma density, the plasma frequency decreases as x increases. On the right side of the boundary, within the magnetised part, the perturbations are associated with the electron cyclotron frequency and the dispersion relation diagram (Fig. 3, right panel) reveals the presence of electron Bernstein waves. We remind the reader that the wave analysis is limited to the electrons here. Although quantities are recorded every 100 time-steps, the plasma frequency is not uniform and decreases for increasing x such that the highest frequency resolved near the cavity (where it is n_i ≈ 0.15n₀) is around the plasma frequency. Indeed, although the plasma frequency decreases along x, the ‘perceived’ plasma frequency oscillates between zero and the sampling rate through the simulation box.

Finally, regarding $E_z^2$ as seen in Fig. 2 (bottom panel), nothing noticeable appears except for the propagation of light waves at the beginning of the simulation. Due to the symmetry of the simulation, $${\partial }_t\left({\epsilon }_0{E}_z^2+B_y^2/{\mu }_0\right)=-2J_z{E}_z+2{\epsilon }_0{c}^2{\partial }_x{E}_z{B}_y,$$(17)

where E_z and B_y are expected to be zero or very small.

Figure 4 shows the energy density stored in the different components of the magnetic field, $B_y^2$ and $B_z^2$, with B_x being null. Like $E_z^2$, the diagram of $B_y^2$, (top panel) shows the presence of light waves during the first few milliseconds of the run. However, unlike $E_y^2$, the fluctuations of $B_y^2$, are larger in the unmagnetised region than in the magnetised region dominated by a z-component. In the bottom panel, $B_z^2$ evolves slowly though time and does not fluctuate. It is the only electromagnetic field initialised to a non-zero value with a spatial profile. However, a slight and weak increase in the field emerges from 0.025 to 0.03 s and moves away from the cavity (see inset in Fig. 4). Its origin will be tentatively addressed in Sect. 3.3.

Fig. 2 Colour plot (position vs. time) of the energy stored in the different electric field components. As an indication, lines with squares represent the propagation of structures at different speeds. From the most horizontal line to the most vertical one (with squares): speed of light in vacuum, thermal speed of the electrons $\sqrt{k_{\text{B}}{T}_e/m_e}$, and ion acoustic speed $\sqrt{k_{\text{B}}{T}_e/m_i}$. The vertical line with circles is located at 0.5 Lbox, where Bz is originally B0/2. The colour bar is in logscale (different for each plot) and SMILEI units (here n0mec2).

Fig. 3 Two-dimensional FFT of Ey between 227 and 257 km (left panel, DCBL location) and between 348 and 408 km (right panel, magnetised part), both between 0 and 0.005 s. Sampling rates are Δx in space and 100Δt in time. For the right panel, as we are in the magnetised part, the pulsation is given in terms of electron cyclotron pulsation.

Fig. 4 Colour plot (position vs. time) of the energy stored in the different magnetic field components. Due to the symmetry, Bx is null. As an indication, lines with squares represent the propagation of structures at different speeds. The inset is a zoom of the red box along the boundary. The colour bar is in logscale and SMILEI units (here n0mec2). Figure 2 provides further details of the labels.

3.2 Thermodynamics quantities

Figure 5 shows the xx component of the dynamic and thermal pressures tensors of ions and electrons, that is to say, $$\begin{array}{l}{\text{P}}_{\text{s,dyn}}=n_s{m}_s{V}_s\otimes V_s\hfill \\ {\text{P}}_{\text{s,th}}={\displaystyle \int \left(v-V_s\right){\left(v-V_s\right)}^{\text{T}}}{f}_s\left(r,v,t\right){\text{d}}^{3}v.\hfill \end{array}$$(18)

The ion dynamic pressure in the magnetised and unmagnetised regions (top left panel) increases because of the ion acceleration by the ambipolar electric field. As the initial ion speed (10⁻⁴c) is below the ion acoustic speed (10⁻³ c), ions are accelerated up to 2 × 10⁻⁴c at the end of the simulation. At the DCBL, we observe an increase in the ion dynamic pressure induced by ions going backward (V_i,x < 0). Indeed, because of the steep increase in B_z, ions are repelled by the magnetic barrier and accelerated by the Hall term in Ohm’s law. At the beginning of the simulation, if we apply Ohm’s law for the x component of the electric field, we find: $$E_x\left(x,0\right)\approx -\frac{1}{q{n}_e}{\partial }_x\left({\text{P}}_{e,xx}\left(x,0\right)+\frac{B_z^2\left(x,0\right)}{2{\mu }_0}\right)\approx \frac{k_{\text{B}}{T}_e}{q{H}_n}-\frac{{\partial }_x{B}_z^2\left(x,0\right)}{2q{n}_e{\mu }_0}.$$(19)

The Hall term generates a strong electric field leftwards, locally overcoming the ambipolar field at the DCBL. This field prevents ions from penetrating too far into the magnetised region while electrons cannot because of their lower momentum (i.e. m_eV_e,x < m_iV_i,x) and lower gyroradius. In addition, there are two positions where the electric field is zero, one on each side of the layer, with one being the minimum of the electric potential and the other its maximum (see Fig. 7 for the electric potential and the associated discussion in Sect. 3.3).

The minimum of the electric potential being in the unmagnetised part is a stagnation point where ions accumulate: they cannot propagate back to the left because of the ambipolar field and cannot propagate to the right because of the steep increase in magnetic pressure. The ion dynamic pressure along x in this simulation is probably the clearest tracer of the evolution of the DCBL. As expected, the ion thermal pressure (bottom left panel) is much lower than the dynamic pressure, though ion heating is observed on each side. This is the opposite for the electron pressure and thermal pressure (bottom right panel). Far from the DCBL, the electron dynamic pressure does not vary through time and follows the initial value given by P_e,xx = 0.01 exp(−x/H_n)[n₀m_ec²] between 10⁻²n₀m_ec² (x = 0) and 10⁻⁴n₀m_ec² (x = L_box). In the DCBL, the electron dynamic pressure changes: it increases on the left flank whereas it decreases on the right flank. These variations are associated with variations in both number density and temperature.

Figure 6 shows the electron temperature component along x as well as the electron number density. We see that the electron pressure gradient is driven by both gradients, that is, in temperature and number density. The depletion of electrons on the magnetised side (down 20% between ~245 km and ~285 km) equates to the accumulation of electrons on the unmagnetised side (up to +80% between ~233 km and ~245 km) such that the number of electrons is conserved within the DCBL over time. These variations are coincident with those of the electron temperature along x: the electron temperature increases by 50% below 245 km while it decreases by 20% above 245 km. It is interesting to note that, in a collisionless simulation, electrons can be both heated and above-all cooled down in the absence of collisions.

Fig. 5 Position vs. time of the xx component of the dynamic (Ps,dyn,xx, top panel) and thermal (Ps,th,xx, bottom panel) pressures, ions (left panel), and electrons (right).

Fig. 6 Electron number density (blue) and temperature along the x component (red) as a function of position at t = 0.027 s. Initial values are represented by the black dashed lines. The initial position of the DCBL where Bz = 0.5B0 is indicated by the vertical dashed line.

3.3 Electron dynamics

One of the benefits of an explicit PIC code is its ability to resolve and grasp physics down to the electron scale, spatial and temporal alike. As mentioned above, different models explore different scales, relying on some assumptions as to the scales they cannot resolve, which could make comparisons between models and their respective results difficult. However, a disadvantage of our simulations is that spatial and temporal averages are required, in particular for the electric field, because of the numerical noise, in part due to the representation of ions and electrons by macroparticles. In addition, one must keep in mind that, because of the electron number density gradient, the plasma characteristics of the plasma (electron Debye length and plasma frequency) are differently resolved through the simulation box. As an example we consider two locations, one around x ≈ 0 and the other around the DCBL at x ≈ L_box/2. On the one hand, at x ≈ 0, n_e ≈ n₀ which means that ∆x/L_e,De(0) ≈ 0.9 and ω_pe(0)∆t = 0.081. On the other, around x ≈ L_box/2, n_e ≈ n₀ exp(−2) which means that ∆x/L_e,De(L_box/2) ≈ 0.9 exp(−1) ≈ 0.33 and ω_pe(0)∆t = 0.03.

In this simulation, we explore the equilibrium at play on the electrons regarding the momentum equation which is along x: $$\begin{array}{lll}{m}_e{\partial }_t\stackrel{-\frac{J_e}{q}}{\overbrace{\left(n_e{V}_{e,x}\right)}}\hfill & =\hfill & -{\partial }_x\stackrel{{\text{P}}_{\text{e,dyn,xx}}}{\overbrace{\left(m_e{n}_e{V}_{e,x}^2\right)}}-{\partial }_x{\text{P}}_{\text{e,th,xx}}\hfill \\ \hfill & \hfill & \underset{\text{Lorentz\hspace{0.17em}force}}{\underbrace{-q{n}_e{E}_x+J_{e,y}{B}_z-J_{e,z}{B}_y}}+\text{residuals,}\hfill \end{array}$$(20)

where the residuals are fluctuations associated with numerical effects which are hard to quantify precisely. Figure 7 shows the balance of the different forces at play on the electrons, namely the electron thermal pressure gradient, the electric and the magnetic forces acting upon electrons (the electron dynamic pressure is ignored), and the associated changes of n_e and B_z. The quantities shown have been smoothed temporally and spatially as follows. We performed a spatial moving average over 21 points, specifically [x₀ − 10∆x, x₀ − 9∆x, …, x₀ + 10∆x]. The subset must not be too large as in this case the sharp structure would be flattened or even disappear. Here, 20 ∆x corresponds roughly to three to four electron Debye lengths at the DCBL. We performed a temporal average over the time interval [t₀; t₀ + 10 000∆t] with one sample every 100 ∆t.

The top left panel of Fig. 7 represents the initial equilibrium of our system on the electrons. We ignore the electron dynamic pressure gradient as it is negligible with respect to other displayed quantities. On each side of the layer, the electric field is balanced by the electron thermal pressure gradient (qn_eE_x ≈ −∂_xP_e,th,xx) and is therefore ambipolar in nature. The problem is that this field is very weak as −∂_xP_e,th,xx ≈ n_e(k_BT_e/qH_n) ≈ 1.36 × 10⁻⁵n_e. At the boundary, the electric field is balanced by the J_e × B_z field. This is responsible for the huge drop in the electric potential by almost 0.04 m_ec², which is four times higher than the initial thermal energy of the electrons. As the electric potential is defined up to a constant, we set max(−qΦ) = 0.01 m_ec². The reader must not be mislead by the potential electric field profile: although it decreases through the DCBL, electrons cannot move into the magnetised part as the magnetic field induces an effective potential barrier. Around 9 ms (top right panel), the electric field separates into two structures, each one moving in opposite direction. The ambipolar field increases following the change in the electron number density. Electrons cumulate on the left side, stopped by the magnetic barrier. The ambipolar field builds up in order to limit the accumulation of electrons. In addition, the electron thermal pressure gradient changes sign approximately where n_e is maximal. The sum of these three forces remains close to zero meaning that Ohm’s law is valid at the scale we are looking at. At around 18 ms (bottom left panel), the DCBL splits into three distinctive regions with different equilibria. Around 240 km, the electric field is purely ambipolar (qn_eE_x + ∂_xP_e,th,xx ≈ 0). The electric field is associated with a steep drop of −qϕ ~ k_BT_e, which falls within the definition of a double layer (Block 1978) as electrons stagnate and accumulate at the foot of the magnetic barrier. Further to the right, the electric field is approximately zero with a flat electric potential, meaning that ∂_xP_e,th,xx ≈ J_e,yB_z. The electron current is therefore diamagnetic in nature. Within this region: $${\text{P}}_{\text{e,th,xx}}+\frac{B_z^2}{2{\mu }_0}~\text{constant,}$$(21)

which corresponds to the balance at the ionopause (Cravens 1986). Finally, −qϕ drops again with a smoother slope associated with the increase in the magnetic field. Here, the electric field is dominated by the Hall term (qn_eE_x ≈ J_e,yB_z), meaning that ions and electrons are decoupled, and the slope is advected with the electrons at V_e,x. In this region, we have J_y,i ≈ −(m_e/m_i)J_y,e.

More interestingly, although the electric field associated with the electron pressure gradient moves away, dissociating itself from the magnetic wall, a little bump in the magnetic field forms at the same exact location (cf. bottom right panel for t ≈ 27 ms). Between the magnetic bump and the magnetic wall, electrons are trapped even if they are heated. At the bump, the electric potential prevents electrons <0.01 eV from moving leftwards. The magnetic bump increases this barrier by a few 0.001 eVs preventing heated electrons above 0.01 eV seen in Fig. 6 from escaping. On the right side, electrons are stopped by the magnetic field only.

Figure 8 shows the pressure profiles and combinations of them at the DCBL for t = 27.24 ms. Except for electromagnetic pressures, the total (dynamic plus thermal) pressures of the ions and electrons exhibit a steep increase at the location of the magnetic ‘bump’ and the electric field peak (see inset) ahead of the magnetic barrier where E² > c² B². From this bump at ~235 km and up to ~257 km, ${\text{P}}_{\text{e,xx}}+B_z^2/\left(2{\mu }_0\right)$ ≈ constant even in the almost unmagnetised region where hotter electrons are found. This suggests that electrons and ions are ‘feeling’ the boundary as the flow is subsonic in the diamagnetic cavity. Even if steady state is not reached, such steepness in the ion and electron pressure spanning over 20–30 grid cells (about ten local Debye lengths and of the order of the electron skin depth) might suggest that both regions are almost disconnected, and particles cannot pass through. If the steady state is reached, and in the absence of collisions, sources, and losses, the pressures must obey: $$\frac{\text{d}}{\text{d}x}\left({\text{P}}_{\text{i,dyn,xx}}+{\text{P}}_{\text{e,dyn,xx}}+{\text{P}}_{\text{i,th,xx}}+{\text{P}}_{\text{e,th,xx}}-\frac{{\epsilon }_0{E}_x^2}{2}+\frac{B_z^2}{2{\mu }_0}\right)=0,$$(22)

as derived from the conservation form of the momentum equation for both species. This form allows the derivation of jump conditions at collisionless shocks when particles cross the surface, bridging their properties upstream and downstream (Rankine-Hugoniot relations). However, even if not applicable here because of non-stationarity, it should be remembered that such an equation holds only if all particles pass through the specific surface, meaning that the equation may hold on each side of the shock but not necessarily at the shock itself, if impermeable. Indeed, the ion pressure is mainly dynamic and drops to zero (P_i,xx, see inset) as the ion current (not shown) is zero going from positive values inside the cavity to negative ones (∂_xn_s V_s,x < 0 and therefore ∂_tn_s > 0). As the simulation neither starts at nor reaches a steady state, ions initially starting where the magnetic field increases were not trapped. Instead, they are quickly pushed away towards the cavity boundary, accelerated by the Hall term present at the beginning of the simulation because their initial kinetic energy is too small (10⁻⁴ m_ec² at the start of the simulation). This causes the negative ion current along x within the DCBL. Similarly, another point is also present further away around 280 km where the ion current goes from negative to positive. However, it is not a stagnation point as ∂_xn_sV_x,s > 0 and therefore ∂_t n_s < 0.

The last properties we investigate are the adiabaticity and gyrotropy of the electrons at the DCBL. Figure 9 shows the diagonal components of the electron thermal pressure tensor as a function of the electron number density. Such a plot allows us to visualise whether or not electrons are adiabatic. The profiles exhibit different slopes through the DCBL. For P_e,th,xx, as we move with increasing x (or from higher to lower electron densities), assuming P_e,th,xx ∝ n^α, we obtain, in succession: α ≈ 1, 1.46 (between blue and red points), 0.5 (between red and yellow points), 1.8 (between yellow and green points), 4.85 (between green and cyan points), and again 1. None of these values appear to be close to the common ones for adiabatic indexes (e.g. 5/3). However, our plasma is not realistic, as ions and electrons are represented by macroparticles. In addition, our simulation is 1D3V and there is no evidence that it does not affect the index. Adiabatic indexes are directly related to the degrees of freedom, which are limited in space to one dimension. Nevertheless, the profiles are evidence that non-adiabatic heating and cooling might occur at the DCBL, in particular where E_x ≈ 0 and the current is diamagnetic (between red and green points).

Therefore, the electron thermal tensor has to be resolved. Furthermore, P_e,th,xx and P_e,th,yy, both components perpendicular to the local magnetic field, behave differently between the blue and yellow locations, around the low magnetised region of the DCBL. Non-diagonal terms of P_e have been ignored as their values are below 10⁻⁵, two orders of magnitude lower than those on the diagonal. In this thin, weakly magnetised region, both temperatures are different and the reason for this remains unknown.

Fig. 7 Evolution over time of the electron number density and magnetic field along z (upper panel) as well as the different forces acting on the electrons: electric field, electron pressure gradient, and x component of the (Je × B) force. The electric potential φ is also displayed and values should be read with respect to the right y-axis. An animated version is available online.

Fig. 8 Pressure profiles at the DCBL. Ps,xx represents the xx component of the total pressure (dynamic plus thermal) for each species at t = 27.24 ms. The inset shows a zoom in to the foot of the ramp.

Fig. 9 Smoothed diagonal components of the electron thermal pressure tensor as a function of the electron number density at the DCBL in logscale, at t = 27.24 ms. Pe,th,yy and Pe,th,zz have been purposely scaled so that they do not to overlap with Pe,th,xx. Five inflexion points (blue, red, yellow, green, and cyan) have been identified, plus one (violet) for a later purpose. Their locations within the DCBL are indicated within the inset where Bz is plotted as a function of x.

3.4 Electron distribution function

Figure 10 shows the electron velocity distribution function (EVDF) at several locations of the DCBL for four distinctive times. At the beginning of the simulation, the EVDF is homogeneous and Maxwellian through the DCBL. After a few tens of milliseconds, an electron beam forms at the foot of the magnetic ramp, deforming the initially isotropic VDF, and then propagates towards the unmagnetised region and the bottom of the magnetic ramp (red point). Its onset is seen in the yellow panel (upper right) at 9.08 ms and 19.16 ms and then passes through the red point (upper middle panel) around 27.24 ms. The most interesting feature is the EVDF at 27.24 ms in the upper middle panel (red axis frame). This is located between the double layer (to the left of the second circle) and the magnetic ramp itself. Electrons are trapped here by the electric potential on their left and the magnetic ramp on their right but this region remains unmagnetised. The EVDF is highly anisotropic and not even gyrotropic. The EVDF is relatively flat along v_x. In addition, the distribution function is more likely noisier in the unmagnetised region. As we averaged over a few time steps, in the magnetised region, the EVDF rotates in a plane that is perpendicular to the magnetic field, smoothing the initial irregularities.

Figure 11 shows normalised phase plots along x for v_x and v_y components at t = 27.24 ms. The magnetic field has been replotted for visualisation purposes. We normalised distributions with respect to the maximum value at a given x to highlight changes in the temperature for example. The double layer is noticeable in the x vs v_x panel with a steep increase in the temperature along x. This is also associated with a negative speed along x for the ions: ions are blocked and reflected by the electric potential of the double layer. Simultaneously, at the double layer, the ion velocity along y stops decreasing and increases again. Nothing appears in the second panel up to this point.

Between the double layer and the yellow point, the magnetic field goes from a plateau around zero and increases to a first step around 240 km. In this region, the phase space x − v_y evolves: the electron velocity along y increases to sustain the magnetic field gradient. Interestingly, the y component of ions increases and becomes positive although weakly compared with that of electrons. This might be caused by ions falling from the magnetic ramp, pushed away by the Hall term, going towards negative x and starting to gyrate clockwise such that they are accelerated towards positive y. In this simulation, ions do not have time to perform a full gyration. As the magnetic field is around 0.1 in SMILEI units, this means that ω_ce = 0.1ω_pe and therefore ω_ci = 10⁻⁵ω_pe. A gyroperiod for an ion is T_ci = 2π/ω_ci = 2 × 10⁵π/ω_pe ≈ 350 ms. The transition between the two magnetic ramps is seen in both plots, with less noisy distributions.

There is a clear change in x − v_y phase space plot. Between the yellow and green dots, both phase-space plots clearly exhibit a decreasing spread in velocity and confirm the electron cooling up to the green dot, before the last magnetic ramp where the Hall term dominates the electric field. Here, the electron temperature increases again to initial values.

Fig. 10 Evolution over time of the EVDF (the colour scale is in arbitrary units and therefore not indicated) in the plane perpendicular to the magnetic field in several locations through the DCBL; same as Fig. 9 but for different times. For each panel, the two top rows are EVDFs which have been sampled between −0.2c and 0.2c every 0.004c. No scales are provided as the purpose is only qualitative. Axis frames are of the same colour as markers indicating a specific location given in the bottom row overlayed on the magnetic field plot. Velocities (x and y axes of EVDFs) are given in terms of c. We performed a temporal average over the time interval [t0; t0 + 10 000∆t] with one sample every 100 ∆t. An animated version is available online.

Fig. 11 Normalised phase-space plots vx vs. x (top panel) and vy vs. x (middle panel) for electrons at t = 27.24 ms in arbitrary units. The figure can be animated online.

4 Discussion

4.1 The electric field and its reliability

Away from the sharp increase in the magnetic field and within the major part of the simulation box, the electric field is ambipolar. It is not clear from the figures whether the electric field is very small or drowned within the fluctuations. This effect is due to the coarse-graining of the simulation, that is, space discretisation, and the relatively low number of particles used to simulate the plasma, which is not able to contain all of the noise. In real plasmas, the fluctuations of the electric field energy from ‘two-particle’ correlations are given by (e.g. Krall & Trivelpiece 1973; Callen 2006): $$\epsilon =\frac{<{ϵ}_0{E}_0^2>}{n{k}_{\text{B}}T}~\frac{1}{\text{Λ}},$$(23)

where Λ is the plasma parameter. Fluctuations of the electric field arise from the finite number of particles within a Debye sphere: these are increasingly damped as the plasma parameter increases. Close to the nucleus, the Debye length is of the order of 1 m or less and the plasma number density is typically between 100 and 1000 m⁻³, such that Λ ≥ 10⁸.

In our simulation, the ambipolar electric field is given by $$E_{\text{amb}}=-\frac{\nabla \cdot {\text{P}}_{\text{e}}}{q{n}_e}\approx \frac{k_{\text{B}}{T}_e}{q{H}_p}{e}_x,$$(24)

where H_p stands for the plasma density scale height. It is then relevant to compare the fluctuations with the mean electric field so that: $$\frac{<E^2>}{{<E>}^2}\approx \epsilon \frac{H_p^2}{{\lambda }_{e,\text{De}}^2}.$$(25)

Observations and theoretical works have shown that H_p ≈ r near the nucleus, with r being the distance from the cometary nucleus (e.g. Balsiger et al. 1986; Edberg et al. 2015; Beth et al. 2019). Therefore, the fluctuations of the electric field are of the order of or lower than the mean electric field where the ambipolar field is the main contribution. The ‘two-particle’ correlations described above are not included in Vlasov theory (Krall & Trivelpiece 1973) and by consequence are not in PIC models either. Nevertheless, the limited number of macro-particles per grid cell in a PIC simulation induces fluctuations that are larger than in real plasma (Melzani et al. 2014). In addition, reducing the electron temperature and therefore the Debye length may aggravate the problem, as shown by Eq. (25). Only substantially increasing the number of particles (from tenfold to hundredfold) may efficiently reduce numerical fluctuations. Another strategy would be using implicit PIC codes such as iPIC3D instead of the explicit SMILEI. Indeed, as mentioned in Markidis et al. (2010), the fast transverse motion of the electrons is damped, which may be seen as a Darwin-like approximation (i.e. neglecting the time derivative of the transverse component of the electric field) though it is not clear how this would work in the unmagnetised part.

Finally, there is one last issue to raise concerning the electric field. As shown by our simulations, the change of the electric potential fulfils the definition of a double layer (ϕ ≳ k_BT_e). Block (1978) showed that the typical width of the double layer is at least of the order of $\sqrt{m_i/m_e}{\lambda }_D$ which is consistent with our simulations. The electric potential drops over a typical spatial scale of 100–200 electron Debye lengths while theory predicts 100 electron Debye lengths (our ion-to-electron mass ratio is 10 000). This result will be important for any future attempt to model this transition region with other approaches or codes. In PIC simulations, it is common to reduce the ion-to-electron mass ratio (see Deca et al. 2017, 2019; Divin et al. 2020, who use a realistic proton mass and an increased electron mass) which might prevent the accurate modelling of this structure. In addition, not all PIC simulations have to resolve the Debye length, only explicit PIC simulations. If the spatial grid is too large (considered as an advantage of the implicit PIC approach over explicit PIC ones) or a spatial average includes too many cells, the double-layer structure may not be resolved as well.

4.2 Agyrotropy and non-adiabaticity of the electrons

As mentioned in Sect. 3.3, electrons are not adiabatic. Understanding the reasons for this would require investigating the electron entropy. Unfortunately, the latter is not currently part of the diagnostic available in SMILEI. Nevertheless, there are already some hints regarding the role of electron adiabaticity on the DCBL. In Koenders (2015), two simulations were carried out: one assuming an adiabatic equation of state for the electrons, the other solving the electron pressure equation. Even if these simulations are hybrid (and thus electrons are considered as a massless fluid), solving the electron pressure adds a second-order correction to the distribution function. The simulations exhibit high discrepancies in the region between the bow shock and the DCBL. In particular, in the adiabatic simulation, the electron temperature monotonically decreases as we go from the solar wind to the DCBL. However, in the non-adiabatic case, electrons are heated with higher temperatures by one to two orders of magnitude compared with the ‘adiabatic’ simulation. This may critically affect the plasma density: electron-ion dissociative recombination plays an important role in terms of loss. Its efficiency depends on the electron temperature. The lower the electron temperature, the higher the plasma loss. Koenders (2015) also showed that the size of the diamagnetic cavity reduces in the non-adiabatic case.

With PIC simulations, we have the ability to verify and check if any assumptions and closures made regarding the electron behaviour at small scales are a posteriori true. Even if our simulations do not incorporate all the required physics to properly model the ionised cometary environment, these first results provide valuable hints. The first closure relation used in MHD and hybrid simulations is the generalised Ohm’s law. By assuming strict neutrality of the plasma and in the low-frequency limit, the electric field cannot be solved from Maxwell’s equations. Instead, the electrons are assumed massless, without inertia, instantaneously adapting to the ambient electromagnetic fields. Such a formulation seems to hold here once quantities are averaged over time and space, still at small spatial scales such as a few electron inertial lengths. This was suggested by Deca et al. (2019). However, within Ohm’s law, another approximation is made: the electron pressure tensor. It can be assumed either isotropic, that is to say, $${\text{P}}_{\text{e}}=n_e{k}_{\text{B}}{T}_e{\text{I}}_3,$$(26)

or gyrotropic $${\text{P}}_{\text{e}}=n_e{k}_{\text{B}}\left[T_{e\Vert }b{b}^{\text{T}}+T_{e\perp }\left(I-b{b}^{\text{T}}\right)\right],$$(27)

where one of the eigenvectors of P_e is along the local magnetic field and eigenvalues associated with the perpendicular directions are equal. The latter means that the shape of the distribution function can be seen as a prolate or oblate spheroid along the local magnetic field direction in velocity space, centred on the electron mean velocity. Our results show that the electron distribution is not isotropic but also not perfectly gyrotropic near the DCBL, although the eigenvectors of P_e are along the x, y and z axes.

Figure 12 shows the different electron temperatures and ‘invariants’ described in Appendix A which highlights the anisotropy of the electron distribution function. The so-called invariants are values derived from the analysis of the characteristic polynomial of P_e. They are called invariants because they do not depend on the frame in which P_e is represented. If the three temperatures of the electron distribution function are distinct, the discriminant ∆ is positive (see Appendix A). If two of them are equal and distinct from the third temperature, ∆ = 0 but −p³ > 0. If all temperatures are equal p³ = q² = ∆ = 0.

As electron anisotropies develop at the DCBL, we check for two possible instabilities: mirror and firehose (see e.g. Rosenbluth 1956; Chandrasekhar et al. 1958; Rudakov & Sagdeev 1958). Firehose instabilities would occur in regions where $P_{\perp }-P_{\Vert }>B_z^2/{\mu }_0$. Even if we have a small region where P_⊥ − P_‖ > 0, the difference is not large enough to overcome $B_z^2/{\mu }_0$ and trigger the instability in our simulation where our magnetic field strength is very large compared with reality. However, in the weakly magnetised region, where both perpendicular temperatures are larger than the parallel temperature, the instability condition for mirror mode generation may be fulfilled. Either assuming T_⊥ = P_e,th,xx/n_e or P_e,th,yy/n_e, the condition $$2n_e{k}_{\text{B}}{T}_{\perp }\left(\frac{T_{\perp }}{T_{\parallel }}-1\right)>\frac{B^2}{{\mu }_0},$$(28)

holds for x ≲ 241 km. We might also include the ions but their temperature is negligible compared to that of the electrons. Nevertheless, these criteria were derived for bi-Maxwellian distributions. In contrast, EVDFs exhibit structures far from Maxwellian or bi-Maxwellian distributions characterised, sometimes, by three different temperatures (cf. Fig. 10, top middle panel at 27.27 ms). Previous works investigated the stability of this region at 1P/Halley. In particular, Ershkovich et al. (1989) performed a stability analysis of both ionospheric structures and balances proposed by Cravens (1986) and Ip & Axford (1987). The latter were shown to be unstable. However, taking mass loading (e.g. photoionisation) and unloading (e.g. dissociative recombination) effects into account helps to stabilise the system, except near the boundary.

Fig. 12 Profiles of the invariants of Pe/ne described in Appendix A (upper panel) and of the different electron temperatures. As the non-diagonal terms of Pe/ne have values of close to zero, the eigenvalues of Pe/ne are its diagonal terms. The vertical black line delimits the region where the mirror instability might occur (on the left side) as the magnetic field strength is low.

4.3 Are any of these plasma characteristics observable with past and future cometary missions?

Definitely one of the most important questions to address is whether or not any of these plasma characteristics are observable with past and future cometary missions. As we show in our simulation, several features are associated with the DCBL. Firstly, we show that a double layer forms near the DCBL with a drop in the electric potential of the order of the electron temperature. This structure is actually too sharp (less than one kilometre and even smaller in reality) to be probed by Rosetta or Comet Interceptor, as electric fields are quite difficult to measure in general. Secondly, we observed a clear evolution of the electron distribution through the boundary. Evidence for a difference in the electron distribution function between inside and outside the diamagnetic cavity was found at 67P (Nemeth et al. 2016). However, one has to keep in mind that these observations of the electron distribution function are limited by several factors: the time resolution of the instrument, the relative speed of the spacecraft, and the spacecraft potential for instance. At this stage, it is difficult to know if probing the EVDF with great spatial and temporal resolution would be possible with the future ESA mission Comet Interceptor. However, it appears critical that in order to probe these structures in the future, we need multi-point measurements and orbiting spacecraft for an extended period of time as suggested by Goetz et al. (2021).

5 Conclusions

In this paper, we present the first fully kinetic PIC simulation of the so-called DCBL that encases the diamagnetic cavity, in the collisionless limit. Although the initial setup cannot be close to reality as some aspects regarding cometary physics have been ignored, simulations provide valuable insights into this structure and the role of the electrons. In particular, in this paper, we focus on the evolution over time of the electromagnetic fields and thermodynamic quantities associated with the electrons and ions as well as the pressure balance at the DCBL. In addition, we investigated the time and spatial dependence of the electron distribution function across the DCBL.

Among the features observed in our PIC simulation, several are of great interest and therefore must be considered and investigated in detail in the future. Firstly, a double layer forms at the DCBL and propagates towards the unmagnetised region. This is of major importance as such a structure may be revealed only through fully kinetic modelling. Secondly, within the DCBL, electrons do not have a Maxwellian distribution and depart from a gyrotropic distribution. In addition, their behaviour does not seem adiabatic, in contrast with assumptions made for MHD and hybrid models. For a better understanding and more accurate modelling of the DCBL, we are planing to perform 2D3V simulations.

Movies

Movie 1 associated with Fig. 7 (fig7-animated) Access here

Movie 2 associated with Fig. 10 (fig10-animated) Access here

Movie 3 associated with Fig. 11 (fig11-animated) Access here

Acknowledgements

We warmly thank the referee for his/her positive comments and promptitude. Work at Umeå University was supported by the Swedish National Space Agency (SNSA) grant 108/18. CSW thanks the Austrian Science Fund (FWF) P32035-N36. The authors would like to acknowledge ISSI for the opportunity it offered for very valuable discussions on this topic as part of the International Team 499 “Similarities and Differences in the Plasma at Comets and Mars”. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the High Performance Computing Center North (HPC2N) at Umeå University in Sweden. The colour scale used in Figs. 2–4, 5, 10, and 11 which is colour-blind friendly has been generated with the MATLAB code provided by Matthias Geissbuehler (see Geissbuehler & Lasser 2013).

Appendix A Anisotropy in the electron distribution function

One of the strengths of PIC simulations resides in the ability to assess the moments of the ion and electron distributions compared with MHD and hybrid simulations. In particular, the second-order moment of the electron pressure tensor is of interest. The commonality between MHD and hybrid is the assumption that the electron pressure tensor is diagonal in the frame of the local magnetic field, with two temperatures, one parallel to the magnetic field and one perpendicular, that is, P_e is given by: $${\text{P}}_{\text{e}}=n_e{k}_{\text{B}}\left[T_{e\Vert }b{b}^{\text{T}}+T_{e\perp }\left(I-b{b}^{\text{T}}\right)\right].$$(A.1)

Hence, the distribution function might be perceived as an ellipsoid prolate or oblate along B centred on the mean electron velocity (first-order moment). In order to compare results from PIC, hybrid, and MHD simulations, we consider P_e to be a good starting point: it is the highest order modelled by hybrid and MHD models. P_e is a symmetric matrix such that its eigenvalues are positive and its eigenvectors are orthogonal to each other. Therefore, P_e may be written as follows: $${\text{P}}_{\text{e}}=n_e{k}_{\text{B}}\left[T_{e1}{\mathcal{V}}_1{\mathcal{V}}_1{}^{\text{T}}+T_{e2}{\mathcal{V}}_2{\mathcal{V}}_2{}^{\text{T}}+T_{e3}{\mathcal{V}}_3{\mathcal{V}}_3{}^{\text{T}}\right],$$(A.2)

where V_j stands for the normalised eigenvectors associated to the temperature T_ej.

In order to assess whether or not the distribution function is isotropic, gyrotropic, or anisotropic, the pressure tensor and its invariants should be scrutinised. The invariants of the pressure tensor are independent of the frame in which the pressure is represented. For instance, the eigenvalues are invariants of P_e. Moreover, there exists different sets of three invariants. One of the common sets is: $$I_1=\text{tr}\left({\text{P}}_{\text{e}}\right)$$(A.3) $$I_2=\frac{1}{2}\left[\text{tr}{\left({\text{P}}_{\text{e}}\right)}^2-\text{tr}\left({\text{P}}_{\text{e}}^2\right)\right]$$(A.4) $$I_3=\det \left({\text{P}}_{\text{e}}\right),$$(A.5)

such that the eigenvalues are solutions of the polynomial equation: $${\lambda }^3-I_1{\lambda }^2+I_2\lambda -I_3=0.$$

As P_e is a symmetric tensor, all λ_i are real. In addition, because of the way the different components of P_e are calculated (see Eq. 18), the eigenvalues are positive. Indeed, it can be shown that ∀_i, j, $P_{i,i}{P}_{j,j}\ge P_{i,j}^2$ (Cauchy-Schwarz inequality) such that I₁, I₂, and I₃ are positive. Using Descartes’ rule, it follows that the polynomial should have three positive roots. Dividing λ_i by the local electron density gives the three different temperatures of the distribution function. Three different cases exist: one single root (multiplicity 3; isotropic distribution), two distinct roots (one with multiplicity 2, the second with multiplicity 1; e.g. gyrotropic distribution), and three distinct roots (anisotropic distribution). The number of positive roots depends on the discriminant ∆ from Cardano’s method, given by: $$\text{Δ}=-p^3-q^2={\left({\lambda }_1-{\lambda }_2\right)}^2{\left({\lambda }_2-{\lambda }_3\right)}^2{\left({\lambda }_3-{\lambda }_1\right)}^2$$(A.6) $$p=-\frac{I_1^2}{9}+\frac{I_2}{3}$$(A.7) $$q=-\frac{I_1}{54}\left(2I_1^2-9I_2\right)-\frac{I_3}{2}.$$(A.8)

As the roots are real, ∆ ≥ 0. If ∆ = 0, only two roots exist. In addition, if p = q = 0, there is one single root (isotropic distribution). p has a mathematical meaning: −2p is the variance of eigenvalues (p < 0); it allows to track any departure of a distribution function from the isotropic case.

References

All Tables

All Figures

	Fig. 1 Initial setup and profiles for ions, electrons, and magnetic field in our simulation box. Ion temperature is not indicated as it is set to 0.
In the text

Fig. 2 Colour plot (position vs. time) of the energy stored in the different electric field components. As an indication, lines with squares represent the propagation of structures at different speeds. From the most horizontal line to the most vertical one (with squares): speed of light in vacuum, thermal speed of the electrons $\sqrt{k_{\text{B}}{T}_e/m_e}$, and ion acoustic speed $\sqrt{k_{\text{B}}{T}_e/m_i}$. The vertical line with circles is located at 0.5 Lbox, where Bz is originally B0/2. The colour bar is in logscale (different for each plot) and SMILEI units (here n0mec2).

In the text

	Fig. 3 Two-dimensional FFT of Ey between 227 and 257 km (left panel, DCBL location) and between 348 and 408 km (right panel, magnetised part), both between 0 and 0.005 s. Sampling rates are Δx in space and 100Δt in time. For the right panel, as we are in the magnetised part, the pulsation is given in terms of electron cyclotron pulsation.
In the text

Fig. 4 Colour plot (position vs. time) of the energy stored in the different magnetic field components. Due to the symmetry, Bx is null. As an indication, lines with squares represent the propagation of structures at different speeds. The inset is a zoom of the red box along the boundary. The colour bar is in logscale and SMILEI units (here n0mec2). Figure 2 provides further details of the labels.

In the text

	Fig. 5 Position vs. time of the xx component of the dynamic (Ps,dyn,xx, top panel) and thermal (Ps,th,xx, bottom panel) pressures, ions (left panel), and electrons (right).
In the text

	Fig. 6 Electron number density (blue) and temperature along the x component (red) as a function of position at t = 0.027 s. Initial values are represented by the black dashed lines. The initial position of the DCBL where Bz = 0.5B0 is indicated by the vertical dashed line.
In the text

	Fig. 7 Evolution over time of the electron number density and magnetic field along z (upper panel) as well as the different forces acting on the electrons: electric field, electron pressure gradient, and x component of the (Je × B) force. The electric potential φ is also displayed and values should be read with respect to the right y-axis. An animated version is available online.
In the text

	Fig. 8 Pressure profiles at the DCBL. Ps,xx represents the xx component of the total pressure (dynamic plus thermal) for each species at t = 27.24 ms. The inset shows a zoom in to the foot of the ramp.
In the text

Fig. 9 Smoothed diagonal components of the electron thermal pressure tensor as a function of the electron number density at the DCBL in logscale, at t = 27.24 ms. Pe,th,yy and Pe,th,zz have been purposely scaled so that they do not to overlap with Pe,th,xx. Five inflexion points (blue, red, yellow, green, and cyan) have been identified, plus one (violet) for a later purpose. Their locations within the DCBL are indicated within the inset where Bz is plotted as a function of x.

In the text

Fig. 10 Evolution over time of the EVDF (the colour scale is in arbitrary units and therefore not indicated) in the plane perpendicular to the magnetic field in several locations through the DCBL; same as Fig. 9 but for different times. For each panel, the two top rows are EVDFs which have been sampled between −0.2c and 0.2c every 0.004c. No scales are provided as the purpose is only qualitative. Axis frames are of the same colour as markers indicating a specific location given in the bottom row overlayed on the magnetic field plot. Velocities (x and y axes of EVDFs) are given in terms of c. We performed a temporal average over the time interval [t0; t0 + 10 000∆t] with one sample every 100 ∆t. An animated version is available online.

In the text

	Fig. 11 Normalised phase-space plots vx vs. x (top panel) and vy vs. x (middle panel) for electrons at t = 27.24 ms in arbitrary units. The figure can be animated online.
In the text

	Fig. 12 Profiles of the invariants of Pe/ne described in Appendix A (upper panel) and of the different electron temperatures. As the non-diagonal terms of Pe/ne have values of close to zero, the eigenvalues of Pe/ne are its diagonal terms. The vertical black line delimits the region where the mirror instability might occur (on the left side) as the magnetic field strength is low.
In the text