Particle transport in hybrid PIC shock simulations: A comparison of diagnostics
D. Trotta D. Burgess G. Prete S. Perri G. Zimbardo


firstpage–lastpage
 

Quasinormal modes of Dirac field in the Einstein-dilaton-Gauss-Bonnet and Einstein-Weyl gravities
    A. F. Zinhailo
    December 30, 2023
=================================================================================================





Recent in-situ and remote observations suggest that the transport regime associated with shock accelerated particles may be anomalous i.e., the Mean Square Displacement (MSD) of such particles scales non-linearly with time. We use self-consistent, hybrid PIC plasma simulations to simulate a quasi-parallel shock with parameters compatible with heliospheric shocks, and gain insights about the particle transport in such a system. For suprathermal particles interacting with the shock we compute the MSD separately in the upstream and downstream regions. Tracking suprathermal particles for sufficiently long times up and/or downstream of the shock poses problems in particle plasma simulations, such as statistically poor particle ensembles and trajectory fragments of variable length in time. Therefore, we introduce the use of time-averaged mean square displacement (TAMSD), which is based on single particle trajectories, as an additional technique to address the transport regime for the upstream and downstream regions. MSD and TAMSD are in agreement for the upstream energetic particle population, and both give a strong indication of superdiffusive transport, consistent with interplanetary shock observations. MSD and TAMSD are also in reasonable agreement downstream, where indications of anomalous transport are also found. TAMSD shows evidence of heterogeneity in the diffusion properties of the downstream particle population, ranging from subdiffusive behaviour of particles trapped in the strong magnetic field fluctuations generated at the shock, to superdiffusive behaviour of particles transmitted and moving away from the shock.



Diffusion – shock waves – acceleration of particles 







§ INTRODUCTION


High energy particle acceleration is a common feature of many astrophysical systems, and characterizing the energetic particle transport is a key challenge for understanding the underlying physical acceleration processes. In most cases, when particle trajectories are available, the method of mean square displacement (MSD) is the most widely used diagnostic to study particle transport. MSD uses ensemble averages, and hence can be unreliable with poor particle statistics. On the other hand, the method of time-averaged mean square displacement (TAMSD) relies on time  rather than ensemble averages, and is a popular diagnostic in research communities concerned with particle transport in biological systems. When dealing with particle tracking in self-consistent plasma simulations of acceleration it can be challenging to obtain sufficiently good particle statistics over long times, due to the non-periodicity and limited extent of such simulations and consequently the varying lifetimes of different particles in the simulation. For this reason, we investigate the use of TAMSD diagnostics to study the transport of accelerated particles in self-consistent hybrid PIC shock simulations, and make comparisons with the results from the MSD method.

It is widely accepted that shocks are efficient particle accelerators in astrophysical systems <cit.>. Generally speaking, shocks convert directed flow energy (upstream) to thermal energy (downstream). In collisionless shocks, a small fraction of the upstream energy can be channeled towards particle acceleration. The mechanisms by which particles are accelerated out of the thermal pool and to high energies are still debated. One of the most important parameters that controls shock structure and behaviour in a collisionless plasma is the angle between the normal to the shock surface and the upstream magnetic field, θ_Bn. Shocks with θ_Bn close to zero are called quasi-parallel, and they are believed to be more efficient at particle acceleration than quasi-perpendicular ones <cit.>. Other important parameters are the shock Alfvénic and sonic Mach numbers, defined as M_A = v_sh/v_A and M_S = v_sh/c_s, respectively,  and the plasma β = v_th^2/v_A^2. Here, v_A is the Alfvén speed in the region upstream from the shock, v_sh is the shock speed in the upstream flow frame, c_s is the sound speed of the upstream flow and v_th is the thermal speed of the upstream flow.

Particle acceleration at shocks is invoked to explain a wide variety of direct and indirect observations. Direct evidence of shock accelerated particles has been found by means of spacecraft measurements since the 1960s <cit.> for interplanetary shocks. At the solar wind termination shock, the suprathermal particle population was found enhanced by Voyager spacecraft observations <cit.>. Indirect evidence of accelerated particles is commonly found in the X-ray signature of supernova remnants <cit.>. Recently, particle acceleration at shocks has also been found to be relevant in the galaxy intracluster medium <cit.>.

A widely used theory for particle acceleration at quasi-parallel shocks is Diffusive Shock Acceleration theory (DSA) <cit.>. The idea underlying DSA theory is that particles with Larmor radii larger than the shock thickness interact with the shock by diffusing upstream and downstream, and thus are continuously returned to the shock, gaining energy at each crossing. This process is often called first order Fermi acceleration. The main prediction of DSA is a power-law distribution in energy for the accelerated particles, with a slope that depends only on the shock compression ratio, at least in the simple time stationary limit. DSA theory also predicts that, in the case of a spatially constant diffusion coefficient, the density profiles of energetic particles should decrease exponentially upstream of the shock, and to be constant downstream. These predictions have been found to be not inconsistent with some observations of interplanetary shocks <cit.>, even though many in-situ spacecraft crossings exhibit more complex energetic particle profiles <cit.>.

Recently, some energetic particle observations have suggested that the transport regime associated with such particles can be anomalous. In fact, energetic particle density profiles at interplanetary shocks are often found to show a power-law decay rather than an exponential one, as one would expect from normal diffusion <cit.>. In addition,  it has been proposed, among other possible mechanisms,  that the discrepancy between Mach numbers obtained from radio and from X-ray observations at galaxy cluster merger shocks may be due to an extension of DSA to superdiffusive transport <cit.>, which leads to superdiffusive shock acceleration <cit.>. Therefore, it is now timely to understand under what conditions the transport of energetic particles at shocks could be considered to be superdiffusive.

