Instability-enhanced transport in low temperature magnetized plasma

TL;DR

This paper demonstrates that in low temperature, collisional magnetized plasma, instabilities at high magnetic fields significantly enhance transport, overriding magnetic confinement and making transport independent of magnetic field strength.

Contribution

It introduces an analytical model for instability-enhanced collision frequency based on a magnetic field independent density ratio, revealing the mechanism behind transport enhancement.

Findings

Transport becomes independent of magnetic field at high magnetization.

Instabilities destroy electron confinement in highly magnetized plasma.

Analytical expression for instability-enhanced collision frequency is proposed.

Abstract

It is shown that the transport in low temperature, collisional, bounded plasma is enhanced by instabilities at high magnetic field. While the magnetic field confines the electrons in a stable plasma, the instability completely destroys the confinement such that the transport becomes independent of the magnetic field in the highly magnetized limit. An analytical expression of the instability-enhanced collision frequency is proposed, based on a magnetic field independent edge-to-center density ratio.

Tables

Table 1. Table 1: Numerical quantities of the simulation used to solve the dispersion relation.

Magnetic field	20	mT
Pressure	3	mTorr
Plasma density	$1.8\times {10}^{16}$	m-3
Electron temperature	3.84	eV
Diamagnetic drift	17.3	km/s
$E\times B$ drift	-4.3	km/s

Equations

$$\frac{d\mathbf{v}}{dt}=\nabla \varphi -\mathbf{v}\times \mathbf{b}-\frac{\nabla n}{n}-\nu \mathbf{v}$$

$$v_y$$

$$\partial n/\partial t+\nabla \cdot (n\mathbf{v})=n{\nu }_{iz}$$

$$v_x^{\prime }+\frac{n^{\prime }}{n}{v}_x={\nu }_{iz}.$$

$$-{\varphi }^{\prime }=-\frac{n^{\prime }}{n}-v_x\left(\frac{1}{\nu }+\nu \right)$$

$$\frac{v_x^2}{2\mu }=-\varphi \Rightarrow \frac{v_x{v}_x^{\prime }}{\mu }=-{\varphi }^{\prime }.$$

$$\left(1-\frac{v_x^2}{\mu }\right)v_x^{\prime }={\nu }_{iz}+\alpha v_x^2.$$

$$(\alpha \mu +{\nu }_{iz})\arctan \left(\frac{v_x{\alpha }^{\frac{1}{2}}}{{\nu }_{iz}^{\frac{1}{2}}}\right)-(\alpha {\nu }_{iz}{)}^{\frac{1}{2}}{v}_x=\mu {\alpha }^{\frac{3}{2}}{\nu }_{iz}^{\frac{1}{2}}x$$

$$v_x={\mu }^{1/2}$$

$$\left[{\left(\frac{\alpha \mu }{{\nu }_{iz}}\right)}^{\frac{1}{2}}+{\left(\frac{{\nu }_{iz}}{\alpha \mu }\right)}^{\frac{1}{2}}\right]\arctan \left[{\left(\frac{\alpha \mu }{{\nu }_{iz}}\right)}^{\frac{1}{2}}\right]-1=\frac{{\mu }^{\frac{1}{2}}\alpha l}{2}$$

$$-{\varphi }^{\prime }=\frac{v_x\left({\nu }_{iz}+\alpha v_x^2\right)}{\mu -v_x^2}$$

$$\frac{n}{n_0}={\left(1+\frac{\alpha v_x^2}{{\nu }_{iz}}\right)}^{-\frac{1}{2}\left(1+\frac{{\nu }_{iz}}{\alpha \mu }\right)}$$

$${\nu }_{iz}=\frac{{\pi }^2}{\alpha l^2}$$

$$v_x={\left(\frac{{\nu }_{iz}}{\alpha }\right)}^{\frac{1}{2}}\tan \left({\nu }_{iz}^{\frac{1}{2}}{\alpha }^{\frac{1}{2}}x\right)\approx \frac{\pi }{\alpha l}\tan \left(\frac{\pi x}{l}\right)$$

$$\frac{{\pi }^2/(\alpha l{)}^2+v_x^2}{\mu -v_x^2}\underset{\alpha \to \infty }{⟶}0$$

