Mach number and plasma beta dependence of the ion temperature perpendicular to the external magnetic field in the transition region of perpendicular collisionless shocks

In various kinds of solar-terrestrial, astrophysical, and laboratory plasmas, collisionless shocks are ubiquitous, at which the upstream kinetic energy of the supersonic plasma flow dissipates into downstream energy of thermal ions and electrons, waves (turbulence), and nonthermal particles.^1,2 Despite various kinds of studies, detailed processes of shock dissipation remain to be clarified. For example, we do not fully understand how energies are partitioned between downstream thermal electrons and ions although their total pressure can be simply predicted by the fluid Rankine-Hugoniot relation.

For supercritical (quasi-)perpendicular shocks, a fraction of the incoming ions can be specularly reflected toward the upstream region but gyrates back to the shock front.^3–10 Such reflected-gyrating ions can gain energy from the motional electric field of the upstream plasma flow and contribute to the increase of the ion temperature component perpendicular to the local magnetic field. Consequently, large temperature anisotropy arises at the shock foot, exciting waves through ion temperature anisotropy instability, which is responsible for the shock ripples.^11,12 Electron preheating at the foot also takes place under some conditions.^3,13,14 The ripple further dissipates ions, increasing ion parallel temperature, and even electron acceleration occurs.¹⁵ In the downstream region, ion distribution is no longer nongyropropic, and its structures are smoothed out by collisionless gyrophase mixing, resulting in downstream ion heating.^16–21 In order to understand such a multistep dissipation process across the shock front, it is important to estimate the initial ion temperature component perpendicular to the shock magnetic field at the foot region. In the present study, using the two-dimensional full particle simulation of low-Mach-number, perpendicular, rippled, and collisionless shocks, we study the ion perpendicular temperature at the shock foot region. We show, for the first time, that it is proportional to the square of the Alfvén Mach number divided by the plasma beta or the square of the sonic Mach number, which is consistent with the analytical scaling relation.⁵ 

We perform two-dimensional (2D) simulations of perpendicular (⁠$θBn=90○$⁠) collisionless shocks by using a standard particle-in-cell code.²² As in our previous works,^15,23,24 the shock is excited by the “relaxation” between a supersonic and a subsonic plasma flow moving in the same direction. The initial state consists of two regions separated by a discontinuity. Both regions have spatially uniform distributions of electrons and ions with different bulk flow velocities, temperatures, and densities, and they have a uniform perpendicular magnetic field with a different strength. The simulation domain is taken in the x-y plane and an in-plane shock magnetic field (B_y0) is assumed. We apply a uniform external electric field E_z0 = u_x1B_y01/c (=u_x2B_y02/c) in both upstream and downstream regions, so that both electrons and ions drift along the x axis. Here, u_x is the bulk flow velocity, and subscripts “1” and “2” denote “upstream” and “downstream,” respectively. At the left (right) boundary of the simulation domain in the x direction, we inject plasmas with the same quantities as those in the initial upstream (downstream) region. We use absorbing boundaries to suppress nonphysical reflection of electromagnetic waves at both ends of the simulation domain in the x direction,²⁵ while the periodic boundaries are imposed in the y direction.

In the present study, we show results of six simulation runs (A–F) with different upstream conditions. We summarize in Table I the upstream plasma parameters, such as the bulk flow velocity u_x1 and the ratio of the electron plasma frequency to the electron cyclotron frequency ω_pe1/ω_ce1. For all runs, we fix the ion-to-electron mass ratio as m_i/m_e = 256, and set v_te1/c = 0.1, where v_te1 and c are the electron thermal velocity upstream and the speed of light, respectively, and subscripts “i” and “e” represent “ion” and “electron,” respectively. Then, we obtain the plasma beta $β1=2(vte1/c)2(ωpe1/ωce1)2$ and $ux1/VA1=(mi/me)1/2(vte1/c)(ωpe1/ωce1)(ux1/vte1)$⁠, where V_A1 is the upstream Alfvén velocity. It is assumed that the electrons and ions have the same plasma beta, β_e1 = β_i1 = β₁, and the same isotropic temperature, T_e1 = T_i1. Here, temperatures and thermal velocities are related as $Te1=mevte12$ and $Ti1=mivti12$⁠. For given upstream frequencies ω_pe1 and ω_ce1, the initial upstream number density $n1≡meωpe12/4πe2$ and the initial magnetic field strength B_y01 = m_eω_ce1/e are derived. Then, the initial downstream parameters are determined by solving the shock jump conditions (Rankine-Hugoniot conditions) for a magnetized two-fluid isotropic plasma consisting of electrons and ions,²⁶ assuming T_i2/T_e2 = 8.0.