One of the main assumptions underlying DSA is that energetic particles undergo normal (i.e., Brownian) diffusion when scattering back and forth across the shock front. Brownian diffusion is characterized by Gaussian statistics. For an ensemble of particles undergoing normal diffusion, the mean square displacement for particle trajectories grows linearly in time, as ⟨Δ s^2⟩ = 2 D t, where D identifies the so-called diffusion coefficient. 

If the time evolution of the MSD is nonlinear, for example, such that ⟨Δ s^2⟩ = 2 D_α t^α, then the system exhibits anomalous transport. The exponent α is known as the anomalous diffusion exponent and its value identifies the transport regime of the system. Anomalous transport has been observed in an enormous variety of systems, ranging from the foraging movements of spider monkeys <cit.> to stock options transactions in finance <cit.>. It is possible to predict, upstream of shock waves, power-law density profiles for shock accelerated particles assuming for them an anomalous diffusion regime <cit.>. In fact, an important feature of anomalous transport is that the systems are characterized by non-Gaussian statistics, which introduces probability distributions with power-law tails <cit.>. In DSA, the time scale for acceleration of particles of given energy is dominated by the slowest portion of the shock acceleration cycle <cit.>. An interesting implication of anomalous, superdiffusive transport is that it may possibly lead to shorter acceleration times than those usually obtained for normal diffusive transport, thus helping to reach higher maximum energies than those predicted by DSA <cit.>.

Simulations of particle acceleration at shocks have been used extensively to bridge theory and observations. Numerical and analytical solutions for the transport equation obtained from DSA have been studied under many approximations <cit.>. A popular approach to address energisation efficiency and particle transport is the use of Monte Carlo simulations <cit.>. A drawback of the usual Monte Carlo approach is the neglect of any self-consistent plasma effects (e.g., internal shock structure) which are believed to be important, particularly at low energies <cit.>. Efficient particle acceleration has been studied both in fully kinetic Particle-In-Cell (PIC) and hybrid (particle ions and fluid electrons) simulations <cit.>. In the fully kinetic model, the spatial scales are resolved down to the Debye length, making it possible to address the complete picture of particle acceleration down to electron dynamics. However, fully kinetic PIC simulations are computationally intensive. On the other hand, in the hybrid model, the electrons are modelled as a fluid and just the ion kinetic scales are resolved, with the advantage that considerably larger domains and time scales can be simulated. 

Particle diffusion has also been studied in the framework of self-consistent plasma simulations. <cit.>, for example, studied the transport of protons in different turbulent regions by means of hybrid simulations. The majority of plasma turbulence simulations are performed on a periodic domain and in the plasma rest frame, making it relatively easy to trace particle trajectories throughout the simulations and compute the MSD for them. However, for studying energetic particles in shock simulations the situation is more complex since only a small fraction of the particles are accelerated to high energies, they may exit the simulation domain because it is nonperiodic and there is no single rest frame due to the velocity change at the shock.

Although these particular features make it harder to track energetic particles in the simulations for sufficiently long times, early studies attempted to study particle transport at shocks using self-consistent simulations. For instance, <cit.> studied particle diffusion in turbulence generated by an ion-ion beam instability and invoked it as a model for the Earth's foreshock. <cit.> used hybrid PIC simulations of turbulence to model the immediate downstream of a perpendicular shock (by using a particular initial condition, consisting of a homogeneous plasma core distribution and a non-gyrotropic ion population) and addressed the particle transport in such a system. <cit.> considered the shock transition layer of a hybrid PIC shock simulation at a single time, imposed boundary conditions on it and evolved test particles in the resulting electromagnetic fields to study diffusion across the shock layer. 

Our goal is to study the particle diffusion regimes at quasi-parallel shocks using self-consistent, hybrid PIC simulations, with the aim of gaining extra information that might advance our understanding of particle transport and so improve Monte Carlo modelling of shock acceleration.We focus on a set of shock parameters relevant for heliospheric shocks and intermediate particle energies. As we mentioned above, it is not trivial to follow many particles for long times in a shock simulation, and this has motivated us to employ a novel diagnostic, namely the time-averaged mean square displacement method (TAMSD). TAMSD relies on time-averaging rather than ensemble-averaging as in MSD, making it possible to extract single particle diffusion exponents from their trajectories <cit.>. TAMSD has been used in a variety of systems where it can be difficult to obtain many particle trajectories to construct an ensemble average; this is the case for many biological systems, especially those for which in-vivo tracking of particles is required <cit.>.

In this paper the trajectories for an ensemble of protons self-consistently interacting with a 1D, quasi-parallel shock are studied. Upstream and downstream diffusion properties are analysed separately, and all the computations are performed in the respective local rest frame. We obtain upstream and downstream diffusion exponents using both MSD and TAMSD. A thorough comparison between the two diagnostics is presented, in order to test the applicability of TAMSD. We have found the two techniques to be consistent, although both MSD and TAMSD have their own limitations. Upstream, strong evidence for superdiffusion is found. In the shock downstream, we observe a richer scenario, encompassing superdiffusion and subdiffusion, and, based on the TAMSD diagnostics we suggest that this is related to different subsets of particles co-existing in different diffusion regimes.
The paper is organised as follows: in Section <ref> we describe both the simulation method and the diagnostics employed. In Section <ref> results for upstream and downstream transport are described, and summarised in Section <ref>. In Appendix <ref> we present some independent tests of the TAMSD diagnostic using Monte Carlo simulations.