$$v_x^{\prime }={\nu }_{iz}/{\cos }^2\left({\nu }_{iz}^{\frac{1}{2}}{\alpha }^{\frac{1}{2}}x\right)=\alpha v_x^2/{\sin }^2\left({\nu }_{iz}^{\frac{1}{2}}{\alpha }^{\frac{1}{2}}x\right)$$

$$-{\varphi }^{\prime }=\alpha {\mu }^{1/2}$$

$$-n^{\prime }/n=2\alpha {\mu }^{1/2}$$

$${\nu }_{iz}=\pi h{\mu }^{1/2}/l$$

$$h=\frac{\pi }{\alpha l{\mu }^{1/2}}$$

$${\nu }_B=h{\mu }^{1/2}l/\pi .$$

$${\nu }_B=\frac{h{m}_e{\omega }_{ce}^{2}L}{\pi (m_i{k}_B{T}_e{)}^{1/2}}.$$

$${\chi }_e\equiv \frac{{\omega }_p^2}{k^2}\frac{n_{e1}}{{\varphi }_1}=\frac{{\omega }_p^2}{k^2}\frac{{\omega }_{\ast }+k^2(\omega +{\omega }_0+i\nu )}{(\omega +{\omega }_0)+k^2(\omega +{\omega }_0+i\nu )}$$

$${\omega }_0=-k{\varphi }^{\prime },{\omega }_{\ast }=-k{n}^{\prime }/n$$

$${\chi }_i\equiv -\frac{{\omega }_p^2}{k^2}\frac{n_{i1}}{{\varphi }_1}=-\frac{\mu {\omega }_p^2}{{\omega }^2}.$$

$$1+{\chi }_e+{\chi }_i=0$$

$$P_0(\omega )=$$

$$+{\omega }^2{\omega }_p^2{\omega }_{\ast }/k^2$$

$$Q(\omega )=$$

$$P(\omega )=P_0(\omega )+i\nu Q(\omega )=0.$$

$$P({\omega }_r+i\nu {\omega }_i)=0\iff {\omega }_i=-Q({\omega }_r)/P_0^{\prime }({\omega }_r)$$

Full text

††preprint: AIP/123-QED

In low temperature plasma discharges, the transport from the ionization region to the walls determines the electron temperature and hence the global discharge properties. The classical drift-diffusion theory of low density magnetized plasma transport Simon (1955); Fruchtman et al. (2005); Sternberg et al. (2006) predicts that the plasma should reach a steady-state equilibrium qualitatively analogous to that of a non-magnetized discharge Tonks and Langmuir (1929); Schottky (1924); Godyak (1986) and that the confinement increases with the magnetic field. However, when the value of the magnetic field is high enough, strong instabilities develop that deconfine the electrons, which enhance the macroscopic transport Curreli and Chen (2014). These so-called universal instabilities develop in bounded plasma, without any external source of energy. The wave energy comes from the natural drift-diffusion motion of the particles.Chen (1984) They can be either Simon-Hoh instabilities (similar to a Rayleigh-Taylor instability with an electrostatic potential),Boeuf and Chaudhury (2013) or electron drift instabilities (comparable to Kelvin-Helmholtz instabilities).Kent et al. (1969) Finally, collisions can also destabilize some particular modes, such as the lower hybrid instability.Smolyakov et al. (2017)

In this Letter, a one-dimensional (1D) model of the cross field equilibrium plasma transport is developed. The model is linearized including time-dependence in order to describe the linear growth of the instability. An electrostatic particle-in-cell / Monte Carlo collision (PIC / MCC) code Birdsall and Langdon (2004) LPPic Croes (2017); Croes et al. (2017a, b); Lucken et al. (2018, 2017) is used to simulate a magnetized plasma column in two-dimensional (2D) Cartesian coordinates. The spectral analysis of the PIC simulation results is compared with the solution of the theoretical dispersion relation, and the effective collision frequency measured in the PIC simulation is compared with the analytical model.

The results are of importance for the modeling of industrial plasma devices such as neutral beam injectors,Hagelaar and Oudini (2011) and electric space propulsion systems.Thomas et al. (2016)

