Large-scale numerical simulation of ionospheric Langmuir turbulence excited by a radio frequency electromagnetic wave

A numerical simulation is presented concerning an L/O mode electromagnetic wave propagating normally into an overdense magnetised plasma with a smooth density gradient leading to excitation of Langmuir turbulence in the vicinity of the reflection point. The simulation parameters are chosen to represent an ionospheric radio frequency heating experiment but may have relevance to other situations. The simulation model is one-dimensional for large-scale electromagnetic waves and two-dimensional for short-scale electrostatic waves. This allows consideration of local modulational and parametric-decay instabilities as well as the larger scale depletion of the driver electromagnetic wave by anomalous absorption due to the excited turbulence. Simulated growth rates are shown to be in broad agreement with expected values and the evolution of the spatial distribution of the turbulence and driver field profile are presented along with simulated scatter radar spectra.


Introduction
Ionospheric radio heating experiments performed over the past decades have utilised the ionosphere as a 'natural plasma laboratory' to study the interaction of electromagnetic waves with plasma in parameter regimes not easily accessible in laboratory experiments [1][2][3][4][5]. As transmitter powers have increased various effects have been observed including secondary electromagnetic emissions (SEE), optical emissions and more recently at the HAARP facility descending artificial ionospheric layers involving enhancement of the ionospheric plasma density [3].
In the initial phase of heating experiments it is believed that the main nonlinear plasma process is that of Langmuir turbulence excited near the reflection layer where the matching condition that the input frequency be near the local plasma frequency is met. This permits the parametric excitation of low frequency ion waves and high frequency Langmuir waves. It is also understood that at later times fieldaligned striations may form leading to upper-hybrid effects at somewhat lower altitudes [4][5][6]. The focus of this work is on the first few 10 s of milliseconds during which Langmuir turbulence may be dominant.
Langmuir turbulence is now a well-understood process [7,8] with many detailed simulations having been performed using fluid [9,10], Vlasov [11], PIC [12,13] and quasilinearkinetic [13] models. However, a common feature of many models is the use of periodic boundary conditions (sometimes with an ad hoc partially periodic/open boundary) and the assumption of a given spatially homogeneous background plasma density and driver field estimated from the initial electromagnetic wave profile. This work attempts to address this limitation with a more satisfactory electromagnetic driver model and a realistic background plasma density gradient. A one-dimensional model [14,15] is generalised to two dimensions, allowing modelling of electrostatic waves excited in directions other than the direction of the background density gradient, but still restricting the electromagnetic wave variation to be one-dimensional parallel to the background density gradient (however we note that in [15] a scheme was developed to allow for off-vertical electromagnetic incidence at a fixed angle). Large-scale two-dimensional electromagnetic simulations of ionospheric heating scenarios have been performed in other works [6,16,17] but not with the aim of simultaneously resolving electrostatic Langmuir turbulence excited by the electromagnetic driver.
Limitations of the simulation include the omission of certain plasma kinetic effects and higher order nonlinearities not handled by the modified Zakharov treatment as well as the still limited dimensionality. The coupling to lower-hybrid waves is also neglected but may be important at the upperhybrid layer on longer timescales. However, the simulation does allow the modelling of driver depletion of the electromagnetic wave which may be important at large input powers and conditions of highly developed turbulent modulation of the plasma.

Numerical method
A modified Zakharov model is used with a complex envelope representation for the high frequency variables: where Ã is a complex amplitude and ω 0 is the transmitted angular frequency. Tildes will henceforth denote complex envelope variables. A direct representation for the low frequency variables is used. The simulation domain is twodimensional with a background density gradient in the z-direction and a background magnetic field lying in the y/z plane. The ionospheric heating scenario and simulation domain are illustrated in figure 1. In the z-direction the main part of the simulation box lies between z=0 and z L , BODY = with absorbing perfectly matched layers (PMLs) of length L PML at either end of the simulation. The y-direction uses periodic boundary conditions. The generalised Zakharov model linearises most of the system retaining only the nonlinear refraction and ponderomotive force terms to allow coupling of high frequency oscillations to low frequency density perturbations. In the following sections the model equations and their spatial discretisation are presented; appendix A describes the timestepping scheme.