§ METHODS




 §.§ Shock simulation


We simulate a shock using a hybrid PIC simulation. In the simulation, protons are modelled as macroparticles and advanced using the standard PIC method. The electrons, on the other hand, are modelled as a massless, charge-neutralizing fluid with an adiabatic equation of state. The hybrid code used throughout this work (HYPSI) is based on the CAM-CL algorithm <cit.>. 

The shock is initiated by the injection method, in which the plasma flows in the x-direction with a defined (super-Alfvénic) velocity V_in. The right-hand boundary of the simulation domain acts as a reflecting wall, and at the left-hand boundary plasma is continuously injected. The simulation is periodic in the y direction. A shock is created as a consequence of reflection at the wall, and it propagates in the negative x-direction. In the simulation frame, the upstream flow is along the shock normal.

In the hybrid simulations distances are normalised to the ion inertial length d_i ≡ c/ω_pi, times to the inverse cyclotron frequency Ω_ci^-1, velocity to the Alfvén speed v_A (all referred to the upstream state), and the magnetic field and density to their upstream values, B_0 and n_0, respectively. The angle between the shock normal and the upstream magnetic field, θ_Bn is 15^∘, with the upstream magnetic field  in the x-y plane. For the upstream flow velocity, a value of  V_in = 6 v_A  has been chosen, and the resulting Alfvénic Mach number of the shock is approximately M_A = 8. The upstream ion distribution function is an isotropic Maxwellian and the ion β_i is 0.5.  The simulation x-y domain  is 3000 × 10 d_i, and hence we describe the shock as quasi-1D. The spatial resolution used is Δ x = Δ y = 0.5 d_i. The final time for the simulation is 750 Ω_ci^-1, the time step for particle (ion) advance is Δ t_  = 0.005 Ω_ci^-1. Substepping is used for the magnetic field advance, with an effective time step of Δ t_B = Δ t_/10. A small, nonzero  resistivity is introduced in the magnetic induction equation. The value of the resistivity is chosen so that there are not excessive fluctuations at the grid scale. The number of particles per cell used is always greater than 300 (upstream), in order to keep the statistical noise characteristic of PIC simulations to a reasonable level.

The number of particles in the simulation is about 10^7, and we store the trajectories of about 5 × 10^5 particles. The time cadence for the particle trajectory outputs is chosen such that the time resolution obtained in the particle dataset is Δ t = 0.1 Ω_ci^-1. The protons are tracked in time starting at 250 Ω_ci^-1,  when the shock is already well-developed, with an upstream region with developed fluctuations and a well-defined downstream region separated from the right hand wall.



An overview of the shock simulation is presented in Figure <ref>, which shows the time evolution of |B| as a function of x position. Three typical trajectories for protons interacting with the shock are also shown (red solid lines). Throughout the simulation, there are upstream waves self-consistently induced by the reflected ion population, that have been extensively studied in previous literature <cit.>. The downstream region exhibits strong fluctuations at several wavelengths, and some of them are relatively coherent in time, an effect related to the quasi-1D nature of the simulation.

As can be seen in Figure <ref>, the upstream waves are artificially truncated. This effect, particularly evident at later times, is due to the finite size of the simulation domain, and is known to affect shock simulations. It can be commented that this is not just a spatial truncation, but also an amplitude truncation, because the extent of the current and the growth (advection) time is truncated. Also, if ions with a given energy escape before properly producing waves, it is likely that particles with the same energy will see a non-fully-developed set of fluctuations, and this could affect their transport properties. However, the results presented here have been confirmed using larger simulations (in spatial extent and time duration), in order to ensure their robustness.

We define the nominal shock position x_sh in the simulation frame as the position at which the magnitude of the magnetic field exceeds the upstream value by at least a factor of 3. The nominal shock position is then fitted in time to define an average shock frame (dashed magenta line in Figure <ref>). The nominal shock position is used to systematically distinguish between upstream and downstream of the shock. The shock speed in the simulation frame (which is also the average downstream frame) is around 2 v_A.



Figure <ref> shows the magnetic field magnitude at the final time of the simulation (top panel). The upstream region is defined as the simulation domain portion with x < x_sh - 25   d_i, and it reduces in size as the simulation evolves. The downstream region is defined as the simulation domain portion with x > x_sh + 80   d_i, which grows with the simulation time. The shock transition region, where there is the main change of flow velocity, is not considered in this work, because there is no well-defined, single flow frame for the analysis of particle diffusion regime.

Energy spectra for upstream and downstream populations are shown in the bottom panels of Figure <ref>. In both spectra, high energy, non-Maxwellian tails are present, extending up to energies corresponding to velocities of 80 v_A. Although this result is well known and has been widely investigated in previous literature <cit.>, it is important to underline that hybrid simulations of quasi-parallel shocks show significant particle acceleration, with energy gains of up to two orders of magnitude, even though the typical simulation times are shorter than the ones used in Monte Carlo simulations and other frameworks. The focus of this work is to find out more about the transport regime for particles in the high energy tail of the distribution, that is relevant to study the initial energisation of ions through their interaction with a quasi-parallel shock.