In this Letter, the frequencies are normalized to the electron cyclotron frequency $\omega_{ce}=eB/m_{e}$, where $B$ is a uniform magnetic field, the length dimensions are normalized to the Larmor radius of thermal electrons $\rho_{L}=\frac{1}{eB}(k_{B}T_{e}m_{e})^{1/2}$, where $T_{e}$ is the electron temperature, that is assumed uniform, the velocities are normalized to the electron thermal velocity $v_{Te}=(k_{B}T_{e}/m_{e})^{1/2}$, and the electrostatic potential $\phi$ is normalized to the electron temperature in volts. Under these conditions, the electron momentum balance equation is:

[TABLE]

where $d/dt$ is the Lagrangian time derivative, $\mathbf{v}$ is the electron fluid velocity, $n$ is the electron density, $\nu$ is the electron collision frequency, and $\mathbf{b}$ is the unit vector in the direction of the magnetic field. It is assumed that the steady-state plasma properties depend only on the variable $x$ and that the magnetic field is along $z$ ($\mathbf{b}=\mathbf{e}_{z}$). We investigate a symmetrical single-ion plasma where the plasma density drops monotonically from the discharge center at $x=0$ to a floating sheath region, near $x=l/2$. The ions have a mass $m_{i}$, their temperature is much lower than the electron temperature, such that ion diffusion is neglected compared to the mobility driven flux, and the plasma is quasineutral.

At steady-state, and making the common assumption that electron inertia is negligible due to the small electron mass Chabert and Braithwaite (2011); Lieberman and Lichtenberg (2005); Allen and Holgate (2016); Riemann (1991); Sternberg et al. (2006); Fruchtman et al. (2005), Eq. (1) is

[TABLE]

In Eq. (2), the first term of the right-hand side represents the diamagnetic drift, the second term is the $E\times B$ drift, and the third term is the drift associated with collisions.

The normalized electron continuity equation is

[TABLE]

where $\nu_{iz}$ is a normalized electron impact ionization frequency. The continuity equation is the same for ions and leads to equal electron and ion fluxes along $x$, under the assumption that the fluid velocity is zero for both species at the discharge center. This allows to use the same fluid velocity $v_{x}$ for both species. We assume that the ion Larmor radius is larger than the discharge dimensions so the ions are not magnetized. In 1D, Eq. (4) becomes

[TABLE]

Eliminating $v_{y}$ from Eqs. (2) and (3), we obtain

[TABLE]

A very simple collisionless momentum conservation equation is used here for the ions:iiiIon collisions could be included here, with a discussion on how to compute the ion momentum transfer collision frequency, Godyak (1986); Lucken et al. (2018) but with no influence on the strongly magnetized regime.

[TABLE]

where $\mu=m_{e}/m_{i}$. The electric field and the pressure gradient terms are eliminated using Eqs. (7) and (5) respectively, and inserted in Eq. (6), which yields

[TABLE]

where $\alpha=1/\nu+\nu$. This first order differential equation is integrated to obtain

[TABLE]

The quasineutral assumption, that is valid in the plasma bulk, breaks at the sheath edge. The boundary condition for the quasineutral plasma is chosen to be the Bohm criterion Riemann (1991)

[TABLE]

that needs to be satisfied at the sheath edge $x=l/2$, even for magnetized plasmas Allen and Holgate (2016); Yankun et al. (2011). Equation 9 yields the equation for the electron temperature

[TABLE]

that can easily be solved numerically for $\nu_{iz}(T_{e})$. From Eq. (7), the electric field is derived:

[TABLE]

and so is the plasma density profile from Eq. (5)

[TABLE]

We now consider the high magnetic field limit in which $\alpha,x,l\gg 1$, for $x/l\neq 0$. Eq. (11) yields

[TABLE]

such that Eq. (9) is also

[TABLE]

At a given relative position $x/l<1/2$, $v_{x}\underset{\alpha\rightarrow\infty}{\longrightarrow}0$.

The ratio between the electric field and the velocity term in Eq. (6) is hence

[TABLE]

everywhere in the plasma bulk ($v_{x}<\mu^{1/2}$). This means that at high magnetic, the electron drift is dominated by the diamagnetic drift, except at the sheath edge, as will be shown below.

In the bulk, Eq. (15) yields

[TABLE]

Using Eq. (7), the electric field has a finite value at the sheath edge

[TABLE]

This is used in Eq. (6) to estimate the diamagnetic drift term at the sheath edge

[TABLE]

