Particle acceleration by sub-proton cyclotron frequency spectrum of dispersive Alfven waves in inhomogeneous solar coronal plasmas

The problem of explaining observed soft X-ray fluxes during solar flares, which invokes acceleration of large fraction of electrons, if the acceleration takes places at the solar coronal loop-top, can potentially be solved by postulating that flare at loop-top creates dispersive Alfven waves (DAWs) which propagate towards the foot-points. As DAWs move in progressively denser parts of the loop (due to gravitational stratification) the large fraction of electrons is no longer needed. Here we extend our previous results by considering $f ^{-1}$ frequency spectrum of DAWs and add ${\rm He^{++}}$ ions using fully kinetic particle-in-cell (PIC) simulations. We consider cases when transverse density gradient is in the range ${ 4-40} c/\omega_{\rm { pe}}$ and DAW driving frequency is $0.3-0.6\omega_{\rm { cp}}$. We find that (i) The frequency spectrum case does not affect electron acceleration fraction in the like-to-like cases, but few times larger percentage of ${\rm He^{++}}$ heating is seen due to ion cyclotron resonance; (ii) In cases when counter propagating DAWs collide multiple-times, much larger electron and ion acceleration fractions are found, but the process is intermittent in time. This is because intensive heating (temperature increase) makes the-above-thermal-fraction smaller; Also more isotropic velocity distributions are seen; (iii) Development of kink oscillations occurs when DAWs collide; (iv) Scaling of the magnetic fluctuations power spectrum steepening in the higher-density regions is seen, due to wave refraction. Our PIC runs produce much steeper slopes than the orginal spectrum, indicating that the electron-scale physics has a notable effect of DAW spectrum evolution.