Electromagnetic coarse grid system
The electromagnetic part of the simulation, represented by the transverse fields E E cB cB , , , x y x y trans, trans,˜˜˜i s one-dimensional such that these variables depend only on the z-coordinate and are assumed to be spatially constant in the y-direction, 14 y y y In the boundary PMLs the equations are modified as in equations (7)- (14) where auxiliary variables FE x y , and FB x y , are introduced according to standard principles [18][19][20], with the damping factor PML s depending quadratically on the distance into the PML: where z PML is the distance into the PML, R is the design amplitude reflection coefficient and λ 0 is an estimate of the plasma-modified L/O mode electromagnetic wavelength at the start of the PML. CUT,PML w is a PML cut-off frequency as in the complex-shifted PML formulation [19] but in practice this term seems to have had little effect on the simulation. The choice of a PML as an absorbing boundary may not be the best choice in that for a dispersive plasma medium the reflection coefficient will depend on frequency and in particular waves which are close to the plasma cut-off frequencies will be poorly absorbed. Possible alternatives include ramping up plasma collisionality or adding a diffusive term but too rapid a gradient would also lead to reflections. In a 2D electromagnetic model of magnetised plasma, waves with oppositely directed phase and group velocities in the direction of the absorber can give rise to numerical instabilities, but for the 1D electromagnetic model used here this is not an issue.
A nested grid scheme is used where the electromagnetic variables, in addition to depending only on the z-coordinate, are represented on a coarser grid in the z-direction with larger spatial step than variables on the fine grid. This alleviates the restrictive Courant explicit timestepping limit that would exist were all variables represented on the same fine grid required to represent the smallest scales excited in the simulation [14] and is justified by the assumption that the electric field of short-scale waves in the simulation is primarily electrostatic and therefore longitudinally polarised. The use of a onedimensional electromagnetic model may be similarly justified. The electromagnetic variables are coupled to the fine grid variables via the current densities J x y , with local averaging performed in the y and z directions as described in section 2.4. External current densities J x y , ,drivẽ are added to launch the input electromagnetic wave, with the source placed slightly in front of the bottom PML zone, spread out over 20 coarse gridpoints (about 1/3 wavelength) with a smoothed profile, and J x y , ,drivẽ chosen to give an input L/O mode with the desired amplitude.
is used, to which is added the slow timescale density perturbation n . slow This is the only nonlinearity included in the fast timescale system, allowing high frequency waves to scatter off low frequency density perturbations.
The fast timescale electrostatic potential is defined by: The Von-Neumann boundary condition in equation (27) used at the top of the simulation is appropriate for the 0th y-Fourier component of the electrostatic potential but is less justified for the nonzero Fourier modes; a condition based on the ansatz k y k z exp i y y static f~-( | | )might be more appropriate but no ill-effects have been observed using the Von-Neumann condition as stated.
The slow timescale response is represented by the variables n slow and n slow  as follows: n t n 28  and is driven by a ponderomotive force term calculated from the high frequency current density. The response is modelled as a single mode ion acoustic wave, assuming the electrons to have the same number density as the ions on the slow timescale, with an ion Landau damping operator iL n used as described in appendix B which makes the ion acoustic response heavily damped. It is known that this is a limited model of the low frequency response due to more complex ion kinetic effects [21]. The values of n slow were initialised with small (of order 4 × 10 −5 n 0 ) random numbers to seed the instabilities in the simulation.

Spatial discretization
On the fine grid a Fourier representation in the y-direction and a finite-difference representation in the z-direction are used. A staggered grid scheme is used where all variables except for J z are regarded as being located on z-gridpoints, numbered by a z-index i and Fourier y-index j, while J z is located in between z-gridpoints, numbered by a half-integer index i + 1/2. All z-derivatives are approximated using centred stencils as in equation (33) whilst the y-derivatives are computed using exact Fourier differentiation as in equation (36). Where required for the magnetic field term and the evaluation of the nonlinear terms, variables are transferred between the staggered grids using spatial averaging as in equation (35). The solution is stored primarily in Fourier y-space but is transformed to and from real y-space to evaluate the nonlinear terms. and FE x y , are stored on half-integer coarse gridpoints. Derivatives on the coarse grid are evaluated in a manner analogous to equation (33).

