Ion temperature anisotropy is a common feature for (quasi-)perpendicular collisionless shocks. By using two-dimensional full particle simulations, it is shown that the ion temperature component perpendicular to the shock magnetic field at the shock foot region is proportional to the square of the Alfvén Mach number divided by the plasma beta. This result is also explained by a simple analytical argument in which the reflected ions get energy from an upstream plasma flow. By comparing our analytic and numerical results, it is also confirmed that the fraction of the reflected ions hardly depends on the plasma beta and the Alfvén Mach number when the square of the Alfvén Mach number divided by the plasma beta is larger than about 20.

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.^{15} 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.^{5}

We perform two-dimensional (2D) simulations of perpendicular ($\theta Bn=90\u25cb$) collisionless shocks by using a standard particle-in-cell code.^{22} 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*_{x1}*B*_{y01}/*c* (=*u*_{x2}*B*_{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,^{25} 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 $\beta 1=2(vte1/c)2(\omega pe1/\omega ce1)2$ and $ux1/VA1=(mi/me)1/2(vte1/c)(\omega pe1/\omega 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} = *β*_{1}, 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\u2261me\omega pe12/4\pi 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,^{26} assuming *T*_{i2}/*T*_{e2} = 8.0.

. | Upstream parameters . | Results . | ||||||
---|---|---|---|---|---|---|---|---|

Run . | u_{x1}/v_{te1}
. | ω_{pe1}/ω_{ce1}
. | β_{1}
. | u_{x1}/V_{A1}
. | M_{A}
. | M_{A}^{2}/β_{1}
. | $\u27e8Ti\u22a5max\u27e9/Ti1$ . | $Ti\u22a5max/Ti1$ . |

A | 1.875 | 2 | 0.08 | 6 | 5.2 | 334 | 183 | 160–205 |

B | 0.937 5 | 4 | 0.32 | 6 | 6.5 | 132 | 63.9 | 53.9–75.4 |

C | 0.468 75 | 8 | 1.28 | 6 | 6.1 | 28.6 | 14.9 | 13.4–17.1 |

D | 1.25 | 2 | 0.08 | 4 | 4.6 | 265 | 101 | 86.5–123 |

E | 0.625 | 4 | 0.32 | 4 | 4.7 | 70.5 | 30.0 | 26.2–35.5 |

F | 0.312 5 | 8 | 1.28 | 4 | 3.9 | 12.5 | 5.01 | 4.80–5.20 |

. | Upstream parameters . | Results . | ||||||
---|---|---|---|---|---|---|---|---|

Run . | u_{x1}/v_{te1}
. | ω_{pe1}/ω_{ce1}
. | β_{1}
. | u_{x1}/V_{A1}
. | M_{A}
. | M_{A}^{2}/β_{1}
. | $\u27e8Ti\u22a5max\u27e9/Ti1$ . | $Ti\u22a5max/Ti1$ . |

A | 1.875 | 2 | 0.08 | 6 | 5.2 | 334 | 183 | 160–205 |

B | 0.937 5 | 4 | 0.32 | 6 | 6.5 | 132 | 63.9 | 53.9–75.4 |

C | 0.468 75 | 8 | 1.28 | 6 | 6.1 | 28.6 | 14.9 | 13.4–17.1 |

D | 1.25 | 2 | 0.08 | 4 | 4.6 | 265 | 101 | 86.5–123 |

E | 0.625 | 4 | 0.32 | 4 | 4.7 | 70.5 | 30.0 | 26.2–35.5 |

F | 0.312 5 | 8 | 1.28 | 4 | 3.9 | 12.5 | 5.01 | 4.80–5.20 |

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 32*l*_{i1} × 6*l*_{i1}, where $li1=c/\omega pi1=(mi/me)1/2(c/vte1)\lambda 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 *ω*_{ci1}*t* = 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 *ω*_{ci1}*t* = 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\u22a5max$, 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\u22a5max$ changes with time. Then, the time mean of $Ti\u22a5max$ for 7 ≤ *ω*_{ci1}*t* ≤ 12 is obtained. For example, we obtain the average value $Ti\u22a5max/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\u22a5max$ are obtained for other runs. The results are summarized in Table I. In Fig. 3, $Ti\u22a5max/Ti1$ is shown as a function of *M*_{A}^{2}/*β*_{1}. The simulation results seem to lie on the line, $Ti\u22a5max/Ti1\u22480.5\xd7MA2/\beta 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\gamma \u22121/(\gamma \u22121)]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\u2032$/*C*, where $vsh\u2032=ux1\u2212vsh=MAVA1$. Next, a part of them is reflected. The bulk velocity of the reflected ions in the shock rest frame is *u*_{(ref)} ≈ −$vsh\u2032$/*C*.^{4} 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*|^{2}, 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 $\u016b$ are given by

respectively. When we introduce the fraction of the reflected ions as *r* = *N*_{(ref)}/*N*_{tot}, then we get

and $\u016b$ = (1 − *r*)*u*_{(in)} + *ru*_{(ref)}. Assuming *T*_{(in)} = *T*_{(ref)} as in one of our previous simulation studies^{27} and eliminating *ū*, we can rewrite Eq. (9) as

Using Eq. (1) together with *T*_{(in)} = *CT*_{i1} and *m*_{i}$vsh\u2032$^{2}/*T*_{i1} = 2*M*_{A}^{2}/*β*_{1}, we finally obtain

The first term in the rhs of Eq. (11) is important for the case of low *M*_{A}^{2}/*β*_{1} only. The compression factor *C* is slightly larger than unity for such shocks. For large *M*_{A}^{2}/*β*_{1}, the second term dominates the rhs of Eq. (11); hence, we can explain our numerical result, $Ti\u22a5max/Ti1\u221dMA2/\beta 1$, shown in Fig. 3, if the factor 8*r*(1 − *r*)/*C*^{2} in Eq. (11) hardly depends on *M*_{A}^{2}/*β*_{1}. 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 8*r*(1 − *r*)/*C*^{2} = 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*)^{2} goes not only to *T*_{i⊥} but also to the thermal energy of reflected ions and waves excited in the shock transition layer.^{28} 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 = (2*T*_{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}^{2}/*β*_{1} 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.^{7} 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/\beta $ 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.).