The grid spacing and the time step of the present simulation runs are set to be Δx = Δy ≡ Δ = λ_De1 and cΔt/Δ = 0.5, respectively, where λ_De is the electron Debye length. The total size of the simulation domain is 32l_i1 × 6l_i1, where $li1=c/ωpi1=(mi/me)1/2(c/vte1)λDe1$ is the ion inertial length of the upstream plasma. We used 25 pairs of electrons and ions per cell in the upstream region and 64 pairs of electrons and ions per cell in the downstream region at the initial state.

The ion temperature component perpendicular to the shock magnetic field is approximated by the arithmetic mean of ion temperatures in the x and z directions, that is, T_i⊥ = (T_x + T_z)/2. In Fig. 1, we show spatiotemporal evolution of B_y (left panel) and T_i⊥ (right panel) for Run D. Here, both of them are averaged over the y direction. The initial unphysical disturbance disappears, and the growth of the shock ripples ceases until ω_ci1t = 7.0 after which the shock structure seems to be in the quasi-steady state. Unlike one-dimensional simulation, quasiperiodic reformation seems unclear although we can still see the front oscillation. Hence, we analyze the data after ω_ci1t = 7.0 for all runs.

As shown in Fig. 1, the shock front moves leftward in our simulation frame. Using the spatiotemporal diagram of B_y, we measured the shock velocity $v$_sh, and obtained the Alfvén Mach number, M_A = (u_x1 − v_sh)/V_A1, in the shock-rest frame. The results for all six runs are shown in Table I.

We extract the representative value of T_i⊥ in the transition layer for each run as follows. First, we take a snapshot of T_i⊥ which is averaged over the y direction and find its maximum value, $Ti⊥max$⁠, for each time step. As shown in Fig. 2, one can find that T_i⊥ has a maximum at the foot region. This fact is also true in an arbitrary epoch. The value of $Ti⊥max$ changes with time. Then, the time mean of $Ti⊥max$ for 7 ≤ ω_ci1t ≤ 12 is obtained. For example, we obtain the average value $Ti⊥max/Ti1=101$⁠, and the maximum and minimum values are 123 and 86.5, respectively. In the same way, the average, maximum, and minimum values of $Ti⊥max$ are obtained for other runs. The results are summarized in Table I. In Fig. 3, $Ti⊥max/Ti1$ is shown as a function of M_A²/β₁. The simulation results seem to lie on the line, $Ti⊥max/Ti1≈0.5×MA2/β1$⁠.

We derive a simple analytical formula to have ion perpendicular temperature at the shock foot. Although similar formulae have been already derived,^5–7 our final equation, Eq. (11), is in an excellent agreement with our simulation results. In our simulation frame, the shock front moves at a velocity v_sh, and upstream and downstream bulk velocities are typically u_x1 and u_x2, respectively. The incoming ions are adiabatically heated at the shock foot. They have a perpendicular temperature, $[Cγ−1/(γ−1)]mivti12$⁠, where C ∼ (n_i,f/n_i1) is the compression factor, and γ and n_i,f are the adiabatic index (γ = 2 in our case) and the typical value of the density at the foot, respectively. Since they are also decelerated due to mass flux conservation, their bulk velocity measured in the rest frame of the shock front becomes u_(in) ≈ $vsh′$/C, where $vsh′=ux1−vsh=MAVA1$⁠. Next, a part of them is reflected. The bulk velocity of the reflected ions in the shock rest frame is u_(ref) ≈ −$vsh′$/C.⁴ Hence, the velocity difference at the shock foot between the incoming and the reflected ions is given below:

A large fraction of energy (per ion), (m_i/2)|Δu|², is consumed for increasing the ion perpendicular temperature.

Here, we consider a simple analytical model to estimate T_i⊥ at the foot region, whose ion distribution function is written as

where the first and the second terms in the rhs describe the incoming and reflected components, respectively, and the number density N_(k), bulk velocity u_(k), and temperature T_(k) are

respectively. A subscript (k) = (in), (ref) denotes each component. Then, it is natural to approximate T_i⊥ as

where N_tot and $ū$ are given by

respectively. When we introduce the fraction of the reflected ions as r = N_(ref)/N_tot, then we get

and $ū$ = (1 − r)u_(in) + ru_(ref). Assuming T_(in) = T_(ref) as in one of our previous simulation studies²⁷ and eliminating ū, we can rewrite Eq. (9) as

Using Eq. (1) together with T_(in) = CT_i1 and m_i$vsh′$²/T_i1 = 2M_A²/β₁, we finally obtain

The first term in the rhs of Eq. (11) is important for the case of low M_A²/β₁ only. The compression factor C is slightly larger than unity for such shocks. For large M_A²/β₁, the second term dominates the rhs of Eq. (11); hence, we can explain our numerical result, $Ti⊥max/Ti1∝MA2/β1$⁠, shown in Fig. 3, if the factor 8r(1 − r)/C² in Eq. (11) hardly depends on M_A²/β₁. The fraction of the reflected ions is typically r ≈ 0.3 and varies from 0.2 to 0.4 during nonstationary processes at the shock front. On the other hand, the value of C is less variable and ranges between 1.0 and 1.1. Then, one can see 8r(1 − r)/C² = 1.1–1.9. This is a factor of a few larger than estimated from Fig. 3. Indeed, our analytical formula, shown in Eq. (11), gives the upper bound of T_i⊥ because the free energy m_i(Δu)² goes not only to T_i⊥ but also to the thermal energy of reflected ions and waves excited in the shock transition layer.²⁸ In practice, one can see that T_i⊥ becomes larger if T_(ref) > T_(in) [see Eq. (9)]. Note that the fraction of the reflected ions, r, estimated in the previous analytical works^4,8,29 is slightly smaller than that obtained from our simulations.

In the present study, we focus on T_i⊥ only. On the other hand, the parallel component of ion temperature (T_i∥ ≈ T_y) is much smaller at the foot region. Therefore, the total ion temperature, T_i = (T_x + T_y + T_z)/3 = (2T_i⊥ + T_i∥)/3, is approximated as T_i ≈ (2/3)T_i⊥.

Using two-dimensional full particle simulations, we have shown that ion perpendicular temperature at the foot of supercritical perpendicular collisionless shocks is proportional to M_A²/β₁ or the square of the sonic Mach number. This fact will give us a simple estimate of the energy partition between downstream thermal ions and electrons although further study is necessary. The ion heating at the foot region of (quasi-)perpendicular shocks has been extensively investigated by many authors mainly using spacecraft observations,^5,9 (semi-)analytical studies,^4,5 and one-dimensional hybrid simulations.⁷ In this paper, we have extended such studies by using two-dimensional full particle simulations which can better capture various kinetic effects including wave excitation and plasma heating in the direction tangential to the shock front. More specifically, we have demonstrated in this paper that our analytical scaling relation, as shown in Eq. (11), is in excellent agreement with two-dimensional full particle simulations of rippled shocks. This result also indicates that the dependence of the fraction of ion reflection on the plasma beta and the Alfvén Mach number is small when $MA2/β$ is larger than about 20, which is consistent with our simulation results.

The authors would like to thank Yutaka Ohira for helpful comments. Computer simulations were performed on the CIDAS supercomputer system at the Institute for Space-Earth Environmental Research in Nagoya University under the joint research program. This work was partly supported by JSPS KAKENHI, Grant Nos. 15K05088, 18H01232 (R.Y.), 2628704, and 19H01868 (T.U.).