Coarse/fine grid coupling
A fixed ratio between the coarse and fine grids is used in the z-direction. The fine gridpoints are regarded as being located at positions given by: z z , coarse + ¢ where z coarse is the location of a coarse grid point and z is the ratio between the coarse and fine grid sizes. To transfer E trans from coarse to fine, it is assumed that the coarse grid variables represent variation that is independent of y, so that they are assigned to the 0th y-Fourier mode. To achieve interpolation in the z-direction cubic b-splines are used, such that a strict interpolation is not used but some smoothing is applied in the process. To perform the averaging of J x y , from fine to coarse the transpose of the interpolation operator is used with an appropriate weighting factor. At the simulation boundaries the stencils are not modified but it is simply assumed that all variables outside the computational domain are 0. This has little effect since the variables are strongly damped in the PML zones.

Boundary conditions
In addition to the PML zones described in section 2.1 boundary conditions must be specified at the ends of the PMLs. Ideally if the PMLs perform well this should not be a major issue, and Dirichlet (perfect conductor) conditions E 0 x y trans, , = at z L PML =and z L L BODY PML = + are used on the coarse grid. On the fine grid Dirichlet conditions are used for all variables except for the electrostatic potential as described in section 2.2. In the context of the staggering scheme this means that a Dirichlet condition is obtained for J z at the bottom boundary and for n fast at the top boundary, but this is modified so that a Dirichlet condition on J z at both boundaries is obtained with a view to keeping the total high frequency charge density constant.
In practice these simple boundary conditions do allow reflections of short-scale waves and low amplitude spatial oscillations accumulate near the boundaries as the simulation proceeds, particularly in n fast and J .
The electron collision frequency was ramped up in the PML zones as in equation (39) in an attempt to mitigate this effect, but the increase in the collision frequency was probably not chosen large enough to make much difference. One may speculate that either the extra collisionality could be increased or perhaps a variable spatial diffusion could be added to try to damp short-scale waves reaching the boundary, with care taken not to introduce further reflections by too fast a gradient. Formally correct outflow boundary conditions instead of a 'sponge' region are challenging to design for this problem due to the magnetised 2D plasma behaviour. In any case the amplitude of the fine grid boundary oscillations has been deemed low enough to ignore so this problem has not been addressed more carefully.

Simulation parameters
The simulation is designed to represent an ionospheric heating experiment by simulating a stretch of plasma at the L/O mode cut-off near the F2 peak of the Ionosphere, typically located between 200 and 300 km above the Earth's surface. Parameters relevant to experiments at the HAARP facility have been chosen as displayed in table 1. The choice of the simulation domain size L BODY is a trade-off between managing computer resources and ensuring a realistic simulation which models the large difference in spatial scale between the density profile and the free-space electromagnetic wavelength. In practice for many experiments the ratio can be significantly larger than that simulated in this work, but we believe that the length simulated is still large enough to be of relevance. The width of the PML layer L PML was chosen large enough to ensure good absorption of electromagnetic waves without excessive computational cost. Note that the input power flux represents an equivalent radiated power (ERP) of only 35 MW assuming that the distance between the transmitter and F2 peak is around 250 km, which is considerably less than the maximum ERP of HAARP which is several GW. Even with this reduced level of input power plasma density perturbations ultimately exceeding the validity of the present model are observed in the cavitation region near the incident O-mode cut-off with density depletions of order 50% observed in collapsing wavepackets. As described in appendix B, the numerical electron Landau damping operator had to be enhanced to prevent the occurrence of negative densities.
The driver current density is applied to a small region of the simulation in front of the bottom PML, with polarisation calculated for an input L/O mode with wavevector in the z-direction: The driver current density is determined from the wave magnetic field amplitudes as J Ac z 2 where the delta function is numerically approximated as a localised function with a width of 20 coarse gridpoints and A is a correction factor of order unity. The source launches waves in both directions but the bottom PML absorbs the downwards wave. The amplitude of the source is ramped up from 0 to its target value with a characteristic time of 40 input wave periods.
Numerical parameters are presented in table 2. A longer spatial step is used for the y-coordinate because electrostatic waves which are excited in the simulation are directed mainly parallel to the background magnetic field (which points nearly perpendicular to the y-direction) and thus have a relatively long wavelength in the y-direction.