INTRODUCTION
Alfven waves play an important role in solar atmosphere as a means of transporting energy from lower to higher layers (Ruderman et al. 1998).On large spatial scales, these waves may contribute to solving solar coronal heating problem in magnetohydrodynamic (MHD) approximation (Heyvaerts & Priest 1983;Boocock & Tsiklauri 2022a,b).On small spatial scales, i.e. in the kinetic approximation such waves may play a major role in resolving problems arising in quantifying observed energetic, super-thermal particles (Tsiklauri et al. 2005;Fletcher & Hudson 2008;McClements & Fletcher 2009).
Some of relevant works on kinetic Alfven waves include the following: Starting with books Wu (2012) and Wu, D. J. & Chen, L. (2020) discuss Kinetic Alfven wave aspects such as theory, experiments, observations and practical applications.A review of Earth magnetospheric aspects of kinetic Alfven waves is presented in Stasiewicz et al. (2000).Some recent works include: Bacchini et al. (2022) have studied kinetic-scale heating by Alfven Waves in magnetic shears.With first-principles kinetic simulations, they show that a large-scale Alfven wave (AW) propagating in an inhomogeneous background decays into kinetic Alfven waves (KAWs), triggering ion and electron energization.They demonstrated that the two species can access unequal amounts of the initial AW energy, experiencing differential heating.Maiorano et al. (2020) study Kinetic Alfven wave generation by velocity shear in collisionless plasmas.The evolution of ★ E-mail: d.tsiklauri@qmul.ac.uk a linearly polarized, long-wavelength Alfven wave propagating in a collisionless magnetized plasma with a sheared parallel-directed velocity flow has been studied by means of two-dimensional hybrid Vlasov-Maxwell (HVM) simulations.By analyzing both the polarization and group velocity of perturbations in the shear regions, they identify them as KAWs.In the moderate amplitude run, kinetic effects distort the proton distribution function in the shear region.This leads to the formation of a proton beam, at the Alfven speed and parallel to the magnetic field.Xiang et al. (2019) studied resonant mode conversion of AWs into KAWs when there is an arbitrary angle between background magnetic field and plasma density gradient in plasma with two different temperature species.Xiang et al. (2019) find that the efficiency of mode conversion depends on the angle, and spatial scales such as the density inhomogeneous gradient scale, inverse of parallel wavenumber, and electron/ion temperature ratio.Malara et al. (2019) studied electron heating by KAWs in solar coronal loop turbulence.A test-particle model describing the energization of electrons in a turbulent plasma was presented.A fluctuating electric field component parallel to the background magnetic field, with properties similar to those of KAWs, is assumed to be present at scales of the order of the proton Larmor radius.Electrons were found to be stochastically accelerated by multiple interactions with such fluctuations, reaching energies of the order of 10 2 eV within tens to hundreds of seconds, depending on the turbulence amplitude.Pezzi et al. (2017b) considered Turbulence generation during the head-on collision of Alfvenic wave packets in the context of Moffatt and Parker problem.Pezzi et al. (2017a) have extended the work in the same context by including Hall magnetohydrodynamics and two hybrid kinetic Vlasov-Maxwell numerical models.Pezzi et al. (2017c) considered the interaction of two colliding Alfven wave packets and found that the extension to include compressive and kinetic effects, while maintaining the gross characteristics of the simpler classic formulation, also reveals interesting features beyond MHD treatment.Valentini et al. (2017) studied a transition to kinetic turbulence at proton scales driven by large-amplitude kinetic Alfven fluctuations.The transition from quasi-linear to turbulent regimes was investigated, focusing in particular on the development of important non-Maxwellian features in the proton distribution function driven by KAW fluctuations.Pucci et al. (2016) studied a model where an initial Alfven wave propagates inside an equilibrium structure which is inhomogeneous in the direction perpendicular to the equilibrium magnetic field.In a previous paper this situation has been considered in a particular configuration where the initial wave vector is parallel to the magnetic field and the wave is polarized perpendicular to the inhomogeneity direction.In Pucci et al. (2016) they consider other configurations, with a different polarization and possible initial oblique propagation.Vasconez et al. (2015) considered kinetic Alfven wave generation by largescale phase-mixing and found that kinetic simulations show that KAWs modify the ion distribution function, generating temperature anisotropy of both parallel and perpendicular to the local magnetic field as well as particle beams aligned along the local magnetic field.Wu & Chao (2003, 2004) investigated the magnetic field-aligned acceleration of auroral super-thermal electrons by KAWs.In the solar physics context, Wu & Fang (1999) studied the nonuniform heating of magnetized plasmas in the solar atmosphere with KAWs using the drift kinetic equation.They investigated electron heating mechanism and proposed that it can in principle explain X-ray brightness distributions observed in solar coronal loops.Wu & Fang (2003) investigated dissipation of KAWs in a solar polar plume placed in so-called coronal hole.They found that electron heating produced by KAW dissipation can in principle balance radiative losses of the plume.Wu & Fang (2007) tackled the problem of explaining the observations that sunspots have a higher temperature than the ambient quiet Sun in the upper chromosphere, although they appear dark in the photosphere.Wu & Fang (2007) showed above 850 km, in the upper chromosphere, KAW dissipation is the main heating mechanism at such heights.Observations of the solar corona find that the heavy ion species experience anisotropic heating that is primarily across the magnetic field.Wu & Yang (2006, 2007) investigated the nonlinear interaction of heavy ions with KAWs.Using three-component plasma model with electrons, protons, and heavy ions.Results show that heavy ions are energized across the magnetic field.Wu et al. (1996Wu et al. ( , 1997) ) discuss a new type of density soliton that is discovered in data from the Freja satellite.A good case is made for it to be KAW based explanation.
Another interesting feature related to KAWs if found in Earth's outer radiation belt.The term time domain structures (TDSs) refers to packets of greater than 1-ms duration intense electric field spikes detected by Van Allen Probes in the Earth's outer radiation belt.Phase space holes, double layers and other solitary electric field structures, referred to as time domain structures (TDSs), often occur around dipolarization fronts in the Earth's inner magnetosphere.In An et al. (2021), authors demonstrate that TDSs can be excited by electrons in nonlinear Landau resonance with kinetic Alfven waves.Hughes et al. (2017) considered kinetic Alfven turbulence and electron and ion heating by PIC simulations.They find that in contrast to dissipation by whistler turbulence, the maximum ion heating rate due to kinetic Alfven turbulence is substantially greater than the maximum electron heating rate.It is also found that the Landau wave-particle resonance is a likely heating mechanism for the electrons and may also contribute to ion heating.
There are also studies that examine causal connection between KAWs and magnetic reconnection is 2D (Huang et al. 2018) and 3D (Liang et al. 2016).Specific effects related to KAW eigenmode in magnetosphere magnetic reconnection has been studied by Dai & Wang (2022).Dai et al. (2017) provided KAW explanation of the Hall fields in magnetic reconnection.
Coming back to solar flares, the observed soft X-ray flux during flares is produced by electron bremsstrahlung, when accelerated electrons that move from magnetic loop top to the foot-points are slowed down by dense layers of the sun.In order to explain the observed soft X-ray flux during solar flares, if electron acceleration happens at loop top, nearly 100% electrons need to be accelerated (Fleishman et al. 2020(Fleishman et al. , 2022)).No acceleration mechanism is known with such high efficiency.This problem can potentially be solved by postulating (Fletcher & Hudson 2008) that flare at loop top creates dispersive Alfven waves (DAWs) which then propagate towards the foot-points, and as they move in progressively denser parts of the loop (due to natural gravitational stratification) the aforementioned high percentage is no longer needed.It has been known that, in homogeneous plasma, when perpendicular wavelength of Alfven wave (AW) approaches kinetic scales such as e.g.ion-inertial length, it acquires magnetic-field-aligned (parallel) electric field, which can efficiently accelerate electrons (Stasiewicz et al. 2000).Further, Tsiklauri et al. (2005) have shown that if DAW propagates in plasma with transverse (with respect to external magnetic field) density inhomogeneity, the generated parallel electric field is orders of magnitude higher than (i) homogeneous plasma case and (ii) Dreicer electric field (one that triggers electron run-away acceleration).Subsequently Tsiklauri (2012) has revisited the problem with full 3D particle-in-cell approach.Ofman (2010) considered similar set up as in Tsiklauri et al. (2005) but instead of considering one DAW harmonic with 0.3 cp he considered  −1 AW spectrum and added He ++ ions and used Hybrid simulation model.Here  cp = / p is proton cyclotron frequency.Note that our approach uses PIC code so it can resolve electron-scale physics contrary to Ofman who used a Hybrid code, which can resolve only ion-scale physics.Now in the present work we essentially revisit Ofman's set up run it for two cases: 1. when transverse density gradient is ≈ 4/ pe (as in Tsiklauri et al. (2005); Tsiklauri (2012)), i.e. on "electron"-scale; 2. when the gradient is on proton scale circa ≈ 40/ pe ; In this paper novel numerical simulation results are presented.Including the scaling of the magnetic fluctuations power spectrum steepening in the higher density regions, and the heating channelled to these regions from the surrounding lower density plasma due to wave refraction.We also present runs where DAW collide multiple times in a situation similar to Howes & Nielson (2013); Daiffallah (2022).We stress that despite the similarities between works of Ofman (2010) and Tsiklauri et al. (2005); Tsiklauri (2012), including the current work, the former used physical parameters for solar wind, while the latter use physical parameters for solar coronal flaring loop tops.Therefore direct comparison between results and magnetic fluctuation spectra scaling laws is not appropriate.In this context, also we note that subsequent studies exist: (i) (Maneva et al. 2015) which performed 2.5D hybrid simulations to investigate the role of initially imposed broad-band wave spectra in a drifting and expanding solar wind plasma and (ii) (Ofman et al. 2017) which study the effects of inhomogeneous (across the magnetic field) background streaming focusing on the fast solar wind.They explore the effects of an initial relative, inhomogeneous ion drift on the perpendicular ion heating and cooling and consider the effects of solar wind expansion.
It should be noted the above mentioned electron fraction/number problem is not without a debate.Kontar et al. (2023) argue that for a solar flare on September 10, 2017, their observations show 100 times smaller fraction of accelerated electrons, i.e. just mere ≈ 10 −2 , compared to the same flare observation in a different observing wavelength by (Fleishman et al. 2020(Fleishman et al. , 2022)), who on contrary, claim that nearly 100% electrons need to be accelerated to match the observations.Kontar et al. (2023)'s explanation for the small fraction is a steady resupply of electrons to the acceleration site.Since the present paper is purely theoretical and is based on numerical simulations, we do not make a judgement which of the aforementioned observations is correct, while insisting that electron acceleration by DAWs is a viable mechanism to accelerate electrons from loop top to foot-points.Nor we aim to match or model particular observational data set, staying within bounds of a theoretical, "the first principles" approach.
To summarize the background work and emphasize the novelty of this work, compared to our previous results, (Tsiklauri et al. 2005;Tsiklauri 2011;Tsiklauri 2012)) is 4-fold: (i) instead of exciting DAWs at the left edge of the domain which mimics loop top we drive DAWs in the middle of the domain, which more realistically represents DAWs generated at the loop top; (ii) we consider wide spectrum of DAWs as specified in Table 1; (iii) we added ion-cyclotron resonant He ++ species; (iv) we considered cases when DAW collide multiple times.

