Interaction of rotational discontinuities with energetic ions in the precursor of the Earth’s bow shock

We combined in-situ solar wind observations by ARTEMIS and MMS missions with kinetic hybrid simulations to study the interaction of solar wind rotational discontinuities (RDs) with the foreshock of the Earth’s bow shock. We found that whistler modes excited by diffuse energetic particles were strongly coupled with RDs and lead to their temporary dissociation. At the same time, RDs trigger the steepening of whistler waves and the generation of ’shocklets’ - small localised shock-like structures, capable of trapping energetic particles and growing up by absorbing the particles energy.


Introduction
Discontinuities embedded into the solar wind constantly encounter the Earth's bow shock and may impact the magnetosphere [1,2,3]. In our recent paper, [4], we studied the interaction of the bow shock with rotational discontinuities (RDs), which are the most common type of discontinuities. One of the peculiar features found in both simulations and observations was the key role of the interaction of RDs with the precursor of the bow shock. This area is filled with shock-reflected particles and electromagnetic waves generated by them. These particles and waves substantially distort an RD and temporarily transform it into a wave train with a shock-like structure at the trailing edge. This process changes the RD's dispersive properties and efficiently thins the RD during the shock crossing. Thus the RD may become less stable to magnetic reconnection in the magnetosheath. In this paper, we focus on the RD-foreshock interaction and show that it can be associated with the nonlinear steepening of kinetic magnetosonic waves produced by beams of energetic particles [5].

Observations
We use observations by two spacecraft missions in the solar wind: Magnetospheric Multiscale Mission (MMS) around the Earth bow shock and ARTEMIS in the pristine solar wind (see details of spacecraft instruments in [6,7], and references therein). We selected an interval when MMS and ARTEMIS observed the same solar wind discontinuity that travelled from ARTEMIS to the  bow shock. In Fig. 1, this discontinuity manifests itself as a sharp rotation of the magnetic field at 07:35 (ARTEMIS) and 7:50 (MMS). This magnetic field rotation is accompanied by solar wind velocity jumps (second panels), i.e. we likely deal with a rotational discontinuity [8,9]. In the pristine solar wind (ARTEMIS), the solar wind discontinuity does not show any variations in the ion spectra (the third panel) and plasma density (the bottom panel). Near the bow shock, however, the discontinuity interacts with the energetic ions reflected from the shock (see ion fluxes in the right column in Fig. 1). This population is observed only at one side of the discontinuity and is accompanied by a strong wave activity seen both in magnetic field and plasma density. These waves are weakly compressional: fluctuations of the density and magnetic field magnitude are much smaller than fluctuations of the magnetic field and plasma flow components. Thus, we likely deal with weakly oblique Alfven or ion cyclotron (whistler) waves driven by cyclotron resonance with the energetic ion beam [10]. A distinctive feature is sharp density and |B| gradients at the trailing edge of the RD. This feature is very similar to so-called 'shocklets', which were associated with steepened waves on magnetosonic-whistler branch, [5]. This 'shocklet' clearly borders the region with energetic ions.

Simulation technique
Collisionless discontinuities are mostly governed by processes on ion scales, so we use hybrid kinetic-ions/fluid-electrons code "Maximus" [11] to simulate passing of an RD through the foreshock occupied by energetic ions. Such codes treat electrons as a massless fluid and explicitly describe ions kinetics. Hybrid codes are widely used for simulations of collisionless plasma [12,13,14] because they are much faster than full particle-in-cell codes. "Maximus" is fully 3d, second order accurate in space and time, and exactly satisfies div B = 0. The code is parallelized with MPI (Message Passing Interface communication protocol) and shows nearly linear speed-up with increasing number of processes.
In the code, the unit length is the far upstream protons inertial length l i , the unit time is 3 the inverse proton gyrofrequencies Ω −1 , the unit density is the far upstream density. All other units can be derived self-consistently, e.g. velocities are normalized by the far upstream Alfven velocity.