Hence, for $\alpha>1/(2\mu^{1/2})$, the diamagnetic drift is larger than the electron thermal velocity. The diamagnetic drift is a purely fluid drift term and should typically not be higher than the thermal velocity. This intuition was confirmed by the PIC simulation and gives a hint that $\alpha$ cannot be arbitrarily large due to the magnetic field. If this is the case, then there should be an effective electron momentum transfer collision frequency $\nu_{B}$, that depends on the magnetic field, and that should satisfy roughly (in normalized dimensions) $\nu_{B}>2\mu^{1/2}$, such that with $\alpha=1/\nu_{B}$, the right hand side of Eq. (19) is less than unity. In general, this effective collision frequency may depend on $x$. However, it will be shown by the PIC simulation that the uniform effective collision model provides relevant corrections to the classical description.

In a 1D system, the $h$ factor that characterizes the ion losses to the walls is $h=n_{s}/n_{0}$, where $n_{s}$ is the plasma density at the sheath edge. Integrating Eq. (5) leads to

[TABLE]

where

[TABLE]

because of Eq. (13). At high magnetic field,

[TABLE]

In dimensional units,

[TABLE]

where $L=\rho_{L}l$ is the denormalized system size. At this point, both the $h$ factor and the effective collision frequency may depend on the magnetic field, which is consistent with instability enhanced transport regimes.

We now examine the plasma unstable behavior. It is assumed that a steady-state solution is perturbed by a harmonic wave propagating in the $y$ direction such that the first order densities $n_{e1}$ and $n_{i1}$ for electrons and ions respectively, and the first order potential $\phi_{1}$ are proportional to $\exp(-i\omega t+iky)$, where $\omega$ is the complex wave frequency (normalized to $\omega_{ce}$), and $k$ is the wavenumber (normalized to $\rho_{L}^{-1}$). Collisions were added to the classical inhomogeneous fluid plasma theoryMikhailovskii (1974) (see Appendix), and the perturbed electron density response to the potential is given by the susceptibility

[TABLE]

where

[TABLE]

and $\omega_{p}^{2}=n_{0}e^{2}/(\epsilon_{0}m_{e}\omega_{ce}^{2})$ is the normalized electron plasma frequency. Equation (24) is also introduced by the gyroviscosity formalism in Smolyakov et al..Smolyakov et al. (2017) The cold, collisionless, non-magnetized ions do not drift in the $y$ direction. Their susceptibility is Bellan (2006)

[TABLE]

The first order Poisson’s equation Bellan (2006)

[TABLE]

provides the wave dispersion relation. Introducing the polynomials

[TABLE]

the dispersion relation is

[TABLE]

Equation (30) is valid as long as the frequency of the instability is smaller than the electron cyclotron harmonics, where a kinetic description of the Bernstein modes is required.Allis et al. (2003); Schmidt (1966); Bellan (2006); Lafleur et al. (2016) Inspecting the following table, $P_{0}$ has always two negative roots and one positive root:

[TABLE]

Let $\omega_{r}$ be the positive root. Since $\nu\ll 1$, collisions can be treated as a perturbation term. The perturbed solution being $\omega_{r}+i\nu\omega_{i}$, to the first order in $\nu$,

[TABLE]

Since $P_{0}^{\prime}(\omega_{r})>0$, the mode stability is determined by the sign of $Q(\omega_{r})$. If $\omega_{*}=\omega_{0}=0$, the positive root is

[TABLE]

a mode transiting from the lower hybrid frequency at low $k$’s to the ion plasma frequency at high $k$’s, and damped by collisions. If $\omega_{0}=0$, but $\omega_{*}\neq 0$, the solution is

[TABLE]

The general instability criterion is found by solving jointly $P_{0}(\omega)=0$ (the mode exists) with $Q(\omega)=0$ (the mode is at stability limit):

[TABLE]

This means that the plasma is unstable when the fluid electron drift is higher than the ion sound speed. The numerical resolution of the two other roots of $P$ shows that they are all stable in the regime of interest here. Using Eqs. (25) and (6) for low wavenumbers, the plasma is unstable if:

[TABLE]