A pre-selection on the dataset of tracked particles is performed: we track particles having final energies above a threshold corresponding to a speed of 10 v_A in the downstream frame. Fragments of upstream trajectories tracked before the first shock encounter are ignored in order to exclude the initial convection of the particle in the upstream flow towards the shock. This has the effect that the initial position for all trajectory fragments will be at the origin in the shock frame.



Figure <ref> shows some examples of particle trajectories in the shock frame. Even at a first glance, it is possible to note that trajectories with completely different dynamics coexist in the simulation: some particles are transmitted downstream and moving away from the shock, other particles seem to move at the same speed of the shock in the region immediately downstream of it, some other ones (the most energetic) are reflected upstream and have multiple encounters with the shock transition.

Diffusion diagnostics can be biased when performed in a drifting frame. Hence, we will treat the upstream and downstream regions of the simulation as independent, with all the calculations performed in the local rest frame. However, the particles are free to travel in between regions so that a particle crossing from one region to another is treated as a new particle in the new region. This ultimately enables us to deal with two different trajectory datasets (or, more accurately, with two different datasets of trajectory fragments).

It is likely that in limiting the study of the particle trajectories to one side only, either upstream or downstream, a bias is introduced in the statistics. We note, however, that for quantities like the cycle time for shock acceleration what matters are the transport regimes and the sum of times spent on each side of the shock before re-crossing the shock <cit.>. 



 §.§ Anomalous diffusion diagnostics


The particle mean square displacement (MSD) is the most common tool to address the particle diffusion regime in a system. For one-dimensional motion, it is defined as 

    ⟨Δ s^2 (t) ⟩  = 1/N_p∑_n=0^N_p [x_n(t) - x_n(0)]^2 ,

where N_p is the total number of particles in the ensemble, x_n(t) and x_n(0) are the n-th particle positions at time t and at the initial time t=0, respectively. From now on, we will use the symbols ⟨ a⟩ to denote ensemble averages, and a to indicate time averages. Note that the MSD is an ensemble average, telling us information about the overall behaviour of the system. Obviously, the accuracy of MSD has a strong dependence on the number of particles used to build the ensemble. Also, in a system with open boundaries, N_p is not necessarily constant, and reduces in time.

In the case of anomalous diffusion the MSD is known to scale with time as

    ⟨Δ s^2 (t)⟩ = D_α (E)   t^α ,

where D_α(E) is the anomalous diffusion coefficient, which, in general, depends on the particle energy. The anomalous diffusion exponent α characterizes the diffusion regime, such that when α is equal to 1 Brownian diffusion is at play, whereas if α < 1 or α > 1 the system is subdiffusive or superdiffusive, respectively. It is possible, given an ensemble of particles for which the diffusion regime is not known, to infer a diffusion exponent, by calculating  ⟨Δ s^2 (t) ⟩  for the ensemble of particles, and then α can be estimated from a linear fit in log-log space.

Often, datasets of trajectories are not ideal for calculating the MSD, since difficulties can arise when dealing with poor particle statistics, leading to noisy ensemble averages. Furthermore, if the system is heterogeneous (e.g., different particle species present, non-uniform background environment), the MSD cannot account for this and simply gives an ensemble-averaged perspective on the overall system diffusion regime.

The time-averaged mean square displacement (TAMSD) method tries to address such problems by using single-particle trajectories and it is commonly used in many areas of biophysics <cit.>. In TAMSD, the ensemble average is replaced by an averaging time window:

    Δ s^2 (τ) = 1/T_f - τ∑_m=1^T_f -τ [x(m+τ) - x(m)]^2 ,

where τ is the averaging window, T_f is the final time for which the particle is tracked and x is the particle position in time. All times are considered to be discretised with uniform sampling. This leads to an estimator of the anomalous exponent based on single particle trajectories if Δ s^2 scales with time in the same way described by Equation <ref>. In the ideal case, for a dataset with infinite time records and for Brownian motion, TAMSD and MSD would produce the same result because of the central limit theorem.

The general procedure to apply TAMSD to a sample of trajectories is to fit an anomalous exponent α to each particle trajectory (as for MSD), and then take an average of the obtained exponents ⟨α⟩. This is known to have some problems that often lead to an underestimation of ⟨α⟩, deriving from the fact that new free parameters are introduced in the diagnostic <cit.>. For example, a time interval needs to be chosen for the fitting of the exponents, and this is nontrivial particularly when long-memory processes are present in the system. However, a set of particles interacting with a shock is intrinsically heterogeneous, and we suggest TAMSD as a tool to identify different ensembles of particles with different diffusion properties.

To have direct comparison between MSD and TAMSD, it is possible to consider an ensemble average of the sample's TAMSD <cit.>. We employ this diagnostic, called ensemble-averaged time-averaged mean square displacement (EA TAMSD), to the particle trajectories obtained in the shock simulation. The EA TAMSD is defined by 

    ⟨Δ s^2  (τ)⟩ = 1/N_p∑_p=1^N_p1/T_f - τ∑_m=1^T_f -τ [x_p(m+τ) - x_p(m)]^2 ,

that is the average of the TAMSD in Equation <ref> over the ensemble of tracked particles. As we will show in the results section, due to its double-averaging feature, the EA TAMSD is not only a way to validate the results found using MSD, but is also a very efficient diagnostic to investigate the overall (i.e., ensemble averaged) properties of the sample. Comparing MSD and EA TAMSD is also important because previous studies have shown that  MSD and EA TAMSD may differ when there is an ergodicity breaking in the sample <cit.>.

