Kinetic Simulations of Nonrelativistic High-mach-number Perpendicular Shocks Propagating in a Turbulent Medium

Strong nonrelativistic shocks are known to accelerate particles up to relativistic energies. However, for diffusive shock acceleration, electrons must have a highly suprathermal energy, implying the need for very efficient preacceleration. Most published studies consider shocks propagating through homogeneous plasma, which is an unrealistic assumption for astrophysical environments. Using 2D3V particle-in-cell simulations, we investigate electron acceleration and heating processes at nonrelativistic high-Mach-number shocks in electron-ion plasma with a turbulent upstream medium. For this purpose, slabs of plasma with compressive turbulence are simulated separately and then inserted into shock simulations, which require matching of the plasma slabs at the interface. Using a novel procedure of matching electromagnetic fields and currents, we perform simulations of perpendicular shocks setting different intensities of density fluctuations (≲10%) in the upstream region. The new simulation technique provides a framework for studying shocks propagating in turbulent media. We explore the impact of the fluctuations on electron heating, the dynamics of upstream electrons, and the driving of plasma instabilities. Our results indicate that while the presence of turbulence enhances variations in the upstream magnetic field, their levels remain too low to significantly influence the behavior of electrons at perpendicular shocks.


INTRODUCTION
Understanding how particles are accelerated at shock waves remains a fundamental challenge for astrophysics and space physics.The description provided by the basic process, diffusive shock acceleration (DSA; e.g., Bell 1978;Blandford & Ostriker 1978;Drury 1983;Blandford & Eichler 1987), is not complete and does not always describe the behavior at interplanetary shocks (see e.g.Lario et al. 2003).The influence of pre-existing turbulence upstream of a shock is an important question in the acceleration theory, yet is still poorly understood.It is known, that fluctuations enhance the trapping of particles near a shock, providing favorable conditions for effective acceleration (see Guo et al. 2021;Perri et al. 2022).Turbulence is ubiquitous in astrophysical environments, thus it is essential to investigate its interplay with shock waves.In this paper we study collisionless shocks in the non-relativistic high-Mach-number regime, with the sonic and Alfvénic Mach numbers M S , M A ≳ 20.Such conditions are observed mainly in supernova remnants (SNRs, Raymond et al. 2023), but also in the bow shocks of Jupiter (Slavin et al. 1985), Saturn (Sulaiman et al. 2016), and Uranus (Bagenal et al. 1987), and occasionally at the Earth's bow shock (Sundberg et al. 2017).
Since it is widely accepted that supernova remnants are plausible candidates for the production of the galactic cosmic rays, we identify these sources as primary objects where our simulations may find application.The fast shocks, v sh ≈ 1, 000 − 10, 000 km/s, form during the interaction of supernova ejecta with the ambient medium, which is turbulent.The impact of the medium inhomogeneities on the evolution of SNRs is an important aspect for understanding the observed emission morphology.Previous numerical studies have demonstrated that magnetohydrodynamic (MHD) turbulence modify the dynamical properties of these systems (Balsara et al. 2001;Zhang & Chevalier 2019;Peng et al. 2020;Bao et al. 2021) and amplify the magnetic field (Giacalone & Jokipii 2007;Guo et al. 2012).
Observations of radio synchrotron emission (e.g., Shklovsky 1954;Dubner & Giacani 2015) and nonthermal X-ray emission (e.g., Koyama et al. 1995;Vink 2012), as well as in the γ-ray waveband (e.g., Aharonian et al. 2007), indicate the presence of accelerated electrons in SNRs.The straightforward explanation for their energy gain is participation in DSA.However, thermal electrons do not satisfy the main requirement of the mechanism: their Larmor radii are much smaller than the shock width.For this reason they must undergo some pre-acceleration mechanism before being injected into DSA, which is known as the electron injection problem.Possible mechanisms and the underlying microphysics of electron acceleration have been extensively investigated (see Amano et al. 2022;Bohdan 2023, for a review).To interpret the broadband emission from supernova remnants, MHD simulations coupled with a prescription of particle acceleration have been utilized (see e.g., Orlando et al. 2021, for a review).Such studies, e.g., Ferrand et al. (2010); Orlando et al. (2012); Ferrand et al. (2014); Pavlović (2017); Pavlović et al. (2018); Brose et al. (2021), usually employing Blasi's semi-analytical model of non-linear DSA (Blasi 2002(Blasi , 2004;;Blasi et al. 2005), which treats the electron injection efficiency as a free parameter.
Previous studies of electron acceleration at the nonrelativistic shocks assumed a homogeneous upstream medium, hence all the turbulence in this region was driven only by shock-reflected particles.It is important to investigate the effect of pre-existing small-scale turbulence in the high-Mach-number regime and compare the previous simulations with new ones, where the upstream medium initially carries fluctuations.Shock waves propagate there through a turbulent medium, which apart from magnetic fluctuations carries density inhomogeneities.While the magnitude of the density fluctuations at SNRs remains unknown, measurements in the heliosphere (Carbone et al. 2021) and in the local interstellar medium (Lee & Lee 2020;Ocker et al. 2021;Fraternale et al. 2022) show amplitudes of δn e /n e ≲ 10% on minute timescales in the spacecraft frame, which roughly corresponds to the kinetic scales of electrons.
The interaction of a shock with pre-existing inhomogeneities has been intensively studied in the magnetohydrodynamics (MHD) regime (e.g.Inoue et al. 2013;Mizuno et al. 2014;Hu et al. 2022).The results show that turbulence is able to distort the shock structure: for instance, it changes the compression ratio and induces rippling of the shock surface.Furthermore, it amplifies the magnetic field via turbulent dynamo.Recent studies on ion scales via hybrid kinetic simulations demonstrate the strong influence of upstream turbulence on particle transport at shocks, as well as efficient proton acceleration (Trotta et al. 2021;Nakanotani et al. 2022).In the hybrid simulations electrons are treated as fluid.Guo & Giacalone (2010, 2012) investigated their acceleration processes with test-particle simulations and found, that pre-existing magnetic turbulence significantly improves the electron-acceleration efficiency.
To fully examine the physics of shock waves propagating in turbulent plasmas, particularly on electron kinetic scales, particle-in-cell (PIC) simulations are needed.This method describes collisionless plasma from first principles, hence it is widely used in studies of shock waves in astrophysics (e.g.Pohl et al. 2020).So far, first attempts to investigate the effects of a nonhomogeneous medium on plasma shocks via kinetic simulations have been made for relativistic shocks propagating in pair plasmas.Static large-scale variations of the upstream density significantly modify the conditions of the downstream plasma (Tomita et al. 2019) and corrugate the shock front, which leads to efficient particle acceleration (Demidem et al. 2023).With a more realistic setup, in which the upstream plasma contains turbulence instead of static density structures, Bresci et al. (2023) confirmed an increased acceleration efficiency at relativistic shocks in the case of strong electromagnetic fluctuations.However, the statistical properties of the turbulence varied throughout the simulation.
In this study, we investigate the effect of pre-existing turbulence at kinetic scales on electron acceleration at non-relativistic high-Mach-number shocks in electronion plasma using a novel technique.In our simulations the upstream medium is fully filled with turbulence, which has stable statistical properties, meaning the intensity and spectrum of the fluctuations do not significantly vary during the evolution of the system.
Besides the Mach numbers the main factor determining the physics of non-relativistic high-Mach-number shocks is the inclination of the large-scale magnetic field.Using the obliquity angle θ Bn between the magnetic field and the shock normal, shocks can be divided into per-pendicular (θ Bn = 90 • ), oblique (θ Bn ≈ 50 • − 75 • ), and quasi-parallel (θ Bn = 0 • ).In this work, we focus on wellstudied perpendicular shocks (Matsumoto et al. 2012(Matsumoto et al. , 2013(Matsumoto et al. , 2015;;Wieland et al. 2016;Bohdan et al. 2017Bohdan et al. , 2019aBohdan et al. ,b, 2020a,b),b).The perpendicular direction of the magnetic field confines the shock front to a thickness of roughly one ion gyroradius, and the turbulence driven by reflected particles is negligible.The global structure of the perpendicular shocks is highly influenced by the electromagnetic Weibel instability (Fried 1959), growing from the interaction between incoming and reflected, in fact gyrating, ions (Kato & Takabe 2010).The Weibel instability strongly amplifies the magnetic field.The model of this process provided by Bohdan et al. (2021) is consistent with in-situ measurements at Saturn's bow shock.The nonlinear evolution of the Weibel modes creates current sheets, which decay at the shock ramp through magnetic reconnection (Matsumoto et al. 2015;Bohdan et al. 2020b).The turbulence introduced by this process creates favorable conditions for electron acceleration.At the leading edge of the shock foot, the Buneman instability is driven by the interaction between cold incoming electrons and reflected ions (Buneman 1958), and electrons can undergo shock-surfing acceleration (SSA; Shimada & Hoshino 2000;Hoshino & Shimada 2002).However, simulations show that most of the non-thermal particles found in the downstream region gained their energy via stochastic Fermi acceleration (SFA; Bohdan et al. 2017Bohdan et al. , 2019b)).In this process, the maximum energy that electrons can achieve is 10-100 times higher than the thermal energy of the downstream electrons.This is insufficient for injection into DSA.Studies involving hybrid simulations (Caprioli & Spitkovsky 2014) and PIC-MHD (van Marle et al. 2022) have also indicated that perpendicular shocks are not efficient particle accelerators, except when the magnetic field is low, or seed cosmic rays (Caprioli et al. 2018) or neutral particles (Ohira 2016) are present.Magnetic reconnection, SSA and SFA, together with the shock potential and adiabatic compression, heat up electrons.
The paper begins with a description of the simulations and the used setup (Section 2).It follows with a presentation of the results in Section 3. Finally, a summary and discussion can be found in Section 4.

METHODS
In the study presented here, we use the MPIparallelized code THATMPI (Niemiec et al. 2008;Bohdan et al. 2017) that tracks 2-spatial and all 3 velocity components of individual particles in a plasma (2D3V configuration).This approximation saves computational resources in comparison with the full 3D simulations, but reasonably accurately reproduces the physics of shock waves (see, e.g., Pohl et al. 2020).
We establish the shock by reflecting an incoming plasma beam off a conducting wall.The interaction of the incoming and the reflected beam creates the shock.Since the simulation is performed in the downstream frame, the upstream plasma which flows toward the shock has to be continuously replenished.Homogeneous plasma can be added in a thin layer.However if that plasma is meant to carry pre-existing fluctuations, the turbulence has to be established separately and then injected into the shock simulation.In the following paragraphs, the method of turbulence generation is described, as well as its injection procedure into a shock simulation.Details of both are presented in the Appendix.

Turbulence driving
In the standard models of MHD turbulence the energy cascades from large scales to small scales (Sridhar & Goldreich 1994;Goldreich & Sridhar 1995).The separation of these scales makes it computationally infeasible to follow the self-consistent evolution of plasma turbulence down to kinetic scales.The Langevin antenna, a method of driving turbulence presented in Tenbarge et al. (2014), allows mimicking the MHD cascade at scales resolved by PIC simulations (Zhdankin et al. 2017;Grošelj 2019).Despite this advantage, the method does not provide a precise control of plasma conditions (they vary in time).Furthermore, turbulence has to be driven continuously in some part of a simulation box, which enforces the use of elongated boxes for simulations in the downstream frame (see Bresci et al. 2023).For these reasons, we choose the so-called decaying turbulence and inject a set of modes at the beginning of a simulation that evolve into turbulent structures.The following paragraph shows why the decaying setup is more suitable than the Langevin antenna for our case.
Our study focuses on compressive turbulence.We prefabricate slabs of turbulent plasma in dedicated simulations with a square simulation box and periodic boundaries.Their size corresponds to the width of the shock simulation into which they are to be injected.In the pre-fabrication stage density fluctuations are obtained by superposition of wave-like disturbances of the local bulk velocity (Giacalone & Jokipii 1994;Grošelj 2019); their exact form is included in Appendix A. Fluctuations of the magnetic field evolve self-consistently, but are weaker than those in density.Even though the decaying turbulence cannot achieve a steady state, its statistical properties vary very slowly during the shock simulation.After a short period of rapid evolution the magnitude of the density fluctuations is stable for at least two ion Larmor times, which means that they are sufficiently long-lived to be inserted into the shock simulation.
On kinetic scales, damping transfers wave energy into plasma heat.Thus, we have to launch the prefabrication of plasma slabs with a low initial temperature and make sure that the temperature at the time of insertion into the shock simulation corresponds to the desired high Mach number.The initial heating constrains the achievable amplitude of density fluctuations to δn/n ≈ 10%, where δn is the root-mean-square of the particle density fluctuations (on the scale of a quarter of ion skin length), and n is the mean density.

Turbulence injection
For efficiency and for maintaining a similar evolution state of the turbulent upstream plasma we regularly add the pre-fabricated plasma slabs to our shock simulation.Each new slab must be matched to the adjacent plasma (at the end of the shock simulation box) to prevent artefacts and transients.Here we briefly present our novel matching technique; a more detailed description can be found in Appendix B.
The matching is done for each computational cell in a selected region called the interpolation zone.We interpolate the magnetic and the electric fields, and then we re-initiate the particles to reproduce the first and second moments of the interpolated distribution function.The zeroth moment of the interpolated distribution must correspond to an integer number of simulated particles.This leads to noise in the charge density that we balance by additional electric field to satisfy Gauss' law.A nonzero divergence of the magnetic field is cleaned by the projection method (see Zhang & Feng 2016).After these two corrections the plasma self-consistently evolves without any numerical transients.
Table 1.The density fluctuation amplitudes and the measured shock properties in our simulations.
Note-All simulations have the same sonic and Alfvénic Mach number, MS ≈ 36 and MA ≈ 32, likewise the plasma beta, βp ≈ 1.The width of the simulation box for each case is roughly Ly ≈ 6λsi.The shock velocities v sh are measured in the upstream frame, and the electron temperatures Te are measured in the downstream region.

Simulation setup
To examine the influence of compressive pre-existing turbulence on perpendicular shocks we performed three simulations, one with homogeneous upstream medium (run H), and two with different amplitudes of density fluctuations, δn/n = 3.5% (run T1) and δn/n = 10% (run T2).For easy comparison with previous kinetic simulations of SNR shocks we used the same plasma parameters as for run B2 of Bohdan et al. (2019a), except for a smaller width of the simulation box, L y ≈ 6λ si instead of L y = 24λ si .This choice is computationally cheaper, but some features like shock front corrugations may not be resolved.Nevertheless, the size is large enough to probe the influence of the turbulence on the electron scale physics, which is the main aim of this work.Table 1 shows the turbulence amplitude and other important plasma parameters for each simulation.The amplitude of the density fluctuations, δn, is calculated as the RMS of the particle density measured in small tiles that are a quarter of ion skin length in size.
We initialize the upstream plasma with twenty particles per cell for both ions and electrons, n 0 = n i = n e = 20, where the subscript "i" represents ions and "e" electrons.The plasma flows toward the reflecting wall with velocity v 0 = v 0 x, where v 0 = −0.2c.The large scale magnetic field is perpendicular to the shock normal, θ 0 = 90 • , and has the in-plane configuration, B 0 = b 0 ŷ.This magnetic field orientation leads to particle gyration in the xz plane, giving particles three degrees of freedom and an adiabatic index Γ = 5/3.The expected shock speed in the upstream frame is v sh ≃ 0.264c.
In the shock simulation the particle species are in quasi-equilibrium T e ≈ T i ≈ 1.6 • 10 −3 m e c 2 /k B ≈ 10 7 K.For such a non-relativistic thermal distribution of particles the thermal speed of electrons is defined as v th,e = 2k B T e /m e ≃ 0.057c.On account of the ion-to-electron mass ratio m i /m e = 100, the thermal speed of ions is ten times smaller.The plasma beta, defined as the ratio of the thermal pressure to the magnetic pressure, is β p ≈ 1 for all runs.The sonic Mach number then follows as M S = v sh /c S ≈ 36 for all runs, where c S = Γk B (T e + T i )/m i ≈ 0.0074c.Likewise, the Alfvén speed v A = B 0 c/ n(m e + m i ) = 0.00829c leads to an Alfvénic Mach number M A ≈ 32.The expected compression ratio at the shock is r ≈ 3.97.
The time step, δt = 1/40ω −1 pe , scales with the electron plasma frequency, ω pe = q 2 e n e /(ϵ 0 m e ).We evolve the system for ten ion cyclotron times, Ω −1 i = m i /(q i b 0 ) ≈ 48, 000δt.The electron skin depth is λ se = c/ω pe = 20∆, where ∆ represents the cell size; for ions it is higher by a factor m i /m e , which gives λ si = 200∆.
The chosen values for our simulation parameters are a compromise between the capabilities of modern supercomputers and the accurate modelling of the shocks.By using a shock speed that is higher than typically observed in Galactic sources, but still non-relativistic, we are able to extend the duration of our simulations.In addition, these parameter choices allow for a direct comparison with previous studies of high-Mach-number shocks.

SIMULATION RESULTS
Here, we compare simulations with different intensities of upstream turbulence, see Table 1.Firstly, we discuss the global structure of a shock wave propagating in pre-existing turbulence, including the shock-reformation period, the shock speed, and the magnetic-field amplification.Then we focus on the shock foot, where we investigate how the density fluctuations influence the Buneman instability.Afterwards, we examine the Weibel instability together with magnetic reconnection.Further, we discuss the vorticity in the context of magnetic-field amplification via the turbulent dynamo and how the orientation of the magnetic field is modified by pre-existing fluctuations.Finally, we present the electron energy spectra along with the temperature of the particles.

Global shock structure
The structure of high-Mach-number supercritical shocks is determined by the reflection of upstream ions at the shock front, which for perpendicular shocks simply is a gyration in the shock-compressed magnetic field and leads to the formation of the shock foot, the shock ramp, and the overshoot-undershoot structure downstream of the shock.Figure 1 shows the structure of a shock propagating in turbulent upstream medium (run T2), at time t ≈ 10Ω −1 i .The upstream plasma at x − x sh ≳ 12λ si carries density fluctuations and significantly weaker fluctuations of the magnetic field.The shock foot spans from 10λ si to 12λ si ahead of the shock front.It contains electrostatic waves excited by the electrostatic Buneman instability, apparent in the E x map.The foot region is followed by the ramp, which extends to the overshoot at x−x sh ≈ −1λ si .This region is dominated by the Weibel instability, which forms filamentary structures in the density and in the magnetic field.The overshoot-undershoot pattern about 10λ si behind the shock front marks the transition to the downstream region.
The top panel of Figure 2 compares the ion-density profiles of all runs.Regardless of the intensity of the upstream density fluctuations, neither the characteristic shock structure nor the compression ratio are affected by the upstream turbulence.The compression is always consistent with the expected ratio of r = 3.97.
Figure 3 shows the evolution of the density profile averaged over the y-direction for a simulation with turbulent upstream plasma (run T2).For visualization purposes, the false-color representation is truncated 20λ si ahead of the shock.The pre-existing turbulent structures are advected to the shock front, thus they appear as dark blue stripes on the plot.The cyclic shock reformation seen as modulations of the shock front is a common feature of high-Mach-number shocks (see e.g., Wieland et al. 2016).They are caused by the nonstationary reflection of ions: the shock re-develops with an average period of 1.64Ω −1 i for all the runs, similar to the 1.50Ω −1 i reported in Wieland et al. (2016) and the 1.55Ω −1 i found in Bohdan et al. (2017).A corresponding quasi-periodic modulation is seen in the shock speed and in the plasma density.Reflected particles reach at most 12λ si ahead of the shock front, and the extent of the shock foot is the same for homogeneous and for turbulent upstream plasma.We do not observe any of the shock distortions that were seen in previous hybrid and kinetic simulations (e.g.Nakanotani et al. 2022;Demidem et al. 2023).There the density structures were much larger than the shock width, whereas here the largest density clumps are with roughly 3λ si considerably smaller than the width of the shock transition, about 15λ si .
In Table 1 we present the shock speed in the upstream frame, calculated over three fully developed reformation cycles.Within the uncertainties, the shock speed does not vary with the upstream density fluctuation level and is consistent with the expected value of 0.263c.
The bottom panel in Figure 2 presents the magneticfield strength normalized by the particle density to remove the effect of compression.The amplification of the magnetic field reaches a factor 3.3 at the overshoot regardless of the intensity of the pre-existing fluctuations, which is consistent with the values reported in Bohdan et al. (2021).The findings suggest that upstream density fluctuations with amplitudes below 10% do not modify the amplification of the magnetic field.

Buneman instability
Interaction of the shock-reflected ions with the incoming electrons causes the growth of electrostatic Buneman waves in the shock foot.They play an important role in electron heating as trapped particles undergo shock surfing acceleration.We investigate how these modes are modified by the pre-existing upstream density fluctuations.We performed a Fourier analysis in the region where the electrostatic field has the largest amplitude during the latest reformation cycle.
Figure 4 compares the power spectrum of the electric field parallel to the wave vector, E ∥ ≡ E ∥ k, for runs with homogeneous (run H, top) and turbulent upstream plasma (run T2, bottom).The evolution of the Buneman waves is non-stationary due to the shock reformation, and to compare the activity of the modes in different simulations, we chose a region where they have the largest magnitude during the last reformation cycle.Before calculation of the power spectrum, we remove the large-scale gradients from the electric field maps.The spectra have a well-defined peak in k x and are broader in the k y -direction, which is in agreement with the multidimensional analysis of the Buneman instability (see e.g., Amano & Hoshino 2009).There are no significant differences between the spectra and the electrostatic energy density between our runs, when we average over reformation cycles: 8.0 • 10 −4 n e m e c 2 for run H,  The relative flow of ions and electrons can be characterized by the velocity difference of the respective beams ∆v.Here, we consider the cold beam of incoming electrons moving along the x-axis.The shock-reflected ions gyrate in the xz plane, since in our simulations the per-pendicular magnetic field has the in-plane configuration.From the linear theory of the Buneman instability, the wave spectrum is dominated by modes parallel to ∆v: where ∆v is the magnitude of the velocity of reflected ions relative to the incoming electrons (Lampe et al. 1974).In 2D3V simulations, however, the z-component of the streaming motion does not drive waves, because the wavevectors must lie in the simulation plane.As a result, the relevant streaming motion of ions is primarily in the x-direction, yet the particles are warm (they gain energy at the shock): the y-component of ∆v is not negligible.Therefore, in our simulations, we expect the waves to be slightly oblique with respect to the xaxis.Moreover, Equation 1 applies to the fastest growing mode, meaning that signal over a range of wavenumbers perpendicular to the relative flow is present in the spectrum.
In this paragraph we investigate whether the wave power spectra from Figure 4 are consistent with the predictions of linear theory.To calculate the expected wavenumbers from Equation 1we estimate ∆v directly from the distribution functions of both plasma species.The relative speed for run H is approximately equal to ∆v ≈ 0.25c, thus we expect the strongest modes with kλ se ≈ 4.0.The position of the peak in the E ∥ power spectrum shown in Figure 4a is kλ se = (3.8,1.5), so the magnitude |k|λ se ≈ 4.1 matches the wavenumber expected from Equation 1.The relative speed computed for run T1 approximately equals ∆v ≈ 0.24c, which leads to the following wavenumber: kλ se ≈ 4.2.The strongest signal in Figure 4b is found at kλ se = (4.2,0.8), hence the magnitude is roughly equal to |k|λ se ≈ 4.3.This again corresponds to the predicted value.For run T2 the relative speed is roughly ∆v ≈ 0.24c, which gives kλ se ≈ 4.2. Figure 4c shows the peak at kλ se = (3.8,1.2), which leads to |k|λ se ≈ 4.0.In this case, the linear theory predictions are also roughly in agreement with the power spectrum obtained from the simulation.

Weibel instability
Interaction between the incoming and the shockreflected ions triggers the Weibel instability.It strongly amplifies the magnetic field and changes its structure at the shock foot and at the ramp.Moreover, the Weibel filaments may become tearing-mode unstable and decay through magnetic reconnection (Matsumoto et al. 2015;Bohdan et al. 2020b).This creates magnetic vortices, which may collate into larger structures.The energy density is normalized by the factor UB 0 = |B0| 2 /(2µ0).The vortex number corresponds to the number of local maxima in z-component of the magnetic vector potential, for visualization purposes the data was smoothed.
The top panel in Figure 5 compares the evolution of the energy density of the Weibel-generated magnetic field for all runs.To calculate this quantity, we consider the magnetic field in the region where the Weibel instability operates: −5λ si < x − x sh < 15λ si .Then, the large-scale gradients are removed to exclude the effect of compression, and the Fourier power spectrum is computed.The energy density is obtained directly from the spectrum, and it is normalized to the energy density of the initial magnetic field U B0 = |B 0 | 2 /(2µ 0 ).To be noted are the cycle-to-cycle variations in the quasiperiodic behavior caused by the shock reformation.On average there is no significant difference in the energy density with and without upstream turbulence, and likewise in the magnetic reconnection activity, since it depends on the strength of the Weibel instability.Indeed, the bottom panel in Figure 5 indicates a similar number of magnetic vortices for all runs, traced here as local maxima in z-component of the magnetic vector potential.

Vorticity
In the MHD regime, the postshock magnetic field can be amplified via the turbulent dynamo (Giacalone & Jokipii 2007;Fraschetti 2013;Xu & Lazarian 2016;Hu et al. 2022).In this process, the interaction of a shock with pre-existing density inhomogeneities causes corrugations of the shock front in the form of ripples, which generate strong rotation and vorticity of the downstream fluid.The frozen-in magnetic field lines are bent and distorted, and that leads to amplification of the field.For this reason, we investigate vorticity, specifically the curl of the particle current density, on the kinetic scales.
Figure 6 shows maps of the magnitude of the curl of the current density and the magnetic field, calculated for the run with turbulence (T2), at the stage of the reformation cycle when the Weibel filaments are well developed (t ≈ 8.5Ω −1 i ).The subscripts "e" and "i" refer to electrons and the ions, respectively.The pre-existing density fluctuations do not carry significant vorticity, hence in the upstream the curl of the currents and the magnetic field is weak.At the shock foot, the vorticity is dominated by the Weibel filaments.The dense filaments found between x − x sh = (0 − 5)λ si are composed of plasma moving towards the shock.They are surrounded by the stream of the shock-reflected particles, inducing shear flows that cause strong curl of both J i and J e at the filament contours.The current density inside the filaments produces the strong curl of the magnetic field, which is seen in the bottom panel of Figure 6 in that region.The Weibel filaments are unstable and dissipate via magnetic reconnection.This is visible as circular structures both in ion and electron current, as well as local maxima of the curl of the magnetic field.The magnetic vortices converge into larger ones and travel through the shock further downstream.A particularly large magnetic vortex lies between −15λ si < x − x sh < −11λ si .The corresponding current comes only from electrons since the dynamic scale of the shock-transmitted ions is much larger than the scale of the magnetic field structure.
The vorticity strongly varies between one reformation cycle and the next.We observed a large vortex in the downstream region only for run T2 during the third reformation cycle.We find no correlation between the presence of the large vortex and the magnetic reconnection activity, which indicates serendipity of this event.On average we do not observe a significant difference between the vorticity level in simulations with and without density fluctuations in the upstream plasma.The vorticity seems to arise only from magnetic reconnection at the Weibel filaments.This may be different if the scale of the density clumps is larger than the shock width or the scale of rippling instabilities.Study of the latter would require a correspondingly wide simulation box.

Obliquity angle
Variations of the magnetic field ahead of the shock front may lead to local changes of the shock obliquity angle.This means that reflection conditions of particles are modified.In this section, we examine whether the variance of the obliquity angle in the foot region depends on the presence of pre-existing density fluctuations.
Figure 7 shows 2D maps of the difference between the local obliquity angle and the inclination angle of the external magnetic field: ∆θ = θ − θ 0 , for run H (the middle panel) and T2 (the bottom panel), at the time t ≈ 8.8Ω −1 i .Additionally, the root mean square of this quantity calculated over the y-direction is included.The RMS value for run H ahead of the shock foot, x ≥ 15λ si , is very much smaller than that in the runs with turbulence it approximately equals: 1.1 • (T1) and 2.3 • (T2).Even so, the variation of the obliquity angle becomes similar for all runs in the regions closer to the shock front: at the foot, at the ramp and at the overshoot.The structure of the magnetic field is strongly modified there by the Weibel instability.The similar values of the RMS for all runs indicates that the instability is not affected by the upstream turbulence (see Section 3.3).

Transmission of fluctuations
Having discussed the influence of the upstream density fluctuations on the shock waves, we now explore how the turbulence is transmitted across the shock, focusing on the amplitude of the velocity and density fluctuations.
The upper panel in Figure 8 shows the amplitude of fluctuations of the ion bulk speed, at scales of roughly 3λ se ×3λ se .The RMS of the fluctuations, δv i , was calculated in regions of the box width, ∆y ≈ 60λ se , and combining two neighbouring cells in the x-direction.The pre-existing turbulence in run T2 carries velocity fluctuations of roughly 0.2% of the shock speed.For run T1 the amplitude is approximately smaller by factor 3. Nevertheless, in all simulations, the amplitude of the fluctuations at the shock ramp and further downstream of the shock reaches roughly 10% and 1% of the shock speed, respectively.
The bottom panel in Figure 8 presents the amplitude of density fluctuations for all simulations.The amplitude δn i is calculated in the same way as for δv i .In the upstream region, x − x sh > 12λ si , the initial fluctuation level carried by the pre-existing turbulence are equal to δn i /n 0 ≃ 3.5% and δn i /n 0 ≃ 10% for runs T1 and T2, respectively.This corresponds to the turbulence amplitudes from Table 1.Furthermore, there is no significant spatial variation of the turbulence level along the x-axis ahead of the shock foot.This sustains that the statistical properties of the fluctuations are stable.For runs T1 and T2 the amplitude of the density fluctuations increases gradually from x − x sh ≃ 10λ si and x − x sh ≃ 9λ si , respectively, whereas for the homoge- The lower panel shows the RMS of the difference between the local obliquity angle and the inclination angle of the external magnetic field for each run.The profiles are obtained from ∆θ maps by calculating the RMS over the y-direction for each time step in the reformation cycle.Then the quantity is time-averaged and normalized by the factor θ0.
neous setup, it happens at the beginning of the shock foot, x − x sh ≃ 12λ si .This amplitude growth is caused by the formation of the Weibel filaments, which at some distance surpass the pre-existing turbulence.For all cases, the fluctuation amplitude reaches δn i /n 0 ≈ 250% at the overshoot and then decreases to δn i /n 0 ≈ 10% further in the downstream.
Recent 2D hybrid simulations of a shock wave propagating in a high-β turbulent medium show significant amplification of the density and the velocity fluctuations (Nakanotani et al. 2022) for density structures of a few tens of ion skin lengths and amplitudes of unity of the upstream density, whereas we see little effect of upstream density variations that are about ten times smaller.We might expect analogous physical features at both scales, but with a smaller size any forcing is applied for a correspondingly shorter time, hence the smaller impact.In our simulations, the magnitude of fluctuations in the shock transition and in the downstream regions is independent of the strength of the pre-existing turbulence.This suggests that the level of the fluctuations is dominated by processes occurring in the shock transition (e.g., the Weibel instability, compression).

Downstream electron temperature
Figure 9 compares the electron energy spectra in the far downstream region, x − x sh = −(32 -64)λ si , for the runs with homogeneous plasma (H) and with δn/n = 10% (T2).The energy distribution is calculated for the late stage of the simulation and in the local plasma rest frame.The electron spectra are represented by a kappa distribution with κ = 13, denoted by the black solid line, demonstrating a weak high-energy tail that is virtually the same in the simulations with and without upstream turbulence.In general, the shape of the downstream electron spectra is a consequence of particle interactions such as shock-surfing acceleration and magnetic reconnection.The similarity of the spectra reflects the absence of significant modifications of the main instabilities at the shock.The temperature of electrons is calculated from the fitted thermal distribution and is presented in Table 1.Within the uncertainties, and averaging over the spatial variation of the downstream temperature, we find the same temperature in all runs.The obtained electron temperatures are in agreement with the value 0.183 ± 0.004 reported for run B2 in Bohdan et al. (2019b).The electron heating model of Bohdan et al. (2020a) includes a combination of different processes operating in the shock transition, modification of which by small-scale density fluctuations was not observed, suggesting that there is no impact on electron heating.
We do not discuss the temperature of ions because the duration of our simulations is too short for these particles to reach equilibrium.However, findings from 1D PIC simulations suggest that spectra of thermal and suprathermal ion populations in the downstream region of quasi-parallel shocks can be represented by a kappa distribution (Arbutina & Zeković 2021).

SUMMARY AND DISCUSSION
Understanding the physics of a plasma shock wave propagating in an inhomogeneous medium is an important challenge in astrophysics and space physics.We present results of kinetic simulations of high-Machnumber shocks with pre-existing density fluctuations, which model conditions at SNRs.Our main goal was to investigate the influence of the turbulence on electronscale phenomena: the Buneman and the Weibel instabilities, electron acceleration, and their heating.The second aim was to examine the global shock structure, including the properties of the magnetic field.Furthermore, we analysed how the upstream fluctuations are transmitted across the shock.Our results can be summarised as follows: 1. Our novel simulation technique establishes a framework for studying shocks propagating in inhomogeneous media.It may be used to model various physical environments and does not bear significant computational costs.
2. The achievable level of pre-existing upstream turbulence on kinetic scales is limited by particle heating.Our empirical test indicate that to maintain M S ≳ 30 the maximum amplitude of density fluctuations should be of the order of δn/n ∼ 10% (on the scale of a quarter of the ion skin length).The limit should be even lower for shock speeds closer to those observed in SNRs.
3. We performed two simulations of a perpendicular high-Mach-number shock with different RMS amplitude of pre-existing density fluctuations, δn/n = 3.5% and δn/n = 10%.The amplitudes were chosen to avoid excessive heating and in concordance with in-situ measurements.We examined the impact of the turbulence on the shock physics and the electron dynamics by direct comparison to a simulation with homogeneous upstream medium.
4. Density fluctuations up to a few λ si do not change the global structure of a perpendicular shock wave, as well as its reformation period and propagation speed.No significant difference was found in the amplification of the magnetic field, on account of weak turbulence-driven vorticity on kinetic scales.
5. Fourier analysis of the Buneman waves in the shock foot shows that their properties remain unchanged under interaction with the density inhomogeneities.Likewise for the Weibel modes and magnetic reconnection at them.The magnetic field fluctuations carried by the pre-existing turbulence locally modify the magnetic field obliquity at the shock foot.At the ramp the topology of the magnetic field is dominated by the Weibel instability.
6.The velocity and the density fluctuations in the downstream region reach a comparable amplitude for all levels of pre-existing density fluctuations that we studied.
7. Energy spectra of downstream electrons are thermal with a weak tail, but the electron preacceleration is not sufficient for injection into DSA.
The temperature and non-thermal fraction are similar for all simulations, and electron acceleration and heating are unmodified.
Taken together, these findings suggest that preexisting density fluctuations of a few λ si and with realistic intensity do not considerably affect the physics of high-Mach-number shocks.In particular, we did not observe any changes in the properties of plasma instabilities or in the efficiency of electron acceleration.This may be different for turbulence on larger scales, as hybrid simulations suggest (Trotta et al. 2021;Nakanotani et al. 2022;Trotta et al. 2023).There turbulent structures introduced by shock corrugation enhance particle acceleration efficiency, as well as amplify the magnetic field via the turbulent dynamo.Analogous studies in kinetic regime require using a wider simulation box and thus more computational resources.
The density fluctuations are also too weak to impact the ion reflection dynamics, which also could modify the instability conditions.On the other hand, the density structures are comparable in size to the Weibel filaments.Although, they have a much smaller amplitude, so they do not modify the Weibel instability.At lower Alfvénic Mach numbers the Weibel modes should be weaker, and correspondingly the impact of pre-existing turbulence larger.
Further studies should explore the effect of the Alfvénic turbulence, where the magnetic field fluctuations dominate.
In contrast to perpendicular shocks, from which electrons cannot escape to the upstream region, at oblique shocks the particles form extended foreshock, where they drive instabilities (Bohdan et al. 2022;Morris et al. 2022).This adds room for particle-turbulence interactions, and hence changes of the reflection conditions and the electron dynamics, as well as modifications of the instabilities.Moreover, 1D3V (Xu et al. 2020;Kumar & Reville 2021) and 2D3V (Morris et al. 2023) simulations show that a significant non-thermal electron population can be formed downstream of such shocks that is well described by a power-law.Therefore, at oblique shocks it is possible to directly investigate the influence of preexisting turbulence on energetic particles.
Full 3D simulations are currently overwhelmingly challenging, and so they must be run with lower resolution and smaller boxes and can only reach the very early stages in the evolution of the system (Matsumoto et al. 2017), in contrast to 2D3V simulations (Bohdan et al. 2017).Nevertheless, including all dimensions may influence the interplay and nonlinear development of the instabilities driven at the shock, and cross-field diffusion can be adequately captured, as recently argued by Orusa & Caprioli (2023).

A. GENERATION OF COMPRESSIVE TURBULENCE
In this work, we consider decaying and compressive turbulence, which is established by initially imposing a local disturbance of the bulk velocity.The velocity of the plasma species at the arbitrary location v(r), contains the bulk component, its disturbance, and the thermal fluctuations: v(r) = v 0 + δv 0 (r) + v th .The disturbance is in the form of the superposition of wave-like modes, with both the transversal and the longitudinal configuration.This is similar to techniques from Giacalone & Jokipii (1994) and Grošelj (2019), but we only impose the velocity disturbances without any magnetic field associated with them.The form of the transversal disturbance is given by the following equation: and for the longitudinal one: The relation between the primed system and the unprimed is given by the rotation matrix For each randomly chosen wavenumber k we chose a random orientation with respect to the x axis: 0 < ϕ k < 2π, and a random phase 0 < β k < 2π.Letter A k denotes the amplitude of each wave.Periodic boundaries require that the wavevector components have to fulfill the following identities: k x L = nπ and k y L = mπ, where L is the box size and n, m are random integers.We note that the initial spectrum strongly depends on combinations of the n, m pairs, but it changes as the system evolve.As we mention above, we do not impose any magnetic field disturbance, but since initially the magnetic field has a constant value, we must input the correct motional electric field E 0 = −v × B 0 .In our simulations, we generate 50 waves (longitudinal or transversal), with randomly chosen n, m such that the wavevector component is in the range 2π L < k x < 6 • 2π L (same for k y ).These pairs of integer numbers give us the orientation of each wave ϕ k .Phase of the wave β k is chosen randomly, furthermore, we assume that all amplitudes are equal: A k = A = const.This setup efficiently generates compressive turbulence, which intensity can be adjusted by the parameter A. The well-developed density structure of turbulence with the intensity δn/n ≈ 10% is presented in Figure 10 (the left panel).The decay time of the fluctuations exceeds 2Ω −1 i , so there are sufficiently long-lived to be inserted into a shock simulation.The right panel in this figure shows the time evolution of the particle kinetic energy E k , and the sum of the particle thermal energy and the energy of the electromagnetic field E th + E EM .The energies are averaged over the simulation box, their initial values are subtracted, and they are normalized to the ion rest energy.Our kinetic simulations cover the dissipation range of turbulence, so the energy of the initial fluctuations is partially transferred into plasma heat.The sonic Mach number of a shock is defined by the end state of the particle temperature, hence the heating constraints the maximum achievable turbulence level for a given Mach number.

B. MATCHING TURBULENT PLASMAS
Injection of a pre-fabricated turbulent plasma slab into a shock simulation requires matching of the plasma at the interface.Otherwise, any mismatch in the electromagnetic field and the current density may drive artificial transients.Here, we discuss the matching of two arbitrary plasma slabs, called the left and the right, pre-fabricated in simulations with periodic boundaries.We chose a region where the electromagnetic fields and the current densities are interpolated: the interpolation zone.The values of the fields in each grid cell of the zone are replaced by the interpolated ones, considering the magnetic field: where the subscripts "L" and "R" denote the left and the right slab, i, j are the cell numbers in x and y-direction respectively, and B is the interpolated value of the magnetic field.The weights w(i) depend only on the i-th location on the simulation grid, and they are imposed in the form of a linear function.The same procedure applies to the electric field.Simple calculations show that the field interpolation does not contribute to large factors in both Ampère's and Faraday's law.

Figure 1 .Figure 2 .
Figure 1.Maps of the electron number density (top panel) and the x-components of the magnetic (middle panel) and the electric fields (bottom panel) at a shock with pre-existing upstream density fluctuations (run T2).The scaling of the fields is logarithmic and sign-preserving, e.g. for Bx/B0 it is sgn(Bx) • [2 + log{max(10 −2 , |Bx|/B0)}].

Figure 3 .
Figure 3.The evolution of the ion number density averaged over the y-direction for run T2.Dark blue stripes represent upstream density fluctuations, which are advected to the shock front.For visualization purposes, the presentation is truncated 20λsi ahead of the shock.

Figure 4 .
Figure 4.The power spectrum of the electric field parallel to the wave vector, E ∥ k, for all runs, in the region where Buneman waves are the strongest during the latest reformation cycle.The power spectrum is shown in units of (ωpemec/e) 2 .

Figure 5 .
Figure5.The evolution of: the energy density of the Weibel-generated magnetic field (top panel) and number of the magnetic vortices (bottom panel), for all simulations.The energy density is normalized by the factor UB 0 = |B0| 2 /(2µ0).The vortex number corresponds to the number of local maxima in z-component of the magnetic vector potential, for visualization purposes the data was smoothed.

Figure 6 .
Figure 6.Maps of the magnitude of the curl of the electron current density (top panel), the ion current density (middle panel), and the magnetic field (bottom panel), for run T2 in the middle of the third reformation cycle (∼ 8.5Ω −1 i ).The current densities are normalised by the factor J0 = n0qec.

Figure 7 .
Figure 7. Variations of the obliquity angle in the region extending from the upstream to the overshoot.The upper and middle panels show maps of the difference between the local obliquity angle and the inclination angle of the external magnetic field for the runs H and T2, respectively, and the time t ≈ 8.8Ω −1 i .The maps are presented in logarithmic and signpreserving scale: sgn(∆θ) • [2.5 + log{max(10 −2.5 , |∆θ|/θ0)}].The lower panel shows the RMS of the difference between the local obliquity angle and the inclination angle of the external magnetic field for each run.The profiles are obtained from ∆θ maps by calculating the RMS over the y-direction for each time step in the reformation cycle.Then the quantity is time-averaged and normalized by the factor θ0.

Figure 8 .
Figure 8. Level of the speed and density fluctuations for all simulations in the region extending from upstream to downstream of the shock.The upper panel shows the RMS of the ion bulk speed fluctuations normalised to the shock speed, and the lower panel shows the RMS of the ion density fluctuations normalised to the initial upstream density.The profiles are obtained from the velocity and the density maps by calculating the RMS over the y-direction and combining two neighboring cells in the x-direction.The profiles are timeaveraged over the latest reformation cycle.

Figure 9 .
Figure9.The electron energy spectrum in the downstream region for runs H (the blue dashed line) and T2 (the red dotted line).The black solid line denotes a fitted kappa distribution.
Figure 10.Left: Normalized ion number density map showing well-developed structure of the compressible turbulence, after 1Ω −1 i .Right: Time evolution of the mean kinetic energy, as well as sum of the thermal energy and the electromagnetic field energy.The initial values are subtracted, the energies are normalized to the ion rest energy.

Figure 11 .
Figure 11.Left: The matching procedure presented for one grid cell.Right: Maps of the magnetic field after around half an ion Larmor time after merging.The interpolation zone spans from x = 5.3λsi to x = 5.8λsi (middle of the panel).