At the sheath edge, $v_{x}=\mu^{1/2}$, and $1/\nu+\nu$ is greater than 1 for all $\nu>0$. Therefore, the plasma is always unstable, at least at the sheath edge, as long as the electrons are magnetized ($l\gg 1$). The destabilization of similar modes by collisions was first found experimentally Birdsall et al. (1953); Hendel et al. (1967) and explained theoretically by Chen Chen (1965a) as a particular type of resistive drift mode.Chen (1965b) More recently, several experiments highlighted drift wave instabilities in magnetized plasma columns.Brochard et al. (2005); Gravier et al. (2004) As shown by Lakhin et al. taking into account finite Larmor radius effects,Lakhin et al. (2018a, b) the condition $-n^{\prime}\phi^{\prime}/n<\mu/4$ is sufficient to ensure electrostatic modes stability in a collisionless plasma, and only collisions can destabilize it in this case.

In order to validate the 1D model, the PIC/MCC simulation needs to be 2D to capture both the direction of the instability and the gradients of the equilibrium solution. Simulations of a square argon discharge were performed at gas pressures of 3, 6, and 12 mTorr and uniform magnetic fields between 0 and 40 mT oriented along $z$, perpendicularly to the simulation plane. The discharge is sustained by a uniform radio-frequency (RF) electric field at 13.56 MHz along $z$. The amplitude of the heating electric field is adjusted at each RF cycle to control the power absorbed by the plasma, fixed to 9.6 kW/m3. The $x$ and $y$ components of the electric field come from the solution of the Poisson’s equation, with the potential set to zero at the conducting walls.

Super-particles lost at the walls are discarded and new particles are generated through self-consistent ionization with a uniform background gas at 300 K. The newly created ions and electrons are initialized with Maxwellian distribution functions at 0.026 eV and 4 eV respectively. The ion temperature remains lower than 0.2 eV in all simulation conditions. The steady-state electron temperature is nearly uniform in the bulk plasma and varies between 3.6 and 5.3 eV depending on the cases. A realistic kinetic scheme for argon is used Lucken et al. (2017) using cross sections from the LXCat database.Biagi (2004); www.lxcat.net ; Phelps (1994) The simulation time step and the cell size satisfy the classical stability criteria (Birdsall and Langdon Birdsall and Langdon (2004)). The number of particles per cell was always greater than 100 in the center (150 in most cases). A simulation at 20 mT and 3 mTorr was performed with 4 times more particles per cell (600 part./cell in the discharge center) with a change of less than 3% for the discharge global properties.

Plasma sheaths form in about one microsecond, and at high magnetic field, instabilities rotating in the electron drift direction develop rapidly. After a transient of about 20 µs (depending on the pressure and the magnetic field), the volume-averaged plasma properties reach a steady-state and the instability saturates with a continuous spectrum featuring a maximum around 3 MHz. The instability wavevector is directed mainly in the azimuthal direction.

Figure 1 shows the electron density and the electron current at 3 mTorr and 20 mT. When the plasma properties are averaged over several tens of microseconds, the discharge aspect is the same as in the non-magnetized case,Lucken et al. (2018) the only qualitative difference being that the electrons are rotating in the direction of the diamagnetic drift. On a shorter time-scale, the instabilities are visible and the waves seem to break against the sheath.

In Figure 2, the electron azimuthal drifts measured in the PIC simulations are compared with simple asymptotic formulae coming from the model equations. Neglecting curvature effects, the azimuthal direction is treated as $y$, and the radial direction as $x$. The colored dashed lines in Figure 2 correspond to the weakly-magnetized limit of Eq. (15),

[TABLE]

for each value of the pressure. The following formula coming from 2D theories of non-magnetized plasma transport was used for the $h$ factor Lucken et al. (2018)

[TABLE]

where $\lambda_{i}$ is the (dimensional) ion mean free path and $T_{i}$ is the ion temperature. The solid black line in Figure 2 corresponds to the diamagnetic drift

[TABLE]

For simplicity, a reference temperature of 3.6 eV was considered in all cases.

Figure 3(a) shows the azimuthal electric field as a function of time and $\theta$ in a typical case where the instability drives the transport. The 2D Fourier transform (FT) of this signal is shown in Figure 3(b). The dashed and solid black lines are the approximate solutions of Eqs. (32) and (33) respectively, while the full numerical solution of Eq. (30) using the parameters of Table 1 is displayed in cyan. The electron density, temperature, and drifts are extracted from the time-averaged solution of the PIC simulation in saturated state.

An effective collision frequency is extracted from the simulation data using Eq. (3)

[TABLE]