Simulation results
This section contains an analysis of the simulation predictions with respect to the evolution of the electromagnetic wave, the small-scale electrostatic turbulence and simulated radar spectra.

Large-scale electromagnetic field evolution
At the early stages of the simulation the electromagnetic wave reflects from the plasma density gradient at the O-mode cutoff and forms an Airy-like standing wave pattern displayed in figures 2(a) and (b) [22]. In figure 2 the electric fields have been smoothed in space to filter out the short-scale electrostatic turbulence and extract the electromagnetic waves. Near to the cut-off the wave is polarised almost parallel to the background magnetic field. Following the development of turbulence and resulting anomalous absorption the driver profile is depleted as shown in figures 2(c) and (d). The initial peak field of E  figure 3, and there is consequently a small amplitude Z-mode wave seen above z 2200 m = in figures 2(c) and (d) which may propagate to the top boundary and be absorbed by the PML located there. Figure 4 displays the z-component of the electromagnetic at the start of the simulation and following the development of turbulence. The decomposition into wave modes was performed by projecting the amplitudes of the 4 transverse electromagnetic variables E E cB , , x y x˜˜a nd cB ỹ onto a basis formed by the 4 cold plasma electromagnetic wave modes with frequency 0 w and wavevector in the z-direction. In this we assume that most of the coarse grid wave energy is in fact at frequency 0 w (which is 0 frequency in the envelope representation) which will not be a good approximation after the development of strongly nonlinear effects, and furthermore the decomposition into different wave modes will be invalid near cut-offs and resonances even for monochromatic waves. At early times there is approximately 90% reflection; the 10% absorption is due to the electron collision term. Following the development of turbulence there is nearly total absorption of the incident wave, and this is also visible in figure 2 where the electromagnetic wave pattern transforms from a near perfect standing wave to a propagating wave pattern. One sees in figure 4(c) that at time 4 ms the anomalous absorption is highly concentrated near the cut-off altitude at z 1860 m.
= Figure 5   Downwards propagating O-mode waves with frequency above this cut-off will propagate into less dense plasma and can potentially escape into vacuum. In figure 5(a) one sees a tightly peaked feature around the incident driver wave at f 0 kHz D = which may correspond to waves generated at higher altitudes due to turbulence which have then propagated downwards to this measurement point in the O-mode.
Signals with frequencies below the local O-mode cut-off cannot couple into the escaping O-mode but can potentially propagate in the Z-mode, which in the cold plasma approximation has a resonance at f 52 kHz. D = -However, as this mode propagates downwards it will eventually reach a point where its frequency equals the local resonance frequency and it should then become an oblique upper-hybrid/Langmuir wave and be absorbed by Landau damping. In figure 5(a) one may identify a small amplitude downshifted spectrum in the local Z-mode although the power is weak when the logarithmic scale is considered. One must consider that at the measurement point there may be short-scale electrostatic waves present which might confuse the interpretation of the spectrum in terms of the cold plasma dispersion relation (in a full treatment the Z-mode becomes an oblique upper-hybrid/ Langmuir wave above its cold plasma resonance frequency), but by analysing E y,trans this should have been mitigated to some extent.
If one assumes that the downwards L/O mode signal at the bottom of the simulation is a reasonable measure of the radiation that would escape from the ionosphere during a heating experiment then one may extract an SEE spectrum as in figure 6. This is a highly simplistic view due to the restriction to 1D vertical up/down propagation imposed by the simulation and the possibility of contamination by other wave modes The Z-mode (or slow X mode) has its electrostatic resonance at z=1346 m based on the 194°angle to the background magnetic field. The R/X mode cut-off occurs below the bottom boundary of the simulation and the Z-mode cut-off above the top boundary for the density profile used. which cannot escape from the ionospheric plasma into vacuum. The spectrum is displayed when measured at two different heights and it is seen that the curves overlap almost exactly, suggesting that most of the secondary O-mode radiation comes from turbulence in a small region near the incident O-mode cutoff altitude. Whilst local signals with large frequency downshifts may be present in regions exhibiting SLT only those in a small frequency band around the local plasma frequency will be able to escape downwards in the O-mode, whilst all upshifted signals should be able to escape [23].