Even though this work is concerned with computing diffusion exponents from particle trajectories obtained from hybrid PIC shock simulations, MSD, TAMSD and EA TAMSD diagnostics have also been tested on a  Monte Carlo model, to ensure the robustness of the inferences made for the PIC simulations. Two examples of such tests are given in Appendix  <ref>. 




§ RESULTS




 §.§ Upstream particle transport


Results regarding the diffusion regimes for particles upstream of the shock are presented in this section. All the calculations shown here are done in the upstream rest frame. The particle sample considered consists only of particles initially reflected at the shock, as the trajectory fragments of particles before the first interaction with the shock are neglected. Throughout this analysis, particles are treated regardless of their energies even though the particle diffusion coefficient is known to be energy dependent, and this may affect the estimation of transport properties. However, this choice is motivated by the fact that the discussion below is focused on investigating the particle transport regime using different diagnostics (MSD, TAMSD, EA TAMSD) rather than to provide a discussion about particles' diffusion coefficients. However, an analysis of particle MSD for different energies is reported in Section <ref>.
It is worth noting that the upstream sample is the poorest from the statistical point of view, as the particles interacting with the shock have a relatively high probability to be transmitted downstream. On the other hand, the particles in the upstream region are likely to be efficiently accelerated, as they can gain energy through the multiple-shock encounter scenario discussed above and seen in the trajectories of Figure <ref>.

MSD and EA TAMSD are shown in Figure <ref>, which gives an overall view of the ensemble behaviour upstream of the shock.

Firstly, we observe that the overall behaviour of the two diagnostics is remarkably similar, validating the use of EA TAMSD. Note that the time τ shown when comparing MSD and EA TAMSD in Figure  <ref> has two different meanings: for the MSD it is the proper simulation time (see Equation <ref>) and for the EA TAMSD it is the length of the averaging time window, defined in Equation <ref>. 

At extremely short times (less than one Ω_ci^-1), a ballistic regime, with the diffusion exponent α =2 is identified in the MSD. This is expected since at short times scattering is still not taking place and, in principle, both diffusion and superdiffusion should be studied over sufficiently long times (i.e., for times longer than a few interaction times). Since the ballistic regime lasts for a very short time, the EA TAMSD does not capture well this ballistic transient, and this is related to the time average operation included in TAMSD.

At moderate times (from 6 to 45 Ω_ci^-1), there is a slight change of behaviour, and the MSD show a quasi-ballistic trend. A fit to the displacement in time gives values for the anomalous diffusion exponent of 1.86 and 1.93 when calculated from MSD and EA TAMSD, respectively. These values are both consistent with superdiffusion. At later times (from 50 to 450 Ω_ci^-1), another behaviour can be seen where MSD and EA TAMSD have exactly the same slope, and enter another superdiffusive regime, with an anomalous diffusion exponent α of 1.4. The reason for this change of slope may be that the most superdiffusive particles (i.e., the ones undergoing a very small number of scattering episodes) leave the box on short timescales.

It is important to underline that, for the upstream sample, EA TAMSD and MSD have a very similar behaviour and are almost indistinguishable for more than a decade. On one hand, this is a good test for the TAMSD diagnostic, which is a novel technique for particle simulations in collisionless plasma. On the other hand, this result is a validation for what was found using MSD.

TAMSD has the advantage that it gives an insight into the diffusive behaviour of single particles, in particular, diffusion exponents for each single particle trajectory.

Figure <ref> shows an example of how single particle TAMSD can be used for the upstream particles. A subset of around 4000 trajectory fragments of particles surviving for at least 200 Ω_ci^-1 in the upstream region has been selected. Single particle TAMSDs have been calculated using Equation <ref> and then single particle diffusion exponents have been fitted to obtain a distribution of single particle diffusion exponents α.  The histogram in Figure <ref> (top) shows the values of α obtained with TAMSD on these trajectory fragments. The fit excludes the last few points of the single particle TAMSD in order to reduce the statistical noise present when using large averaging windows (see TAMSD fluctuations at large τ in Figures <ref> & <ref>) and to obtain a robust estimation of the values of α.

The fluctuations that TAMSD shows at long time scales are due to poor statistics in the trajectory fragments. From Equation <ref>, for larger times fewer points contribute to the time averaging window, so if a particle returns close to its initial position the value of Δ s^2 drops drastically, inducing the fluctuation seen in the single particle TAMSD.

The mean value for the single particle α distribution  shown in Figure <ref> is approximately 1.75, i.e., between the transient with very strong superdiffusion (almost ballistic) and the long term superdiffusive behaviour of 1.4 found using MSD and EA TAMSD and shown in Figure <ref>. This value for ⟨α⟩ is still consistent with superdiffusion. We note that a transition between different anomalous transport regimes is also found in laboratory plasmas <cit.>. It is also worth noting that every single particle α value is consistent with superdiffusion.








 §.§ Downstream particle transport


In this section we present results obtained for the region downstream of the shock. All the calculations are done in the downstream rest frame. The downstream region is very different from that upstream since the plasma has been thermalised by the shock and the level of fluctuations is much higher. From the point of view of particle transport, particles with very different histories are found downstream. Some particles are transmitted and moving away from the shock, others are returning upstream (gaining energy in the process), and some others stay trapped in the strong magnetic field fluctuations present immediately after the shock transition. Some traces of this intrinsic heterogeneity are expected to be found when addressing the diffusion regime.

