Mildly relativistic magnetized shocks in electron-ion plasmas I. Electromagnetic shock structure

Mildly relativistic shocks in magnetized electron-ion plasmas are investigated with 2D kinetic particle-in-cell simulations of unprecedentedly high resolution and large scale for conditions that may be found at internal shocks in blazar cores. Ion-scale effects cause corrugations along the shock surface whose properties somewhat depend on the configuration of the mean perpendicular magnetic field, that is either in or out of the simulation plane. We show that the synchrotron maser instability persists to operate in mildly relativistic shocks in agreement with theoretical predictions and produces coherent emission of upstream-propagating electromagnetic waves. Shock front ripples are excited in both mean-field configurations and they engender effective wave amplification. The interaction of these waves with upstream plasma generates electrostatic wakefields.


INTRODUCTION
The origin of energetic particles is a long-standing problem of major importance in astrophysics. While it is widely assumed that cosmic rays (CRs) with energies up to ∼ 10 15 eV are produced at non-relativistic shocks of Galactic supernova remnants, higher-energy particles, in particular the so-called ultra-high-energy cosmic rays (UHECRs) with energies above ∼ 10 18 eV, are presumably generated in extragalactic systems with relativistic plasma outflows -active galactic nuclei (AGN) and/or gamma-ray bursts (GRBs). Non-thermal synchrotron and inverse Compton emission in blazar jets extends in broad energy range from radio up to TeV γ rays, indicating the presence of ultrarelativistic electrons. Recently established possible association of one of the high-energy neutrino sources with a flaring blazar TXS 0506+056 (Aartsen et al. 2018) shows that also CR hadrons can be produced in AGN.
High-energy particles in AGN and GRBs are often assumed to be accelerated at shock waves associated with the jets. These shocks have Lorentz factors, γ sh , ranging from mildlyrelativistic to ultrarelativistic values. Many such systems are magnetized, exhibiting inherently quasi-perpendicular and E-mail: jacek.niemiec@ifj.edu.pl superluminal conditions. Superluminal shocks are mediated by magnetic reflection of the incoming flow off the shockcompressed magnetic field (e.g. Langdon, Arons & Max 1988;Gallant et al. 1992;Hoshino et al. 1992). Coherent gyration of particles at the shock front breaks up in bunches of charge and triggers the synchrotron maser instability (SMI), which excites large-amplitude electromagnetic waves of the extraordinary mode (X-mode) that can escape towards the upstream region. This precursor wave emission has been confirmed through one-dimensional (1D) (e.g. Langdon, Arons & Max 1988;Hoshino & Arons 1991;Gallant et al. 1992;Hoshino et al. 1992;Amato & Arons 2006;Plotnikov & Sironi 2019) and two-dimensional (2D) (e.g. Sironi & Spitkovsky 2009Iwamoto et al. 2017Iwamoto et al. , 2018Plotnikov, Grassi & Grech 2018;Iwamoto et al. 2019) PIC simulations. In the electronion plasmas, interactions of the incoming electrons with the precursor waves can also generate large-amplitude longitudinal electrostatic oscillations, so-called wakefield (Lyubarsky 2006). As demonstrated by Hoshino (2008), a large-amplitude coherent electromagnetic wave propagating in the plasma can expel electrons in front of the wave packet and so induces a longitudinal polarization electric field. Electron expulsion results because the so-called ponderomotive force is proportional to the gradient of the wave pressure and acts much stronger on electrons than ions. The electric field excites longitudinal electron motions that lead to the electrostatic Langmuir waves. The formation of large-amplitude wakefields results from the parametric decay instability (PDI; e.g. Kruer 1988). In this wave-wave interaction the large-amplitude electromagnetic (pump) wave decays into a Langmuir wave and a scattered electromagnetic (light) wave. If the pump-wave frequency is much larger than the plasma frequency, Forward Raman Scattering (FRS) is triggered, in which the scattered electromagnetic wave and the Langmuir wave propagate in the same direction as the pump wave. The wavelength of the Langmuir wave is close to the electron inertial length, and its phase velocity approaches the group velocity of the pump wave, that is close to the speed of light. Electrons and ions can be energized to very high energies in a manner analogous to wakefield acceleration (WFA) during the nonlinear collapse of the Langmuir waves (Hoshino 2008). WFA was first proposed in laboratory plasmas (Tajima & Dawson 1979) and later applied to UHECR acceleration (e.g. Chen, Tajima & Takahashi 2002). It was then demonstrated through laser plasma experiments and simulations (e.g. Kuramitsu et al. 2008) that the WFA produces power-law energy spectra with a spectral index of 2.
Relativistic magnetized shocks have recently been studied with 2D PIC simulations for the case of pair plasmas (Iwamoto et al. 2017(Iwamoto et al. , 2018Sironi & Spitkovsky 2009;Plotnikov, Grassi & Grech 2018), electron-ion plasmas (Sironi & Spitkovsky 2011;Stockem et al. 2012;Iwamoto et al. 2019) and also mixed-composition plasmas (Stockem et al. 2012). Iwamoto et al. (2017) demonstrated that simulations need to have high numerical resolution to capture the precursor waves, in which case coherent waves persist even in weakly magnetized plasmas, dominated by the relativistic Weibel instability (e.g. Kato & Takabe 2010;Sironi & Spitkovsky 2011). In pair plasmas, the precursor wave amplitudes were found to be systematically smaller in 2D simulations than in the 1D case, but are still sufficient to disturb the upstream medium. 2D simulations with magnetic field in the simulation plane showed that also ordinary mode (O-mode) waves are excited, which at low magnetizations are amplified by the Weibel instability (Iwamoto et al. 2018). The amplitudes in pair plasmas are in general much smaller than at ion-electron shocks (Iwamoto et al. 2019). In conditions of high electron magnetization the wave energy exceeds that in pair plasmas by almost two orders of magnitude, and the 2D amplitude is close to the 1D level. This amplification at high-γ sh shocks is attributed to a positive feedback process associated with the ion-electron coupling through the induced wakefields. In the turbulent wakefields close to the shock the electrons can be efficiently heated so that the energy equipartition between electrons and ions may be achieved before the flow arrives at the shock front. At the same time non-thermal electrons and ions can be generated.
Most published studies address ultra-relativistic shocks with Lorentz factors γ sh ≥ 10. The mildly relativistic regime, γ sh ≈ 3, has been explored only with low-resolution studies which for superluminal shocks show very weak (Sironi & Spitkovsky 2011) or no wakefield (Lyubarsky 2006). It has been estimated that only for electron-ion shocks with γ sh (mi/me) 1/3 , where mi/me is the ion-to-electron mass ratio, the electrons will form ring-like phase-space distribution unstable to SMI. If so, one would expect little elec- tron energization upstream of the shock in blazar jets, which has important consequences for their synchrotron and inverse Compton emission (e.g. Sikora et al. 2013).
Here we revisit the efficiency of WFA and the level of the electron-proton coupling at mildly-relativistic magnetized shocks in electron-ion plasma with unprecedentedly highresolution 2D PIC simulations. We also account for ion-scale corrugations of the shock surface by employing a very large computational box. We study strictly perpendicular shocks, in which the strength of the precursor wave is expected to be largest for all superluminal obliquities (Lyubarsky 2006;Sironi & Spitkovsky 2011). We assume a shock Lorentz factor of γ sh 3 and the plasma magnetization (the ratio of the Poynting flux to the kinetic energy flux) σ = 0.1. These values are in the range of those expected for internal shocks in AGN jets. In this first paper we discuss the shock structure and the generation of plasma instabilities and waves. In a forthcoming publication (Ligorini et al., in preparation, Paper II) we present the particle acceleration and heating mechanisms and discuss the energy transfer from ions to electrons downstream of the shock. Section 2 presents the simulation setup. Section 3 shows results for the out-of-plane field orientation, which are compared to the case with the in-plane magnetic field in Section 4. Section 5 presents a summary and conclusions of this first part of our study.

SIMULATION SETUP
We use a modified version of the relativistic electromagnetic PIC code TRISTAN (Buneman 1993) with MPI-based parallelization (Niemiec et al. 2008) and the option to trace individual particles. The simulation setup is shown in Fig. 1. An electron-ion beam flows with speed v0 in the negative xdirection. It bounces off a reflective wall at the left side of the simulation box and collides with the incoming flow to form a shock propagating in the positive x-direction.
In our 2D3V simulations we use a two-dimensional spatial grid but follow three components of particle momenta and electromagnetic fields. The beam carries a large-scale homogeneous magnetic field, B0, oriented perpendicular to the shock normal, and the associated motional electric field, E0 = −v0 × B0. We study two configurations of the largescale field with respect to the simulation plane, described by the angle ϕB (see Fig. 1): the out-of-plane magnetic field, B0 = B0zẑ and ϕB = 90 o and the in-plane setup with B0 = B0yŷ and ϕB = 0 o .
The beam Lorentz factor, γ0 = 2.03, results in the shock Lorentz factor γ sh 3.3 in the upstream rest frame. The total plasma magnetization, σ = 0.1, is written with simulation-frame magnetic-field strength, B0, and ion density, Ni, as σ = B 2 0 /(µ0Ni(me + mi)γ0c 2 ), where c is the speed of light, µ0 is the permeability of free space, me and mi are the electron and ion mass, respectively . The reduced ionto-electron mass ratio, mi/me = 50, determines the electron and ion magnetizations, σe 5.1 and σi σ, through 1/σ = 1/σe + 1/σi. We verified that our results do not change if a higher mass ratio of mi/me = 100 is used.
The unit of length used here is the ion skin depth, λsi = c/ωpi, where ωpi = e 2 Ni/γ0 0mi is the relativistic ion plasma frequency. Here, e is the electron charge, and 0 is the vacuum permittivity. Time is expressed in units of the upstream ion cyclotron frequency Ωci = (eB0)/(miγ0). We ran our 2D simulations up to tmax = 84.3 Ω −1 ci and complementary 1D simulations reach tmax = 163.1 Ω −1 ci . The time-step is δt = 1/1131 ω −1 pi = 1/3556.8Ω −1 ci . Iwamoto et al. (2017) noted that numerical investigations of magnetized shocks require high resolution, otherwise the precursor waves may be artificially damped. Based on extensive tests described in Appendix A, we set the grid resolution to λse = 80∆, where λse = me/miλsi is the electron skin depth and ∆ is the size of the grid cells. The corresponding ion skin depth is λsi 566∆. This resolution is twice larger than that adopted in Iwamoto et al. (2017Iwamoto et al. ( , 2018. This unprecedentedly high resolution allowed us to detect precursor waves in the mildly-relativistic regime, that were invisible with lower resolution (e.g. Sironi & Spitkovsky 2011). Since our convergence tests show no dependence of the results on the number of particles per cell, Nppc, here we use Nppc = 10 per particle species.
Relativistic shock simulations are extremely prone to the numerical Cherenkov instability that artificially heats and slows down the plasma beam (Yee 1966;Birdsall & Langdon 1991;Hockney & Eastwood 1981). We minimize these unphysical effects, by using Friedman filters, a fourth-order accurate FTFD field-pusher (Greenwood et al. 2004), and also by injecting cold plasma. This numerical model also stabilizes the beam against the so-called finite-grid instability arising from an unresolved Debye length. The performance of the model has been extensively verified through test simulations. Our moving injection layer makes sure that the simulated plasma contains all particles and waves propagating upstream, while we minimize the propagation time of the unperturbed beam.
The transverse size of our simulation box is Ly = 5760∆ 10λsi, enough to capture structures in the shock surface with a characteristic length of several λsi. The box length, Lx, increases during simulations and reaches a final size of Lx = 160000∆ 283λsi.

SHOCKS WITH OUT-OF-PLANE MAGNETIC FIELD
A mildly relativistic strictly perpendicular shock in ionelectron plasma with large-scale magnetic field pointing out of the simulation plane is followed up to tmaxΩci = 84.3. Unlike for highly relativistic flows (e.g. Sironi & Spitkovsky 2011;Iwamoto et al. 2019), the mildly relativistic shock is not laminar. At tΩci ∼ 6.5 it develops corrugations visible in both the density and the electromagnetic field, that are fully developed at tΩci ∼ 8.5. In Section 3.1 we first present the , and the magnetic field fluctuations, δBz (c). Logarithmic scaling is applied, which is sign-preserving for electromagnetic fields (e.g. sgn(Bz) · {2 + log[max(|Bz|/B 0 , 10 −2 )]}), so that field amplitudes below 10 −2 B 0 are not resolved. Panel d) shows the transversely averaged profile of the electric field, Ex , and panel e) displays the profile of δBz taken along y/λ si = 6.
structure of the semi-laminar shock at tΩci 7.5 to demonstrate that the SMI already operates at this early stage in line with theoretical expectations. In Section 3.2 we discuss the fully-evolved rippled shock.

Semi-laminar shock stage
Fig. 2 displays the initial development of the shock front, then located at x 11λsi. The density compression by a factor around 3 is the theoretically expected value, κ = 3.18, for relativistic plasma with adiabatic index Γ = 2 . Upstream of the shock, at x 12λsi, one can see Xmode waves as plane-wave fluctuations in Bz that move with the speed of light and have a wave vector kBz = kBz,xx. The tip of the waves is at x ≈ 23λsi, which is the light travel distance ct from the reflective wall for tΩci = 7.5. One can also Figure 3. Fourier power spectra for Bz and Ex at t = 7.5Ω −1 ci , calculated upstream of the shock in the region x/λ si = 13 − 18 (compare Fig. 2). The solid white line represents the theoretical cutoff of the precursor waves. see co-moving longitudinal fluctuations in Ex at longer wavelength (Figs. 2(b) and 2(d)) with the same phase velocity. The normalized amplitude of these electrostatic waves averaged over the three oscillations observed, Ex/B0c 1.8·10 −2 , is a factor ten smaller than that of the X-mode waves. Note, that already at this very early stage the shock surface is perturbed.
The emission of X-mode waves indicates the operation of SMI at the shock (e.g. Hoshino 2008;Iwamoto et al. 2017Iwamoto et al. , 2018. We calculated Fourier spectra of Bz and Ex for a region upstream of the shock at x/λsi = 13 − 18 (see Fig. 2). The waves localized in the x/λsi ∼ 20 − 23 region were emitted during the initial beam reflection off the conducting wall, when the shock had not yet formed. They are heavily affected by the initial conditions and hence not considered in our analysis. The X-mode waves can reach the precursor only, if the x-component of their group velocity is faster than that of the shock, which imposes a limit to the wave vector. Fig. 3 demonstrates that most of the wave power is indeed observed at a kx larger than that limit, supporting the association of the X-mode waves with SMI at the shock.
The interaction of the precursor waves with the magnetized electron-proton plasma upstream should lead to electrostatic wakefield fluctuations that are evident at kEx,xλsi 2.5 − 3.0 in the power spectrum of Ex in Fig. 3(b). The signal at kEx,xλsi ∼ 1 − 4 and kEx,yλsi ∼ 3 − 20 is due to filamentation and discussed in detail in Section 3.2.3.
In Appendix B we derive the expected frequency of the Xmode waves as ω /Ωce 1, where Ωce 2.25ωpe is the electron cyclotron frequency and the prime denotes a quantity measured in the upstream frame. Then, wakefield generation by Raman scattering should yield ω L ωpe and k L 1/λse (Kruer 1988;Hoshino 2008). In the downstream (simulation) frame kL,y = k L,y and where we inserted our parameters to derive the last expression. Despite the poor wavenumber sampling of the signal in Fig. 2, the match is reasonable.

Precursor waves
Fig. 4 demonstrates that at time tΩci = 56.2 the magneticfield fluctuations extend to the far upstream. There, at x/λsi 148, one can only find waves that have been emitted very early when the shock was semi-laminar shock and have a very large x-component of their group speed. Hence the precursor waves retain their plane-wave character. Behind this region, closer to the shock, the waves also have an oblique component. Near the shock and up to x/λsi ∼ 120 (see also Fig. 7) these oblique waves form a quasi-regular pattern of oblique stripes. The emergence of the oblique wave component and formation of oblique stripes are related to the ripples in the shock surface, as we demonstrate below. A Fourier-Laplace analysis in two selected regions of the shock precursor confirms that the waves in the upstream region are X modes, presumably generated through SMI. The regions are stationary in the upstream plasma rest frame, and their initial location in the simulation box is marked in Fig. 4. The electromagnetic field is correspondingly transformed to the upstream frame. The Fourier-Laplace power spectra are shown in Figs. 5 and 6. They can be compared with the theoretical dispersion relation for the electron SMI discussed in Appendix B. Since the group velocity of the waves emitted at the shock is small except for the wave numbers near the light mode, the majority of the waves cannot outrun the shock and propagate upstream, and hence the dispersion relation is indicated with white dots only along ω = k c in Figs. 5 and 6. The low-frequency modes of SMI generated by ions are mostly subluminous and do not propagate ahead of the shock. In any case, they would not be detectable with our time window, tΩci 56.23 to tΩci 56.79, and the sampling every 10 time steps. Fig. 5 demonstrates that in Region 1 far ahead of the shock the observed signal very well matches the theoretical dispersion relation for the electron SMI. In particular, the wave power is mostly along the light mode, ω = k c, and a few harmonic modes exist for a wide wave vector range. The signal between the harmonics arises from fluctuations in the shockcompressed magnetic field. The power spectrum in Region 2, shown in Fig. 6, is heavily influenced by shock rippling and . Map of normalized magnetic-field fluctuations, δBz, at time tΩ ci = 56.2. The shock is located at x/λ si 84. Sign-preserving logarithmic scaling is applied (see Fig. 2). Regions 1 and 2 highlighted with blue squares mark the initial positions of the regions chosen for the Fourier-Laplace spectra presented in Figs. 5 and 6. Note that we show only part of the precursor, that extends up to x/λ si ≈ 177 ≈ ct.  also by the non-linear evolution of the wave modes, but retains qualitative agreement with the electron maser model. Between Region 1 and 2 one finds a slow transition from parallel to oblique modes that represents a spatial mapping of the temporal development of shock rippling.

Effects of shock rippling on wave properties
The shock ripples visible in Figs. 4 and 7 propagate in −ydirection with average speed v rippl ≈ 0.8c, commensurate with that of ions gyrating at the shock. Their mean separation along the shock surface, λ rippl 3.3λsi 2 rg, and extension in x-direction, 1.6λsi rg, reflect the ion gyroradius, rg. We associate the shock ripples with the modulation of shock-reflected ions along the shock surface first described by Burgess & Scholer (2007) for low-Mach-number nonrelativistic shocks. The instability occurs only in simulations with out-of-plane magnetic field, for which the ions gyrate in the simulation plane. In contrast, parallel-propagating waves driven by ion temperature anisotropy are frequently observed with in-plane field in the regime of low Mach numbers (e.g. Winske & Quest 1988;Umeda et al. 2014), but also in studies of high-Mach-number perpendicular nonrelativistic shocks (Wieland et al. 2016) and ultrarelativistic perpendicular shocks with moderate magnetizations, σ 0.1 (Sironi, Spitkovsky & Arons 2013).
In our simulation with ϕB = 90 • the shock ripples quickly grow from small-scale fluctuations to a long-wave mode visible in Figs. 4 and 7, in particular in the ion density. Their structure is highly dynamic on time-scales shorter than Ω −1 ci and driven by the magnetic-field compression and charge separation induced by the different inertia of electrons and ions. Arcs of increased magnetic field, electron density, and associated electric field are generated. The maps of Bz and Ex in Figs. 7(c-d) suggest that these arcs are the origin of the observed pattern of oblique waves.
The oblique structure of the precursor waves result from relativistic retardation and aberration of light and precursorwave emission in a direction normal to the local front of the arcs. The arcs move with v rippl close to c in the negative y-direction. Retardation of emission in x-direction in the simulation frame gives rise to the stripes of high wave intensity seemingly oriented at an angle of ϑ ≈ 37 • with the y-axis. Aberration provides a large x-component of the wave speed, so that the precursor waves can outpace the shock. The emission normal to the arc front results from phase bunching of the electron distribution (Hoshino & Arons 1991;Sprangle, Granatstein & Drobot 1977) which requires that the fre- quency of the wave be slightly higher than the plasma cyclotron frequency. Then, the particles on average gyrate less than 2π in a wave period and slip behind the waves. After a number of wave periods their distribution is bunched in gyrophase. The wave emission thus will be defined by the structure of the compressed magnetic field at the arcs, in which the electrons gyrate. The combined effects of the gyrophase bunching, retardation, and aberration cause the direction of precursor wave emission to vary with the evolving shape of the ripples, leading to a wide range of angles, as visible in Figs. 4 and 7(c-d) and apparent in 2D Fourier power spectra of fluctuations in Bz and Ex shown in Figs. 8(a-b) for waves in Region B marked in Fig. 7. The dominant emission pattern, however, comes from the average ripple profile and is compatible with the ripple scale-length. Fig. 7f shows the averaged and smoothed profile of the Bz magnetic-field fluctuations taken along an oblique direction, as marked with a blue parallelogram in panel (c). Highintensity patches are spaced every (2−3.5)λsi, consistent with λ rippl . Similar wave profiles are observed in Ex and Ey (not shown). Short-wavelength fluctuations in Bz have associated electric field in the x − y plane, approximately perpendicular to the wave vector, which suggests they are X-mode waves. Their spectra have cutoffs, indicated as white lines in Fig. 8(ab), that arise from the requirement that the waves be faster than the shock. While the peculiar structure of the precursor waves in our mildly relativistic shock results from shock rippling, the emission mechanism appears to be generic and corresponds to the well-known electron synchrotron maser.
In addition to the dominant oblique component, the power spectrum of Bz oscillations in Fig. 8(a) shows parallel waves with kBz,xλsi (10 − 30) that have an electric counterpart in Ey, not in Ex. The wakefield in Ex has a wave number kEx,xλsi 2 and wide distribution in kEx,y reflecting the entire range of obliquity of the precursor waves. The wave number of the wakefield agrees with that expected in the standard electron SMI scenario, as estimated in equation 1. It non-linearly couples to magnetic-field and density perturbations in the same wave band (Figs. 8a and 8c).

Filamentation via parametric instability
The Fourier spectrum of electron density in Fig. 8(c) shows significant wave power at kyλsi ∼ 10 − 30 and kxλsi ∼ 2 − 4. The corresponding density perturbations are highlighted in Fig. 7(e). They form oblique filamentary structures whose transverse scale is a few λse. We interpret these perturbations as result of the parametric filamentation instability (Kaw, Schmidt & Wilcox 1973;Drake et al. 1974) triggered when intense electromagnetic waves interact with the incoming upstream plasma. Similar filaments in density and magnetic field have been recently identified in high-resolution studies of ultrarelativistic magnetized pair shocks (Iwamoto et al. 2017(Iwamoto et al. , 2018Plotnikov, Grassi & Grech 2018) and electronion shocks (Iwamoto et al. 2019). Their presence indicates coherence and self-focusing of the precursor waves in 2D systems. The filaments observed in pair plasmas largely retain their structure during advection toward the shock. In electron-ion plasma, instead, the filaments quickly merge to form long, ion-scale turbulent structures ahead of the shock. At our mildly-relativistic shock the filaments resemble those at pair shocks. They are observable very far upstream, up to x/λsi ∼ 200, but their structure is disrupted by the oblique waves, and their amplitude is only δNe/N0 ≈ 0.1. The corresponding spectral signature in the magnetic field is very weak (compare Fig. 8(a)). Although at mildly relativistic shocks the precursor waves are less prominent than in the ultrarelativistic case, and the parametric instability is weakly driven, our high-resolution simulations can still detect coherent precursor-wave emission. Fig. 9 shows profiles in the precursor region of δBz, taken at y/λsi = 6, and Ex y-averaged to filter out the contribution of the oblique large-amplitude precursor waves. These profiles can be compared to corresponding profiles in 1D simulation that we describe in Appendix C. The magnetic-field fluctuation amplitude normalized to the upstream field strength, δB/B0, and energy density normalized to the upstream electron kinetic energy, pe = δB 2 /(2µ0γ0Nemec 2 ), are listed in Table 1 in comparison with results of other studies. For out-of-plane magnetic field the only relevant polarization is δB = δBz. The amplitudes are averaged in the region of x/λsi = 129 − 134, located about 2λsi from the shock. Our 1D test simulation is described in detail in Appendix C. The amplitude of the precursor waves is comparable in 2D and 1D simulations. In 2D, the average wave amplitude out to x/λsi ≈ 230 is even larger than the strongly coherent oscillations at the tip of the precursor, that were generated in the early linear phase (see Fig. 4). Fig. B2 demonstrates that in 1D the waves closer to the shock are weaker, possibly due to heating of the electrons that suppresses higher harmonics in the electron SMI (Amato & Arons 2006). In a 2D simulation the same should happen, and inhomogeneities at the shock may further cause a reduction of wave coherency. But we observe that shock rippling is highly organized and produces a semi-coherent, modulated train of oblique precursor waves. Thus, instead of being destructive, the ripples amplify the precursor-wave amplitude.

Precursor wave amplitudes
The magnetic-field amplitude can be compared to that observed at ultra-relativistic shocks. Since electron magnetization is a relevant parameter, in Table 1 we list available results for shocks with σe = 5. As expected, both δB/B0 and pe are smaller than those in 1D simulations of pair shocks ), but they are much larger than those obtained in the high-resolution 1D and 2D simulations of pair shocks by Iwamoto et al. (2019). The wave energy is much smaller than that at ion-electron shocks with γ0 = 40 (Iwamoto et al. 2019), which both in 1D and 2D exceeds that in pair plasmas by almost two orders of magnitude, and the 2D amplitude is close to that in 1D. The high wave intensity at high-γ sh ionelectron shocks was attributed to the so-called positive feedback, in which incoming electrons accelerated by the wakefield cause enhanced precursor-wave emission, that in turn amplifies the wakefield. At the mildly relativistic shocks described here, the wakefield does not reach a very high amplitude (see Section 3.2.5), and the positive feedback is not effective. However, electromagnetic precursor wave amplification up to the 1D level is achieved through shock rippling.

Wakefield waves
As noted before, the ponderomotive force on the upstream plasma leads to longitudinal plasma motions and associated electrostatic wakefield, whose average amplitude does not exceed Ex /(B0c) ∼ 0.01. The weakness of the wakefield reflects the relatively low amplitude of the precursor waves, that can be expressed through the so-called strength parameter a that can be estimated as (Iwamoto et al. 2017) a γ0 √ σe ωpe ω δB B0 .
The average magnetic-field amplitude is δBz/B0 0.19 (see Table 1), and all together we find for the strength parameter a ≈ (0.14 − 0.35). (3) The corresponding amplitude of the wakefield can be estimated following Hoshino (2008), using ξ = 1/2 as for a linearly polarized wave: Ex B0c in agreement with our simulation result. Shock rippling causes enhanced emission of the precursor waves at oblique angles which in turn produces oblique Langmuir waves. They are averaged out of the wakefield profile shown in Fig. 9(a), and so the local wave amplitude may be much larger than the mean amplitude. In fact, from time tΩci ∼ 40 on, when the oblique precursor-wave structure is well established, episodes of stronger semi-coherent wave emission from the shock lead to stronger wakefield in the near-upstream region, that non-linearly evolves. Fig. 10 shows a stack plot of averaged wakefield profiles for a time period 28.1 Ω −1 ci , starting at t0Ωci = 56.2. The profiles are given in shock-centered coordinates, x − x shock . Far upstream of the shock, the electrostatic waves propagate away from the shock, but within 70 λsi of the shock the wakefield on average moves back toward the shock. We interpret this downstream-directed motion of the wakefield as result of Forward Raman Scattering (FRS) operating at our mildly relativistic shock.
The ponderomotive force is proportional to the gradient of the wave pressure and can also act inside the precursor, if the electromagnetic waves are modulated (Hoshino 2008). The enhanced emission of the precursor waves through shock rippling amplifies the waves and triggers the nonlinear FRS. In this process the scattered electromagnetic waves successively decay into other electromagnetic waves and Langmuir waves. As the frequency of the scattered wave is lower than that of the pump wave, broadband precursor wave spectra extending from the initial ω Ωce down to ω > ωpe are generated. Similarly, broadband Langmuir waves are produced with k L < ωpe/c and ω L ωpe (Hoshino 2008). In the upstream plasma rest frame the electromagnetic and Langmuir waves all have phase velocities in the upstream direction, but in the simulation frame part of these waves move toward downstream.
The presence of large-amplitude wakefield propagating toward the shock is of importance for electron energization upstream of the shock. We will show in Paper II that the waves can scatter electrons and boost them toward the shock, contributing to ion-to-electron energy transfer in the precursor region.

COMPARISON WITH THE IN-PLANE SETUP
In this section we compare the electromagnetic structure of a mildly relativistic shock with upstream magnetic field lying in the plane of the simulation, B0 = Byŷ, thus ϕ = 0 o , with that for out-of-plane field discussed in Section 3. With in-plane magnetic field the shock quickly acquires its steadystate form, and so we show the structure and discuss its properties only at time tΩci = 84.3.

Shock front and downstream turbulence
To be noted from Fig. 11 is the lower shock speed, v = 0.41c, compared to the out-of-plane case, that is implied by its location at x/λsi 108. The density compression is κ 3.2. Both are consistent with a shock in plasma with three degrees of freedom and adiabatic index Γ = 4/3 (Plotnikov, Grassi & Grech 2018). Ion gyration at the shock happens in the x−z plane, which suppresses the gyration-driven rippling mode seen with out-of-plane magnetic field.
Fluctuations in density and electromagnetic field can be observed together with corrugations in the shock surface, that develop very early in the simulation and quickly evolve into a large-scale rippling mode with λ rippl 5λsi and amplitude λsi. They propagate along the mean magnetic field and are most likely driven by the anisotropy in the ion temperature that results from ion reflection from the shock and was geometrically suppressed in the out-of-plane simulation. Fig. 11(e) demonstrates that at the shock T i ⊥ T i with respect to the mean magnetic field direction, which triggers the Alfvén Ion Cyclotron (AIC) instability that is known to produce ripples in low-Mach-number shocks (e.g. Winske & Quest 1988;Umeda et al. 2014) and can generate magnetic-field fluctuations at the front of relativistic pair shocks (Iwamoto et al. 2018). Shock-front corrugations are also a source of downstream turbulence, through the mechanism of the vorticity generation via a process similar to the Richtmyer-Meshkov instability (e.g. Mizuno et al. 2011Mizuno et al. , 2014.

Upstream waves
Short-scale precursor waves and large-scale electrostatic wakefield are evident in Fig. 11 and in the profiles covering the entire upstream region at time tΩci = 84.3, that we plot in Fig. 12 in the same manner as earlier for 2D out-ofplane simulation and in Appendix C for a 1D test run. Fig. 13 shows an enlarged view of the region x/λsi = (120−125), presenting also electric-field fluctuations. Electromagnetic waves with By fluctuations parallel to the large-scale magnetic field, Figure 13. Enlarged view of the region x/λ si = 120 − 125 with magnetic-field and the electric-field profiles plotted in black and red, respectively, and normalized as in Fig. 12. The phase (anti-)correlation between y and z field components indicate X-mode (a) and O-mode (b) waves. B0 = Byŷ, and anticorrelated with Ez fluctuations are Xmode waves. Correlated Bz and Ey fluctuations are O-mode waves. Oscillations in Bx result from an oblique propagation of the X-mode wave, as we discuss below.
The magnetic-field fluctuation amplitudes are listed in Table 1. The time evolution of the field amplitudes is shown in Fig. C1 in Appendix C. We list the total amplitude of the magnetic field oscillations, δB = δB 2 x + δB 2 y + δB 2 z , not differentiating between the X-mode and O-mode waves. The amplitude of the X-mode wave, δBX /B0 = δB 2 x + δB 2 y /B0 0.12, is larger than the amplitude of the O-mode wave, δBO/B0 = δBz/B0 0.08. The total precursor wave amplitude, δB/B0 0.15, is slightly smaller than the amplitudes obtained in our 2D run with ϕB = 90 o and in the 1D simulation.
With out-of-plane magnetic field strong shock ripples increase the precursor-wave amplitude to the level observed in 1D simulation. With in-plane field we see a similar amplification. Fig. 11(b) shows the emission of precursor waves in bunches that correspond to rippling features at the shock. The rippling driven by the AIC instability is relatively weak and cannot fully compensate losses in coherency due to random fluctuations in shock structure at the shock surface and the thermal damping of the waves. The precursor wave amplitude is large enough, though, to induce the wakefield and thus accelerate and heat particles. Fig. C1 shows a roughly constant wave amplitude throughout the 2D simulation with the in-plane setup, which is due to an early formation of the shock ripples.
The linear theory of the SMI predicts X-mode emission at a level above that of O-mode waves (Melrose, Hewitt & Dulk 1984;Lee, Kan & Wu 1980;Wu & Lee 1979). In earlier 2D simulations of ultrarelativistic pair-plasma shocks with in-plane magnetic field the O-mode energy was observed to exceed that in X modes at very small electron magnetizations, σe 10 −2 (Iwamoto et al. 2018), which was attributed to mode conversion from X to O modes at the turbulent shock transition. Early in the shock evolution, charged particles gyrate in the x − z plane and emit X-mode waves with δBy. As plasma instabilities generate fluctuations in Bz of amplitude comparable to B0, the net magnetic field undulates in the y − z plane, and the X-mode waves carry both δBy and δBz fluctuations. Upon transmission to the upstream region the δBz components may be converted into O-mode waves. Indeed, Fig. 12 shows the tip of the Bz turbulence behind that of By, indicating that the O-mode waves are produced a bit later than the X-mode waves, after the shock front has developed substantial turbulence. The slightly smaller amplitude of the O-mode waves with respect to the X-mode waves reflects the moderate level of Bz fluctuations in the shock front, δBz/B0 ≈ 1 (compare results for σe 10 −2 in Iwamoto et al. (2018)). Fig. 14 shows Fourier power spectra of By, Bz, Ex, and the electron density in Region B marked in Fig. 11. Most of the wave power in δBy and δBz is located to the right of the theoretical cutoff calculated in Appendix B, and the waves are faster than the shock. The wave vector range of the precursor waves is similar to that in the out-of-plane setup, and an oblique component is likewise present. Coherent precursor waves are persistent in mildly relativistic shocks, albeit with smaller amplitude than in the ultrarelativistic regime.
Also with B0 in the simulation plane we observe transverse density filaments upstream of the shock, whose amplitude, δNe/N0 ≈ 0.5, is much larger than in the out-of-plane simulation. They are better aligned with the x-direction, though. It appears that the parametric filamentation instability is not affected by the weak ripples at the shock.
The wakefield is evident upstream of the shock, as for ϕB = 90 o . Fig. 14(c) shows a signal at kEx,x 2.5λ −1 si that is marginally consistent with equation 1. The wave power is slightly less than that observed with out-of plane magnetic field. As in Section 3.2.5 we take the typical wavenumber of the precursor waves, k ≈ 2.8λ −1 se , and the typical frequency of the dominant X-mode waves, ω ≈ 3.2 ωpe, and add the contributions from the X-and O-modes to estimate the strength parameter and calculate the amplitude of the wakefield (equatio 4), Therefore, the wakefield should exert similar effects on the upstream plasma in simulations with in-plane and out-ofplane magnetic field.

SUMMARY AND CONCLUSIONS
This is the first of two articles in which we investigate mildlyrelativistic magnetized shocks in electron-ion plasma. In this paper we explore with PIC simulations the electromagnetic shock structure and the production of plasma instabilities and waves. Paper II shall be devoted to particle acceleration, heating, and the energy transfer from ions to electrons downstream of the shock.
Our high-resolution studies show that the SMI operates at mildly relativistic shocks as theoretically predicted and produces coherent emission of upstream-propagating electromagnetic waves. The waves are substantially weaker than at ultrarelativistic shocks (Iwamoto et al. 2017(Iwamoto et al. , 2018), but reach a persistent level that is similar in 2D and 1D simulations. In 2D shock corrugation provides wave amplification that compensates other destructive 2D effects. Shock ripples appear for both in-plane and out-of-plane mean magnetic field, but their generation mechanism differs -modulation of ion gyration (Burgess & Scholer 2007) for ϕB = 90 o and the AIC temperature-anisotropy instability with ϕB = 0 o . In both cases the ripples heavily influence the upstream plasma and the structure of downstream turbulence.
For out-of-plane mean magnetic field the precursor waves are of the X-mode type. Both the emission direction about 40 deg off the shock normal and the wave amplification are caused by the shock ripples. With in-plane magnetic field the AIC-generated shock-front corrugations have a slightly lower amplitude, and the waves propagate mostly along the shock normal. Magnetic turbulence at the shock causes part of the precursor waves to be O-mode waves. For both magnetic-field orientations we observe in the upstream plasma the electrostatic Langmuir modes -the wakefields -and the density filaments that result from the parametric filamentation instability. Except for brief periods, their average amplitude is moderate and smaller than at ultrarelativistic shocks.
The important role of shock rippling has not been demonstrated so far for relativistic shocks. At perpendicular high-Lorentz-factor shocks Sironi, Spitkovsky & Arons (2013) found shock corrugations consistent with Burgess & Scholer (2007) only in a limited range of plasma magnetization, 3×10 −4 σ 10 −1 , probably on account of electron heating at Weibel filaments for σ ≤ 10 −4 and in the SMI-mediated precursor for σ > 0.1. The ripples at mildly relativistic shocks may be similarly suppressed at low magnetizations due to the Weibel instability that still operates in this regime (e.g., Kato & Takabe 2010), but at σ > 0.1 one does not expect precursor wave emission stronger than for the σ = 0.1 analysed here (see Iwamoto et al. 2019), and rippling may persistently develop. The same should apply to AIC-instability-generated corrugations.
Our 2D simulations show intense coherent precursor waves generated by the SMI irrespective of the magnetic-field configuration. One should expect strong precursor waves also in 3D, that are a mixture of X-mode and O-mode waves. The intensity of the ordinary mode is difficult to estimate, because these waves arise from local variations in the gyration direction at the shock front that may be less coherent in 3D than in the in-plane 2D simulations, possibly leading to weaker O-mode emission (Iwamoto et al. 2018). Plotnikov & Sironi (2019) showed for ultra-relativistic shocks in pair plasma that at high σe 1 (σe = 5.1 in this study) the physics of precursor-wave emission in 3D is better represented with the out-of-plane 2D model. This result will likely hold in electron-ion plasma, since the SMI mechanism operates as in pair plasma. As we demonstrated, at mildly relativistic shocks the precursor-wave strength and structure is significantly affected by shock rippling. Rippling along the lines of Burgess & Scholer (2007) requires a suppression of fluctuations parallel to the mean magnetic field that might be difficult to achieve in 3D, but ripples generated through the AIC instability amplify the precursor waves to comparable amplitudes, and so precursor-wave amplification may be expected in 3D as well.
2019/33/B/ST9/02569 (J.N.). This research was supported by PLGrid Infrastructure. Numerical experiments were conducted on the Prometheus system at ACC Cyfronet AGH. This work was supported by JSPS-PAN Bilateral Joint Research Project Grant Number 180500000671. Part of the numerical work was conducted on resources provided by the North-German Supercomputing Alliance (HLRN) under projects bbp00003, bbp00014, and bbp00033.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author. Figure B1. Dispersion relation for X-mode waves produced by cold gyrating electrons and ions in plasma with electron magnetization σe = 5.1. The oscillation frequency, ωr, and growth rate, ω i , for the electron SMI are shown with red and pink dotted lines, respectively. Blue and cyan dots correspondingly represent the ion maser instability. The black line shows the lightwave dispersion relation, ω = kc. Fig. B1 shows the X-mode dispersion relation for relativistic magnetized electron-ion plasmas with parameters used in this study. The dispersion relation was derived following Hoshino & Arons (1991), but assuming that both electrons and ions form cold rings in momentum space while they gyrate about the magnetic-field lines with γ 3.3. Hoshino & Arons (1991) treated electrons as a hot background fluid with temperature γemec 2 when deriving the dispersion relation of the ion-generated SMI. The ring distribution that we use provides the same effective perpendicular temperature. In Fig. B1 we show the real part of frequency, ωr, and the imaginary part, ωi, for gyrating electrons and ions with zero bulk (drift) velocity. The k vector is perpendicular to the magnetic field. The emission spectrum has a clear harmonic structure for both the electron-generated SMI (ωr/Ωce 1) and the ion-generated SMI (ωr/Ωce 0.1). We considered the first ten harmonics of the Bessel function when calculating the dispersion relation.

APPENDIX B: LINEAR THEORY OF THE SYNCHROTRON MASER INSTABILITY IN ELECTRON-ION PLASMAS
For the electron SMI, the growth rate of the fundamental mode is comparable to that of the higher harmonics. The phase velocity of the waves is about c at the wavenumber of maximum growth. Note that the group velocity is small except near ω = kc. The frequencies and growth rates of the unstable ion-modes are lower than the electron ones by approximately the ion-to-electron mass ratio, and the growth rate slightly increases with higher harmonic number. One should thus expect two stages of the SMI to occur -first the electron SMI, then the ion maser instability. However, the ions emit mainly left-handed elliptically polarized magnetosonic waves that are mostly sub-luminous, the phase velocities of the modes are ω/kc v sh , so that the waves emitted at the shock cannot outrun it to reach the precursor.
The polarization of the electron SMI waves is that of X Figure B2. Upstream wave profiles for Ex (a) and δBz (b) in a 1D simulation at time tΩ ci = 84.7.
modes, and after transmission to the upstream medium the waves would propagate with the standard dispersion relation of X-mode waves that reads where the prime denotes quantities measured in the upstream plasma frame and η ≡ ck /ω is the refraction index. Equation B1 includes only the electron contribution and can be used for ω ≥ ωpe √ 1 + σe (Hoshino & Arons 1991). Following Iwamoto et al. (2017) we estimate the theoretical cutoff wave number, above which the waves emitted at the shock can escape upstream. Lorentz transforming equation B1 one obtains the dispersion relation in the simulation rest frame: Here, η 2 = c 2 k 2 /ω 2 is the refraction index in the simulation frame and θ is the angle between the x−axis and the wave propagation direction. The first term on the right hand side is relevant only for η 1, in which case it is approximately σ/γ 2 . We can therefore always write it in that form and derive the simplified dispersion relation ω ≈ (1 + σe γ 2 )ω 2 pe + c 2 k 2 .
Precursor waves propagate toward the shock upstream. The cutoff wave number is then determined by equating the xcomponent of the wave group velocity, vgr,x = vgr cos θ with the shock velocity, cβ sh , yielding: The cutoff wave number for O-mode waves, that we use in Fig. 14, can be estimated in an analogous way. The dispersion relation in the simulation frame is the same as for the electromagnetic wave in unmagnetized plasma (σ = 0), ω 2 = k 2 c 2 + ω 2 pe (B5) which leads to the cutoff wave number for the ordinary mode: kx = β sh γ sh ω 2 pe c 2 + k 2 y .

APPENDIX C: RESULTS FOR 1D SIMULATION
The 2D simulations presented in Sections 3 and 4 can be compared with a 1D run, in which shock rippling and obliquely Figure C1. Temporal evolution of normalized precursor wave amplitudes, δB/B 0 (red lines), and precursor wave energy normalized to the upstream electron kinetic energy, pe (green lines), for 2D out-of-plane (thick solid lines), 2D in-plane (dash-dotted-lines) and 1D simulations (dashed lines).
emitted waves are suppressed. In the ultra relativistic regime the amplitude of SMI-generated precursor waves are systematically larger in 1D than in 2D, on account of the loss in phase coherence imposed by inhomogeneity in the shock surface (Iwamoto et al. 2017(Iwamoto et al. , 2018. There is a positive feedback in electron-ion plasma, however, through which electrons accelerated in the shock upstream enhance the precursor wave emission that in turn induce stronger wakefield, accelerating the incoming electrons even more efficiently, leading to energy equipartition between electrons and ions (Lyubarsky 2006;Hoshino 2008). At ultrarelativistic shocks with high electron magnetizations, σe 1, the amplitude of the precursor waves in 2D is comparable to that in 1D (Iwamoto et al. 2019). In our 2D simulations of moderately relativistic shocks the positive feedback process is not operative, and the wave amplification may be attributed to shock rippling. A 1D test simulation is performed to evaluate these issues.
The setup of the 1D simulation is the same as that for 2D simulations, but the transverse dimension of the computational box is only 5 cells wide, making it effectively 1D. Fig. B2 shows the wave profile for the 1D simulation at time of tΩci = 84.7, which is tmax in 2D simulations. The shock speed is similar to that measured in the 2D out-of plane simulation. The fluctuations in Ey and Bz are anti-correlated, indicating that the waves are of X-mode type.
The precursor wave profiles in the 1D run can be directly compared with those obtained in 2D (see Fig. 9 and also Fig. 12). The wave amplitude is also listed in Table 1. Whereas in 2D wave amplification by the shock ripples causes a high wave amplitude near the shock, in the 1D case the electrons heated upstream of the shock reduce the later emission of the precursor waves whose amplitude is then low (Amato & Arons 2006). This behavior is also evident in the time evolution of the normalized wave amplitude, δB/B0, and wave energy, pe, that we show in Fig. C1.
One can note that the wave evolution is similar in both 2D runs, and in the 2D in-plane case the total precursor waves amplitude is only slightly smaller than the one observed for ϕB = 90 o . This shows that the shock rippling-mediated wave amplification operates with similar efficiency in both 2D setups, despite the different mechanisms. Fig. C2 shows 1D Fourier power spectra upstream of the shock, in the region x/λsi = (130 − 140). The signal band in Figure C2. 1D power spectrum for Bz (black) and Ex (red) in the region x/λ si = 130 − 140 at time tΩ ci = 84.7 (cf. Fig. B2).
the magnetic field oscillations, Bz, is consistent with the SMI precursor waves observed in 2D simulations. The electrostatic component in Ex has a wavenumber of kE x ≈ 2, consistent with the theoretical wave number for SMI-generated wakefield. This paper has been typeset from a T E X/L A T E X file prepared by the author.