where $\langle\Gamma_{er}\rangle$ and $\langle\Gamma_{e\theta}\rangle$ are the electron fluxes in the radial and azimuthal directions respectively, averaged along $\theta$. The effective collision frequency of Eq. (39) is normalized to the classical collision frequency known from the literature Lieberman and Lichtenberg (2005) and plotted in Figure 4. The error bars correspond to all the data collected at various radial positions, from 1 to 14 mm from the discharge center. At low magnetic field, the effective collision frequency is equal to the classical collision frequency. In the strongly magnetized regime, the effective collision frequency is correctly described by Eq. (23), assuming an $h$ factor that does not depend on the magnetic field (Eq. (37)).

A 1D model of the cross field plasma transport was derived and validated using 2D PIC simulations. It was shown that if the diamagnetic drift at the sheath edge remains below the electron thermal velocity, then the effective electron collision frequency should depend on the magnetic field. This result is consistent with the plasma unstable behavior predicted by the linear theory of wave perturbations, and by PIC simulations. The effective collision frequency scales as $\omega_{ce}^{2}$ at high magnetic field. The transport across a magnetic field in a weakly ionized plasma is much better described by assuming that the $h$ factor does not depend on the magnetic field than by assuming a classical collision frequency.iiiiiiThe dependence of the $h$ factor on the magnetic field, from the classical regime to the instability-enhanced regime, will be the focus of another paper. ,Sternberg et al. (2006) The plasma transport becomes completely insensitive to the magnetic field due to the electron drift resistive instability.

*Acknowledgments – * The authors are grateful to A. J. Lichtenberg, A. Tavant, R. Martorelli, and O. Gurcan for many enlightening discussions. This work was granted access to the HPC resources of CINES under the allocation 2017-A0020510092 made by GENCI, and it was partially funded by the Agence Nationale de la Recherche (ANR-16-CHIN-0003-01) and the U.S. Department of Energy Office of Fusion Science Contract DE-SC0001939.

*Appendix: Electron response to a perturbed electrostatic potential

The normalized electron momentum conservation equation is

[TABLE]

Assuming

[TABLE]

where $\epsilon$ is a small parameter, we can expand $\mathbf{v}$ with respect to $\epsilon$. To the 0-th order,

[TABLE]

which yields the $E\times B$ and diamagnetic drifts, respectively

[TABLE]

Eq. (40) is then to the first order

[TABLE]

where $\mathbf{v}^{(1)}=\mathbf{v}-\mathbf{v}_{E}-\mathbf{v}_{d}$ is the first order term of the electron velocity. Hence,

[TABLE]

where

[TABLE]

and

[TABLE]

are the polarization drift terms due to the $E\times B$ and the diamagnetic drifts, respectively (including the friction force). It is useful to derive the fluxes corresponding to each of the drift terms.

[TABLE]

The second term is zero due to the formula

[TABLE]

Thus,

[TABLE]

Using again Eq. (50),

[TABLE]

The relationship is known as the gyroviscous cancellation.Smolyakov et al. (2017) By neglecting all particle source and loss terms ($\nu_{iz}\ll\omega$), the electron continuity equation is

[TABLE]

So to the first order in $\epsilon$,

[TABLE]

or

[TABLE]

The only term contributing to the motion of the electron guiding centers is the $E\times B$ drift Weiland (2012)

[TABLE]

Using the property

[TABLE]

[TABLE]

For the density gradient polarization drift,

[TABLE]

We assume

[TABLE]

with $n_{1}\ll n_{0}$ and $\phi_{1}\ll\phi_{0}$, and $n_{1}$ and $\phi_{1}$ proportional to $\exp(-i\omega t+iky)$. To the first order in $\phi_{1}$ and $n_{1}$:

[TABLE]

where $\omega_{*}=-kn_{0}^{\prime}/n_{0}$ and $\omega_{0}=-k\phi_{0}^{\prime}$. Similarly,

[TABLE]

Finally, Eq. (51) yields

[TABLE]

To the first order in $\omega_{*}$, and neglecting second and third order derivatives of $n_{0}$ and $\phi_{0}$ (no shear),

[TABLE]

Eq. (53) is therefore:

[TABLE]

Hence,

[TABLE]

This result is valid within the assumptions of the model derived here (Eqs. (41) and (42)), but can be generalized by the Padé approximation Pegoraro and Schep (1986); Schep et al. (1994) to all wavenumbers.