With the same motivation discussed for the upstream region, particles are treated regardless of their energy. A discussion about the downstream MSD behaviour at different energies is reported in Section <ref>. 

An ensemble of downstream particle trajectory fragments is built in the same fashion as described earlier. A lower boundary is set just after the shock transition zone (80 d_i downstream of the nominal shock transition). Another boundary is placed close to the reflecting wall (20 d_i in front) to remove any artificial scattering that the wall might produce, and when a particle crosses this “wall boundary” it is discarded. Every particle that crosses the left boundary to enter the downstream region is treated as a new particle trajectory fragment, and particles escaping the left boundary are discarded and tracked as new particles in the transition zone.

Figure <ref> shows a comparison between MSD and EA TAMSD for the downstream sample. The scenario obtained is clearly different from the upstream case. At very short times, the ballistic behaviour is recovered with MSD. Once again, the EA TAMSD does not capture well this regime because of the time average included in the diagnostic. At intermediate times (from 3 to 30 Ω_ci^-1), a temporarily superdiffusive regime is observed once again. The behaviour at intermediate times is different from the one found in the upstream particle sample. First of all, the diffusion exponents found with MSD and EA TAMSD, 1.65 and 1.85 respectively, are smaller than the typical value of 1.9 found upstream. This result can be interpreted considering that scattering for downstream particles starts to happen, statistically, at smaller time scales than the upstream case, and this is due to the higher level of fluctuations downstream of the shock. It is important to note that, for intermediate times, the diffusion exponents obtained with MSD and EA TAMSD differ by 0.2, which may be due to the heterogeneity of the sample and the process of ensemble averaging, as discussed later. 

For longer times (30 to 150 Ω_ci^-1), both diagnostics have the same slope, with an α very close to unity, indicating normal diffusion. The recovery of normal diffusion yields to two important considerations. On one hand, we find a diffusion regime that is compatible with the theoretical assumptions of DSA, but on the other hand, the result must be interpreted keeping in mind that the ensemble average involved in both the diagnostics averages out the behaviour of single particles or smaller populations of particles that may be in different transport regimes.

With single particle TAMSD, we can identify subsets of particles having different behaviours. Before discussing single particle TAMSD, it is worth  noting that, for times larger than 150 Ω_ci^-1,  Δ s^2 exhibits a large scale oscillation. This is an effect due to the increase of the finite size of the downstream region, confirmed by performing different simulations with different box sizes.     

Results for single particle TAMSD are shown in Figure <ref>. Again, a subset of particles that survive for at least 200 Ω_ci^-1 has been selected to perform this diagnostic. The single particle diffusion exponents have been fitted excluding the last few points of Δ s^2, where TAMSD become noisy (as for the upstream region). The histogram representing the α values obtained for the particle ensemble has a mean value of 1.27, representing superdiffusion. Although this value of ⟨α⟩ is smaller than that found upstream (c.f. ⟨α⟩ = 1.74), it is still indicative of superdiffusion with an exponent consistent with those obtained from data analysis <cit.>. Like the upstream case, the ⟨α⟩ obtained with TAMSD is in between the two regimes found with MSD and EA TAMSD (see Figure <ref>). 

When confronting the results obtained with different diagnostics, it is useful to note that EA TAMSD (defined in Equation <ref>) leads to an estimation of α over the ensemble, that may depend on time, and is robust against the noise of MSD. TAMSD, on the other hand, leads to an estimation of single particle α values without α varying in time, so it introduces the possibility of examining the heterogeneity of the ensemble. To this extent, when analysing the downstream particle transport properties, the key result is not the mean value of the TAMSD histogram, but rather its large spread over values covering the whole domain of possible anomalous transport exponents, that demonstrates that particles with completely different diffusion regimes coexist in the downstream region.

It can be seen in the bottom-right plot of Figure <ref> that single particle TAMSD have a remarkable spread, much larger than those observed in the upstream samples, ranging from subdiffusive regimes (α < 1) to superdiffusive (α > 1). We suggest that subdiffusive diffusion exponents represent particles trapped in the strong magnetic fluctuations present downstream of the shock (and in particular immediately after the shock transition). Superdiffusive single particle anomalous exponents can, on the other hand, correspond to particles transmitted downstream which move away from the shock without undergoing much scattering.




 §.§ TAMSD advantages and limitations



To demonstrate our interpretation of the many possible transport realizations seen in the downstream region, we show in Figure <ref> three examples of particles in different transport regimes.  

 On the left, a particle trapped downstream is shown. The α value obtained with single particle TAMSD is 0.7, indicating a subdiffusive behaviour. The same analysis is carried out in the middle plot, and gives a normal diffusive value of α, confirmed from the particle trajectory observed, that undergoes many scattering episodes but still samples more space than for the particle shown in the left panel. Finally, a particle with superdiffusive α value is shown (right panel), and it is an example from the population of particles moving away from the shock without much scattering. 
The diffusion properties of particle ensembles are, however, statistical properties, i.e., the single trajectories may behave very differently from the average over many particles. The examples presented in Figure <ref> demonstrates that TAMSD can be very powerful diagnostic when dealing with poor particle statistics and/or heterogenous particle samples.  



 §.§ MSD and energy dependence





