Electron acceleration driven by the lower-hybrid-drift instability: an extended quasilinear model

Density inhomogeneities are ubiquitous in space and astrophysical plasmas, in particular at contact boundaries between different media. They often correspond to regions that exhibits strong dynamics on a wide range of spatial and temporal scales. Indeed, density inhomogeneities are a source of free energy that can drive various instabilities such as, for instance, the lower-hybrid-drift instability which in turn transfers energy to the particles through wave-particle interactions and eventually heat the plasma. We aim at quantifying the efficiency of the lower-hybrid-drift instability to accelerate and/or heat electrons parallel to the ambient magnetic field. We combine two complementary methods: full-kinetic and quasilinear models. We report self-consistent evidence of electron acceleration driven by the development of the lower-hybrid-drift instability using 3D-3V full-kinetic numerical simulations. The efficiency of the observed acceleration cannot be explained by standard quasilinear theory. For this reason, we develop an extended quasilinear model able to quantitatively predict the interaction between lower-hybrid fluctuations and electrons on long time scales, now in agreement with full-kinetic simulations results. Finally, we apply this new, extended quasilinear model to a specific inhomogeneous space plasma boundary: the magnetopause of Mercury, and we discuss our quantitative predictions of electron acceleration in support to future BepiColombo observations.


Introduction
Inhomogeneities in magnetic field, velocity, density, temperature, etc. from fluid down to kinetic scales are commonly encountered in space and astrophysical plasmas (Amatucci 1999). The gradient associated to such inhomogeneous plasma regions is the source of "free" energy that can drive various kind of plasma instabilities. For instance, in the case of density gradient on scale close to the ion gyroradius, a situation commonly encountered in many plasma environments, the plasma is unstable against the so-called drift instabilities. These instabilities arise from the relative motion between ions and electrons, and turn out to be of paramount importance in shaping plasma boundaries found in space, allowing for strong anomalous mass and energy transport not achievable by standard collision-like diffusion.
In situ measurements show that lower-hybrid waves (hereafter LHW) with frequency close to the lower-hybrid frequency f LH ≈ √ ω ci ω ce are ubiquitous in magnetized space plasma environments. Such waves are commonly observed at Earth's magnetotail (Huba et al. 1978;Retinò et al. 2008;Zhou et al. 2009Zhou et al. , 2014Khotyaintsev et al. 2011;Norgren et al. 2012;Le Contel et al. 2017) and Earth's magnetopause (André et al. 2001;Bale et al. 2002;Vaivads et al. 2004;Graham et al. 2017Graham et al. , 2019Tang et al. 2020). In these two regions LHW are commonly observed in the vicinity of magnetic reconnection sites where strong density gradients do form. Their role on the onset (and/or relaxfederico.lavorenti@oca.eu ation) of magnetic reconnection has been addressed in the past and still represents a key point in the context of reconnection research (Daughton 2003;Lapenta et al. 2003;Lapenta et al. 2018;Yoo et al. 2020).
Moreover, LHW are also observed at plasma shock fronts such as the terrestrial bow shock (Walker et al. 2008), interplanetary shocks in the solar wind (Krasnoselskikh et al. 1985;Zhang & Matsumoto 1998;Wilson et al. 2013), and supernova remnants (Laming 2001). Finally, LHW have been observed in induced ionosphere of comet 67P (André et al. 2017;Karlsson et al. 2017;Goldstein et al. 2019), of the planets Venus (Scarf et al. 1980;Shapiro et al. 1995), Mars (Sagdeev et al. 1990) and at Earth's ionosphere (Reiniusson et al. 2006). In this context, space observations of supra-thermal electron populations in conjunction with LHW represents one of the basic points motivating the interest in the study of the interaction of these waves with electrons (Norgren et al. 2012;Zhou et al. 2014;Le Contel et al. 2017;Broiles et al. 2016;Goldstein et al. 2019).
More than being just a mechanism at work in space plasma environments, electron acceleration by LHW is a mechanism commonly used in tokamak experiments to heat electrons and hence the plasma along the toroidal magnetic field lines (Bécoulet et al. 2011;Pericoli-Ridolfini et al. 1999). This naturally suggests that LHW generated in space are an efficient driver for electron acceleration in space plasma environments (Broiles et al. 2016). However, the mechanism for the generation of LHW in laboratory differs much from the one at play in space plasmas.
Article number, page 1 of 12 arXiv:2104.05011v1 [physics.plasm-ph] 11 Apr 2021 A&A proofs: manuscript no. main In plasma fusion experiments the LHW are commonly excited by an external pump enabling to sustain the waves and so the electron acceleration process on long time scales with parameters controlled by the experimenter himself. In inhomogeneous space plasma instead LHW are generated by the development of plasma instabilities, which nonlinear saturation might reduce the efficiency of the acceleration process when compared to plasma fusion experiments.
In a natural plasma environment the two instabilities at play for the generation of LHW are (i) the modified-two-stream instability (hereafter MTSI) driven by a supra-thermal ion beam  and (ii) the lower-hybrid-drift instability (hereafter LHDI) driven by the relative drift between ions and electrons (Krall & Liewer 1971;Krall & Trivelpiece 1973;Gary 1993). The electron acceleration driven by LHW generated by the MTSI has been widely addressed in the literature using quasilinear theory Shapiro et al. 1999), full-kinetic simulations (Mc-Clements et al. 1993;Bingham et al. 2002) and experiments (Rigby et al. 2018). The MTSI is considered to be the typical source for the above mentioned LHW observations at plasma shock fronts (Krasnoselskikh et al. 1985;Shapiro et al. 1995;Shapiro et al. 1999) due to the reflection of a large fraction of solar wind ions by the shock front; however, LHW are routinely observed also in the absence of such beams of reflected ions, a condition that prevents to invoke the MTSI as the underlying mechanism. In such cases, the LHW are instead generated by the LHDI. The driver for the development of the LHDI has to be found in the "strong" density gradients reaching length scales of the order of, or even shorter than, the ion gyro-radius. This is the case, for instance, for the above-mentioned observations at Earth's magnetosphere and in cometary plasmas. From now on, we shall focus on LHW generated by the LHDI and their interaction with the electron population.
The LHDI fastest growing modes propagate perpendicular to both the density gradient, let say the x-direction, and the ambient magnetic field direction, the z-direction. The phase velocity is of the order of the ion thermal speed (Gary 1993). However, the LHDI modes are unstable over a narrow cone angle around this direction (in this case the y-direction) proportional to the square root of the ion-to-electron mass ratio k z /k y ≈ √ m e /m i (Gary & Sanderson 1978). This means that the oblique LHDI modes have a component of the phase velocity parallel (resp. perpendicular) to the ambient magnetic field of the order of the electron (resp. ion) thermal speed. As a result, electrons can be resonantly accelerated by LHDI fluctuations in the direction parallel to the ambient magnetic field through wave-particle interactions, mechanism hereafter called LHDI electron acceleration.
The goal of this paper is to investigate the efficiency of the LHDI electron acceleration. As of today, the state-of-the-art on this mechanism is provided by the analytical quasilinear model proposed by Cairns & McMillan (2005), hereafter called the QL model. Such QL model is well suited to study the early stage of the electron acceleration, but eventually it breaks down when the nonlinear feedback from the particles distribution to the wave becomes important. More precisely the QL model breaks down as soon as nonlinear effects locally modify the distribution function shape (i.e. around the resonant velocity), hereafter called nonlinear LD-like effects because analogous to the well known nonlinear Landau damping effects (see Brunetti et al. (2000) and references therein).
Other past works have addressed the problem of the interaction between pump-generated LHW and electrons including such nonlinear LD-like effects (Singh et al. 1996(Singh et al. , 1998Zacharegkas et al. 2016). However, the configuration adopted in these works is not well suited for space plasma configurations where LHW are typically driven by a plasma instability, like the LHDI, and not by an external pump. To our best knowledge, the LHDI electron acceleration mechanism has not yet been studied using a self-consistent full-kinetic model.
In the past, the efficiency of LHDI electron acceleration has been addressed by means of reduced analytical models like the QL model, mostly because of the limitations on computational resources. Indeed, the computational power required to solve the Vlasov equation for LHDI electron acceleration is challenging because (i) ion and electron kinetic physics must be included self-consistently in the model, meaning that a full-kinetic numerical simulation is required and (ii) the wave propagation occurs over an angle, in the y−z plane, perpendicular to the inhomogenity direction, the x-direction, meaning that three-dimensional numerical simulations are required. All in all, investigating the electron acceleration generated by the LHDI requires full-kinetic 3D-3V simulations. This is one of the methods used in this work.
In this paper we investigate the electron acceleration associated to the LHDI through a comparison of complementary numerical simulations. We use the quasilinear approach and direct full-kinetic 3D-3V simulations. This enables us to assess the intrinsic limits of the QL model and to investigate the consequences of nonlinear LD-like effects on the LHDI electron acceleration. We present first direct numerical evidence of LHDI electron acceleration from full-kinetic 3D-3V simulations, and we build up an extended quasilinear (eQL) model that takes into account the effect of such nonlinear saturation to quantitatively estimate electron acceleration under realistic space plasma parameters.
The paper is organised as follows: Sec. 2 describes the QL model and the full-kinetic 3D-3V simulation model we use here. Sec. 3 presents the results of both models using two common sets of plasma parameters ("strong" and "weak" gradient configurations). Sec. 4 compares the results of the two models, shows the limitations of the QL model as compared to the full-kinetic one, and presents a novel eQL model. Finally, we apply this new eQL model to the magnetopause of Mercury in view of the future observations of the BepiColombo space mission. Sec. 5 summarizes and concludes this work.

Models and methods
In this study, we use two different models of plasma evolution. First, a QL model of LHDI electron interaction based on the work of Cairns & McMillan (2005). Second, a full-kinetic 3D-3V plasma simulations of a plasma boundary initially unstable to the LHDI. The former is a simplified model of the plasma dynamics that does not account for the full response of the plasma itself to nonlinear interactions, and therefore is considered a reduced model. The latter instead is fully self-consistent, even if constrained by a specific parameter choice, and therefore it is considered an ab initio model. The full-kinetic model, being more general than the QL one, is used to assess the limits of the latter and eventually, to build an extended description that properly models the LHDI electron acceleration still within a quasilinear framework.

Quasilinear analytical model
The wave-particle interaction between LHW (generated from LHDI) and electrons can be modelled using a powerful analytical tool: quasilinear theory (Bernstein & Engelmann 1966;Alexandrov et al. 1984). Quasilinear theory is based on a second order perturbative expansion of the Vlasov equation averaged over the spatial variables. The system of quasilinear equations describes (i) the diffusion in velocity space of the electron distribution function trough a diffusion coefficient proportional to the electric field energy (Eqs. 1-2), and (ii) the time evolution of the electric field energy (Eqs. 3-4). The state-of-the-art QL model for LHDI electron interaction is the one developed in Cairns & McMillan (2005) and summarized here: There, f e (v , t) is the electron distribution function, k ⊥ (resp. k ) is the wavevector perpendicular (resp. parallel) to the ambient magnetic field, S k = E 2 k /8π is the electric field energy density in wavevector-space, S k,max is the maximum value of S k attained at saturation, n 0 is the plasma density, ω LH , ω ci are the lowerhybrid and ion cyclotron frequencies, ω(k ⊥ , k ) is the spectrum of the wave, γ LHDI is twice the linear growth rate of the LHDI and δ(x) is the Dirac delta function. In the following, ρ i is the ions gyroradius, and v thi = ρ i ω ci is the ion thermal speed.
After normalization, this nonlinear system of coupled partial differential equations (Eqs.1-4) is solved by numerical integration using a time staggered leapfrog scheme. In k-space, we limit the computation to the region of LHDI fastest growing modes, i.e. 0.7 < k ⊥ ρ e < 1 and 0 < k ρ i < 1, as done in Cairns & McMillan (2005). Since the wavevectors are limited to this region, the wave spectrum turns out to be more or less flat with a corresponding frequency for the fastest growing mode given by the lower-hybrid frequency, ω(k ⊥ , k ) ω LH . The grid in velocity space is chosen by testing the convergence of the solution.
The QL model, Eqs.
(1-4), depends on four parameters: m i /m e , ω pe /ω ce , γ LHDI and S k,max . The first two parameters m i /m e , ω pe /ω ce define the plasma itself. The last two parameters γ LHDI and S k,max (or analogously S max = d 3 kS k,max ) define the linear growth and saturation level of the instability, and they only depend on the initial plasma configuration. The analytical expressions for these two quantities are given by: where β i is the ion plasma beta, and the inverse gradient scale length n is defined as The growth rate γ LHDI in Eq. (5) has been obtained using a linearized kinetic model by Davidson et al. (1977).The electric energy at saturation S max in Eq. (6) has been obtained using an analytic quasilinear approach by Davidson (1978), and it has been later tested numerically by Brackbill et al. (1984) using 2D fullkinetic simulations. The LHDI can saturate through two different processes depending on the initial value of the density gradient: ion trapping (resp. current relaxation) for high (resp. low) values of the density gradient (and therefore of the drift speed) (Brackbill et al. 1984). The derivation of Eqs. (5-6) is based on the assumption that the only source of drift in the plasma is the density gradient. Thus, all particle drifts, except the diamagnetic drift v Di , are considered negligible. As a consequence, v Di /v thi = n ρ i . In the fullkinetic simulations, presented in the next Section, this assumption is substantially verified.
Quasilinear models are inherently limited because they do not include the nonlinear feedback from the modified plasma dispersion function on the electromagnetic fields. To overcome this limitation, we present in the next Section a full-kinetic 3D-3V numerical plasma model.

Setup full-kinetic 3D-3V simulations
The full-kinetic model of the plasma is based on a direct solution of the Vlasov-Maxwell system of equations using a Lagrangian PIC (particle-in-cell) approach.
To run the simulations of a plasma boundary unstable to the LHDI, we have used the explicit, electromagnetic, relativistic, PIC code SMILEI (Derouillat et al. 2018). The ambient magnetic field is directed along the z-axis and the density gradient along the x-axis. In order to model both wave propagation (predominantly along the y-axis, perpendicular to both the magnetic field and the density gradient direction) and the electrons waveparticle interaction (predominantly along the z-axis, parallel to the magnetic field), we consider a three-dimensional (3D) numerical box. Compared to previous numerical investigations of the LHDI (Brackbill et al. 1984;Gary & Sgro 1990;Hoshino et al. 2001;Shinohara & Hoshino 1999;Lapenta & Brackbill 2002) that focused on the wave generation mechanism only through 2D-3V simulations in the equivalent of our (x, y) plane, we also include in this study the out-of-plane direction that hosts the electron acceleration resonant processes. An overview of the numerical setup is shown in Fig. 1, the right panel shows the 3D numerical box used here, highlighting (i) a slice of the ion density field in the (x, y) plane, and (ii) the plane (y, z) most unstable to the LHDI in yellow.
The simulations are initialized ensuring pressure balance and using the Vlasov equilibrium proposed in Alpers (1969); Pu et al. (1981), which shapes a plasma boundary with density and magnetic field asymmetries, uniform temperature, no electric field and no velocity nor magnetic field shear. Hereafter, we shall refer to the side I (resp. side II) of this boundary as the high (resp. low) density side. The expressions for the initialization profiles are: Where, T = T i + T e is the uniform temperature, PBC is the constant of pressure balance (set to have the plasma beta β I = 10/3 and β II = 1/12), and w is a constant that defines the width of the layer -set to 0.98 (resp. 0.5) in the "strong" (resp. "weak") gradient case. We use a density and magnetic field asymmetry of n I /n II = 10 and B I /B II = 0.5. The shape of these initialization profiles is shown in the left panel of Fig. 1. In the following, all quantities are normalized to ion quantities in side I: the ion gyrofrequency ω I ci =eB I /m i c, the ion skin depth d I i =c m i /4πn I e 2 and the Alfvén speed V I A =B I / 4πm i n I . The numerical box dimensions are L x =L y =2π and L z =20π with a number of cells N x =N y =288 and N z =1472. The box is elongated -by a factor √ m i /m e =10 -in the magnetic field direction in order to reliably include the narrow cone angle over which unstable LHDI modes grow. The timestep used in the simulations is dt = 3.4 × 10 −4 to satisfy the CFL stability condition. The simulations are ran for several tens of the ion gyro-periods until electron acceleration saturates. We use 100 macro-particles per cell and a second order spline interpolation for the macro-particles. We use a reduced ion-to-electron mass ratio m i /m e = 100 for the sake of computational reasons. The implications of using such a reduced mass ratio are discussed in Sec. 4. We use an ion-to-electron temperature ratio T i /T e = 1, and a plasma-to-cyclotron frequency ratio of ω I pe /ω I ce = 4. The magnitude of the density asymmetry used in these simulations encompasses the typical parameters observed in small planetary magnetospheres such as Mercury (Gershman et al. 2015).
Two different simulations are investigated using two different layer widths: (i) a steeper boundary case (with the inverse gradient scale length n =1, defined in Eq. 7), hereafter called "strong gradient" simulation, and (ii) a smoother boundary case (with n =0.5) hereafter called "weak gradient" simulation. These parameters are chosen, first, in order to ensure that the LHDI fluctuations amplitude are well above the PIC noise level, and, second, in order to saturate through the two different mechanisms introduced in Sec. 2.1: ion trapping or current relaxation.

Evolution of LHDI and associated electron acceleration: simulations results
First, we show the results of the full-kinetic model. Then, we show those obtained with the QL model using the same parameters as for the full-kinetic simulations. Finally, we compare the results of the two models showing the range of validity and the limitations of the QL model.

Results from the full-kinetic 3D-3V simulations
In both "strong" and "weak" gradient, full-kinetic simulations the layer is unstable to the LHDI due to the presence of a density gradient on ion kinetic scales. The LHDI fluctuations grow exponentially in the layer for times t < t sat as predicted by kinetic linear theory. The fastest growing mode (hereafter FGM) is electrostatic and directed along the y-axis with wavevector k y ρ e ≈ 1, frequency ω ω LH and a growth rate γ ω LH in agreement with the linear estimation (Eq. 5) for both simulations. The growth of the electric field energy -normalized to the ion thermal energy and integrated over the unstable layer -is shown in Fig. 2, green curves, for both simulations. At t ≈ t sat (corresponding to the first vertical dashed lines in each panel of Fig. 2), the electric field fluctuation's growth saturates. Using the growth rate and the saturation level from Eqs. (5-6), we compute the saturation time analytically as The saturation mechanism of the LHDI have been extensively studied in the literature in the past, see for instance Brackbill et al. (1984); Chen & Birdsall (1983); Davidson (1978). In the strong gradient simulation the LHDI is observed to saturate due to ion trapping, while in the weak gradient simulation, the LHDI saturates due to current relaxation, as expected. In both cases, the saturation levels are comparable with the ones obtained by Eq. (6). Subsequently, in the strong gradient simulation, we observe for t > t sat a rapid decrease of the electric fluctuations amplitude, characteristic of an overshoot pattern, shown in the electric-to-thermal energy ratio in Fig. 2 (left panel, green curve) at time t = t sat . Such overshoot of the electric field fluctuations is not observed in the weak gradient simulation. An important point to mention is that, in both simulations, the electric field energy always remains much smaller than the thermal energy of the particles (it never exceeds 1%). Note that a common issue with PIC simulations is the possible occurrence of spurious numerical heating of the macroparticles population during the simulation. This effect, due to the so-called "finite grid instability" (Birdsall & Langdon 1991;Markidis et al. 2010;McMillan 2020), appears when the Debye length is not resolved enough by the numerical grid. In such a case, the electrons would be numerically heated, which would hide the wave-particle interaction we are looking for. We have carefully checked that our PIC simulations are free from such numerical spurious heating, by comparing the evolution of the supra-thermal electron density in the unstable layer to the one outside the layer on both sides, where the plasma is stable, to check that the numerical heating of the electrons is negligible compared to the one arising from wave-particle interaction in the layer.
To quantify the efficiency of the LHDI in accelerating electrons parallel to the ambient magnetic field, we define a suprathermal electrons density as follows: In the following, the supra-thermal electrons density N e,sup is used as a quantitative tracer of the LHDI electron acceleration.
The growth and saturation of N e,sup is shown for both simulations in Fig. 2, red curves. In both simulations, the efficiency of the LHDI electron acceleration process can be described through a three-phase evolution: (i) the linear phase 0 < t < t sat , (ii) the quasilinear phase Article number, page 4 of 12 Lavorenti et al.: LHDI electron acceleration in space plasmas Fig. 2. Full-kinetic simulations: evolution of the electron supra-thermal density N e,sup (t)/N e,sup (0) (red curves) and of the electric field energy normalized to the ion thermal energy S /nT i (t) (green curves), both curves are integrated over the unstable layer (3.5 < x < 4 for strong gradient and 3.5 < x < 4.5 for weak gradient). Note that time axes are different for the two simulations. t sat < t < τ NL and (iii) the strongly nonlinear phase t > τ NL . The characteristic time scales t sat and τ NL are highlighted in Fig. 2 by vertical dashed lines. The latter time scale constitutes the core of the full-kinetic modelling, therefore, an extensive discussion is provided in Sec. 4.2. In the linear phase, the electric field grows exponentially as predicted from linear theory. As expected, no electron acceleration is observed. In the quasilinear phase, electron acceleration starts to take place as observed in both simulations. In this phase the efficiency of the acceleration is directly related to the driver intensity (i.e. the stronger the electric field the stronger the acceleration), as expected from the QL diffusion equations (Eqs. 1-2). Finally, in the strongly nonlinear phase, although the electric field amplitude remains constant and finite, the acceleration stops due to the onset of nonlinear LD-like effects. Such a multi-phase evolution is observed in both strong and weak gradient simulations (see Fig. 2). A main difference between the two simulations is the presence -in the strong gradient simulation -of an overshoot in the electric energy, not observed in the weak gradient simulation. Such overshoot corresponds to a sharp peak in the electric-to-thermal energy ratio at t ≈ t sat , which subsequently leads to a strong electron acceleration for the overshoot duration (see Fig. 2, left panel around t sat ). However, the overall contribution of this phenomenon to the total amount of accelerated electrons is shown not to be significant in Sec. 4. In the end, for both simulations, the supra-thermal electron density in the inhomogeneous layer is increased by around 15% − 20%. However, in the two cases this same final value is attained through a different evolution: in the strong gradient simulation the acceleration is faster but stops sooner (at τ NL ≈ 10), while in the weak gradient simulation the acceleration occurs at a slower rate but remains efficient on longer time scales (up to τ NL ≈ 40). Interestingly, these two effects seems to compensate each other.
The nonlinear time scale τ NL , not accessible to quasilinear models, represents the original result of our full-kinetic simulations; an extensive discussion is left to Sec. 4.

Results of the standard quasilinear model
To enable a quantitative comparison between the QL model described in Sec. 2.1 and the full-kinetic one, the former is solved using the same plasma parameters as for the unstable layer of the full-kinetic simulations. More precisely, the input parameters of the QL model are obtained by averaging the plasma density, and the magnetic field in the layer 3.5 < x < 4 (resp. 3.5 < x < 4.5) for the strong (resp. weak) gradient full-kinetic simulation.
Numerical integration of the QL model (Eqs. 1-4) provides the time evolution of the electron distribution function f e (v z ) shown in Fig. 3 at different time instants. We observe a diffusion in velocity space for both "strong" and "weak" gradient cases, corresponding to electron acceleration by LHW. The characteristic timescale for such acceleration turns out to be of the order of hundreds of ion cyclotron periods. As expected, this time scale is longer in the weak gradient case than in the strong gradient one.
Eventually, electron acceleration described by QL theory slows down after long times. Indeed, as already explained in Cairns & McMillan (2005), the characteristic time scale to accelerate an electron of velocity v by an amount δv scales as ∼ v 5 for δv ∼v . With the parameters used in our two simulations, the electron acceleration obtained from the QL model becomes negligible after time τ Di f f 150 (resp. τ Di f f 400) in the strong (resp. weak) gradient case, as shown in Fig. 3.
Such LHDI electron acceleration appears much weaker than those presented in the work of Cairns & McMillan (2005). This discrepancy is directly associated to the choice of the parameter S max /nT i , the electric-to-thermal energy at the saturation of the LHDI, which is proportional to the quasilinear diffusion coefficient (Eq. 2). In this study, we use S max /nT i = 0.01 (resp. 0.002) for the strong (resp. weak) gradient case, to obtain the results shown in Fig. 3. These values are obtained from Eq. (6) and are consistent with those observed in our full-kinetic simulations. Differently, in the seminal paper of Cairns & McMillan (2005) the authors have used a value S max /nT i = 0.5 that is two orders of magnitude larger than what is expected from the theory of saturation of the LHDI in Davidson (1978) using the configuration considered in that study. Such choice of an overestimated value of the parameter S max /nT i unavoidably leads to an overestimation of the quasilinear diffusion coefficient, and therefore of an higher electron acceleration efficiency. This points out that a consistent choice of the electric-to-thermal energy ratio at saturation is crucial in quasilinear models in order to accurately assess a reliable value of the quasilinear diffusion coefficient.
Despite the different choice of the parameter S max /nT i , the time evolution of the electron distribution function obtained with the QL model agrees qualitatively with the one presented in Cairns & McMillan (2005).
Nonetheless, the comparison with the evolution obtained from the full-kinetic model is left to the next section (Sec. 4.1).

Discussion: towards an extended quasilinear model and beyond
In this section, we first compare the full-kinetic and standard quasilinear results, we show the discrepancies between the two models, and we highlight the physical processes that give rise to such discrepancies. Then, in order to overcome the intrinsic limits of the standard quasilinear theory, we build an extended quasilinear model (eQL) that includes the consequences of nonlinear LD-like effects. Finally, we validate such eQL and extrapolate how it scales with the plasma parameters of interest (e.g. Article number, page 5 of 12 A&A proofs: manuscript no. main ion-to-electron mass ratio, gradient length) in order to enable its use beyond the two specific set of parameters of our full-kinetic simulations.

Comparison between the full-kinetic and the QL models
To quantitatively compare the results of the full-kinetic and QL models on LHDI electron acceleration, we focus on the evolution of the supra-thermal electron density (defined in Eq. 12), shown for both strong and weak gradient cases in Fig. 4, where the two characteristic times t sat and τ NL -previously identified from the full-kinetic simulations -are reminded for sake of clarity (vertical dashed lines). The three-phase evolution of LHDI electron acceleration observed in the full-kinetic simulations (red lines in Fig. 4) is not explicable in terms of the QL model (blue lines in Fig. 4). Two main discrepancies are identified. First, we observe a minor discrepancy between the two models around t = t sat in the strong gradient case, see Fig. 4 left panel. At this time, the supra-thermal electron density obtained from the full-kinetic simulation (red curve) is higher than that obtained from the QL model (blue curve). This enhanced electron acceleration is due to the overshoot of the electric field in the full-kinetic simulation, effect not included in the QL model, and not observed in the weak gradient case. However, this short and sharply peaked phenomenon brings a negligible contribution to the total amount of supra-thermal electrons at the end of the simulation.
Second, we observe a strong discrepancy between the two models for times t > τ NL (phase III in Fig. 4). Indeed, on the one hand, the electron acceleration stops at time t ≈ τ NL in fullkinetic simulations (red curve) due to the onset of nonlinear LDlike effects, while on the other hand, the electron acceleration goes on in the QL model (blue curve) up to time t ≈ τ Di f f . The fact that the value of τ NL (obtained from full-kinetic simulations) is about one order of magnitude lower than τ Di f f (obtained from QL simulations) is the main reason for the strong discrepancy between the two models. Such discrepancy is due to the fact that the QL model does not include the nonlinear LD-like effects that are responsible for the abrupt stop of electron acceleration observed at t = τ NL in the full-kinetic simulations. Thus, even though the two models are in good agreement in the linear and quasilinear phases (phases I and II in Fig. 4), the final value of the supra-thermal electron density is strongly overestimated by the QL model. In the following, we shall use the increase in the supra-thermal electron density, defined as: to quantify the discrepancy between the models. In particular, the value of ∆N e,sup is 50 (resp. 20) times higher in the strong (resp. weak) gradient case for the QL model as compared to the full-kinetic model.

Need for an extended quasilinear model
Quasilinear theory is a powerful tool to study wave-particle interaction in an analytical framework due to its relative simplicity. Therefore, it represents an ideal way to provide quantitative predictions on LHDI electron acceleration/heating. However, due to its derivation via a perturbative approach, QL theory does not include "strong" nonlinear effects (i.e. those effects, arising at sufficiently long times, that alter the wave fluctuations due to the feedback from the modified distribution function to the fields). Although they are not included in the QL model, these effects are well reproduced by the full-kinetic simulations given its ab initio nature. The comparison between both models (Sec. 4.1) highlights the need to build an extended QL model (hereafter called eQL) that overcomes the intrinsic limitations of standard QL model and provide a dynamics consistent with that observed in full-kinetic models. In this purpose, we argue that such eQL model: does not require to include the possible overshoot of the electric field, as we have shown it is not a dominant process in energizing electrons, requires to include the eventual inhibition of LHDI electron acceleration by nonlinear LD-like effects through the parameter τ NL .
The nonlinear time τ NL is the parameter used in the eQL to define when nonlinear LD-like effects become dominant in the layer, inhibiting efficient wave-particle interactions. From such time, the QL diffusion coefficient (Eq. 2) becomes negligible even though the amplitude of the electric field remains constant, as shown in Fig. 5 for both full-kinetic simulations. This is due to its integral dependence on k 2 /k 2 ⊥ . The underlying physical process is the energy transfer, in k-space, from resonant, oblique modes to out-of-resonance, more perpendicular modes. Therefore, we explicitly switch off the QL diffusion coefficient for t > τ NL in the eQL model (Eq. 14-15).
We now focus on the estimation of this new parameter τ NL . Previous works on the LHDI, addressing its nonlinear evolution using 2D full-kinetic simulations, suggest a connection between LHDI modes and drift-kink-instability (hereafter DKI) modes (Pritchett et al. 1996;Shinohara & Hoshino 1999;Lapenta & Brackbill 2002). DKI modes are characterized by (i) frequency and growth rate γ DKI smaller than the ion gyrofrequency and (ii) wavevectors perpendicular to the magnetic field (Daughton 1998(Daughton , 1999. Following these studies, we assume that the mechanism underlying the energy transfer from oblique to strictly perpendicular modes at t > τ NL , observed in our full-kinetic simulations, is a LHDI-DKI coupling. Consequently, we estimate that the time of onset of LHDI nonlinear LD-like effects scales as the inverse linear growth rate of the DKI, i.e. τ NL ∝ 1/γ DKI . The DKI linear growth rate γ DKI scaling with plasma physical parameters has been addressed in Daughton (1998) using kinetic theory, and validated in Daughton (1999) using 2-fluid theory. In our work, we use such results to build our eQL model (Eq. 16).
Taken all those considerations into account, the eQL model reads: where τ 0 = 1.5 is a constant obtained from our two full-kinetic simulations. This value is obtained by fitting the evolution of the supra-thermal electron density of the full-kinetic simulations with the one of the eQL model.

Validation of the extended quasilinear model
The eQL presented in Sec. 4.2 (Eqs. 14-16) retains all the advantages of the standard QL model (namely being analytical and easily integrable using a numerical solver) while extending its range of validity to asymptotically long time scales. This now enables to provide quantitative predictions of LHDI electron acceleration. Indeed, as shown by the time evolution of the suprathermal electron density in Fig. 4, unlike the standard QL predictions (blue), the eQL model (orange) is able to reproduce the predictions of the nonlinear Vlasov-Maxwell theory (red). The relative discrepancy between the two curves does never exceed 5% in both simulations. This validation of the eQL model is of particular interest to quantify the LHDI electron acceleration over long times.
The results of the eQL model strongly depend on the estimation of the nonlinear time τ NL , identified as the inverse linear growth rate of the DKI 1/γ DKI . This choice is justified both (i) a priori, by previous numerical works that showed evidence that the long time evolution of a LHDI-unstable layer gets coupled to a DKI (Pritchett et al. 1996;Shinohara & Hoshino 1999); and (ii) a posteriori, by validating that the scaling γ DKI ∝ ( n ρ i ) 2 in Eq. (16) is consistent with the outputs of our two full-kinetic simulations.
In our two full-kinetic simulations -that have enabled to both define and validate the eQL model -we have used a reduced ion-to-electron mass ratio m i /m e = 100. Using such a reduced mass ratio is a standard procedure in PIC plasma simulations, since it allows to reduce the scale separation among the two species in order to run the numerical simulation on reasonable amount of CPU time (still about one million computational hours per simulation, in our case), while maintaining the necessary temporal and spatial scales separations between both species dynamics. However, the properties of the LHDI actually depend on the ion-to-electron mass ratio, because of the hybrid character of the instability. Therefore, the quantitative predictions of our full-kinetic simulations are not directly applicable to a realistic plasma configuration (with a physical proton-toelectron mass ratio 1836). This is exactly where an eQL model becomes extremely useful.

Using the extended quasilinear model to assess LHDI electron acceleration at physical mass ratio
The eQL model does not suffer from the strong computational constraints of full-PIC simulations and enables us to address, in this section, the question of LHDI electron acceleration using a realistic proton-to-electron mass ratio. The results of the eQL models are summarized in Fig. 6, showing both the electron distribution function evolution averaged in the inhomogeneous layer (top panels) and the resulting supra-thermal electron density (bottom panels), for both strong and weak gradient setups. Here, we use the two setups described in Sec. 2, where only the mass ratio is modified to address the case m i /m e = 1836. We emphasize that full kinetic simulations using such a physical mass ratio would have been extremely challenging computationally, so that we consider the eQL model as a way to extrapolate the results of our full-kinetic simulations (done using a reduced mass-ratio) to a physical plasma environment.
Compared to results shown in Sec. 4.3 for a reduced mass ratio, using a physical mass ratio the nonlinear time increases following Eq. (16). This means that the LHW-particle interaction occurs over longer times, thus more electrons are accelerated. However, on the one hand, this effect can be compensated by the reduction of the electric field energy of order ∼ m e /m i (Eq. 6) if the LHDI saturates through current relaxation, i.e. in the weak gradient case, see Fig. 6 right panel. On the other hand, if the LHDI saturates through ion trapping the saturation level does not depend on the mass-ratio (Eq. 6), therefore, the LHDI electron acceleration is more efficient by a factor ∼ m i /m e due to the increase in the nonlinear time. All in all, the mass ratio does not affect the weak gradient case, while it increases the fraction of LHDI accelerated electrons by a factor ∼ 1836/100 in the strong gradient case. These predictions -solely based on scaling arguments -are well reproduced by the numerical solution of the eQL using a realistic proton-to-electron mass-ratio, as shown in Fig. 6, and are also confirmed by the analytical estimates developed below in this section. Now, we focus on the value of the supra-thermal electron density at asymptotically long times, since, all in all, this is the quantity relevant for space plasma observations of LHDI accelerated electrons. Under some simplifying assumptions we derive an analytical expression for this quantity in the framework of the eQL model.
First, we approximate the evolution of the supra-thermal electron density by a linear interpolation in the time interval t sat < t < τ NL to get which is constant in the interval t sat < t < τ NL since we assume that (i) the distribution function is weakly modified by the interaction with the wave, i.e. f e (t, v ) ≈ f e (0, v ), and (ii) the amplitude of the resonant electric field wave |E(k ⊥ , k )| remains constant after the saturation of the LHDI for t > t sat . Finally, under the previous assumptions, using the expression for the QL diffusion coefficient (Eq. 2), assuming a Maxwellian distribution function for the electrons, the Eqs. (17-18) lead to Before going on, we stress here that the increase in supra-thermal density ∆N e,sup is proportional to the amplitude of the electric field at saturation S max . This emphasizes how the input parameter S max /nT i impacts the output of QL theory and further supports the discussion regarding the difference between our results and the ones by Cairns & McMillan (2005) in Sec. 3.2. Finally, expressing the electric field at saturation S max using Eq. (6), the nonlinear time τ NL using Eq. (16), the saturation time t sat using Eq. (11), and the LHDI growth rate γ LHDI using Eq. (5), and under the assumption ω pe > ω ce , the supra-thermal density increase in Eq.  where χ reads: With the parameters considered in this study, and typically encountered in space plasmas, χ m e m i 1. As a consequence, the supra-thermal electron density increase at long times is well approximated analytically by  (22) depending on the LHDI saturation mechanism at play.
This estimation leads to a total increase in the supra-thermal electron density summarized in Tab. 1 for the different sets of parameters used throughout this work. These values are also shown as horizontal grey dotted lines in both panels of Fig. 4 -for reduced mass ratio -and Fig. 6 for a physical mass ratio. From the values in Tab. 1, we infer that our analytical approximation (Eq. 20) is valid (error of the order of 10%) in the limit of "weak" LHDI electron acceleration (i.e. ∆N e,sup /N 0 e,sup 1). Stronger discrepancies arise for stronger electron acceleration (see third row in Tab. 1), as expected.

Application to Mercury's magnetopause
Mercury's magnetopause represents an excellent "textbook" example of a plasma boundary with a ion kinetic scale density gradient potentially LHDI-unstable. Previous space missions at Mercury -Mariner 10 (Russell et al. 1988) and MESSEN-GER (Solomon et al. 2007) -did not bring an instrumental payload capable of providing simultaneous measurements of electric field in the lower-hybrid frequency range and of the electron distribution function. Therefore, the expected LHW physics at Mercury cannot yet be tackled from past observations. However, the physics at these scales/range is exactly one of the main scientific objectives of the ongoing ESA/JAXA space mission BepiColombo (Benkhoff et al. 2010). Different plasma instruments onboard the Mio spacecraft will provide measurements of the electron distribution function (MPPE/MEA, Mercury Plasma Particle Experiment / Mercury Electron Analyzer) and of the electric field in the lower-hybrid frequency range (PWI, Plasma Wave Investigation) enabling the possibility to directly address the physics of Mercury's strongly inhomogeneous magnetopause. In this section, we quantify the expected efficiency of LHDI electron acceleration at Mercury's magnetopause using the eQL model developed in our work in support for future BepiColombo observations.
As previously pointed out in Sec. 2.2, the parameters used in our strong and weak gradient cases are taken in the range expected at Mercury's magnetopause (Slavin et al. 2008;Gershman et al. 2015). Therefore, we make use of the results of our eQL model with physical mass ratio in Sec. 4.4 to assess the importance of electron acceleration driven by the LHDI at Mercury's magnetopause.
First, we discuss the features of the LHW that are expected to be generated by the development of the LHDI at Mercury's magnetopause. These waves have frequency f ≈ f LH ≈ 50-100 Hz and wavelength λ ≈ 2π/ρ e ≈ 5-15 km (using hereafter typical magnetic field values at Mercury's magnetopause of 10-30 nT, ion temperature 30-70 eV, and density 10-30 cm −3 , see Milillo et al. (2020)). Assuming that those waves are advected by the shocked solar wind flow with a speed V S W ≈ 100 km/s, the Doppler-shifted frequency of the LHW lies in the range f = f ± kV S W ≈ 0-200 Hz. Moreover, the electric field amplitude of these LHW at saturation is E ≈ 10-100 mV/m, as obtained from Eq. (6). These typical frequency and energy range of LHW at Mercury suggests the use of electric field instruments onboard Mio spacecraft to address the LHW physics at Mercury's magnetopause. In particular, the sensors MEFISTO and WPT -part of the PWI consortium -can provide electric field observations in this frequency range with a sensitivity of the order of 10 −3 mV/m, well-below the expected amplitude of these waves (Kasaba et al. 2020).
Second, we discuss the features of electron distribution functions possibly modified by resonant interaction with the previously discussed LHW. This resonant wave-particle interaction accelerates sub-thermal electrons (with speed v z v the ) to suprathermal energies (with speed v z 2v the ) in the direction parallel to the ambient magnetic field, as shown in Fig. 6 top panels. In the following we assume an unperturbed electron temperature of 50 eV (typical values at Mercury's magnetopause being ∼ 20 − 100 eV (Ogilvie et al. 1974)). This acceleration process is well in the range of observations of the electron instrument MPPE/MEA onboard the Mio spacecraft of the BepiColombo mission (Saito et al. 2010), since this instrument includes two electron analyzers that can measure the three-dimensional en-ergy distribution of electrons in the range 3-3000 eV (in solar wind mode) or 3-25500 eV (in magnetospheric mode), with a time resolution of 1 second (Milillo et al. 2020). Here, we simulate the instrumental response of MPPE/MEA when encountering electrons accelerated by a LHDI at Mercury's magnetopause, as shown in Fig. 7, with the unperturbed electron distribution function in blue, and the one resulting from interaction with LHDI at long times in red. In Fig. 7, the uncertainties on the simulated response of the sensor MEA are obtained using the uncertainty on the G-factor of the instrument reported in Saito et al. (2010) (of the order of 10%). So, on the one hand we predict that the modifications in the electron distribution function above 100 eV are observable by MPPE/MEA in the case of a strongly inhomogeneous magnetopause (width around one ion gyro-radius, or less), as shown in the left panel of Fig. 7. On the other hand, since this process is very sensitive to the width of the magnetopause layer, we expect negligible modifications in the electron distribution function for less inhomogeneous magnetopause conditions (width around two ion gyro-radii or more), as shown in the right panel of Fig. 7. With typical magnetic field and temperature values at Mercury's magnetopause the ion gyroradius is around 20 − 80 km.
We also assess that the interaction between LHDI and electrons increases both (i) the density of supra-thermal electrons N e,sup (i.e. electrons with energy higher than 100 eV, see Eq. 12), and (ii) the electron temperature. In the simulated responses with strong (resp. weak) density gradient, the former increases by 20% (resp. 0.3%), due to the interaction with the LHW, and the latter increases by 80% (resp. 1.5%). We stress here how these two scalar quantities could be used as tracers of LHDI-electron interaction events in MPPE/MEA data, and how this response is limited to the direction parallel/anti-parallel to the ambient magnetic field (i.e. around θ = 0 • and 180 • in the electron pitch-angle distributions). Third, we present electron observations made by the NASA spacecraft Mariner 10 during its first Mercury flyby on 29 March 1974 (Ogilvie et al. 1974). In Fig. 8 we show magnetic field (upper panels) and electrons (bottom panels) data for the inbound (left) and outbound (right) magnetopause crossings. In particular, we observe a bimodal character of the electron energy distribution during the crossing (bottom panels), with one ambient electron population around 50-70 eV and a second energized population around 300-500 eV. This is consistent with the signatures expected for LHDI-accelerated electrons. Although these observations suggest that the LHDI might play a role in the electron energization at Mercury's Magnetopause, the Mariner 10 data cannot unambiguously confirm this hypothesis because of the lack of electric field observations. Moreover, the low time resolution of Mariner 10 electron measurements (1/6 s −1 ), the narrow energy range (13.4-687 eV) and the lack of a statistically significant number of crossings make the Mariner 10 observations prevent from a conclusive evidence of LHDI-accelerated electrons at Mercury's magnetopause.
The limits of such Mariner 10 measurements will be overcame by the more advanced and complete payload of the ESA-JAXA BepiColombo space mission, especially the joint electric field observations of LHW (with PWI/MEFISTO and PWI/WPT) and electrons (with MPPE/MEA). We are therefore confident that the Mio plasma and electric field instrumental suites of BepiColombo will soon enable to shade light on the statistical relevance of LHDI and associated electron heating in the global dynamics of Mercury's magnetosphere frontier with the solar wind.

Conclusions
In this work, we have addressed the question of electron acceleration efficiency by lower hybrid waves generated by the lower hybrid drift instability.
For this purpose, we have performed 3D, full-kinetic numerical simulations to provide numerical evidence of electrons acceleration parallel to the ambient magnetic field by resonant waveparticle interaction with LHDI waves. Our self-consistent nonlinear model has also enabled us to address the consequences of the saturation of the instability on the eventual stoppage of the electron acceleration. To our best knowledge, this is the first time that this process is self-consistently observed in full-kinetic simulations. This represents the first original contribution of this work.
Moreover, we have provided quantitative estimates of the efficiency of this resonant acceleration process. To model this process we have (i) used a standard QL model based on the work of Cairns & McMillan (2005), and (ii) developed an extended QL model that includes the consequences of nonlinear LD-like effects on long time scales evolution of the electron distribution function. We have compared the results of these two quasilinear models and the full-kinetic model. Such comparison highlighted the limitations at long time scales of the standard quasilinear theory that paved the way for an extended one, designed to overcome such limitations. Such comparison also enabled to validate the eQL model. This new extended quasilinear (eQL) model successfully captures the electron acceleration properties found on full-kinetic simulation results, and represents the second original contribution of this work.
The eQL model, thanks to its simplicity, has enabled us to explore a range of parameter space that is not accessible to fullkinetic simulations due to computational constraints. In partic- From top to bottom, we show the magnetic field module (first row), the magnetic field components (second row), the electron energy spectra as a function of time (third rows), and the electron energy spectrum at the time of the crossings (fourth row) -time corresponding to the grey vertical dashed lines in the first three panels. ular, we have addressed the efficiency of LHDI electron acceleration at Mercury's magnetopause, using a realistic proton-toelectron mass ratio. In this context, we estimate that LHDI electron acceleration is an efficient mechanism to energize electrons during periods of strong magnetospheric compressions (when the magnetopause boundary steepens on scales of the order or lower than the ion gyroradius). This efficiency strongly depends on the density gradient at the magnetopause. Under such conditions of "steep" magnetopause at Mercury, we expect strong signatures of LHDI electron acceleration in the BepiColombo instrumental suite response (PWI and MPPE/MEA).
We are confident that our extended quasilinear model will enable to quantitatively study the efficiency of electron acceleration in inhomogeneous space plasmas, in support to past and future space plasma observations. More important, this work also provides a validated framework to extend the range of validity and applicability of quasilinear modeling, not only associated to lower hybrid drift instability -as in this study -but also to a broader variety of waves and instabilities in space plasmas.