Simulation setup
We used a 3d Cartesian domain with 6000 × 50 × 50 cells and a spatial resolution 0.5l i in each direction. First we simulated interaction of an RD with a quasi-parallel (inclination angle θ = 15 o ) collisionless shock. The shock was initialized by the 'rigid piston' method: the supersonic flow with V x = −3.15V a was moving in the negative x direction faces a reflecting wall at x = 0 and formed a shock travelling in the positive x direction. An RD with initial width 40l i and rotational angle 240 o was embedded into the upstream flow and is was convected towards the shock front. The RD and shock normals were parallel. Hereafter we will refer to this simulation as Run A.
To examine the role of energetic particles in the RD evolution, we also conducted a test run (Run B) where an isolated standing RD interacts with a streaming population of shockaccelerated particles. This simulation setup resembles a foreshock occupied by diffuse accelerated particles streaming relatively to the inflow. Properties of energetic particles were estimated from Run A. The main difference is that in Run B energetic particles were uniformly distributed in space, while in Run A their density decreases with a distance from the shock, and their velocity distribution shifts to higher velocities (see the bottom panel in Fig. 2).

Simulation results
In Fig. 2 we illustrate a general picture of the RD-shock interaction. The left column shows an RD in the pristine solar wind (just as in the left column of Fig. 1). The shock front is at x f ≈ 150l i , while the RD is at x rd ≈ 1450l i . Shock-reflected beam can be seen in the (x, V x ) phase space (bottom panel) at x < 1200l i . Closer to the shock front, the beam becomes more diffuse in the momentum space, indicating the 1st order Fermi acceleration. The beam has not touched the RD yet, so the RD is not distorted. The right column shows the same RD after it entered the region occupied by diffuse distribution of energetic particles (x rd ≈ 1000l i ). The RD has dissociated into a wave-train and is hardly distinguishable. Similarly to the waves observed by MMS, these waves are weakly compressional: density and magnetic field amplitudes oscillate within 10% of the respective mean values, while variations of the transverse magnetic field are about 30%. These waves are right-hand polarised and propagate away from the shock in the upstream reference frame, as will be also seen from the wave analyses below.
There are sharp density ρ and |B| jumps at the trailing edge of this wave-train. These jumps slow down energetic ions and push them back to the shock (see the bottom panel in Fig. 2). This is consistent with the fact that in observations energetic ions reside only at one side of the 'shocklet' (see the right column in Fig. 1). Similar 'shocklets,' resulting from the steepening of kinetic magnetosonic waves, were studied in [5]. So it is natural to suppose that the waves are magnetosonic in our case as well. To check it we performed additional wave analyses in the vicinity of the RD.
We chose two windows co-moving with the RD: the upstream window [x rd − 200l i , x rd ] (bordered by the red dashed lines in Fig. 2), and the downstream window [x rd , x rd + 200l i ] (bordered by the green dashed lines). A 2d Fourier transform of B y and B z in those two windows yielded the wave power spectrum |B ⊥ | 2 (k x , ω) = |B y | 2 +|B z | 2 , whereB y,z is the Fourier transform of B y,z (x, t). The resulting upstream and downstream spectra are shown in the left and middle panels in Fig. 3 respectively. The total wave power drops downstream of the RD (see the middle panel of Fig. 3), which means that there are much fewer waves there. In observations, the waves completely disappear downstream of an RD. A strong mode with a nearly linear dispersion is present only upstream of the RD and completely disappears downstream. This mode might be associated with cold reflected beam, which is substantially slowed down and diluted at the RD. At the same time, a low-frequency ω < Ω i branch is present on both sides of the RD and is possibly a beam-whistler mode excited by the diffuse energetic population.
To examine the role of the diffuse energetic population in generation of the low-frequency mode we simulated such a population streaming across an isolated standing RD (Run B). We computed a wave power spectrum for this run and found it the same for the upstream and downstream windows. So we show the combined spectrum in the right panel of Fig. 3. In this case, the spectrum is broader, but the similarity with Run A is apparent. The dispersion curve is convex, resembling those of the magnetosonic-whistler branch. The right-hand polarisation, seen in the third row of Figs. 2 and 4 is also in agreement with the whistler mode.
A general picture of RD-energetic particles interaction is shown in Fig. 4 and is very similar to the one in Fig. 2. The RD dissociates into a train of right-hand polarised weakly compressional waves. These waves efficiently trap energetic particles and further steepen, absorbing the energy of these particles. A small 'shocklet' is also seen at x ≈ 300l i . Thus the diffuse population of energetic particles upstream the bow shock strongly interact with solar wind RDs by means of steepened magnetosonic-whistler modes.