Evolution of electrostatic turbulence
In the initial stage of the simulation exponential growth of electrostatic waves occurs due to parametric instabilities driven by the electromagnetic driver wave. The strongest instability is a purely-growing modulational instability (historically also known as an oscillating two stream instability) occurring at the final standing wave maximum where the driver amplitude is strongest and the driver frequency is very close to the local plasma frequency. Figures 7-9 display the evolution of the slow density perturbation evolving between 0.5 and 1.5 ms near z=1800 m. Waves are mostly excited along the direction of the background magnetic field and the simulation growth rate displayed in figure 9(a) is seen to be consistent with a theoretical calculation using local plasma parameters displayed in figure 9(b) (this calculation is outlined in appendix C). The simulation wavenumber resolution in the y-direction is limited due to a small domain size in the y-direction. Of note in figure 9 is that the growth rate parallel to the z-direction is 0 because the region is located above the Z-mode resonance at z=1346 m for propagation in the z-direction such that the magnetic field effectively cuts off Langmuir/upper-hybrid electrostatic wave propagation in the z-direction for wave frequency near the driver frequency (for an unmagnetized plasma a much broader growth characteristic would be expected). This feature has the important consequence that a 1D simulation in which waves are restricted to be directed in the vertical direction might not exhibit growth of electrostatic waves at the highest standing wave maximum and the turbulence may be non-physically restricted to lower altitudes. However a larger input amplitude could overwhelm this effect.   At lower altitudes with greater mismatch between input wave frequency and local plasma frequency the parametricdecay instability (PDI) is dominant [24], and figures 10 and 11 show the slow density perturbation around z=920 m where the PDI is seen to occur, but with a lower theoretical growth rate than the modulational instability at higher altitudes (a simulation growth rate was difficult to extract). The PDI corresponds to the inner ring in figures 9(b) and 11(b) and the outer ring is the modulational instability which in this region has a very narrow region of unstable wavenumbers as seen in figure 11(b). For the PDI case the density perturbation oscillates as well as grows in time.
After some time strong Langmuir turbulence (SLT) due to the modulational instability develops near the cut-off of the input electromagnetic wave, with very large amplitude density depletions formed in spatially localised cavitons or wavepackets. Figure 12 shows density perturbations representative of SLT near the largest initial standing wave pattern maximum. The spatial extent of the cavitons is seen to be much larger perpendicular to the magnetic field than parallel, and there is a tendency for 'strings' of wavepackets aligned along the magnetic field to form which has been observed in previous simulations [9]. Note that the depth of the density depletions created exceeds the model validity (which requires n n slow 0  ) so only qualitative details should be inferred    from the simulation in this region after the development of the strong turbulence.
The temporal and spatial development of the intensity of the turbulence is displayed in figure 13. Initially the strongest electrostatic turbulence in the simulation is located near the cut-off of the input electromagnetic wave, but one sees in figure 13 that the turbulent electrostatic energy density decreases slightly at large z and that turbulence appears at lower altitudes as the simulation progresses. Initially the imprint of the standing wave pattern is seen in the distribution of the turbulence, but at later times this is smoothed out, which may be due to the increased absorption of the driver wave destroying the standing wave pattern or due to the propagation of turbulence between different regions.

Simulated scatter radar spectra
Radars with frequency in the range 100 MHz-1 GHz have been employed in ionospheric heating experiments to measure the spectrum of density perturbations excited by the separate driver wave [25][26][27][28]. Some radars are sensitive enough to measure scattering from the natural ionosphere allowing measurement of electron temperature and density while others are only sensitive enough to measure the enhanced scatter spectrum due to the excited turbulence. The experimental methods vary but we here consider a radar pulse launched and scattered vertically, consistent with the spatial domain used in the simulation. With a frequency much higher than the plasma frequency a radar signal will have a well-defined spatial wavevector k c R R w = and we can consider backscattering off electron density fluctuations having frequency scat w and spatial wavevector To extract a representative spectrum from the simulation we convolve the values of n slow and n fast with a spatial kernel k z exp 2i scat -( ) / where a Hann window is used to represent localisation to different height ranges. Figure 14 displays the spectrum extracted from the ion timescale fluctuation n .
slow The lowest height region in figure 14(a) corresponds to a region approximately at the primary-decay 'matching height' of the radar where the local plasma frequency is such that k k IA scat L UH scat 0 [25], allowing resonant parametric decay of the heating wave along the line of sight of the radar into oblique upper-hybrid/ Langmuir waves and ion acoustic waves with wavevectors equal to k . scat  This gives rise to distinct peaks in the return spectrum at k Ck IA scat s scat w =  ( ) corresponding to scattering off upwards and downwards propagating ion acoustic waves (or more correctly, waves with phase velocities in the upwards and downwards directions). In figure 14(c) a different 'filled-in' spectrum is seen known as the 'caviton continuum' where it is understood that the SLT present at this higher altitude is associated with a much broader spectrum. Due to the extreme density fluctuations in the simulation in the SLT zone the analysis of this region is questionable. Figure 14

