Particle transport in hybrid PIC shock simulations: A comparison of diagnostics
MMNRAS , 1–17 (2019) Preprint 30 September 2019 Compiled using MNRAS L A TEX style file v3.0
Particle transport in hybrid PIC shock simulations: Acomparison of diagnostics
D. Trotta (cid:63) , D. Burgess , G. Prete , S. Perri and G. Zimbardo Queen Mary University of London, School of Physics and Astronomy, London E1 4NS, UK Dipartimento di Fisica, University of Calabria, Italy
Accepted 2019 September 26. Received 2019 September 25; in original form 2019 March 4
ABSTRACT
Recent in-situ and remote observations suggest that the transport regime associatedwith 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, hybridPIC plasma simulations to simulate a quasi-parallel shock with parameters compati-ble with heliospheric shocks, and gain insights about the particle transport in such asystem. For suprathermal particles interacting with the shock we compute the MSDseparately in the upstream and downstream regions. Tracking suprathermal particlesfor sufficiently long times up and/or downstream of the shock poses problems in par-ticle plasma simulations, such as statistically poor particle ensembles and trajectoryfragments of variable length in time. Therefore, we introduce the use of time-averagedmean square displacement (TAMSD), which is based on single particle trajectories, asan additional technique to address the transport regime for the upstream and down-stream regions. MSD and TAMSD are in agreement for the upstream energetic particlepopulation, and both give a strong indication of superdiffusive transport, consistentwith interplanetary shock observations. MSD and TAMSD are also in reasonable agree-ment downstream, where indications of anomalous transport are also found. TAMSDshows evidence of heterogeneity in the diffusion properties of the downstream parti-cle population, ranging from subdiffusive behaviour of particles trapped in the strongmagnetic field fluctuations generated at the shock, to superdiffusive behaviour of par-ticles transmitted and moving away from the shock.
Key words:
Diffusion – shock waves – acceleration of particles
High energy particle acceleration is a common feature ofmany astrophysical systems, and characterizing the ener-getic particle transport is a key challenge for understand-ing the underlying physical acceleration processes. In mostcases, when particle trajectories are available, the method ofmean square displacement (MSD) is the most widely useddiagnostic to study particle transport. MSD uses ensembleaverages, and hence can be unreliable with poor particlestatistics. On the other hand, the method of time-averagedmean square displacement (TAMSD) relies on time ratherthan ensemble averages, and is a popular diagnostic in re-search communities concerned with particle transport in bi-ological systems. When dealing with particle tracking in self-consistent plasma simulations of acceleration it can be chal-lenging to obtain sufficiently good particle statistics over (cid:63)
E-mail: [email protected] long times, due to the non-periodicity and limited extent ofsuch simulations and consequently the varying lifetimes ofdifferent particles in the simulation. For this reason, we in-vestigate the use of TAMSD diagnostics to study the trans-port of accelerated particles in self-consistent hybrid PICshock simulations, and make comparisons with the resultsfrom the MSD method.It is widely accepted that shocks are efficient parti-cle accelerators in astrophysical systems (e.g., Blandford &Eichler 1987; Fisk 2015). Generally speaking, shocks convertdirected flow energy (upstream) to thermal energy (down-stream). In collisionless shocks, a small fraction of the up-stream energy can be channeled towards particle accelera-tion. The mechanisms by which particles are accelerated outof the thermal pool and to high energies are still debated.One of the most important parameters that controls shockstructure and behaviour in a collisionless plasma is the anglebetween the normal to the shock surface and the upstreammagnetic field, θ Bn . Shocks with θ Bn close to zero are called © a r X i v : . [ phy s i c s . s p ace - ph ] S e p Trotta et al. quasi-parallel, and they are believed to be more efficient atparticle acceleration than quasi-perpendicular ones (see, forexample, Burgess & Scholer 2015). Other important param-eters are the shock Alfv´enic and sonic Mach numbers, de-fined as M A = v sh / v A and M S = v sh / c s , respectively, and theplasma β = v / v A . Here, v A is the Alfv´en speed in the re-gion upstream from the shock, v sh is the shock speed in theupstream flow frame, c s is the sound speed of the upstreamflow and v th is the thermal speed of the upstream flow.Particle acceleration at shocks is invoked to explain awide variety of direct and indirect observations. Direct evi-dence of shock accelerated particles has been found by meansof spacecraft measurements since the 1960s (e.g, Bryantet al. 1962; Scholer et al. 1982; van Nes et al. 1984) forinterplanetary shocks. At the solar wind termination shock,the suprathermal particle population was found enhanced by Voyager spacecraft observations (Decker et al. 2005; Deckeret al. 2008). Indirect evidence of accelerated particles is com-monly found in the X-ray signature of supernova remnants(see Reynolds 2008, for a review). Recently, particle accel-eration at shocks has also been found to be relevant in thegalaxy intracluster medium (e.g. van Weeren et al. 2010,2017; Kang et al. 2017).A widely used theory for particle acceleration at quasi-parallel shocks is Diffusive Shock Acceleration theory (DSA)(Krymskii 1977; Axford et al. 1978; Bell 1978; Blandford &Ostriker 1978; Drury 1983; Jokipii 1987). The idea under-lying DSA theory is that particles with Larmor radii largerthan the shock thickness interact with the shock by diffus-ing upstream and downstream, and thus are continuouslyreturned to the shock, gaining energy at each crossing. Thisprocess is often called first order Fermi acceleration. Themain prediction of DSA is a power-law distribution in energyfor the accelerated particles, with a slope that depends onlyon the shock compression ratio, at least in the simple timestationary limit. DSA theory also predicts that, in the caseof a spatially constant diffusion coefficient, the density pro-files of energetic particles should decrease exponentially up-stream of the shock, and to be constant downstream. Thesepredictions have been found to be not inconsistent with someobservations of interplanetary shocks (e.g. Giacalone 2012),even though many in-situ spacecraft crossings exhibit morecomplex energetic particle profiles (see van Nes et al. 1984).Recently, some energetic particle observations have sug-gested that the transport regime associated with such par-ticles can be anomalous. In fact, energetic particle densityprofiles at interplanetary shocks are often found to showa power-law decay rather than an exponential one, as onewould expect from normal diffusion (Perri & Zimbardo 2007,2008, 2009a,b; Sugiyama & Shiota 2011; Perri & Zimbardo2015). In addition, it has been proposed, among other pos-sible mechanisms, that the discrepancy between Mach num-bers obtained from radio and from X-ray observations atgalaxy cluster merger shocks may be due to an extensionof DSA to superdiffusive transport (Zimbardo & Perri 2017,2018), which leads to superdiffusive shock acceleration (Perri& Zimbardo 2012a). Therefore, it is now timely to under-stand under what conditions the transport of energetic par-ticles at shocks could be considered to be superdiffusive.One of the main assumptions underlying DSA is thatenergetic particles undergo normal (i.e., Brownian) diffusionwhen scattering back and forth across the shock front. Brow- nian diffusion is characterized by Gaussian statistics. For anensemble of particles undergoing normal diffusion, the meansquare displacement for particle trajectories grows linearlyin time, as (cid:104) ∆ s (cid:105) = Dt , where D identifies the so-called dif-fusion coefficient.If the time evolution of the MSD is nonlinear, forexample, such that (cid:104) ∆ s (cid:105) = D α t α , then the system ex-hibits anomalous transport. The exponent α is known asthe anomalous diffusion exponent and its value identifies thetransport regime of the system. Anomalous transport hasbeen observed in an enormous variety of systems, rangingfrom the foraging movements of spider monkeys (Ramos-Fern´andez et al. 2004) to stock options transactions in fi-nance (Meerschaert & Scalas 2006). It is possible to pre-dict, upstream of shock waves, power-law density profiles forshock accelerated particles assuming for them an anomalousdiffusion regime (Perri & Zimbardo 2007, 2008). In fact, animportant feature of anomalous transport is that the sys-tems are characterized by non-Gaussian statistics, which in-troduces probability distributions with power-law tails (e.g.Klafter et al. 1987; Metzler & Klafter 2000). In DSA, thetime scale for acceleration of particles of given energy is dom-inated by the slowest portion of the shock acceleration cy-cle (e.g., Drury 1983). An interesting implication of anoma-lous, superdiffusive transport is that it may possibly leadto shorter acceleration times than those usually obtainedfor normal diffusive transport, thus helping to reach highermaximum energies than those predicted by DSA (Perri &Zimbardo 2015).Simulations of particle acceleration at shocks have beenused extensively to bridge theory and observations. Numer-ical and analytical solutions for the transport equation ob-tained from DSA have been studied under many approx-imations (e.g. Malkov 1997; Kang & Jones 1997; Caprioli2012). A popular approach to address energisation efficiencyand particle transport is the use of Monte Carlo simula-tions (e.g. Ellison et al. 1990; Jones & Ellison 1991; Baringet al. 1995; Wolff & Tautz 2015; Bykov et al. 2017). A draw-back of the usual Monte Carlo approach is the neglect ofany self-consistent plasma effects (e.g., internal shock struc-ture) which are believed to be important, particularly atlow energies (e.g. Sundberg et al. 2016). Efficient particleacceleration has been studied both in fully kinetic Particle-In-Cell (PIC) and hybrid (particle ions and fluid electrons)simulations (e.g. Amano & Hoshino 2007, 2009; Riquelme& Spitkovsky 2011; Giacalone et al. 1992; Giacalone & Elli-son 2000; Caprioli & Spitkovsky 2014a). In the fully kineticmodel, the spatial scales are resolved down to the Debyelength, making it possible to address the complete pictureof particle acceleration down to electron dynamics. However,fully kinetic PIC simulations are computationally intensive.On the other hand, in the hybrid model, the electrons aremodelled as a fluid and just the ion kinetic scales are re-solved, with the advantage that considerably larger domainsand time scales can be simulated.Particle diffusion has also been studied in the frame-work of self-consistent plasma simulations. Servidio et al.(2016), for example, studied the transport of protons in dif-ferent turbulent regions by means of hybrid simulations. Themajority of plasma turbulence simulations are performed ona periodic domain and in the plasma rest frame, makingit relatively easy to trace particle trajectories throughout MNRAS000
E-mail: [email protected] long times, due to the non-periodicity and limited extent ofsuch simulations and consequently the varying lifetimes ofdifferent particles in the simulation. For this reason, we in-vestigate the use of TAMSD diagnostics to study the trans-port of accelerated particles in self-consistent hybrid PICshock simulations, and make comparisons with the resultsfrom the MSD method.It is widely accepted that shocks are efficient parti-cle accelerators in astrophysical systems (e.g., Blandford &Eichler 1987; Fisk 2015). Generally speaking, shocks convertdirected flow energy (upstream) to thermal energy (down-stream). In collisionless shocks, a small fraction of the up-stream energy can be channeled towards particle accelera-tion. The mechanisms by which particles are accelerated outof the thermal pool and to high energies are still debated.One of the most important parameters that controls shockstructure and behaviour in a collisionless plasma is the anglebetween the normal to the shock surface and the upstreammagnetic field, θ Bn . Shocks with θ Bn close to zero are called © a r X i v : . [ phy s i c s . s p ace - ph ] S e p Trotta et al. quasi-parallel, and they are believed to be more efficient atparticle acceleration than quasi-perpendicular ones (see, forexample, Burgess & Scholer 2015). Other important param-eters are the shock Alfv´enic and sonic Mach numbers, de-fined as M A = v sh / v A and M S = v sh / c s , respectively, and theplasma β = v / v A . Here, v A is the Alfv´en speed in the re-gion upstream from the shock, v sh is the shock speed in theupstream flow frame, c s is the sound speed of the upstreamflow and v th is the thermal speed of the upstream flow.Particle acceleration at shocks is invoked to explain awide variety of direct and indirect observations. Direct evi-dence of shock accelerated particles has been found by meansof spacecraft measurements since the 1960s (e.g, Bryantet al. 1962; Scholer et al. 1982; van Nes et al. 1984) forinterplanetary shocks. At the solar wind termination shock,the suprathermal particle population was found enhanced by Voyager spacecraft observations (Decker et al. 2005; Deckeret al. 2008). Indirect evidence of accelerated particles is com-monly found in the X-ray signature of supernova remnants(see Reynolds 2008, for a review). Recently, particle accel-eration at shocks has also been found to be relevant in thegalaxy intracluster medium (e.g. van Weeren et al. 2010,2017; Kang et al. 2017).A widely used theory for particle acceleration at quasi-parallel shocks is Diffusive Shock Acceleration theory (DSA)(Krymskii 1977; Axford et al. 1978; Bell 1978; Blandford &Ostriker 1978; Drury 1983; Jokipii 1987). The idea under-lying DSA theory is that particles with Larmor radii largerthan the shock thickness interact with the shock by diffus-ing upstream and downstream, and thus are continuouslyreturned to the shock, gaining energy at each crossing. Thisprocess is often called first order Fermi acceleration. Themain prediction of DSA is a power-law distribution in energyfor the accelerated particles, with a slope that depends onlyon the shock compression ratio, at least in the simple timestationary limit. DSA theory also predicts that, in the caseof a spatially constant diffusion coefficient, the density pro-files of energetic particles should decrease exponentially up-stream of the shock, and to be constant downstream. Thesepredictions have been found to be not inconsistent with someobservations of interplanetary shocks (e.g. Giacalone 2012),even though many in-situ spacecraft crossings exhibit morecomplex energetic particle profiles (see van Nes et al. 1984).Recently, some energetic particle observations have sug-gested that the transport regime associated with such par-ticles can be anomalous. In fact, energetic particle densityprofiles at interplanetary shocks are often found to showa power-law decay rather than an exponential one, as onewould expect from normal diffusion (Perri & Zimbardo 2007,2008, 2009a,b; Sugiyama & Shiota 2011; Perri & Zimbardo2015). In addition, it has been proposed, among other pos-sible mechanisms, that the discrepancy between Mach num-bers obtained from radio and from X-ray observations atgalaxy cluster merger shocks may be due to an extensionof DSA to superdiffusive transport (Zimbardo & Perri 2017,2018), which leads to superdiffusive shock acceleration (Perri& Zimbardo 2012a). Therefore, it is now timely to under-stand under what conditions the transport of energetic par-ticles at shocks could be considered to be superdiffusive.One of the main assumptions underlying DSA is thatenergetic particles undergo normal (i.e., Brownian) diffusionwhen scattering back and forth across the shock front. Brow- nian diffusion is characterized by Gaussian statistics. For anensemble of particles undergoing normal diffusion, the meansquare displacement for particle trajectories grows linearlyin time, as (cid:104) ∆ s (cid:105) = Dt , where D identifies the so-called dif-fusion coefficient.If the time evolution of the MSD is nonlinear, forexample, such that (cid:104) ∆ s (cid:105) = D α t α , then the system ex-hibits anomalous transport. The exponent α is known asthe anomalous diffusion exponent and its value identifies thetransport regime of the system. Anomalous transport hasbeen observed in an enormous variety of systems, rangingfrom the foraging movements of spider monkeys (Ramos-Fern´andez et al. 2004) to stock options transactions in fi-nance (Meerschaert & Scalas 2006). It is possible to pre-dict, upstream of shock waves, power-law density profiles forshock accelerated particles assuming for them an anomalousdiffusion regime (Perri & Zimbardo 2007, 2008). In fact, animportant feature of anomalous transport is that the sys-tems are characterized by non-Gaussian statistics, which in-troduces probability distributions with power-law tails (e.g.Klafter et al. 1987; Metzler & Klafter 2000). In DSA, thetime scale for acceleration of particles of given energy is dom-inated by the slowest portion of the shock acceleration cy-cle (e.g., Drury 1983). An interesting implication of anoma-lous, superdiffusive transport is that it may possibly leadto shorter acceleration times than those usually obtainedfor normal diffusive transport, thus helping to reach highermaximum energies than those predicted by DSA (Perri &Zimbardo 2015).Simulations of particle acceleration at shocks have beenused extensively to bridge theory and observations. Numer-ical and analytical solutions for the transport equation ob-tained from DSA have been studied under many approx-imations (e.g. Malkov 1997; Kang & Jones 1997; Caprioli2012). A popular approach to address energisation efficiencyand particle transport is the use of Monte Carlo simula-tions (e.g. Ellison et al. 1990; Jones & Ellison 1991; Baringet al. 1995; Wolff & Tautz 2015; Bykov et al. 2017). A draw-back of the usual Monte Carlo approach is the neglect ofany self-consistent plasma effects (e.g., internal shock struc-ture) which are believed to be important, particularly atlow energies (e.g. Sundberg et al. 2016). Efficient particleacceleration has been studied both in fully kinetic Particle-In-Cell (PIC) and hybrid (particle ions and fluid electrons)simulations (e.g. Amano & Hoshino 2007, 2009; Riquelme& Spitkovsky 2011; Giacalone et al. 1992; Giacalone & Elli-son 2000; Caprioli & Spitkovsky 2014a). In the fully kineticmodel, the spatial scales are resolved down to the Debyelength, making it possible to address the complete pictureof particle acceleration down to electron dynamics. However,fully kinetic PIC simulations are computationally intensive.On the other hand, in the hybrid model, the electrons aremodelled as a fluid and just the ion kinetic scales are re-solved, with the advantage that considerably larger domainsand time scales can be simulated.Particle diffusion has also been studied in the frame-work of self-consistent plasma simulations. Servidio et al.(2016), for example, studied the transport of protons in dif-ferent turbulent regions by means of hybrid simulations. Themajority of plasma turbulence simulations are performed ona periodic domain and in the plasma rest frame, makingit relatively easy to trace particle trajectories throughout MNRAS000 , 1–17 (2019) article transport in shock simulations the simulations and compute the MSD for them. However,for studying energetic particles in shock simulations the sit-uation is more complex since only a small fraction of theparticles are accelerated to high energies, they may exit thesimulation domain because it is nonperiodic and there is nosingle rest frame due to the velocity change at the shock.Although these particular features make it harder totrack energetic particles in the simulations for sufficientlylong times, early studies attempted to study particle trans-port at shocks using self-consistent simulations. For in-stance, Kucharek et al. (2000) studied particle diffusion inturbulence generated by an ion-ion beam instability and in-voked it as a model for the Earth’s foreshock. Scholer et al.(2000) used hybrid PIC simulations of turbulence to modelthe immediate downstream of a perpendicular shock (byusing a particular initial condition, consisting of a homo-geneous plasma core distribution and a non-gyrotropic ionpopulation) and addressed the particle transport in such asystem. Caprioli & Spitkovsky (2014b) considered the shocktransition layer of a hybrid PIC shock simulation at a sin-gle time, imposed boundary conditions on it and evolvedtest particles in the resulting electromagnetic fields to studydiffusion across the shock layer.Our goal is to study the particle diffusion regimes atquasi-parallel shocks using self-consistent, hybrid PIC simu-lations, with the aim of gaining extra information that mightadvance our understanding of particle transport and so im-prove Monte Carlo modelling of shock acceleration.We fo-cus on a set of shock parameters relevant for heliosphericshocks and intermediate particle energies. As we mentionedabove, it is not trivial to follow many particles for longtimes in a shock simulation, and this has motivated us toemploy a novel diagnostic, namely the time-averaged meansquare displacement method (TAMSD). TAMSD relies ontime-averaging rather than ensemble-averaging as in MSD,making it possible to extract single particle diffusion expo-nents from their trajectories (Qian et al. 1991). TAMSD hasbeen used in a variety of systems where it can be difficultto obtain many particle trajectories to construct an ensem-ble average; this is the case for many biological systems,especially those for which in-vivo tracking of particles is re-quired (see, for example H¨ofling & Franosch 2013; Weronet al. 2017; Metzler et al. 2014, and references therein).In this paper the trajectories for an ensemble of pro-tons self-consistently interacting with a 1D, quasi-parallelshock are studied. Upstream and downstream diffusion prop-erties are analysed separately, and all the computations areperformed in the respective local rest frame. We obtainupstream and downstream diffusion exponents using bothMSD and TAMSD. A thorough comparison between thetwo diagnostics is presented, in order to test the applica-bility of TAMSD. We have found the two techniques to beconsistent, although both MSD and TAMSD have their ownlimitations. Upstream, strong evidence for superdiffusion isfound. In the shock downstream, we observe a richer sce-nario, encompassing superdiffusion and subdiffusion, and,based on the TAMSD diagnostics we suggest that this isrelated to different subsets of particles co-existing in dif-ferent diffusion regimes. The paper is organised as follows:in Section 2 we describe both the simulation method andthe diagnostics employed. In Section 3 results for upstreamand downstream transport are described, and summarised in Section 4. In Appendix A we present some independent testsof the TAMSD diagnostic using Monte Carlo simulations. We simulate a shock using a hybrid PIC simulation. Inthe simulation, protons are modelled as macroparticles andadvanced 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. Thehybrid code used throughout this work (HYPSI) is basedon the CAM-CL algorithm (Matthews 1994; Sundberg et al.2016).The shock is initiated by the injection method, in whichthe plasma flows in the x -direction with a defined (super-Alfv´enic) velocity V in . The right-hand boundary of the sim-ulation domain acts as a reflecting wall, and at the left-handboundary plasma is continuously injected. The simulation isperiodic in the y direction. A shock is created as a conse-quence of reflection at the wall, and it propagates in thenegative x -direction. In the simulation frame, the upstreamflow is along the shock normal.In the hybrid simulations distances are normalised tothe ion inertial length d i ≡ c / ω pi , times to the inverse cy-clotron frequency Ω ci − , velocity to the Alfv´en speed v A (allreferred to the upstream state), and the magnetic field anddensity to their upstream values, B and n , respectively. Theangle between the shock normal and the upstream magneticfield, θ Bn is 15 ◦ , with the upstream magnetic field in the x - y plane. For the upstream flow velocity, a value of V in = v A has been chosen, and the resulting Alfv´enic Mach number ofthe shock is approximately M A = . The upstream ion distri-bution function is an isotropic Maxwellian and the ion β i is0.5. The simulation x − y domain is 3000 × d i , and hencewe describe the shock as quasi-1D. The spatial resolutionused is ∆ x = ∆ y = 0.5 d i . The final time for the simulationis 750 Ω − ci , the time step for particle (ion) advance is ∆ t =0.005 Ω − ci . Substepping is used for the magnetic field ad-vance, with an effective time step of ∆ t B = ∆ t / . A small,nonzero resistivity is introduced in the magnetic inductionequation. The value of the resistivity is chosen so that thereare not excessive fluctuations at the grid scale. The numberof particles per cell used is always greater than 300 (up-stream), in order to keep the statistical noise characteristicof PIC simulations to a reasonable level.The number of particles in the simulation is about 10 ,and we store the trajectories of about 5 × particles. Thetime cadence for the particle trajectory outputs is chosensuch that the time resolution obtained in the particle datasetis ∆ t = 0.1 Ω − ci . The protons are tracked in time startingat 250 Ω − ci , when the shock is already well-developed, withan upstream region with developed fluctuations and a well-defined downstream region separated from the right handwall.An overview of the shock simulation is presented in Fig-ure 1, which shows the time evolution of | B | as a function of x position. Three typical trajectories for protons interactingwith the shock are also shown (red solid lines). Throughoutthe simulation, there are upstream waves self-consistently MNRAS , 1–17 (2019)
Trotta et al.
Figure 1.
Total magnetic field as a function of space ( x -axis) and time ( y -axis). The colorbar shows the total magnetic field magnitudein logarithmic scale. The magenta dotted-dashed line identifies the shock nominal position. Red lines represent three examples of protonsinteracting with the shock. induced by the reflected ion population, that have been ex-tensively studied in previous literature (e.g., Sundberg et al.2016). The downstream region exhibits strong fluctuationsat several wavelengths, and some of them are relatively co-herent in time, an effect related to the quasi-1D nature ofthe simulation.As can be seen in Figure 1, the upstream waves are ar-tificially truncated. This effect, particularly evident at latertimes, is due to the finite size of the simulation domain,and is known to affect shock simulations. It can be com-mented that this is not just a spatial truncation, but alsoan amplitude truncation, because the extent of the currentand the growth (advection) time is truncated. Also, if ionswith a given energy escape before properly producing waves,it is likely that particles with the same energy will see anon-fully-developed set of fluctuations, and this could affecttheir transport properties. However, the results presentedhere have been confirmed using larger simulations (in spatialextent and time duration), in order to ensure their robust-ness.We define the nominal shock position x sh in the simu-lation frame as the position at which the magnitude of themagnetic field exceeds the upstream value by at least a fac-tor of 3. The nominal shock position is then fitted in time todefine an average shock frame (dashed magenta line in Fig-ure 1). The nominal shock position is used to systematicallydistinguish between upstream and downstream of the shock.The shock speed in the simulation frame (which is also theaverage downstream frame) is around 2 v A .Figure 2 shows the magnetic field magnitude at the finaltime of the simulation (top panel). The upstream region is defined as the simulation domain portion with x < x sh − d i ,and it reduces in size as the simulation evolves. The down-stream region is defined as the simulation domain portionwith x > x sh + d i , which grows with the simulation time.The shock transition region, where there is the main changeof flow velocity, is not considered in this work, because thereis no well-defined, single flow frame for the analysis of par-ticle diffusion regime.Energy spectra for upstream and downstream popula-tions are shown in the bottom panels of Figure 2. In bothspectra, high energy, non-Maxwellian tails are present, ex-tending up to energies corresponding to velocities of 80 v A .Although this result is well known and has been widely in-vestigated in previous literature (e.g., Scholer & Terasawa1990; Giacalone 2004; Sugiyama 2011; Gargat´e & Spitkovsky2012; Caprioli & Spitkovsky 2014a), it is important to un-derline that hybrid simulations of quasi-parallel shocks showsignificant particle acceleration, with energy gains of up totwo orders of magnitude, even though the typical simulationtimes are shorter than the ones used in Monte Carlo simu-lations and other frameworks. The focus of this work is tofind out more about the transport regime for particles in thehigh energy tail of the distribution, that is relevant to studythe initial energisation of ions through their interaction witha quasi-parallel shock.A pre-selection on the dataset of tracked particles isperformed: we track particles having final energies above athreshold corresponding to a speed of 10 v A in the down-stream frame. Fragments of upstream trajectories trackedbefore the first shock encounter are ignored in order to ex-clude the initial convection of the particle in the upstream MNRAS000
Total magnetic field as a function of space ( x -axis) and time ( y -axis). The colorbar shows the total magnetic field magnitudein logarithmic scale. The magenta dotted-dashed line identifies the shock nominal position. Red lines represent three examples of protonsinteracting with the shock. induced by the reflected ion population, that have been ex-tensively studied in previous literature (e.g., Sundberg et al.2016). The downstream region exhibits strong fluctuationsat several wavelengths, and some of them are relatively co-herent in time, an effect related to the quasi-1D nature ofthe simulation.As can be seen in Figure 1, the upstream waves are ar-tificially truncated. This effect, particularly evident at latertimes, is due to the finite size of the simulation domain,and is known to affect shock simulations. It can be com-mented that this is not just a spatial truncation, but alsoan amplitude truncation, because the extent of the currentand the growth (advection) time is truncated. Also, if ionswith a given energy escape before properly producing waves,it is likely that particles with the same energy will see anon-fully-developed set of fluctuations, and this could affecttheir transport properties. However, the results presentedhere have been confirmed using larger simulations (in spatialextent and time duration), in order to ensure their robust-ness.We define the nominal shock position x sh in the simu-lation frame as the position at which the magnitude of themagnetic field exceeds the upstream value by at least a fac-tor of 3. The nominal shock position is then fitted in time todefine an average shock frame (dashed magenta line in Fig-ure 1). The nominal shock position is used to systematicallydistinguish between upstream and downstream of the shock.The shock speed in the simulation frame (which is also theaverage downstream frame) is around 2 v A .Figure 2 shows the magnetic field magnitude at the finaltime of the simulation (top panel). The upstream region is defined as the simulation domain portion with x < x sh − d i ,and it reduces in size as the simulation evolves. The down-stream region is defined as the simulation domain portionwith x > x sh + d i , which grows with the simulation time.The shock transition region, where there is the main changeof flow velocity, is not considered in this work, because thereis no well-defined, single flow frame for the analysis of par-ticle diffusion regime.Energy spectra for upstream and downstream popula-tions are shown in the bottom panels of Figure 2. In bothspectra, high energy, non-Maxwellian tails are present, ex-tending up to energies corresponding to velocities of 80 v A .Although this result is well known and has been widely in-vestigated in previous literature (e.g., Scholer & Terasawa1990; Giacalone 2004; Sugiyama 2011; Gargat´e & Spitkovsky2012; Caprioli & Spitkovsky 2014a), it is important to un-derline that hybrid simulations of quasi-parallel shocks showsignificant particle acceleration, with energy gains of up totwo orders of magnitude, even though the typical simulationtimes are shorter than the ones used in Monte Carlo simu-lations and other frameworks. The focus of this work is tofind out more about the transport regime for particles in thehigh energy tail of the distribution, that is relevant to studythe initial energisation of ions through their interaction witha quasi-parallel shock.A pre-selection on the dataset of tracked particles isperformed: we track particles having final energies above athreshold corresponding to a speed of 10 v A in the down-stream frame. Fragments of upstream trajectories trackedbefore the first shock encounter are ignored in order to ex-clude the initial convection of the particle in the upstream MNRAS000 , 1–17 (2019) article transport in shock simulations x [ d i ] | B | [ B ] τ = 750 Ω − ci -1 E p [ m p v A ] -6 -5 -4 -3 -2 -1 N ( E ) Upstream -1 E p [ m p v A ] -6 -5 -4 -3 -2 -1 N ( E ) Downstream
Figure 2.
Top: Magnetic field magnitude (averaged over y ) at the final time of the simulation. Bottom: Upstream (left panel) anddownstream (right panel) particle energy spectra calculated in the upstream and downstream rest frames, respectively. flow towards the shock. This has the effect that the initialposition for all trajectory fragments will be at the origin inthe shock frame.Figure 3 shows some examples of particle trajectoriesin the shock frame. Even at a first glance, it is possibleto note that trajectories with completely different dynam-ics coexist in the simulation: some particles are transmitteddownstream and moving away from the shock, other par-ticles seem to move at the same speed of the shock in theregion immediately downstream of it, some other ones (themost energetic) are reflected upstream and have multipleencounters with the shock transition.Diffusion diagnostics can be biased when performed in adrifting frame. Hence, we will treat the upstream and down-stream regions of the simulation as independent, with all thecalculations performed in the local rest frame. However, theparticles are free to travel in between regions so that a par-ticle crossing from one region to another is treated as a newparticle 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 particletrajectories to one side only, either upstream or downstream,a bias is introduced in the statistics. We note, however, thatfor quantities like the cycle time for shock acceleration whatmatters are the transport regimes and the sum of times spenton each side of the shock before re-crossing the shock (Drury1983; Perri & Zimbardo 2012b; Zimbardo & Perri 2013). The particle mean square displacement (MSD) is the mostcommon tool to address the particle diffusion regime in asystem. For one-dimensional motion, it is defined as (cid:104) ∆ s ( t )(cid:105) = N p N p (cid:213) n = [ x n ( t ) − x n ( )] , (1) MNRAS , 1–17 (2019)
Trotta et al.
Figure 3.
Typical particle trajectories extracted from the sample, in the shock frame. The colour bar shows the particle energy inlogarithmic scale and in units of m p v A . The dashed magenta line represents the nominal shock position. where N p is the total number of particles in the ensemble, x n ( t ) and x n ( ) are the n -th particle positions at time t andat the initial time t = , respectively. From now on, we willuse the symbols (cid:104) a (cid:105) to denote ensemble averages, and a toindicate time averages. Note that the MSD is an ensembleaverage, telling us information about the overall behaviourof the system. Obviously, the accuracy of MSD has a strongdependence on the number of particles used to build theensemble. Also, in a system with open boundaries, N p is notnecessarily constant, and reduces in time.In the case of anomalous diffusion the MSD is knownto scale with time as (cid:104) ∆ s ( t )(cid:105) = D α ( E ) t α , (2)where D α ( E ) is the anomalous diffusion coefficient, which,in general, depends on the particle energy. The anomalousdiffusion exponent α characterizes the diffusion regime, suchthat when α is equal to 1 Brownian diffusion is at play,whereas if α < or α > the system is subdiffusive orsuperdiffusive, respectively. It is possible, given an ensembleof particles for which the diffusion regime is not known, toinfer a diffusion exponent, by calculating (cid:104) ∆ s ( t )(cid:105) for the ensemble of particles, and then α can be estimated from alinear fit in log-log space.Often, datasets of trajectories are not ideal for calculat-ing the MSD, since difficulties can arise when dealing withpoor particle statistics, leading to noisy ensemble averages.Furthermore, if the system is heterogeneous (e.g., differentparticle species present, non-uniform background environ-ment), the MSD cannot account for this and simply givesan ensemble-averaged perspective on the overall system dif-fusion regime.The time-averaged mean square displacement(TAMSD) method tries to address such problems byusing single-particle trajectories and it is commonly used inmany areas of biophysics (e.g. Golding & Cox 2004; Wonget al. 2004; Saxton 2012). In TAMSD, the ensemble averageis replaced by an averaging time window: ∆ s ( τ ) = T f − τ T f − τ (cid:213) m = [ x ( m + τ ) − x ( m )] , (3)where τ is the averaging window, T f is the final time forwhich the particle is tracked and x is the particle position in MNRAS000
Typical particle trajectories extracted from the sample, in the shock frame. The colour bar shows the particle energy inlogarithmic scale and in units of m p v A . The dashed magenta line represents the nominal shock position. where N p is the total number of particles in the ensemble, x n ( t ) and x n ( ) are the n -th particle positions at time t andat the initial time t = , respectively. From now on, we willuse the symbols (cid:104) a (cid:105) to denote ensemble averages, and a toindicate time averages. Note that the MSD is an ensembleaverage, telling us information about the overall behaviourof the system. Obviously, the accuracy of MSD has a strongdependence on the number of particles used to build theensemble. Also, in a system with open boundaries, N p is notnecessarily constant, and reduces in time.In the case of anomalous diffusion the MSD is knownto scale with time as (cid:104) ∆ s ( t )(cid:105) = D α ( E ) t α , (2)where D α ( E ) is the anomalous diffusion coefficient, which,in general, depends on the particle energy. The anomalousdiffusion exponent α characterizes the diffusion regime, suchthat when α is equal to 1 Brownian diffusion is at play,whereas if α < or α > the system is subdiffusive orsuperdiffusive, respectively. It is possible, given an ensembleof particles for which the diffusion regime is not known, toinfer a diffusion exponent, by calculating (cid:104) ∆ s ( t )(cid:105) for the ensemble of particles, and then α can be estimated from alinear fit in log-log space.Often, datasets of trajectories are not ideal for calculat-ing the MSD, since difficulties can arise when dealing withpoor particle statistics, leading to noisy ensemble averages.Furthermore, if the system is heterogeneous (e.g., differentparticle species present, non-uniform background environ-ment), the MSD cannot account for this and simply givesan ensemble-averaged perspective on the overall system dif-fusion regime.The time-averaged mean square displacement(TAMSD) method tries to address such problems byusing single-particle trajectories and it is commonly used inmany areas of biophysics (e.g. Golding & Cox 2004; Wonget al. 2004; Saxton 2012). In TAMSD, the ensemble averageis replaced by an averaging time window: ∆ s ( τ ) = T f − τ T f − τ (cid:213) m = [ x ( m + τ ) − x ( m )] , (3)where τ is the averaging window, T f is the final time forwhich the particle is tracked and x is the particle position in MNRAS000 , 1–17 (2019) article transport in shock simulations time. All times are considered to be discretised with uniformsampling. This leads to an estimator of the anomalous expo-nent based on single particle trajectories if ∆ s scales withtime in the same way described by Equation 2. In the idealcase, for a dataset with infinite time records and for Brow-nian motion, TAMSD and MSD would produce the sameresult because of the central limit theorem.The general procedure to apply TAMSD to a sample oftrajectories is to fit an anomalous exponent α to each parti-cle trajectory (as for MSD), and then take an average of theobtained exponents (cid:104) α (cid:105) . This is known to have some prob-lems that often lead to an underestimation of (cid:104) α (cid:105) , derivingfrom the fact that new free parameters are introduced inthe diagnostic (see, for example, Kepten et al. 2015; Bur-necki et al. 2015). For example, a time interval needs to bechosen for the fitting of the exponents, and this is nontrivialparticularly when long-memory processes are present in thesystem. However, a set of particles interacting with a shockis intrinsically heterogeneous, and we suggest TAMSD as atool to identify different ensembles of particles with differentdiffusion properties.To have direct comparison between MSD and TAMSD,it is possible to consider an ensemble average of the sam-ple’s TAMSD (Kepten et al. 2013). We employ this diag-nostic, called ensemble-averaged time-averaged mean squaredisplacement (EA TAMSD), to the particle trajectories ob-tained in the shock simulation. The EA TAMSD is definedby (cid:104) ∆ s ( τ )(cid:105) = N p N p (cid:213) p = T f − τ T f − τ (cid:213) m = [ x p ( m + τ ) − x p ( m )] , (4)that is the average of the TAMSD in Equation 3 over theensemble of tracked particles. As we will show in the resultssection, due to its double-averaging feature, the EA TAMSDis 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. Compar-ing MSD and EA TAMSD is also important because pre-vious studies have shown that MSD and EA TAMSD maydiffer when there is an ergodicity breaking in the sample(Metzler et al. 2014).Even though this work is concerned with computing dif-fusion exponents from particle trajectories obtained from hy-brid PIC shock simulations, MSD, TAMSD and EA TAMSDdiagnostics have also been tested on a Monte Carlo model,to ensure the robustness of the inferences made for the PICsimulations. Two examples of such tests are given in Ap-pendix A. Results regarding the diffusion regimes for particles up-stream of the shock are presented in this section. All thecalculations shown here are done in the upstream rest frame.The particle sample considered consists only of particles ini-tially reflected at the shock, as the trajectory fragments ofparticles before the first interaction with the shock are ne-glected. Throughout this analysis, particles are treated re-gardless of their energies even though the particle diffusion coefficient is known to be energy dependent, and this mayaffect the estimation of transport properties. However, thischoice is motivated by the fact that the discussion below isfocused on investigating the particle transport regime usingdifferent diagnostics (MSD, TAMSD, EA TAMSD) ratherthan to provide a discussion about particles’ diffusion co-efficients. However, an analysis of particle MSD for differ-ent energies is reported in Section 3.4. It is worth notingthat the upstream sample is the poorest from the statisti-cal point of view, as the particles interacting with the shockhave a relatively high probability to be transmitted down-stream. On the other hand, the particles in the upstreamregion are likely to be efficiently accelerated, as they cangain energy through the multiple-shock encounter scenariodiscussed above and seen in the trajectories of Figure 3.MSD and EA TAMSD are shown in Figure 4, whichgives an overall view of the ensemble behaviour upstreamof the shock. Firstly, we observe that the overall behaviourof the two diagnostics is remarkably similar, validating theuse of EA TAMSD. Note that the time τ shown when com-paring MSD and EA TAMSD in Figure 4 has two differentmeanings: for the MSD it is the proper simulation time (seeEquation 1) and for the EA TAMSD it is the length of theaveraging time window, defined in Equation 3.At extremely short times (less than one Ω − ci ), a ballis-tic regime, with the diffusion exponent α =2 is identifiedin the MSD. This is expected since at short times scatter-ing is still not taking place and, in principle, both diffusionand superdiffusion should be studied over sufficiently longtimes (i.e., for times longer than a few interaction times).Since the ballistic regime lasts for a very short time, theEA TAMSD does not capture well this ballistic transient,and this is related to the time average operation included inTAMSD.At moderate times (from 6 to 45 Ω − ci ), there is a slightchange of behaviour, and the MSD show a quasi-ballistictrend. A fit to the displacement in time gives values for theanomalous diffusion exponent of 1.86 and 1.93 when calcu-lated from MSD and EA TAMSD, respectively. These valuesare both consistent with superdiffusion. At later times (from50 to 450 Ω − ci ), another behaviour can be seen where MSDand EA TAMSD have exactly the same slope, and enteranother superdiffusive regime, with an anomalous diffusionexponent α of 1.4. The reason for this change of slope maybe that the most superdiffusive particles (i.e., the ones un-dergoing a very small number of scattering episodes) leavethe box on short timescales.It is important to underline that, for the upstream sam-ple, EA TAMSD and MSD have a very similar behaviour andare almost indistinguishable for more than a decade. On onehand, this is a good test for the TAMSD diagnostic, whichis a novel technique for particle simulations in collisionlessplasma. On the other hand, this result is a validation forwhat was found using MSD.TAMSD has the advantage that it gives an insight intothe diffusive behaviour of single particles, in particular, dif-fusion exponents for each single particle trajectory. Figure 5shows an example of how single particle TAMSD can beused for the upstream particles. A subset of around 4000trajectory fragments of particles surviving for at least 200 Ω − ci in the upstream region has been selected. Single par-ticle TAMSDs have been calculated using Equation 3 and MNRAS , 1–17 (2019)
Trotta et al. -1 N p -1 τ [ Ω − ci ] ∆ s MSDEA TAMSD α = 2.0 α = 1.9 α = 1.4 Figure 4.
Top: Number of particles in the upstream region as a function of time. Bottom: MSD and EA TAMSD for upstream sample(red and blue lines, respectively). The black lines show the diffusion exponents α obtained by fitting the MSD in the indicated ranges. then single particle diffusion exponents have been fitted toobtain a distribution of single particle diffusion exponents α . The histogram in Figure 5 (top) shows the values of α obtained with TAMSD on these trajectory fragments. Thefit excludes the last few points of the single particle TAMSDin order to reduce the statistical noise present when usinglarge averaging windows (see TAMSD fluctuations at large τ in Figures 5 & 7) and to obtain a robust estimation of thevalues of α .The fluctuations that TAMSD shows at long time scalesare due to poor statistics in the trajectory fragments. FromEquation 3, for larger times fewer points contribute to thetime averaging window, so if a particle returns close to itsinitial position the value of ∆ s drops drastically, inducingthe fluctuation seen in the single particle TAMSD.The mean value for the single particle α distributionshown in Figure 5 is approximately 1.75, i.e., between thetransient with very strong superdiffusion (almost ballistic) and the long term superdiffusive behaviour of 1.4 found us-ing MSD and EA TAMSD and shown in Figure 4. This valuefor (cid:104) α (cid:105) is still consistent with superdiffusion. We note thata transition between different anomalous transport regimesis also found in laboratory plasmas (Bovet et al. 2015). It isalso worth noting that every single particle α value is con-sistent with superdiffusion. In this section we present results obtained for the regiondownstream of the shock. All the calculations are done inthe downstream rest frame. The downstream region is verydifferent from that upstream since the plasma has been ther-malised by the shock and the level of fluctuations is muchhigher. From the point of view of particle transport, par-ticles with very different histories are found downstream.Some particles are transmitted and moving away from theshock, others are returning upstream (gaining energy in the
MNRAS000
MNRAS000 , 1–17 (2019) article transport in shock simulations Figure 5.
Top: Histogram of single particle diffusion exponents obtained from TAMSD for the upstream region. Bottom: Examples of asubset of typical upstream particle trajectories (left) and their single particle TAMSDs (right). Lines corresponding to α = and α = are also plotted. process), and some others stay trapped in the strong mag-netic field fluctuations present immediately after the shocktransition. Some traces of this intrinsic heterogeneity are ex-pected to be found when addressing the diffusion regime.With the same motivation discussed for the upstreamregion, particles are treated regardless of their energy. A dis-cussion about the downstream MSD behaviour at differentenergies is reported in Section 3.4.An ensemble of downstream particle trajectory frag-ments is built in the same fashion as described earlier. Alower boundary is set just after the shock transition zone(80 d i downstream of the nominal shock transition). An-other boundary is placed close to the reflecting wall (20 d i infront) to remove any artificial scattering that the wall mightproduce, and when a particle crosses this “wall boundary” itis discarded. Every particle that crosses the left boundary toenter the downstream region is treated as a new particle tra-jectory fragment, and particles escaping the left boundaryare discarded and tracked as new particles in the transition zone. Figure 6 shows a comparison between MSD and EATAMSD for the downstream sample. The scenario obtainedis clearly different from the upstream case. At very shorttimes, the ballistic behaviour is recovered with MSD. Onceagain, the EA TAMSD does not capture well this regimebecause of the time average included in the diagnostic. Atintermediate times (from 3 to 30 Ω − ci ), a temporarily su-perdiffusive regime is observed once again. The behaviourat intermediate times is different from the one found in theupstream particle sample. First of all, the diffusion expo-nents found with MSD and EA TAMSD, 1.65 and 1.85 re-spectively, are smaller than the typical value of 1.9 foundupstream. This result can be interpreted considering thatscattering for downstream particles starts to happen, statis-tically, at smaller time scales than the upstream case, andthis is due to the higher level of fluctuations downstreamof the shock. It is important to note that, for intermediatetimes, the diffusion exponents obtained with MSD and EATAMSD differ by 0.2, which may be due to the heterogene- MNRAS , 1–17 (2019) Trotta et al. -1 N p -1 τ [ Ω − ci ] ∆ s MSDEA TAMSD α = 2 α MSD =1.6574 α EA TAMSD =1.8538 α =1.0425 Figure 6.
Top: Number of particles in the downstream region as a function of time. Bottom: MSD and EA TAMSD for the downstreamparticle sample (red and blue lines, respectively). The black lines show the diffusion exponents α obtained by fitting the MSD in theindicated ranges. ity of the sample and the process of ensemble averaging, asdiscussed later.For longer times (30 to 150 Ω − ci ), both diagnostics havethe same slope, with an α very close to unity, indicating nor-mal diffusion. The recovery of normal diffusion yields to twoimportant considerations. On one hand, we find a diffusionregime that is compatible with the theoretical assumptionsof DSA, but on the other hand, the result must be inter-preted keeping in mind that the ensemble average involvedin both the diagnostics averages out the behaviour of singleparticles or smaller populations of particles that may be indifferent transport regimes.With single particle TAMSD, we can identify subsetsof particles having different behaviours. Before discussingsingle particle TAMSD, it is worth noting that, for timeslarger than 150 Ω − ci , ∆ s exhibits a large scale oscillation.This is an effect due to the increase of the finite size of thedownstream region, confirmed by performing different sim- ulations with different box sizes. Results for single particleTAMSD are shown in Figure 7. Again, a subset of particlesthat survive for at least 200 Ω − ci has been selected to per-form this diagnostic. The single particle diffusion exponentshave been fitted excluding the last few points of ∆ s , whereTAMSD become noisy (as for the upstream region). The his-togram representing the α values obtained for the particleensemble has a mean value of 1.27, representing superdiffu-sion. Although this value of (cid:104) α (cid:105) is smaller than that foundupstream (c.f. (cid:104) α (cid:105) = . ), it is still indicative of superdiffu-sion with an exponent consistent with those obtained fromdata analysis (e.g. Perri et al. 2015; Zimbardo & Perri 2018).Like the upstream case, the (cid:104) α (cid:105) obtained with TAMSD is inbetween the two regimes found with MSD and EA TAMSD(see Figure 6).When confronting the results obtained with different di-agnostics, it is useful to note that EA TAMSD (defined in MNRAS000
Top: Number of particles in the downstream region as a function of time. Bottom: MSD and EA TAMSD for the downstreamparticle sample (red and blue lines, respectively). The black lines show the diffusion exponents α obtained by fitting the MSD in theindicated ranges. ity of the sample and the process of ensemble averaging, asdiscussed later.For longer times (30 to 150 Ω − ci ), both diagnostics havethe same slope, with an α very close to unity, indicating nor-mal diffusion. The recovery of normal diffusion yields to twoimportant considerations. On one hand, we find a diffusionregime that is compatible with the theoretical assumptionsof DSA, but on the other hand, the result must be inter-preted keeping in mind that the ensemble average involvedin both the diagnostics averages out the behaviour of singleparticles or smaller populations of particles that may be indifferent transport regimes.With single particle TAMSD, we can identify subsetsof particles having different behaviours. Before discussingsingle particle TAMSD, it is worth noting that, for timeslarger than 150 Ω − ci , ∆ s exhibits a large scale oscillation.This is an effect due to the increase of the finite size of thedownstream region, confirmed by performing different sim- ulations with different box sizes. Results for single particleTAMSD are shown in Figure 7. Again, a subset of particlesthat survive for at least 200 Ω − ci has been selected to per-form this diagnostic. The single particle diffusion exponentshave been fitted excluding the last few points of ∆ s , whereTAMSD become noisy (as for the upstream region). The his-togram representing the α values obtained for the particleensemble has a mean value of 1.27, representing superdiffu-sion. Although this value of (cid:104) α (cid:105) is smaller than that foundupstream (c.f. (cid:104) α (cid:105) = . ), it is still indicative of superdiffu-sion with an exponent consistent with those obtained fromdata analysis (e.g. Perri et al. 2015; Zimbardo & Perri 2018).Like the upstream case, the (cid:104) α (cid:105) obtained with TAMSD is inbetween the two regimes found with MSD and EA TAMSD(see Figure 6).When confronting the results obtained with different di-agnostics, it is useful to note that EA TAMSD (defined in MNRAS000 , 1–17 (2019) article transport in shock simulations Figure 7.
Top: Histogram of single particle diffusion exponents obtained from TAMSD for the downstream region. Bottom: Examplesof a subset of typical downstream particle trajectories (left) and single particle TAMSDs (right).
Equation 4) leads to an estimation of α over the ensemble,that may depend on time, and is robust against the noiseof MSD. TAMSD, on the other hand, leads to an estima-tion of single particle α values without α varying in time, soit introduces the possibility of examining the heterogeneityof the ensemble. To this extent, when analysing the down-stream particle transport properties, the key result is notthe mean value of the TAMSD histogram, but rather itslarge spread over values covering the whole domain of possi-ble anomalous transport exponents, that demonstrates thatparticles with completely different diffusion regimes coexistin the downstream region.It can be seen in the bottom-right plot of Figure 7that single particle TAMSD have a remarkable spread, muchlarger than those observed in the upstream samples, rangingfrom subdiffusive regimes ( α < ) to superdiffusive ( α > ).We suggest that subdiffusive diffusion exponents representparticles trapped in the strong magnetic fluctuations presentdownstream of the shock (and in particular immediately after the shock transition). Superdiffusive single particleanomalous exponents can, on the other hand, correspondto particles transmitted downstream which move away fromthe shock without undergoing much scattering. To demonstrate our interpretation of the many possibletransport realizations seen in the downstream region, weshow in Figure 8 three examples of particles in differenttransport regimes. On the left, a particle trapped down-stream is shown. The α value obtained with single parti-cle TAMSD is 0.7, indicating a subdiffusive behaviour. Thesame analysis is carried out in the middle plot, and gives anormal diffusive value of α , confirmed from the particle tra-jectory observed, that undergoes many scattering episodesbut still samples more space than for the particle shownin the left panel. Finally, a particle with superdiffusive α value is shown (right panel), and it is an example from the MNRAS , 1–17 (2019) Trotta et al. x [ d i ] t [ Ω − c i ] Subdiffusive τ [ Ω − ci ] ∆ s α =0.76691 x [ d i ] t [ Ω − c i ] Diffusive τ [ Ω − ci ] ∆ s α =1.0085 x [ d i ] t [ Ω − c i ] Superdiffusive τ [ Ω − ci ] ∆ s α =1.671 Figure 8.
Top: Three particle trajectory examples. Bottom: Single particles TAMSDs (blue lines) and their fits (dashed magenta lines). population of particles moving away from the shock with-out much scattering. The diffusion properties of particle en-sembles are, however, statistical properties, i.e., the singletrajectories may behave very differently from the averageover many particles. The examples presented in Figure 8demonstrates that TAMSD can be very powerful diagnosticwhen dealing with poor particle statistics and/or heteroge-nous particle samples.
As discussed previously, the analyses presented so far treatedupstream and downstream particles regardless of their en-ergy, although the diffusion coefficient depends from the par-ticles’ energy (see Equation 2). Figure 9 shows upstream (a)and downstream (b) MSD calculations for particles in threedifferent energy ranges. It should be noted that, since theshock transition zone is neglected in this study, and all thecalculations are done in the local rest frame, the particles’ energy remains roughly constant throughout the trajectoryfragments analysed.The dependence of the diffusion coefficient from the par-ticles’ energy in the upstream region is clearly found in Fig-ure 9 (left panels), where the (cid:104) ∆ s (cid:105) is found to be systemati-cally larger for particles with larger energies. The black line,showing the upstream MSD for the whole particle ensem-ble, represents an ‘average’ of what obtained in the differentenergy ranges.Figure 9 shows that the MSD diagnostic for upstreamparticles is consistent with a superdiffusive scenario for allthe energy ranges considered. It is worth noting that this re-sult was predictable using single particle TAMSD upstream(demonstrating one of the strengths of TAMSD). In par-ticular, the histogram of single particle anomalous diffusionexponents reported in Figure 5 shows that all the upstreamparticles have single particle diffusion exponent greater thanone, indicating superdiffusion. MNRAS000
As discussed previously, the analyses presented so far treatedupstream and downstream particles regardless of their en-ergy, although the diffusion coefficient depends from the par-ticles’ energy (see Equation 2). Figure 9 shows upstream (a)and downstream (b) MSD calculations for particles in threedifferent energy ranges. It should be noted that, since theshock transition zone is neglected in this study, and all thecalculations are done in the local rest frame, the particles’ energy remains roughly constant throughout the trajectoryfragments analysed.The dependence of the diffusion coefficient from the par-ticles’ energy in the upstream region is clearly found in Fig-ure 9 (left panels), where the (cid:104) ∆ s (cid:105) is found to be systemati-cally larger for particles with larger energies. The black line,showing the upstream MSD for the whole particle ensem-ble, represents an ‘average’ of what obtained in the differentenergy ranges.Figure 9 shows that the MSD diagnostic for upstreamparticles is consistent with a superdiffusive scenario for allthe energy ranges considered. It is worth noting that this re-sult was predictable using single particle TAMSD upstream(demonstrating one of the strengths of TAMSD). In par-ticular, the histogram of single particle anomalous diffusionexponents reported in Figure 5 shows that all the upstreamparticles have single particle diffusion exponent greater thanone, indicating superdiffusion. MNRAS000 , 1–17 (2019) article transport in shock simulations Upstream Downstream (a) (b) h s ih s i Figure 9. (a) , Bottom panel : Upstream MSDs for increasing particles’ energy (red, green and blue lines) and for the whole particleensemble (black line). The dashed lines show the α exponents calculated in Section 3.1 and shown in Figure 4. (a) , Top panel : Numberof particles in the upstream region as a function of time for different energies (same colors as bottom panel). (b) : same as (a) , but forthe downstream region. The α exponents shown here (black dashed lines) are the ones obtained in Section 3.2 and presented in Figure 6. The transition from α = . to α = . is seen in allenergy bands, although there is a trend to see this later atlower energies. The fact that the time at which the transitionto α = . happens changes is also influenced by a finite-box size selection effect. At later times, the results for thelowest energy particles becomes less significant because thenumber of particles in the ensemble decreases, as particlesreencounter the shock.In the shock downstream region, the dependence of thediffusion coefficient on the particle energy is also clear fromFigure 9 (right hand panels). Once again, when the en-ergy increases, the square displacement of particles becomeslarger, and the MSD calculated for the entire particle ensem-ble represents an ‘average’ between the behaviours observedat different energies.When the energy filter is applied to the downstreamparticle ensemble, the regime corresponding to α ∼ at latertimes (discussed in Section 3.2), indicating normal diffusion,is not recovered. We suggest that this may be an effect re-lated to the fact that the downstream particle ensemble isintrinsically heterogeneous.Like for the upstream region, this analysis confirmswhat found using other diagnostics, namely that the shockdownstream region has a much richer scenario as far as theparticle transport properties are concerned. Energetic particle transport regimes have been investigatedin the framework of 1D, quasi-parallel shock hybrid PICsimulations. The simulation domain was divided into up-stream and downstream regions, neglecting the transportregimes present within the shock transition. The transportregimes obtained for the two regions were found to have im-portant differences. To the best of our knowledge, this is thefirst time that diffusion exponents have been extracted us-ing particle trajectories advanced in the self-consistent fieldsproduced with hybrid plasma simulations.The usual technique to calculate anomalous diffusionexponents, based on the method of mean square displace-ment, relies on ensemble averages, and so it is very sensitiveto poor particle statistics. It is often hard to track particlesacross the domain for long enough times in the frameworkof hybrid PIC simulations. Also, the dynamics of differentparticle populations are often very different in the shocksimulation, and the ensemble average does not always cap-ture this feature very well. For this reason, we have used anovel application of time-average mean square displacementtogether with the usual MSD method. TAMSD is a commondiagnostic in scientific areas that have to deal with poor par-ticle statistics, such as biophysics. It tries to overcome theproblems of MSD by introducing a time averaging process
MNRAS , 1–17 (2019) Trotta et al. on the single particle trajectories, and enables single particleestimates of the diffusion exponent.In our simulation, we have found strong evidence for up-stream superdiffusion. This evidence primarily arose fromthe calculation of MSD in the upstream frame for the re-flected particles, in which we found extremely strong su-perdiffusion for intermediate times (with α ≈ α of 1.4 for long times. TheEA TAMSD, which is a novel technique for plasma simula-tions, is in extremely good agreement with the more classicMSD, although there are some differences between the twodiagnostics at intermediate times (6 to 45 Ω − ci ), associatedwith transient behaviour. Single particle TAMSD enabledus to look at the statistical distribution of single particleanomalous diffusion exponents, in agreement with the otherdiagnostics. TAMSD is a powerful tool when we deal withpoor particle statistics and with trajectory fragments of vari-able time extension.Downstream of the shock, a much richer scenario ofparticle transport regimes is obtained. Both MSD and EATAMSD show evidence for superdiffusion at intermediatetimes (6 to 30 Ω − ci ), even though the anomalous diffusionexponents inferred using these two diagnostics are different(and equal to 1.65 and 1.85, respectively). A possible expla-nation for this difference in the diagnostics at intermediatetimes lies in the intrinsic heterogeneity of the particle sam-ple. For later times, a diffusion exponent closer to unity, sug-gesting a tendency to Brownian motion is recovered usingboth MSD and EA TAMSD. When TAMSD is performed onsingle particle trajectories for the downstream sample, thesingle particle diffusion exponents are found to have a meanof 1.27 (consistent with moderate superdiffusion). However,the spread of the histogram of single particle α values ismuch larger than the one found for the upstream sample,ranging from highly subdiffusive to highly superdiffusive be-haviours. We consider that this is the indication of differentparticle populations coexisting in the same region: from par-ticles trapped in the shock fluctuations to particles movingaway from the shock without interacting with the surround-ing environment, (see examples presented in Figure 8).We have presented evidence of anomalous diffusion ofenergetic ions in a single 1D simulation of a collisionlessquasi-parallel shock. Such simulations obviously have limi-tations due to finite domain size and simulation time, so theapplicability of our results to other shocks (in terms of shockand/or simulation parameters) remains an open question.Future work will utilise fully 2D and 3D simulations(Trotta & Burgess 2019), addressing diffusion in the mag-netic field perpendicular and parallel directions. Moreover,it is possible to cross-correlate the diffusive regimes of dif-ferent 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 fur-ther investigation. We also intend to examine the trans-port regimes present in the shock transition layer, that havebeen neglected throughout this work, but which are proba-bly crucial for understanding particle acceleration, especiallyat low energies. It is worth underlining that, throughout thiswork, we studied a shock with parameters compatible withthose observed in the heliosphere. Further investigations willbroaden the parameter space, in order to address whether these results are relevant for higher Mach number, astro-physical shocks. ACKNOWLEDGEMENTS
DT acknowledges support of a studentship funded by thePerren Fund of the University of London. This researchwas supported by the UK Science and Technology FacilitiesCouncil (STFC) grant ST/P000622/1. This research utilisedQueen Mary’s Apocrita HPC facility, supported by QMULResearch-IT, http://doi.org/10.5281/zenodo.438045.
REFERENCES
Amano T., Hoshino M., 2007, ApJ, 661, 190Amano T., Hoshino M., 2009, ApJ, 690, 244Axford W. I., Leer E., Skadron G., 1978, in Dergachev V. A.,Kocharov G. E., eds, Cosmophysics. pp 125–134Baring M., Ellison D., Jones F., 1995, Advances in Space Re-search, 15, 397Bell A. R., 1978, MNRAS, 182, 147Blandford R., Eichler D., 1987, Physics Reports, 154, 1Blandford R. D., Ostriker J. P., 1978, ApJ, 221, L29Bovet A., Fasoli A., Ricci P., Furno I., Gustafson K., 2015, Phys.Rev. E, 91, 041101Bryant D. A., Cline T. L., Desai U. D., McDonald F. B., 1962,Journal of Geophysical Research, 67, 4983Burgess D., Scholer M., 2015, Collisionless Shocks in Space Plas-mas. Cambridge University PressBurnecki K., Kepten E., Garini Y., Sikora G., Weron A., 2015,Scientific Reports, 5, 11306Bykov A. M., Ellison D. C., Osipov S. M., 2017, Phys. Rev. E,95, 033207Caprioli D., 2012, J. Cosmology Astropart. Phys., 7, 038Caprioli D., Spitkovsky A., 2014a, ApJ, 783, 91Caprioli D., Spitkovsky A., 2014b, The Astrophysical Journal,794, 47Decker R. B., Krimigis S. M., Roelof E. C., Hill M. E., ArmstrongT. P., Gloeckler G., Hamilton D. C., Lanzerotti L. J., 2005,Science, 309, 2020Decker R. B., Krimigis S. M., Roelof E. C., Hill M. E., ArmstrongT. P., Gloeckler G., Hamilton D. C., Lanzerotti L. J., 2008,Nature, 454, 67Drury L. O., 1983, Reports on Progress in Physics, 46, 973Ellison D. C., Moebius E., Paschmann G., 1990, ApJ, 352, 376Fisk L. A., 2015, Journal of Physics: Conference Series, 642,012009Gargat´e L., Spitkovsky A., 2012, ApJ, 744, 67Giacalone J., 2004, ApJ, 609, 452Giacalone J., 2012, ApJ, 761, 28Giacalone J., Ellison D. C., 2000, J. Geophys. Res., 105, 12541Giacalone J., Burgess D., Schwartz S. J., Ellison D. C., 1992,Geophysical Research Letters, 19, 433Golding I., Cox E. C., 2004, Proceedings of the National Academyof Sciences, 101, 11310H¨ofling F., Franosch T., 2013, Reports on Progress in Physics, 76,046602Jokipii J. R., 1987, ApJ, 313, 842Jones F. C., Ellison D. C., 1991, Space Sci. Rev., 58, 259Kang H., Jones T. W., 1997, ApJ, 476, 875Kang H., Ryu D., Jones T. W., 2017, ApJ, 840, 42Kepten E., Bronshtein I., Garini Y., 2013, Phys. Rev. E, 87,052713Kepten E., Weron A., Sikora G., Burnecki K., Garini Y., 2015,PLOS ONE, 10, 1 MNRAS000
Amano T., Hoshino M., 2007, ApJ, 661, 190Amano T., Hoshino M., 2009, ApJ, 690, 244Axford W. I., Leer E., Skadron G., 1978, in Dergachev V. A.,Kocharov G. E., eds, Cosmophysics. pp 125–134Baring M., Ellison D., Jones F., 1995, Advances in Space Re-search, 15, 397Bell A. R., 1978, MNRAS, 182, 147Blandford R., Eichler D., 1987, Physics Reports, 154, 1Blandford R. D., Ostriker J. P., 1978, ApJ, 221, L29Bovet A., Fasoli A., Ricci P., Furno I., Gustafson K., 2015, Phys.Rev. E, 91, 041101Bryant D. A., Cline T. L., Desai U. D., McDonald F. B., 1962,Journal of Geophysical Research, 67, 4983Burgess D., Scholer M., 2015, Collisionless Shocks in Space Plas-mas. Cambridge University PressBurnecki K., Kepten E., Garini Y., Sikora G., Weron A., 2015,Scientific Reports, 5, 11306Bykov A. M., Ellison D. C., Osipov S. M., 2017, Phys. Rev. E,95, 033207Caprioli D., 2012, J. Cosmology Astropart. Phys., 7, 038Caprioli D., Spitkovsky A., 2014a, ApJ, 783, 91Caprioli D., Spitkovsky A., 2014b, The Astrophysical Journal,794, 47Decker R. B., Krimigis S. M., Roelof E. C., Hill M. E., ArmstrongT. P., Gloeckler G., Hamilton D. C., Lanzerotti L. J., 2005,Science, 309, 2020Decker R. B., Krimigis S. M., Roelof E. C., Hill M. E., ArmstrongT. P., Gloeckler G., Hamilton D. C., Lanzerotti L. J., 2008,Nature, 454, 67Drury L. O., 1983, Reports on Progress in Physics, 46, 973Ellison D. C., Moebius E., Paschmann G., 1990, ApJ, 352, 376Fisk L. A., 2015, Journal of Physics: Conference Series, 642,012009Gargat´e L., Spitkovsky A., 2012, ApJ, 744, 67Giacalone J., 2004, ApJ, 609, 452Giacalone J., 2012, ApJ, 761, 28Giacalone J., Ellison D. C., 2000, J. Geophys. Res., 105, 12541Giacalone J., Burgess D., Schwartz S. J., Ellison D. C., 1992,Geophysical Research Letters, 19, 433Golding I., Cox E. C., 2004, Proceedings of the National Academyof Sciences, 101, 11310H¨ofling F., Franosch T., 2013, Reports on Progress in Physics, 76,046602Jokipii J. R., 1987, ApJ, 313, 842Jones F. C., Ellison D. C., 1991, Space Sci. Rev., 58, 259Kang H., Jones T. W., 1997, ApJ, 476, 875Kang H., Ryu D., Jones T. W., 2017, ApJ, 840, 42Kepten E., Bronshtein I., Garini Y., 2013, Phys. Rev. E, 87,052713Kepten E., Weron A., Sikora G., Burnecki K., Garini Y., 2015,PLOS ONE, 10, 1 MNRAS000 , 1–17 (2019) article transport in shock simulations Klafter J., Blumen A., Shlesinger M. F., 1987, Phys. Rev. A, 35,3081Krymskii G. F., 1977, Akademiia Nauk SSSR Doklady, 234, 1306Kucharek H., Scholder M., Matthews A. P., 2000, Nonlinear Pro-cesses in Geophysics, 7, 167Malkov M. A., 1997, ApJ, 491, 584Matthews A. P., 1994, Journal of Computational Physics, 112,102Meerschaert M. M., Scalas E., 2006, Physica A: Statistical Me-chanics and its Applications, 370, 114Metzler R., Klafter J., 2000, Physics Reports, 339, 1Metzler R., Jeon J.-H., Cherstvy A. G., Barkai E., 2014, Phys.Chem. Chem. Phys., 16, 24128van Nes P., Reinhard R., Sanderson T. R., Wenzel K.-P., ZwicklR. D., 1984, Journal of Geophysical Research: Space Physics,89, 2122Perri S., Zimbardo G., 2007, ApJ, 671, L177Perri S., Zimbardo G., 2008, Journal of Geophysical Research(Space Physics), 113, A03107Perri S., Zimbardo G., 2009a, Advances in Space Research, 44,465Perri S., Zimbardo G., 2009b, ApJ, 693, L118Perri S., Zimbardo G., 2012a, ApJ, 750, 87Perri S., Zimbardo G., 2012b, ApJ, 754, 8Perri S., Zimbardo G., 2015, ApJ, 815, 75Perri S., Zimbardo G., Effenberger F., Fichtner H., 2015, A&A,578, A2Prete G., Perri S., Zimbardo G., 2019, Advances in Space Re-searchQian H., Sheetz M., Elson E., 1991, Biophysical Journal, 60, 910Ramos-Fern´andez G., Mateos J. L., Miramontes O., Cocho G.,Larralde H., Ayala-Orozco B., 2004, Behavioral Ecology andSociobiology, 55, 223Reynolds S. P., 2008, Annual Review of Astronomy and Astro-physics, 46, 89Riquelme M. A., Spitkovsky A., 2011, ApJ, 733, 63Saxton M. J., 2012, Biophysical Journal, 103, 2411Scholer M., Terasawa T., 1990, Geophysical Research Letters, 17,119Scholer M., Ipavich F. M., Gloeckler G., Hovestadt D., 1982, Jour-nal of Geophysical Research: Space Physics, 88, 1977Scholer M., Kucharek H., Giacalone J., 2000, Journal of Geophys-ical Research: Space Physics, 105, 18285Servidio S., Haynes C. T., Matthaeus W. H., Burgess D., CarboneV., Veltri P., 2016, Physical Review Letters, 117, 095101Strauss R. D. T., Effenberger F., 2017, Space Sci. Rev., 212, 151Sugiyama T., 2011, Physics of Plasmas, 18, 022302Sugiyama T., Shiota D., 2011, ApJ, 731, L34Sundberg T., Haynes C. T., Burgess D., Mazelle C. X., 2016, TheAstrophysical Journal, 820, 21Trotta D., Burgess D., 2019, MNRAS, 482, 1154Trotta E. M., Zimbardo G., 2015, Journal of Plasma Physics, 81,325810108van Weeren R. J., R¨ottgering H. J. A., Br¨uggen M., Hoeft M.,2010, Science, 330, 347van Weeren R. J., et al., 2017, Nature Astronomy, 1, 0005Weron A., Burnecki K., Akin E. J., Sol´e L., Balcerek M., TamkunM. M., Krapf D., 2017, Scientific Reports, 7, 5404Wolff M., Tautz R. C., 2015, A&A, 580, A58Wong I. Y., Gardel M. L., Reichman D. R., Weeks E. R., ValentineM. T., Bausch A. R., Weitz D. A., 2004, Phys. Rev. Lett., 92,178101Zimbardo G., Perri S., 2013, ApJ, 778, 35Zimbardo G., Perri S., 2017, Nature Astronomy, 1, 0163Zimbardo G., Perri S., 2018, MNRAS, 478, 4922
APPENDIX A: TESTING TAMSD WITHMONTE CARLO MODELS
Although TAMSD is a common diagnostic used to addressanomalous transport properties of many biological systems,it has never been used in the framework of numerical sim-ulations of collisionless plasmas. On the other hand, self-consistent trajectories from hybrid PIC shock simulationshave rarely been used to address the diffusion regimes atplay in the plasma environment. Therefore, we have testedTAMSD on particle trajectories obtained using a simpler,Monte Carlo model for anomalous transport. In the MonteCarlo model, particles are advanced integrating a Langevintype equation in the form of a stochastic differential equa-tion (e.g. Strauss & Effenberger 2017): dx i = v ran dt i . (A1)The distribution of v ran controls the diffusion regime, thatis, scattering episodes for the particles are modelled throughthe distribution of v ran . The random velocity is changed ran-domly at each scattering time t i . In the case of superdiffu-sion, the scattering times t i have a power-law-tailed prob-ability distribution. The Monte Carlo model employed togenerate the trajectories used for testing the anomalous dif-fusion diagnostics has extensively been used to bridge obser-vations and models of anomalous transport in plasmas bothfor collisionless shocks and plasma turbulence (Prete et al.2019).In this model, it is possible to set the exponent µ of thepower-law distribution of scattering times. The µ exponentis related to the anomalous diffusion exponent through therelation (e.g. Klafter et al. 1987): α = − µ. (A2)The particle scattering times are then defined through thedistribution ψ ( t ) = (cid:40) A t ≤ t A (cid:16) tt (cid:17) − µ t > t , (A3)where A is a normalisation constant and t is a free parame-ter. The details about this model can be found in Zimbardo& Perri (2013); Trotta & Zimbardo (2015).It is possible to compare the theoretical values for α with those obtained from the particle trajectories analysisusing MSD, TAMSD and EA TAMSD. It is worth notingthat for the self-consistent plasma model, the theoretical α is inferred from the analysis of particle trajectories (i.e., itis not a free parameter). MNRAS , 1–17 (2019) Trotta et al.
Figure A1.
Top: Single particle diffusion exponents histogram, calculated using TAMSD. Bottom Left: A sample of single particleTAMSD curves and their comparison with the theoretical diffusion exponent α theory = . Bottom right: MSD (red line) and EA TAMSD(blue line) curves for the particle sample. The dashed magenta line represents the α exponent obtained fitting MSD and EA TAMSD. The first test uses a diffusion exponent α = , corre-sponding to normal (Brownian) diffusion. The equations ofmotion (Equation A1) are integrated for an ensemble of 10 particles, so that the ensemble averages do not suffer fromstatistical noise. The final time for the particle trajectoriesis of 36000 time units with a timestep of unit time. Resultsare presented in Figure A1. Results for MSD, TAMSD andEA TAMSD are consistent and capture the expected theo-retical α used in the simulation. The single particle TAMSDexponents have been calculated by fitting the single particleTAMSD curves (Figure A1, bottom-left) in the range oftimes before they start to become noisy. The slight underes-timation of α using TAMSD ( α TAMSD ∼ . vs α theory = . )is an effect from which the TAMSD diagnostic is known tosuffer. We note that MSD and EA TAMSD capture the valueof the expected α with a smaller error than single particleTAMSD and, moreover, they return the same value when α is fitted ( α MSD = α EA TAMSD ∼ . ). The extremely strongagreement between MSD and EA TAMSD also has a theo-retical explanation: when normal diffusion is achieved, thecentral limit theorem is satisfied, so time and ensemble av-erages must converge for asymptotic times. For comparison, a superdiffusive case has been anal-ysed as well by means of the Monte Carlo method. Here,an anomalous diffusion exponent α of 1.5 (i.e., µ =2.5) hasbeen imposed in Equation A3. The number of trajectoriesfollowed and the final simulation time are the same as theprevious case. Figure A2 presents the results of the test.The single particle diffusion exponents have been fitted fromsingle particle TAMSDs with the criterion described above.The mean single particle diffusion exponent (cid:104) α TAMSD (cid:105) valueis 1.43 and slightly underestimates the theoretical value of1.5. Again, this is a known effect for the TAMSD diagnostic.However, the inferred (cid:104) α TAMSD (cid:105) is in good agreement withthe theoretical value, as can be also seen in the histogramshown in the top panel of Figure A2.The diffusion exponents obtained using MSD and EATAMSD are both consistent with the expected α theory withinreasonable errors, although their slopes have been fitted fora shorter time window. The slight discrepancy between theMSD and the EA TAMSD curves, as well as the ∼
1% differ-ence in the calculated α values may be due to the fact thatin the superdiffusive regime the central limit theorem doesnot hold anymore, and so we would expect some differencewhen taking time and ensemble averages. MNRAS000
1% differ-ence in the calculated α values may be due to the fact thatin the superdiffusive regime the central limit theorem doesnot hold anymore, and so we would expect some differencewhen taking time and ensemble averages. MNRAS000 , 1–17 (2019) article transport in shock simulations Figure A2.
Top: Single particle diffusion exponents histogram, calculated using TAMSD. Bottom Left: A sample of single particleTAMSD curves and their comparison with the theoretical diffusion exponent α theory = . . Bottom right: MSD (red line) and EA TAMSD(blue line) curves for the particle sample. The dashed lines are obtained by fitting the MSD (green) and the EA TAMSD (magenta) inthe time range indicated. The results of the tests are summarised in Table A1. Itis worth noting that the diagnostics used to infer the anoma-lous diffusion exponents are consistent, but show slight dif-ferences. However, TAMSD is a flexible and powerful toolto use when dealing with datasets affected by poor particlestatistics and trajectories with varying duration. The testspresented in this appendix help us understand the robust-ness of the diagnostics presented.
This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS , 1–17 (2019) Trotta et al.
Table A1.
Summary of anomalous diffusion exponents calculated with MSD, TAMSD and EA TAMSD for Monte Carlo simulationswith assigned α α theory α MSD α EA TAMSD (cid:104) α TAMSD (cid:105)000