The need for high-resolution studies

The main problem with performing numerical simulations of shocks in plasma is the constraint by the available computational resources. Present-day simulations employ billions of computational particles but at the same time the size of the simulation grid and its resolution and a final simulation time must be properly chosen for the characteristic physical distances and frequencies of the problem under study. The inclusion of ions in PIC simulations constitutes a further challenge: in fact , if the simulated system is only composed by electrons and positrons with mass me, the characteristic time scale that must be resolved is typicallyω⁻¹_pe ∼ √

me, and it is possible to run a simulation for thousands ofω⁻¹_pe and analyse the system up to a strongly non-linear phase. However, studying protons in the simulations requires much longer time scales. In order to correctly capture the ion behaviour, a PIC simulation must resolve a good number of ion characteristic time scales, ω⁻¹_pi ∼ √

m_i. If one were to perform a simulation for a realistic mass ratio mi/me≈ 1836 this would correspond to running the simulation for a time interval that is

√

m_i/me ≈ 42 times longer than for the case involving only leptons. The spatial scale also constitutes a problem, since ions require correspondingly larger spatial scales. It can be seen that in principle the simulation box size that is necessary to model electron-proton plasmas increases compared to that required by leptonic plasmas by a factor of ∼

√

D, where D is the number of spatial dimensions resolved by the code [Bret and Dieckmann, 2010]. For this reason PIC

simulations that use realistic electron-to-ion mass ratio are very unusual and mostly restricted to 1D regimes [e.g., Dieckmann et al., 2006] and to 2D simulations that need to resolve a very small spatio-temporal domain [e.g., Honda et al., 2000], while 2D PIC simulations that cover a large domain and 3D PIC simulations normally introduce reduced mass ratio of 16 ÷ 100 (e.g. [e.g., Spitkovsky, 2007, Niemiec et al., 2008, Hoshino et al., 1992, Martins et al., 2009, Stockem et al., 2012, Iwamoto et al., 2019]). However, caution must always be exercised since some phenomena such as magnetic reconnection and Buneman instability have been shown to be influenced by the mass ratio [Guo and Lu, 2007, Bohdan et al., 2017].

A second, but very important factor is the simulation the space resolution, i.e., the number of simulation cells used to resolve electron and ion skin depths. High-frequency short-wavelength oscillations characteristic of, e.g., precursor waves emitted by the shock towards the upstream plasma may be easily damped out in case of insufficient resolution. Consequently, simulations which aim at investigating these kind of waves must set this parameter carefully. If, on one hand, applying a lower resolution eases the computational burden, on the other hand it may fully or partially suppress physical phenomena. Very low resolutions ofλse = 10∆x have been used in, e.g., Sironi and Spitkovsky [2011] and λse ≈ 2.3∆x in [Stockem et al., 2012]. Iwamoto et al. [2017] investigated the relation between the value of these simulation parameters and the numerical damping of the SMI shock-emitted waves by performing numerical convergence tests in Nppc and the resolution. They observed that the precursor wave amplitude systematically increases with the resolution until the value at which a convergence is achieved, that occurs for λse = 40∆x in their model. This proves that the precursor wave emission is very sensitive to numerical resolution. As the code used for my simulations is structured differently from the one in [Iwamoto et al., 2017], I performed my own set of test simulations, to find the appropriate simulation parameters.

I have performed two sets of test simulations for a numerical convergence study.

In the first set I used 2D simulations with the setup the same as in our run with the out-of-plane magnetic field, ϕB = 90◦, and I probed different resolutions in terms of the electron skin depth values of λse = 15∆, 40∆ and 80∆, and different number of particles per cell, Nppc= 10,20 and 40. The transverse size of the numerical grid was

Figure 6.2: Profiles of B_zmagnetic field fluctuations upstream of the shock located at x/λsi' 33 at time tΩci' 22.5 and different Nppc= 10 (a), 20 (b) and 40 (c) for the 2D convergence tests.

Profiles are taken in the middle of the box along y/λsi= 1.5 and are normalised to the initial magnetic field strength. Hereλse= 15∆.

Figure 6.3: Profiles of B_zmagnetic field fluctuations upstream of the shock located at x/λsi' 33 at time tΩci ' 22.5 and different λse = 15∆ (a), 40∆ (b) and 80∆ (c) for the 2D convergence tests. Profiles are taken in the middle of the box along y/λsi= 1.5 and are normalised to the initial magnetic field strength. Here Nppc= 10.

adjusted to have Ly' 3λsi in each case. Figure 6.2 shows the profiles of Bz magnetic field fluctuations taken in the middle of the box along y/λsi = 1.5 at time tΩci ' 22.5 for λse = 15∆ and different Nppc. Similarly, Figure 6.3 presents results for the case of Nppc = 10 and different λse. In Figure 6.4 we compare the normalised amplitudes of the waves at different λse and Nppc. The amplitudes are averaged in the region of x/λsi= 35 − 40, located about 2λsi from the shock. One can see that the number of particles per cell does not influence the amplitude of the precursor waves, and so that I choose Nppc= 10 for our 2D large-scale production runs. However, the amplitude increases with the electron skin depth.

Figure 6.4: Normalized amplitudes of δBz at different λse and Nppc obtained in 1D (blue dashed line) and 2D (solid blue, green, and red lines) convergence test runs. The amplitudes are calculated at tΩci' 22.5 and are averaged in the region of x/λsi= 35 − 40, located about 2λsi

from the shock (compare Fig. 6.19). Note that results for 2D tests and different Nppcoverlap.

To verify whetherδBz/B0is saturated at the maximum resolution ofλse= 80∆ tested with 2D simulations, we performed the second set of test runs with 1D simulations with the same λse and also extending the probed resolution range for the cases with λse = 100∆, 120∆, and 140∆. For 1D tests we take Nppc= 10 and use the 2D simulation code but with 5 cells in the y-direction, which is the minimum number of grid points in the transverse direction allowed in our code. 1D tests show that the convergence is indeed obtained for λse = 80∆, which is the value used in our production runs. Note that the amplitudes in 1D simulations are slightly larger than in 2D runs at the same skin depth because the shock has not yet reached a steady state (see discussion in Section 6.2.2(d) and compare Fig. 6.19).