Discussion
A computer model has been developed to study important aspects of the propagation, reflection and absorption of powerful electromagnetic waves as they interact with plasma. The model assumes that the EM signals vary only in the vertical direction, but that they may excite a range of small-scale electrostatic plasma oscillations which are represented in two dimensions. By representing the electromagnetic degrees of freedom on a coarser grid than the plasma fluid degrees of freedom a longer timestep can be taken in light of the inequality v c, Te  allowing considerable computational speedup.
The model has been applied to consider the behaviour of 3.2 MHz EM waves with an intensity of 45 Wm 2 m polarised in the L/O mode and propagating into cut-off at an angle of 14°to the magnetic axis, a scenario relevant to ionospheric heating experiments at facilities like HAARP, EISCAT and SURA. The simulation was conducted over a range of some 1.9 km below the cut-off altitude for the incident L/O wave, and for about 700 m above the cut-off. Initially the incident wave experiences strong reflection and forms a standing wave pattern with a peak field of 3 Vm 1 polarised along the direction of the background magnetic field at the final standing wave peak. The wave turbulence initially develops in accordance with parametric instability theory with growth rates in agreement with those given by asserting the driver amplitude based on the initial standing wave pattern and using the local plasma density in the calculation. At high altitudes with small detuning between the incident wave frequency and the local plasma frequency the modulational/OTSI instability is favoured with theoretical growth rates of order 6000 s −1 but at lower altitudes with larger detuning the PDI is favoured. After about 3 ms SLT develops at high altitudes from the nonlinear phase of the modulational instability giving almost total absorption of the incident wave. The depletion of the input driver wave then leads to a nontrivial distribution of turbulence in the domain which is a new development in understanding the problem as a whole. The simulation may need to be run for longer to obtain a quasi-steady state.
Predictions of the scatter radar spectra one might expect from this perturbed plasma are in broad agreement with the established picture in the literature, with peaked spectra Figure 13. RMS electrostatic electric field amplitude at various times. A high-pass spatial filter was first applied to E z to attenuate the electromagnetic driver wave before evaluating the local squared electric field amplitude and then averaging in the y-direction and smoothing in z before taking the square root. predicted at the low altitude 'matching height' where waves with wavevector equal to the scattering wavevector are resonantly excited by the PDI. Broader spectra are predicted at high altitudes due to SLT wavepackets with a broad k-spectrum. For the mock radar frequency of 167 MHz the ion acoustic sidebands are peaked at a shift of 2.2 kHz.
Analysis of the returned EM signal (performed over a 3-10 ms time range) predicts a primarily downshifted spectrum (having width around 10 kHz for a 30 dB attenutation with respect to the returned driver signal) which appears to be generated at high altitudes before propagating down through the plasma with little modification.
The still limited dimensionality in the simulation may hide certain effects involving the three-dimensional nature of collapsing wavepackets in SLT. The restriction to electromagnetic wave propagation parallel to the background density gradient only allows vertical incidence to be considered and additional effects involving the input electromagnetic beam size might be revealed with a more comprehensive treatment. Strictly speaking the input wave does not 'propagate' in the z-direction as its group velocity is directed in a slightly different direction to its phase velocity which is in the z-direction, but we have ignored this distinction when discussing the 1D electromagnetic model.
Computationally, a patching scheme might be useful when large regions of statistically homogeneous turbulence are created. In such a scheme an anomalous resistivity could be calculated from local simulations of electrostatic turbulence and applied to a global electromagnetic model.
Usually with Zakharov-type simulations a manipulation and approximation is made (neglecting t 2 2 ¶ ¶ compared with t i 0 w - ¶ ¶ ) such that each wave mode is represented by just one complex field [9,10], but in this work every real field has simply been replaced by a complex one on the fast timescale (additionally 4 variables J x y z , , and n fast are used on the fast timescale while a reduced model might require fewer). One advantage of using separate high and low frequency systems is that the electron response can be modelled differently for each. Complex envelope representations can also be useful if they allow longer timesteps to be taken [29], but here this has not been implemented. For the timestep of 8 ns used there are 39 timesteps per wave period which would be sufficient to resolve the waves in full. The timestep is dictated by the electromagnetic Courant stability limit on the coarse grid and the electron Landau damping operator. If implicit schemes were used [16,29] then it should be possible to eliminate the need for the nested grid method (with the cost of extra storage of field variables on the fine grid and the implicit solves). The simplicity of the plasma model used is a major drawback, with no consideration for modification of the electron distribution function by wave-particle interactions, or higher order nonlinearities present in deep density depletions. The correct form of the dissipation and/or extra nonlinear effects required to resolve the collapsing density depletions has not been investigated in this work. Previous work [8,10] has emphasised the importance of transit-time damping of localised coherent electrostatic wave structures instead of ordinary random-phase Landau damping. It would perhaps be more reasonable to restrict the simulation method to smaller input amplitudes to avoid the extreme nonlinearities generated in the SLT zone.
Additionally, the combined effect of the magnetic field and plasma kinetic behaviour is not considered, with Landau damping obliquely to the magnetic field being modelled in the same fashion as in the parallel direction. Similarly the ion response would benefit from a more detailed model.