THE MODEL
We use 2.5D version of EPOCH which is a multi-dimensional, fully electromagnetic, PIC code (Arber et al. 2015).In EPOCH code physical quantities are in SI units.In this paper we use the normalization for the graphical presentation of the results as follows: Time and distance are shown in  −1 pe and / pe .Electric and magnetic fields are shown in units of  0 =  pe   / and  0 =  pe   /, respectively.Here  pe = [ 0  2 /(   0 )] 0.5 is electron plasma frequency. 0 = 10 16 m −3 is electron number density in the lowest density regions of the simulation domain ( = 0 and  =  max ).This sets electron plasma frequency at  pe = 5.64 × 10 9 Hz radian.The size of simulation box is either as  = 5000 and  = 200 grid points when transverse density gradient is ≈ 4/ pe (electron scale); or  = 5000 and  = 2000 grid points when transverse density gradient is 40/ pe (ions scale).The grid size in all numerical simulations presented in this work is electron Debye radius   =  th,e / pe .Here  th,e = √︁   /  is an electron thermal velocity.The plasma beta in this study is fixed at  = 2(  ℎ, /) 2 ( pp / cp ) 2 =  0 (0, 0)  /( 2 0 /(2 0 )) = 0.02.Thus  = 0.006746 <   /  = 1/16 = 0.0625, which means that we are so-called Inertial Alfven wave (IAW) regime.According to Stasiewicz et al. (2000) nomenclature is such that in the different regimes inertial ( <   /  ) and kinetic ( >   /  ).The latter waves are called Kinetic Alfven waves (KAW).We could computationally afford somewhat unrealistic electron to proton mass ratio of   /  = 1/16, because our largest run on  = 5000 and  = 2000 grid points with end-simulation-time of  end = 375.0/cp takes 1 day, 1 hour, 18 minutes on 2048 CPU cores (16 x 128 core nodes) using Dell PowerEdge R6525 compute nodes each with 2 x AMD EPYC 7742 (Rome) 2.25 GHz 64-core processors i.e. 128 cores per node.We admit that indeed the mass ratio   /  = 1/16 is an unrealistic one compared to the actual   /  = 1/1836 for the case of a proton.The mass ratio of   /  = 1/16 sets the ratio of Alfven speed to speed of light of   / =  cp / pp = 0.25, which is larger than the realistic one of   / = 0.023.However, many studies use this simplification to speed up the numerical computation.Our previous study Tsiklauri (2007) has considered the effect of variation of the mass ratio and found that amplitude attained by the generated  ∥ decreases linearly as the inverse of the mass ratio   /  , i.e.  ∝ 1/  .This result means that for a realistic mass ratio of   /  = 1836 our empirical scaling law produces  ∥ = 14 V/m for solar coronal parameters, which is a significant value for electron acceleration.In the present simulation we use 100 electrons, 100 protons and 100 He ++ ions per cell which means that for there are total of 3 × 10 8 particles in the simulation i.e. 1 × 10 8 per species for the 5000 × 200 spatial grid case.When we use 5000 × 2000 grid there are obviously 10 more particles i.e. total of 3 × 10 9 particles and 1 × 10 9 particles per species.
We impose constant background magnetic field of strength  0 = 0.032075 Tesla along -axis.This corresponds to  0 = 1.0( pe   /).This means that with the normalization used for the visualization purposes, normalized background magnetic field is unity.Also such choice sets the ratio of electron cyclotron and plasma frequencies as unity  ce / pe = 1.00.Electron, proton (ion) and He ++ temperature at the simulation box edge i.e. in the homogeneous density parts of the domain are fixed at   (0, 0) =   (0, 0) =  He ++ (0, 0) ≡  0 = 2 × 10 7 K.This temperature along with   (0, 0) ≡  0 = 10 16 m −3 , sets plasma parameters of a typical dense solar coronal flaring loops.In order to keep charge neutrality the relation between charges species is 1.0 0 = 0.9 0 +0.05He ++ 0 , because He is twice ionized.Here subscript 0 refers to location (0, 0) i.e. edge of the domain here density is uniform.
We note that the results presented here are equally applicable to both magnetospheric auroral density cavities (Génot et al. 1999(Génot et al. , 2000;;Thompson & Lysak 1996;Daiffallah & Mottez 2017;Daiffallah 2022) and solar coronal loops (Tsiklauri et al. 2005;Tsiklauri 2007Tsiklauri , 2011;;Tsiklauri 2012;Tsiklauri 2016), as far as there is a density gradient present, be it negative (density depletion) in the case of the cavities or positive (i.e.density enhancement) in the case of solar coronal loops.Discussing magnetosphere analogy here sets an appropriate interdisciplinary context of the current study.We note that only in the Earth magnetosphere DAWs are currently observed in situ, as the Sun is too far away from DAWs to be detected remotely.
We model density variation of all plasma constituent particle species, across the background magnetic field, of a solar coronal loop as Eq.( 1) shows that in the middle part along y-coordinate the number density is increased by a factor of four compared to edges on either side.Also this density profile has two density gradients in the vicinity of  = 0.25 max and  = 0.75 max , both having a width of about 0.2 max = 4/  .The background temperature of positive ions and electrons is varied as This makes sure that the thermal pressure in all plasma particle species is constant across y-coordinate.Therefore, because  0 = , such initial conditions have the total pressure (the sum of thermal and magnetic pressures) being constant.A graphical visualization of the density and temperature profiles according to Equations 1 and 2 have been published before by the author and as a further reference can be found as figure 2 from Tsiklauri (2016).
When we consider a single harmonic, we trigger generation of dispersive Alfven waves by driving the middle part in x-coordinate i.e.  =  max /2, for all y-s, in the following manner Such driving generates left-hand-polarized DAWs which then travel from the middle part of the domain to the left ( = 0) and right ( =  max ) edges.We chose left-hand-polarization for ion cyclotron resonance to occur for He ++ species.Note that such DAW driving is more realistic to solar flare conditions with flare taking place at loop top, i.e. at  =  max /2 than compared to earlier cases (Tsiklauri et al. 2005;Tsiklauri 2007Tsiklauri , 2011;;Tsiklauri 2012), where the driving was taking place at the domain's left edge,  = 0, and the DAWs where appearing on the right edge,  =  max , due to periodic boundary conditions.In our present numerical runs the boundary conditions are periodic for every physical quantity be it particles of fields.In equations ( 3) and ( 4),   is the DAW driving frequency, fixed at   = 0.3 cp .With such driving frequency, there is no significant proton-cyclotron damping and also the generated DAW is Alfvenic in the frequency sense.The onset time of the driver,  0 is fixed at 3.0 −1 cp i.e. 48.000 −1 pe for the case of   /  = 16.This means that the driver onset time is about 3 2  −1 cp (= 9 −1 cp ).This means that after the driver onset time we constantly inject DAWs with amplitude of 0.05 0 , i.e. 5% of the background magnetic field.The initial amplitudes of the  ⊥ are chosen to yield circular polarization   =   = 0.05( 0 ).This corresponds to the relative DAW amplitude 5% of the background magnetic field strength.This value of the amplitude is similar to that considered by Ofman (2010) who select 3% of the background magnetic field.We would like to comment that both studies suffer from the same shortcoming of too large value of / 0 of 3 − 5%, which is otherwise appropriate for many MHD scale waves, as observed extensively in the corona (e.g.Liu & Ofman (2014)).The above MHD waves are 8-10 orders of magnitude longer than the DAW discussed in the present paper or the ones discussed in Ofman (2010).However, this magnetic perturbation ratio and corresponding magnetic energy flux is unrealistically large for kinetic scale DAWs be it on ion or electron scale.Evidence of the observed magnetic fluctuation spectra throughout the heliosphere shows strong power law decrease of the energy with frequency (e.g., Alexandrova et al. (2021), Lotz et al. (2023)).The steep scaling of the power spectrum with  was also demonstrated in the present model (Figure 3) (see below).Therefore, the assumption of / 0 of 3-5% for these high frequency kinetic scale waves or DAWs is unrealistic.The main reason both PIC (this study) and Hybrid (Ofman (2010)) consider few percent amplitude waves is because considering 8-10 order of magnitude smaller amplitude is not possible because any kinetic model suffers from so-called 'shot noise'.Kinetic models typically consider 100 or 100s of particles per cell which is many orders of magnitude smaller than the realistic number of particles per volume.The level of 'shot noise' scales as ∝ 1/ √ , where  is number of particles per cell.In order to reduce amplitude we need to reduce 'shot noise', but that is computationally impossible because it would require many billions of particles per cell.However we remark that in the solar wind the magnetic field is 6 orders of magnitudes weaker that the magnetic field in typical active region loops, and the relative amplitudes of DAWs is expected to be orders of magnitude larger in the solar wind compared to flaring loops.
When we consider a broadband 1/  spectrum, as in Ofman ( 2010), the following DAW driving is imposed at the middle part in xcoordinate,  =  max /2: Here  is the number modes which is varied as  = 128, 512.() are 128 or 512 random numbers between 0 and 1 which stay the same throughout the simulation and they mimic some randomness introduced in the sine and cosine phase of the order of 5%, as in Ofman (2010).See Table 1 for details.DAW driving frequency is in the range of 0.3 − 0.6 cp , i.e. in Eqs.5 and 6  min = 0.3 cp and  max = 0.6 cp .This means that for protons our numerical simulation runs always stay clear from proton cyclotron resonance because in our notation  cp = / p denotes proton cyclotron frequency.Note that for He ++ ions, which has 2 positive charges and has a mass of 4 nucleons (2 protons and 2 neutrons), cyclotron frequency is  He ++ = 2/(4  ) = 0.5/ p = 0.5 cp , which is near the middle (0.45 cp ) of DAW frequency range of 0.3−0.6cp .Therefore we intentionally have the effect of He ++ ion cyclotron resonance heating as in Ofman (2010), which, in turn, was motivated by observations of high perpendicular temperature observations of heavy ion species such as e.g.He ++ (Kasper et al. 2013) and O VI (Telloni, D. et al. 2007).In the near future, there will be significant discoveries in this field with the advent of the Solar Orbiter Heavy Ion Sensor (HIS), which is a time-of-flight ion mass spectrometer dedicated to measuring heavy ions in the solar wind with unprecedented precision (Livi, S. et al. 2023).The factor  in Eqs.5-6 is calculated using For  = 128  = 0.11006650 and  = 512  = 0.053259745.Such choice insures that the broadband spectrum of DAWs delivers the same power as one alters number of harmonics, which will be demonstrated below.
Run 1 from Table 1 uses just one harmonic as per Eqs.3-4 with the end simulation time of  end  cp = 75.This is essentially reproduction of a previous run considered in Tsiklauri (2011) for the reference purposes and its input is used in production of Fig. 4. Run 2 from Table 1 uses Eqs.3-4 with the sum of 512 harmonics, but each with amplitude 512 times smaller than Run 1.Here the end simulation time is  end  cp = 75.This is essentially reproduction of a previous run considered in Tsiklauri (2011) for reference purposes and its input is used in production of Fig. 4. The purpose of this run is to ascertain that using the sum of 512 single frequency harmonics with amplitudes 512 times smaller than Run 1 produces identical results.
Run 3 from Table 1 uses 128 harmonics as per Eqs.5-6 with  from Eq.7 with the end simulation time of  end  cp = 75.The purpose of this this run is to ascertain changing number of harmonic still delivers the same power compared to Run 4 (512 harmonics).
Run 4 from Table 1 uses 512 harmonics as per Eqs.5-6 with  from Eq.7 with the end simulation time of  end  cp = 75.This the regarded as "the main run" in a sense that: (i) it shows the differences brought by wide spectrum of DAWs; (ii) most detailed simulation results will be shown for this run.
Run 5 from Table 1 has the end simulation time of  end  cp = 750.Otherwise it is identical to the main run (Run 4).The purpose of this run is to consider situation when DAWs collide multiple times.In numerical runs with  end  cp = 75 DAWs travel from the domain middle part in x-coordinate i.e.  =  max /2 towards  = 0 and  =  max , nearly reaching the said boundaries.Because of periodic boundary conditions DAWs will re-enter from the domain opposite sides, as the time progresses.This means that for Run 5  end  cp = 750 with 10-times longer simulation time the DAWs will collide at least 9 times, near stopping short of the 10th collision at  end  cp = 750.With Run 5 we wish to study how multiple DAW frontal collisions will affect particle acceleration in a situation similar to Howes & Nielson (2013); Daiffallah (2022).
Run 6 from Table 1 has the length-scale of the density gradient in -direction increased 10 times to nearly "ion"-scales / ≃ 40/ pe .This is achieved by setting 10-times the number of grids in -direction, i.e.  = 5000 and  = 2000 grid points in each direction.Otherwise it is identical to the main run (Run 4).The purpose of this run is to study how particle acceleration is affected by up-scaling the transverse density gradient from electron scales / ≃ 0.2 max = 4/ pe to ion-scales / ≃ 0.2 max = 40/ pe .In Run 6 we set  end  cp = 375, the maximal possible with the computational resources available.This run took 1 day, 34 minutes on 2048 CPU cores (16 x 128 core nodes) using Dell PowerEdge R6525 compute nodes.
Runs 7,8,9 from Table 1 are homogeneous versions of Runs 4,5,6 respectively.Note the different end simulation times used in each run.The purpose of these numerical runs is to quantify the difference in properties of particle acceleration in the homogeneous and inhomogeneous cases respectively.

THE RESULTS
In Fig. 1 we show a snapshot of   −  0 ,   ,   ,   ,   ,  He ++ ,   ,   ,  He ++ , at  = 75/ cp .We gather from this plot that the perturbation of magnetic field along the background magnetic field is generated as shown in panel (a).It closely follows pattern of the generated background magnetic filed-aligned electric field   (not shown for brevity), i.e. the generated field-aligned   and   −  0 have wave-like structure confined to the density gradient across the magnetic field with amplitudes ≈ 2%.Similarly as in previous results (panel (b)) is not phase-mixed while   (panel (c)) shows a clear pattern of phase mixing when the Alfvenic wave (DAW to be precise) travels slower in the over-dense region 5/ pe <  < 15/ pe as Alfven speed is proportional to the inverse square root of the plasma density.In fact, panels (a),(b),(c) look very similar to panels (a),(b),(c) in Fig. 1 from Tsiklauri (2011), in that amplitudes are similar, and phase-mixing pattern looks the same.However, the notable difference is that in Fig. 1 now (i) instead of exciting DAWs at the left edge of the domain now we use more realistic DAWs generated at the solar coronal loop top, i.e.  =  max /2 and (ii) we consider wide spectrum of DAWs as specified in Table 1.From panels (d) and (e) we gather that electron and proton number densities show significant transverse to the background magnetic field perturbations.We see in panels (g) and (h) the electron and proton temperatures also show wake like perturbation confined near density gradient regions and these gradient regions are the locations of a noticeable temperature increase.We note again that panels (d),(e),(g),(h) look very similar to panels (a),(b),(c),(d) in Fig. 2 from Tsiklauri (2011), with the proviso of the middle of the domain driving and wide spectrum.Tsiklauri (2011) did not consider He ++ hence the panels (f) and (i) which show He ++ number density and temperature respectively, are new.He ++ number density shows (panel (f)) significant non-linear, cat-claw like perturbations with the maximum concentrated near the location DAW driving at  =  max /2.He ++ temperature (panel (i) shows significant increase near the location DAW driving.This is due to the fact that for He ++ ion cyclotron condition is met.
Fig. 2 has a purpose to quantify how particle acceleration main features, fractions and the budget, split between different plasma species, is modified by the introduction of  −1 frequency spectrum of DAWs with driving frequency in the range 0.3 − 0.6 cp .Thus, in Fig. 2 we show distribution function dynamics for electron   ,   ,  components (panels (a), (b), (c)); ion (proton)   ,   ,  components (panels (e), (f), (g)); He ++   ,   ,  components (panels (i), (j), (k)); In panels (a), (b), (c), (e), (f), (g), (i), (j), (k), the inner dashed curve show the distribution function (distribution by velocities) at  = 0, while outer solid curve at the final simulation time.Note that for protons and He ++ the inner dashed curves for  = 0 look almost as solid due to large number of data points (each distribution function is resolved with 100,000 points in all panels).By comparing panels (a), (b), (c), (e), (f), (g) with similar run e.g.panels (a),(b),(c),(d),(e),(f) in Fig. 3 from Tsiklauri (2011) we gather many similarities.In particular we see in panel (a) that electron beams are formed near the DAW phase velocity.A study of Tsiklauri (2012) has shown that the electron beam formation along the uniform background magnetic field is due to Landau resonance between electrons and DAW.In fact Tsiklauri (2012) showed if one varies ion(proton) mass, which prescribes DAW phase speed, then location of the beam peak in the distribution function changes accordingly, signalling occurrence of the Landau resonance.In panels (b) and (c) we see that distribution functions with respect to   and   at time zero and the final simulation time do overlap.Hence, we deduce from these panels that no sizable electron heating in the perpendicular to the magnetic field direction is seen.While, for protons we see in corresponding panels (f) and (g) that primarily perpendicular to the magnetic field broadening of the distribution occurs due to DAW energy injection/driving.No change in the proton parallel to the magnetic field distribution function (panel (e)) is seen, suggesting that on this time-scale proton parallel dynamics is not affected by DAW driving.In panel (d) Tsiklauri (2012), make use of the following quantities/indexes: where  ,, are electron, proton and H++ velocity distribution functions and <> brackets denote average over y-coordinate, because temperature and density vary across y-coordinate.  , is the width of each density gradients according to Eq.1. , is the full width in the -direction.The numerator of Eq.( 8) yields number of super-thermal electrons in the parallel to the magnetic field direction, divided by the area (2  , ×  , ) in the region where particle acceleration takes place.Note that in this study because we consider some runs that are 10 longer in time, than previously, and hence expect commensurately a larger amount of energy being injected into the physical system due to continuous input from DAW wave driving, temperature of each plasma species is expected to increase.Therefore, crucial difference from the previous definitions as in Tsiklauri (2011); Tsiklauri (2012),   ℎ, () and   ℎ,, () are now functions of time.Equations 8 and 9 provide fraction (i.e.percentage) of super-thermal plasma species parallel and perpendicular to the regular magnetic field.Thus, in panel (h) (of Figure 2) solid line with triangles is for  ⊥, (); dashed line with diamonds is for  ⊥, (); dash-dotted line with stars is for  ∥, ().
We gather from panel (h) that, within  end  cp = 75,  ⊥, () attains value of 0.63, i.e. 63% of He ++ become super-thermal.As seen in panels (j) and (k) this is mainly due to the ion-cyclotron resonance in the perpendicular to  0 direction.So it is He ++ heating, not acceleration as such.A similar behavior for  ⊥, () attaining values of 23% is seen for protons in panel (h) i.e. by looking at panels (f) and (g) we see perpendicular proton heating, but not as dramatic as for He ++ because cyclotron resonance condition is not met for protons.In panel (h) we note that  ∥, () attains even smaller value so fraction of super-thermal electrons is 17%.Compared to panel (h) from Fig. 3 c)); ion (proton)   ,   ,  components (panels (e), (f), (g)); He ++   ,   ,  components (panels (i), (j), (k)); In panels (a), (b), (c), (e), (f), (g), (i), (j), (k), the inner dashed curve show the distribution function (distribution by velocities) at  = 0, while outer solid curve at the final simulation time.In panel (d) we show max( |E x /E x0 | ) (solid line with triangles) and < |E x (x, y, t)/E x0 | > (dashed line with diamonds).In panel (h) solid line with triangles is for  ⊥, ( ); dashed line with diamonds is for  ⊥, ( ); dash-dotted line with stars is for  ∥, ( ).In panel (l) we show He ++ , proton and electron temperature evolution in time.Solid line with triangles is for < T He ++ >; dashed line with diamonds is for < T i >; dash-dotted line with stars is for < T e >.See averaging definitions in the text.All panels show data for the Run 4. for (, , )/(  ×   ); dash-dotted line with stars is for < T e > = =  , =  =1, =1   (, , )/(  ×   ).We gather that < T He ++ > reaches 2.4 × 10 7 K while < T i > and < T e > stay around their initial average values of 1.2 × 10 7 K.Note at the domain edges initial temperature is 2 × 10 7 K and the temperature varies with y-coordinate according to Eq.2.All panels show data for the Run 4, which has   = 5000,   = 200 and  end  cp = 75.
In Fig. 3 we show scaling of the magnetic fluctuations power spectrum with wavenumber .In particular similar to Ofman (2010) we plot Fast Fourier Transform (FFT) of   (,  =  max /2,  =  max ) (solid line);   (,  =  max /20,  =  max ) (dashed line);  ,  (), according to Equation 10, (dash-dotted line).Thus solid line shows   fluctuations power spectrum in the middle of the high-density region; dash line shows the same but for a homogeneous, edge of the domain; dash-dotted shows the same using the following analytical expression: where   is left hand polarized wave number calculated for   = 0.3   , using the general expression with   = 0.3 cp .We note that the above expression of   is based on the dispersion relation for L-and R-polarized DAWs, which are as following (Tsiklauri 2011) , here upper signs are for the L-polarization and lower signs for the R-polarization.The latter equation indicates that for L-polarization possible physical frequency range is 0 <  <    and that there is a ion-cyclotron resonance at  =    .Protons can resonate because they rotate in the same direction as the electric field vector of the wave.In the case of R-polarization (not considered here) possible physical range is 0 <  <   with an electron-cyclotron resonance present at  =   .In Fig. 3  ,  () with dash-dotted line shows two expected peaks around   = 0.3   (at location  =   ) and   = 0.6   (at location  = 2  ).We also see the steepening of magnetic fluctuations power spectrum in the higher-density regions, possibly due to wave refraction, as in previous Hybrid simulation results (Ofman 2010), but now PIC runs produce much steeper slopes than in Hybrid runs.For example, in Fig. 3 we see the scaling as  −3.5 and  −4.9 for the in low density,   (,  =  max /20,  =  max ), and high density,   (,  =  max /2,  =  max ), regions, respectively.Corresponding power-law indexes for the (Ofman 2010) Hybrid run (see their Fig.15) are −1.7 and −2.5.Thus we see additional steepening and that electron-scale physics, which can only be resolved by a PIC simulation, has a notable effect of DAW spectrum evolution.As mentioned earlier, this comparison should taken with care because the considered physical parameters in the two studies are different: solar wind and flaring coronal loop tops.
In Fig. 4 we put to test Eq.( 7).In particular we would like to verify that Eq.( 7) inputs the same power via DAW driving as a single harmonic case considered in Tsiklauri (2011).In panel (a) we plot the following quantities  ∥, () for Run 1 (solid line);  ∥ , () for Run 2 (diamonds);  ∥, () for Run 3 (dashed line);  ∥ , () for Run 4 (triangles).We make two observations from panel (a) that: (i)  ∥, () for Run 1 and Run 2 overlap to a plotting precision, which means that single harmonic and a sum of 512 single harmonics with amplitudes 512 times smaller than Run 1 produces identical results; (ii)  ∥ , () for Run 3 and Run 4 overlap to a plotting precision, which means using 128 harmonics as per Eqs.5-6 with  from Eq.7, i.e. changing number of harmonic still delivers the same power compared to Run 4 (512 harmonics).A more general conclusion from panel (a) is that as far is electron acceleration along the magnetic field is concerned, the fraction of super-thermal electrons remains broadly similar with one uses single, 128, and/or 512 harmonics with frequencies in the considered range.Broadly similar conclusion can be reached about average temperature <   () > dynamics in panel (d).In other words broadband spectrum of DAWs does not affect much  ∥, () dynamics.Situation is quite different when we look at dynamics of the proton and He ++ hearing across the magnetic field indexes  ⊥, () and  ⊥, () shown in panels (b) and (c) respectively.Although as in panel (a) we see that Run 1 (solid line) and Run 2 (diamonds) data as well as Run 3 (dashed line) and Run 4 (triangles) data overlap to a plotting precision, Run 1/Run 2 attain higher values of  ⊥, () and  ⊥, () at later times compared to Run 3/Run 4. As these indexes capture fraction of accelerated particles, this result supports a conclusion that single harmonic (Run 1) or 128 harmonics with 128 smaller amplitude (Run 2) are more efficient in accelerating protons and He ++ , compared to broadband spectrum of waves case (Run 3 and Run 4).This can be understood by the fact that in Run 1/Run 2 DAW act coherently, while in Run 3/Run 4 there is also effect of 5% random phase present which "washes out" efficiency of particle acceleration which seems to be a coherent process based on Landau resonance.The heating due to DAWs is normally non-resonant, while the heating due to left-hand polarized ion-cyclotron waves is resonant with the resonant condition affected by the Doppler shift.We note that in our case Doppler shift is zero, as no drifts/flows of particles present.In panels (e) and (f) we see reverse situation compared to panels (b) and (c) in that broadband Run 3/Run 4 attain higher average temperatures compared to coherent Run 1/Run 2. This then clearly indicates that broadband spectrum that contains cyclotron resonance produces profuse heating of He ++ and to some extent protons.As summary of Figure 4 we remark that the frequency spectrum case does not seem to affect electron acceleration fraction in the like-to-like cases, but few times larger percentage of He ++ heating is seen due to ion cyclotron resonance; In Fig. 5 we show panels similar to Figure 1 but for Run 5. Purpose of Run 5 is to study how multiple DAW front collisions will affect particle acceleration.Thus in Run 5 we set  end  cp = 750 with 10-times longer simulation time, which insures that the DAWs will collide at least 9 times.We note in panel (a) that perturbation of   is no longer confined the density gradient regions and spread across ycoordinate, hence indicating oblique propagation across the magnetic field at small angle.We note that panel (b) is similar to panel (b) in Figure 1.The notable diffidence can be seen in panels (c)-(i) in that: (i) the whole density/temperature inhomogeneity has shifted its location from  =  max /2 to larger values of y-coordinate; (ii) noticeable bending of the entire background density/temperature inhomogeneity structure can be seen, i.e. edges at  = 0 and  =  max have lower y-coordinates than  =  max /2, i.e. we see development of kink oscillations when driving DAWs collide.When seeing animated version of Fig. 5 (not shown here) the background density/temperature inhomogeneity structure's  = 0 and  =  max edges move up and down in concert, while the middle part  =  max /2 oscillates out of phase, giving rise to kink-like motion.This kinking seems to come from DAW wave collision, which seems reasonable and intuitive from the first physical principles, and yet new according to the author's opinion.In panels (g) and (h) we see significant increase of temperatures in the density gradient regions, while for panel (i) the temperature increase is near the DAW injection location at  =  max /2, we think, this is due to the ion cyclotron resonance for He ++ .
In Fig. 6 we plot panels similar to Figure 2 but for Run 5. Commenting on the differences compared to Figure 2 we note the following new features: (i) now electron distribution function is also broadened in y-(panel (b)) and z-(panel (c)) directions.This is because of 10 times longer driving by DAWs injects much larger amount of energy into the system.Thus, electron population is heated up and we see more isotropic velocity distributions in all three spatial directions; (ii) we note in panel (h) an intermittent in time electron, proton and He ++ acceleration fractions.This is because intensive heating (plasma temperature increase) makes the-above-thermal-fraction smaller; Broadly speaking, when in Figure 2 we see  ⊥, () attains value of 0.63, i.e. 63% of He ++ become super-thermal.But in Figure 6  ⊥, () after attaining is peak value it starts to decease down to 0.25, this is because plasma temperature increase makes the-above-thermal-fraction (i.e denominator increase) smaller; (iii) now much higher temperatures are attained (panel (i)) by all plasma species namely <   >= 1.8 × 10 7 K, <   >= 1.4 × 10 7 K, <   ++ >= 5.3 × 10 7 K.The temperature increase for He ++ is the highest of all three species, because of ion-cyclotron resonance condition is met.In this context, we mention a new study that analyses relativistic electron measurements obtained by the High Energy   Telescope (HET) aboard Solar Orbiter in the energy range from 200 keV to 10 MeV (Fleth et al. 2023).The latter authors compile a list of 21 greater 200 keV energy electron time periods of enhanced flux.Out of this list they detect 3 events with clearly anisotropic electron pitch-angle distributions above 1 MeV.Also 5 other events have been seen within 0.5 AU from the Sun that show no discernible anisotropy.Hence, upcoming data from HET aboard Solar Orbiter will be able to test predictions made from numerical simulation studies such as this one.In our case in Figure 1 electron distribution function is anisotropic, while in Figure 6 it is isotropic for the reasons explained above.
In Fig. 7 we plot similar physical quantities to Figure 3 but for Run 5. We note that now we see further steepening of magnetic fluctuations power spectrum both in the higher-density regions and the domain edges.In particular, Fig. 7 we see the scaling as  −6.2 and  −8.0 for the in low density,   (,  =  max /20,  =  max ), and high density,   (,  =  max /2,  =  max ), regions, respectively.It is difficult to comment on the specific reason for the spectrum slope steepening as a result of multiple DAW collisions, but generally steepening of the slope means increased dissipation.Hence steepened slopes means that multiple DAW collisions seeded turbulence which made dissipation more intense.We admit, this is a speculative argument based on an intuitive judgement.
In Fig. 8 we show panels Figure 2 but for Run 6, which has the length-scale of the density gradient in -direction increased 10 times to approximately "ion"-scales / ≃ 40/ pe .Hence in Run 6 there are  = 5000 and  = 2000 grid points in each direction.The main purpose of this run is to study how particle acceleration is modified by up-scaling the transverse density gradient spatial size 10 times.There are two notable differences compared to Run 4 (Figure 2) in that (i) now He ++ peak heating by ion-cyclotron is higher with the maximum of  ⊥, () now attaining value of 0.79 (79%); (ii) The maximal value of <   ++ >= 5.6 × 10 7 K, which is also a higher value.From these finding we conclude that increasing density gradient 10 times enhances ion cyclotron heating for He ++ ions which is evidenced by the higher super-thermal percentages and also the higher the temperature gained by the He ++ plasma species.Hence the main conclusion of this Run 6 is that increasing density gradient scale across the magnetic field 10 times (almost to ion scales) does not affect the particle acceleration features by DAWs.We stress this point because the usual criticism such numerical simulations attract (Private Communication with various Plasma Physics colleagues) is that when one up-scales density gradient -the effect of parallel electric field generation and particle acceleration by DAWs goes away.Here we demonstrate otherwise.
In Fig. 9 we show similar physical quantities to Figure 2 but for Run 7. We note the following differences compared to Run 4 (Fig. 2) in that in this homogeneous plasma case, i.e. no transverse density/Alfven speed inhomogeneity across the magnetic field case, there is no electron acceleration present at all by the end simulation time, which supports the conclusion that the electron acceleration by DAWs is purely caused by the transverse plasma density inhomogeneity (panel (a)).There is an indication of (i) He ++ beams being formed along the magnetic field (panel (i)); and (ii) the bulk heating still occurring across the magnetic field due to cyclotron resonance in panels (j) and (k), by the end simulation time.The former seems not well under-    stood from physical grounds as to why homogeneity should support ion-cyclotron magnetic field-aligned resonant beams, while the latter is well understood by the said resonance.
In Fig. 10 we show similar physical quantities to Figure 2 but for Run 8. We note the following differences compared to Run 5 (cf.Fig. 6) in that in this homogeneous plasma case, despite 10 times longer end-simulation time of  end  cp = 750 we still not see any Bfield aligned electron acceleration (panel (a)).There is an indication of He ++ distribution functions become isotropic in   ,   ,   directions (panels (i),(j),(k)).From panel (h) we deduce that  ⊥, () attains value of 0.89, i.e. 89% of He ++ become super-thermal.From panel (l) we note that average temperatures of all plasma species became increased: <   >= 2.1 × 10 7 K, <   >= 2.5 × 10 7 K, <   ++ >= 1.3 × 10 8 K. Especially large increase is seen in <   ++ > due to continuous action of ion-cyclotron resonance.
In Fig. 11 we show similar physical quantities to Figure 2 but for Run 9. We note the following differences compared to Run 6 (cf.Fig. 8) in that in this homogeneous plasma case, despite 10 times larger domain-size, hence 10 times larger density gradient in transverse to the magnetic field y-direction, and  end  cp = 375 (maximum time that was computationally feasible), we again note no electron acceleration along x-direction (panel (a)).The other conclusions that we reached for Fig. 10 (Run 7) are applicable here in this Run 9 too.Namely that  ⊥, () is peaking at 0.89; and isotropization of He ++ distribution functions in   ,   ,   directions is seen; as well as average temperatures of all plasma species attaining values of <   >= 2.05 × 10 7 K, <   >= 2.3 × 10 7 K, <   ++ >= 1.01 × 10 8 K.