As discussed previously, the analyses presented so far treated upstream and downstream particles regardless of their energy, although the diffusion coefficient depends from the particles' energy (see Equation <ref>). Figure <ref> shows upstream (a) and downstream (b) MSD calculations for particles in three different energy ranges. It should be noted that, since the shock transition zone is neglected in this study, and all the calculations are done in the local rest frame, the particles' energy remains roughly constant throughout the trajectory fragments analysed.

The dependence of the diffusion coefficient from the particles' energy in the upstream region is clearly found in Figure <ref> (left panels), where the ⟨Δ s^2 ⟩ is found to be systematically larger for particles with larger energies. The black line, showing the upstream MSD for the whole particle ensemble, represents an `average' of what obtained in the different energy ranges.

Figure <ref> shows that the MSD diagnostic for upstream particles is consistent with a superdiffusive scenario for all the energy ranges considered. It is worth noting that this result was predictable using single particle TAMSD upstream (demonstrating one of the strengths of TAMSD). In particular,  the histogram of single particle anomalous diffusion exponents reported in Figure <ref> shows that all the upstream particles have single particle diffusion exponent greater than one, indicating superdiffusion.

The transition from α = 1.9 to α = 1.4 is seen in all energy bands, although there is a trend to see this later at lower energies. The fact that the time at which the transition to α = 1.4  happens changes is also influenced by a finite-box size selection effect. At later times, the results for the lowest energy particles becomes less significant because the number of particles in the ensemble decreases, as particles reencounter the shock.