Conclusion
The two-dimensional simulation presented shows the excitation of Langmuir turbulence by an electromagnetic driver wave injected from the bottom side. The simulation parameters are representative of ionospheric heating experiments and results of electromagnetic emissions and radar scattering agree with theoretical expectations. The two-dimensional model allows for excitation of electrostatic waves in directions other than the background density gradient. By = -simulating both electromagnetic waves and electrostatic waves simultaneously it is possible to obtain a more complete and self-consistent description. One expects this method to be most relevant at high driver powers when the turbulence is sufficiently intense to interact nontrivially with the driver wave as seen in figures 2 and 13. By varying such parameters as the input wave strength, the magnetic field angle and the plasma density gradient it could be possible to explore a range of different effects and perhaps adapt the method to other problems where a large-scale driver transfers energy to much shorter scale turbulence as well as to lower-hybrid turbulence at large angles to the magnetic field lines.
Historically numerical Landau damping operators have been applied as multiplications in Fourier space [9]. In this work the solution is available in y-Fourier space but in the z-direction a different approach is required. An ad hoc scheme has been developed where an operator F k F k * = -( ) ( ) defined in z-Fourier space is applied on a grid in real z-space using a finite difference scheme as: where it is assumed that the operator is symmetric about k=0, and 2N−1 is the number of points included in the stencil. To determine the stencil coefficients j a the cost function C: representing the error in the Fourier approximation in wavenumber space is minimised with respect to . j a In practice H k ( ) corresponds simply to a Fourier cosine series taken over the range of wavenumbers represented on the simulation grid. For both the ion and electron Landau damping operators N=41 is used in the finite difference stencil.
The ion Landau damping operator is written as: , which was used in the simulation, is consistent with previous work [9]. | | ) | | ( ) / It is known that the ion response obtained in this manner has limitations [21], and the effect of the magnetic field on the electron Landau damping operator has been ignored. More seriously however it is seen in figure B2 that the extra damping added to the electron Landau damping operator dominates over the thermal part. The numerical damping operator therefore only qualitatively represents the physically expected high k cut-off given by electron Landau damping and as mentioned earlier it does not capture the true form of the dissipation of collapsing wavepackets.
Fitting over the entire sampled k-space interval up to the Nyquist wavenumber with uniform weighting is unlikely to be an optimal strategy, and the stencils with N = 41 are computationally expensive to evaluate. It is necessary to use a large number of points to ensure that the numerical damping added at small k is not excessive as displayed in figure B2(b). The method involves Fourier analysis of functions with discontinuities which is a questionable procedure. When evaluating the convolutions, boundaries in the simulation are treated by padding with zeros.