CONCLUSIONS
We would like to comment on the feasibility of the DAWs driving considered in this work and more generally what is known what the type of waves solar and stellar flares can generate.In the solar flare context, the first study which initiated investigation of electron acceleration by a single harmonic with sub-proton cyclotron DAWs   = 0.3 cp (to avoid proton cyclotron resonance) in the transversely inhomogeneous plasma is by Tsiklauri et al. (2005).Subsequently, Fletcher & Hudson (2008) explored a possibility when flare-generated reconfiguring coronal field launches a torsional AW pulse.We refer the reader to Figure 1 from Fletcher & Hudson (2008) which shows that AWs are launched at the loop top and propagate towards the footpoints.They study in detail feasibility of electron acceleration by estimating various parameters such as parallel electric field based on two-fluid theory.Although Fletcher & Hudson (2008) refer to AW pulse, estimates are made based on plasma kinetic dispersion relations, which are based on harmonic waves, like the ones considered in the present work.We believe initial AW pulse, rapidly phase-mixes on the transverse inhomogeneity and reaches short, kinetic scales.As can be seen from Fig. 1, panels (b) and (c), by the end simulation time DAWs develop only 3 wavelength by the end simulation time.So it is plausible that inital AW pulse generates L-hand polarized DAWs on such time-scale.Also, it is known that in the sub-proton cyclotron frequency regime, linearly polarized AWs naturally decay into and L-and R-hand circularly polarized DAWs (Krall & Trivelpiece 1986) because their phase speeds are different.Horne & Miyoshi (2016) suggest that electromagnetic ion cyclotron (EMIC) waves are important for electron acceleration and loss from the Earth radiation belts.They suggest that these waves are generated by unstable ion distributions that form during geomagnetically disturbed times.Horne & Miyoshi (2016) show that magnetosonic waves could be a source of EMIC waves as a result of propagation and a process of linear mode conversion.It is feasible the same mechanism also works in solar corona, with solar flares playing the same role as geomagnetic storms, which is essentially magnetic reconnection.The observed MHD wave amplitudes of 0.05 0 , i.e. 5% of the background magnetic field are typically observed.Such waves are routinely observed in the solar corona and also found in MHD simulations e.g.Li et al. (2023) reported the evolution of a traveling kink pulse to a standing kink wave along coronal loops that has been induced by a solar flare.As per their Figure 5 amplitude can be estimated roughly as zero-to-maximal deflection divided by 2 minutes (120 s)  = Δ/Δ = 12 Mm/120 s = 100 km/s which is about 10 % of their estimated Alfven seed of 950 km/s.Above we discussed MHD waves that were observed on spatial scales of 10,000-100,000 km, while the wavelengths of the waves considered in the present study are on order of 10 −2 m, 8-10 orders of magnitude shorter wavelength than the MHD waves.It may seem that the discussion of MHD scale waves is not relevant here, but presently only such scale waves are observed because DAWs are too small scale to be observable presently.Here we refer the reader back to the relevant discussion of unrealsitically high wave amplitudes in the kinetic regime considered in this study, and the objective computational reasons for that, stated in a paragraph just above equation 5.Moreover, there are some other mechanisms that can effectively excite small-scale KAWs in non-uniform media, such as the resonant mode conversion of AWs into KAWs (Xiang et al. 2019) and the density striation instability (Wu & Chen 2013;Chen et al. 2015).Therefore, it is meaningless to argue about the magnitude of the amplitude of the unmeasurable KAWs in the solar corona.Kryshtal et al. (2020) considered generation of low-frequency plasma waves in the lower/middle chromosphere region of magnetic loop foot-points when in the pre-flare state plasma.They showed that the generation of KAWs and kinetic ion-acoustic waves can occur both, in plasma with Coulomb conductivity and in the presence of small-scale Bernstein turbulence in the solar chromospheric setting.Ofman & Wang (2022) model the excitation and dissipation of slow magneto-sonic MHD waves in hot coronal flaring loops according to EUV observations using a 3D magnetohydrodynamic (MHD) visco-resistive simulations.Provornikova et al. (2018) contributes to an understanding of physical properties of observed solar coronal perturbations by investigating the excitation of waves by hot plasma injections from below and the evolution of flows and wave propaga-    tion along the solar coronal loop.Such ejections are can presumably come from solar flares.We remark that DAWs are observed in situ routinely in Earth magnetosphere (Lysak 2023) but the Sun is too far away for DAWs to be detected remotely.The 1/f spectrum of AWs naturally appears due to turbulent scale cascade and is a sign of scale-free (fractal-like) behaviour.
We would like to comment about feasibility of the spatial scale of transverse inhomogeneity considered in this work, namely, 4 − 40/ pe .In SI units for the physical parameters considered this corresponds to 0.21 − 2.1 m.For the real proton (for the mass ratio   /  = 1836) and electron inertial length are / pp = 2.3 m and / pe = 0.05 m, respectively.So in fact when we state that we consider inhomogeneity on "ion scale" 40/ pe = 2.1 m, it is not very different from the realistic / pp = 2.3 m.However, neither Parker Solar Probe nor Solar Dynamics Observatory can resolve magnetic loop transverse threads of such small width.The reason we consider such small scales is because current PIC models cannot consider realistic physical scales.This does not mean that such finescale transverse structuring does not exist and future observations may well detect it.
In this work we extended our previous results (Tsiklauri et al. 2005;Tsiklauri 2011;Tsiklauri 2012)) in five-ways: (i) instead of exciting DAWs at the left edge of the domain which mimics loop top we now drive DAWs in the middle of the domain, which more realistically represents DAWs generated at the solar coronal loop top during a flare; (ii) we now consider wide spectrum of DAWs as specified in Table 1; (iii) we added ion-cyclotron resonant He ++ species; (iv) we considered cases when DAW collide multiple times with now 10 times longer in time than our previous numerical runs; (v) we considered cases when the density gradient scale, and commensurately transverse to the magnetic field domain spatial size, is now increased 10 times compared to our previous results.Hence we studied cases when transverse density gradient is in the range 4 − 40/ pe and DAW driving frequency is 0.3 − 0.6 cp .Detailed discussion of the results established in this work can be found in the results section.
Here we mention only the main findings as following: (i) The considered DAW driving frequency spectrum 0.3 − 0.6 cp does not affect electron acceleration fraction in the like-to-like cases, but few times larger percentage of He ++ heating is seen due to ion cyclotron resonance; (ii) When we consider DAW multiple collisions, then much larger electron and ion acceleration fractions are seen, but the process is intermittent in time and more isotropic velocity distributions are also found.We attribute this to intensive heating (plasma temperature increase) makes the-above-thermal-fraction smaller; (iii) We witness formation of kink oscillations when DAWs collide, which seems a novel effect; We observe noticeable bending of the entire background density/temperature inhomogeneity structure, i.e. edges at  = 0 and  =  max have lower y-coordinates than  =  max /2, and we see development of kink oscillations when driving DAWs collide.We see also that the background density/temperature inhomogeneity structure's  = 0 and  =  max edges move up and down in concert, while the middle part  =  max /2 oscillates out of phase, giving rise to a kink-like motion.Due to the limitation of data storage, we could only output 20 snapshots of the kink oscillations from time zero to the final simulation time.This is not sufficient to infer the period of the kink oscillation, but we can confirm that this is the fundamental harmonic, with the maximal oscillation occurring at  =  max /2.(iv) We study scaling of the magnetic fluctuations power spectrum and find that steepening in the higher-density regions occurs, due to wave refraction.This means that inclusion of electron-scale physics has a notable effect of DAW spectrum evolution.