In the shock downstream region, the dependence of the diffusion coefficient on the particle energy is also clear from Figure <ref> (right hand panels). Once again, when the energy increases, the square displacement of particles becomes larger, and the MSD calculated for the entire particle ensemble represents an `average' between the behaviours observed at different energies.

When the energy filter is applied to the downstream particle ensemble, the regime corresponding to α∼ 1  at later times (discussed in Section <ref>), indicating normal diffusion, is not recovered. We suggest that this may be an effect related to the fact that the downstream particle ensemble is intrinsically heterogeneous.

Like for the upstream region, this analysis confirms what found using other diagnostics, namely that the shock downstream region has a much richer scenario as far as the particle transport properties are concerned.






§ SUMMARY AND CONCLUSIONS


Energetic particle transport regimes have been investigated in the framework of 1D, quasi-parallel shock hybrid PIC simulations. The simulation domain was divided into upstream and downstream regions, neglecting the transport regimes present within the shock transition. The transport regimes obtained for the two regions were found to have important differences. To the best of our knowledge, this is the first time that diffusion exponents have been extracted using particle trajectories advanced in the self-consistent fields produced with hybrid plasma simulations.

The usual technique to calculate anomalous diffusion exponents, based on the method of mean square displacement, relies on ensemble averages, and so it is very sensitive to poor particle statistics. It is often hard to track particles across the domain for long enough times in the framework of hybrid PIC simulations. Also, the dynamics of different particle populations are often very different in the shock simulation, and the ensemble average does not always capture this feature very well. For this reason, we have used a novel application of time-average mean square displacement together with the usual MSD method. TAMSD is a common diagnostic in scientific areas that have to deal with poor particle statistics, such as biophysics. It tries to overcome the problems of MSD by introducing a time averaging process on the single particle trajectories, and enables single particle estimates of the diffusion exponent.

In our simulation, we have found strong evidence for upstream superdiffusion. This evidence primarily arose from the calculation of MSD in the upstream frame for the reflected particles, in which we found extremely strong superdiffusion for intermediate times (with α ≈ 1.9), and a superdiffusive behaviour with an α of 1.4 for long times. The EA TAMSD, which is a novel technique for plasma simulations, is in extremely good agreement with the more classic MSD, although there are some differences between the two diagnostics at intermediate times (6 to 45 Ω_ci^-1), associated with transient behaviour. Single particle TAMSD enabled us to look at the statistical distribution of single particle anomalous diffusion exponents, in agreement with the other diagnostics. TAMSD is a powerful tool when we deal with poor particle statistics and with trajectory fragments of variable time extension.

Downstream of the shock, a much richer scenario of particle transport regimes is obtained. Both MSD and EA TAMSD show evidence for superdiffusion at intermediate times (6 to 30 Ω_ci^-1), even though the anomalous diffusion exponents inferred using these two diagnostics are different (and equal to 1.65 and 1.85, respectively). A possible explanation for this difference in the diagnostics at intermediate times lies in the intrinsic heterogeneity of the particle sample. For later times, a diffusion exponent closer to unity, suggesting a tendency to Brownian motion is recovered using both MSD and EA TAMSD. When TAMSD is performed on single particle trajectories for the downstream sample, the single particle diffusion exponents are found to have a mean of 1.27 (consistent with moderate superdiffusion). However, the spread of the histogram of single particle α values is much larger than the one found for the upstream sample, ranging from highly subdiffusive to highly superdiffusive behaviours. We consider that this is the indication of different particle populations coexisting in the same region: from particles trapped in the shock fluctuations to particles moving away from the shock without interacting with the surrounding environment, (see examples presented in Figure <ref>).

We have presented evidence of anomalous diffusion of energetic ions in a single 1D simulation of a collisionless quasi-parallel shock. Such simulations obviously have limitations due to finite domain size and simulation time, so the applicability of our results to other shocks (in terms of shock and/or simulation parameters) remains an open question.

Future work will utilise fully 2D and 3D simulations <cit.>, addressing diffusion in the magnetic field perpendicular and parallel directions. Moreover, it is possible to cross-correlate the diffusive regimes of different particles subsets with other important parameters, such as the energy gains, fluctuation levels, scattering times, or particle pitch angles and this will be the object of further investigation. We also intend to examine the transport regimes present in the shock transition layer, that have been neglected throughout this work, but which are probably crucial for understanding particle acceleration, especially at low energies. It is worth underlining that, throughout this work, we studied a shock with parameters compatible with those observed in the heliosphere. Further investigations will broaden the parameter space, in order to address whether these results are relevant for higher Mach number, astrophysical shocks. 



§ ACKNOWLEDGEMENTS

DT acknowledges support of a studentship funded by the Perren Fund of the University of London. This research was supported by the UK Science and Technology Facilities Council (STFC) grant ST/P000622/1. This research utilised Queen Mary's Apocrita HPC facility, supported by QMUL Research-IT, http://doi.org/10.5281/zenodo.438045.



mnras
=10000
 





§ TESTING TAMSD WITH MONTE CARLO MODELS



Although TAMSD is a common diagnostic used to address anomalous transport properties of many biological systems, it has never been used in the framework of numerical simulations of collisionless plasmas. On the other hand, self-consistent trajectories from hybrid PIC shock simulations have rarely been used to address the diffusion regimes at play in the plasma environment. Therefore, we have tested TAMSD on particle trajectories obtained using a simpler, Monte Carlo model for anomalous transport. In the Monte Carlo model, particles are advanced integrating a Langevin type equation in the form of a stochastic differential equation <cit.>:

    dx_i = v_ran dt_i .

The distribution of v_ran controls the diffusion regime, that is, scattering episodes for the particles are modelled through the distribution of v_ran. The random velocity is changed randomly at each scattering time t_i. In the case of superdiffusion, the scattering times t_i have a power-law-tailed probability distribution. The Monte Carlo model employed to generate the trajectories used for testing the anomalous diffusion diagnostics has extensively been used to bridge observations and models of anomalous transport in plasmas both for collisionless shocks and plasma turbulence <cit.>.

In this model, it is possible to set the exponent μ of the power-law distribution of scattering times. The μ exponent is related to the anomalous diffusion exponent through the relation <cit.>: 

    α =  4- μ .

The particle scattering times are then defined through the distribution 

    ψ(t) = 
    	
    	A       t ≤ t_0
    
    	A (t/t_0)^-μ    t > t_0
    	  ,

where A is a normalisation constant and t_0 is a free parameter. The details about this model can be found in <cit.>.

It is possible to compare the theoretical values for α with those obtained from the particle trajectories analysis using MSD, TAMSD and EA TAMSD. It is worth noting that for the self-consistent plasma model, the theoretical α is inferred from the analysis of particle trajectories (i.e., it is not a free parameter).

The first test uses a diffusion exponent α = 1, corresponding to normal (Brownian) diffusion. The equations of motion (Equation <ref>) are integrated for an ensemble of 10^4 particles, so that the ensemble averages do not suffer from statistical noise. The final time for the particle trajectories is of 36000 time units with a timestep of unit time.

Results are presented in Figure  <ref>. Results for MSD, TAMSD and EA TAMSD are consistent and capture the expected theoretical α used in the simulation. The single particle TAMSD exponents have been calculated by fitting the single particle TAMSD curves (Figure  <ref>, bottom-left) in the range of times before they start to become noisy. The slight underestimation of α using TAMSD (α_TAMSD∼ 0.96 vs α_theory = 1.0) is an effect from which the TAMSD diagnostic is known to suffer. We note that MSD and EA TAMSD capture the value of the expected α with a smaller error than single particle TAMSD and, moreover, they return the same value when α is fitted (α_MSD = α_EA   TAMSD∼ 0.99). The extremely strong agreement between MSD and EA TAMSD also has a theoretical explanation: when normal diffusion is achieved, the central limit theorem is satisfied, so time and ensemble averages must converge for asymptotic times. 

For comparison, a superdiffusive case has been analysed as well by means of the Monte Carlo method. Here, an anomalous diffusion exponent α of 1.5 (i.e., μ=2.5) has been imposed in Equation <ref>. The number of trajectories followed and the final simulation time are the same as the previous case. Figure  <ref> presents the results of the test. The single particle diffusion exponents have been fitted from single particle TAMSDs with the criterion described above. The mean single particle diffusion exponent ⟨α_TAMSD⟩ value is 1.43 and slightly underestimates the theoretical value of 1.5. Again, this is a known effect for the TAMSD diagnostic. However, the inferred ⟨α_TAMSD⟩ is in good agreement with the theoretical value, as can be also seen in the histogram shown in the top panel of Figure <ref>. 



The diffusion exponents obtained using MSD and EA TAMSD are both consistent with the expected α_theory within reasonable errors, although their slopes have been fitted for a shorter time window. The slight discrepancy between the MSD and the EA TAMSD curves, as well as the ∼ 1% difference in the calculated α values may be due to the fact that in the superdiffusive regime the central limit theorem does not hold anymore, and so we would expect some difference when taking time and ensemble averages.

The results of the tests are summarised in Table <ref>. It is worth noting that the diagnostics used to infer the anomalous diffusion exponents are consistent, but show slight differences. However, TAMSD is a flexible and powerful tool to use when dealing with datasets affected by poor particle statistics and trajectories with varying duration. The tests presented in this appendix help us understand the robustness of the diagnostics presented.


