Figure 1 .
Fig.2has a purpose to quantify how particle acceleration main features, fractions and the budget, split between different plasma species, is modified by the introduction of  −1 frequency spectrum of DAWs with driving frequency in the range 0.3 − 0.6 cp .Thus, in Fig.2we show distribution function dynamics for electron   ,   ,  components (panels (a), (b), (c)); ion (proton)   ,   ,  components (panels (e), (f), (g)); He ++   ,   ,  components (panels (i), (j), (k)); In panels (a), (b), (c), (e), (f), (g), (i), (j), (k), the inner dashed curve show the distribution function (distribution by velocities) at  = 0, while outer solid curve at the final simulation time.Note that for protons and He ++ the inner dashed curves for  = 0 look almost as solid due to large number of data points (each distribution function is resolved with 100,000 points in all panels).By comparing panels (a), (b), (c), (e), (f), (g) with similar run e.g.panels (a),(b),(c),(d),(e),(f) in Fig.3fromTsiklauri (2011) we gather many similarities.In particular we see in panel (a) that electron beams are formed near the DAW phase velocity.A study ofTsiklauri (2012) has shown that the electron beam formation along the uniform background magnetic field is due to Landau resonance between electrons and DAW.In factTsiklauri (2012) showed if one varies ion(proton) mass, which prescribes DAW phase speed, then location of the beam peak in the distribution function changes accordingly, signalling occurrence of the Landau resonance.In panels (b) and (c) we see that distribution functions with respect to   and   at time zero and the final simulation time do overlap.Hence, we deduce from these panels that no sizable electron heating in the perpendicular to the magnetic field direction is seen.While, for protons we see in corresponding panels (f) and (g) that primarily perpendicular to the magnetic field broadening of the distribution occurs due to DAW energy injection/driving.No change in the proton parallel to the magnetic field distribution function (panel (e)) is seen, suggesting that on this time-scale proton parallel dynamics is not affected by DAW driving.In panel(d)  we show max(|E x /E x0 |) (solid line with triangles) and < |E x (x, y, t)/E x0 | > = =  , =  =1, =1|  (, , )/ 0 |/(  ×   ) (dashed line with diamonds), i.e. maximal and average of the abso-

Figure 5 .
Figure 5. Similar to Figure 1 but for Run 5.

Figure 6 .
Figure 6.Similar to Figure 2 but for Run 5.

Figure 7 .
Figure 7. Similar to Figure 3 but for Run 5.

Figure 8 .
Figure 8. Similar to Figure 2 but for Run 6.

Figure 9 .
Figure 9. Similar to Figure 2 but for Run 7.

Figure 10 .
Figure 10.Similar to Figure 2 but for Run 8.

Figure 11 .
Figure 11.Similar to Figure 2 but for Run 9.

Table 1 .
Table of numerical runs considered.