Fast correlation heating in moderately coupled electron–ion plasmas

Abstract

If the electrons in a plasma are suddenly heated, the resulting change in Debye shielding causes the ion kinetic energy to quickly increase. For the first time, this correlation heating, which is much faster than collisional energy exchange, is rigorously derived for a moderately coupled, electron–ion plasma. The electron–ion mass ratio is taken to be the smallest parameter in the Bogoliubov–Born–Green–Kirkwood–Yvon hierarchy, smaller even than the reciprocal of the plasma parameter. This ordering differs from conventional kinetic theory by making the electron collision rates faster than the ion plasma frequency, which allows stronger coupling and makes the ion heating a function only of the total energy supplied to the electrons. The calculation uses known formulae for correlations in a two-temperature plasma, for which a new, elementary derivation is presented. Suprathermal ions may be created more rapidly by this mechanism than by ion–electron Coulomb collisions. This means that the use of a femtosecond laser pulse could potentially help to achieve ignition in certain fast ignition approaches to inertial confinement fusion.

1 Introduction

Plasmas in various experimental contexts are heated via a rapid input of energy primarily to the electrons. Examples include inertial confinement fusion (ICF) plasmas (Lindl 1995), ultracold neutral plasmas (Killian et al. 1999) or photoionised plasmas created by X-rays from a Z-pinch (Rochau et al. 2014). The large mass difference between ions and electrons means that Maxwellianisation of each species occurs much faster than collisional energy exchange between species (Spitzer 1956; Montgomery & Tidman 1964), so the sudden electron heating results in a transitory two-temperature state that persists until collisions eventually restore thermal equilibrium. Since this equilibration is slow, one might ask the following question: in a two-temperature plasma with massless electrons, meaning collisional energy exchange is negligible, does suddenly heating the electrons have any effect on the ion temperature?

Intriguingly, the answer to this question is yes: the ion temperature rapidly jumps upwards. After the electrons are heated, the ions are shielded less effectively and nearby ions repel each other more strongly. Therefore, the ions tend to push away from each other, which converts potential energy into ion kinetic energy (see figure 1). This mechanism is an example of correlation heating, in which changing the spatial correlations between particles leads to a rapid conversion of potential energy into kinetic energy.

Figure 1. Physical explanation of correlation heating. The black circles represent ions, and the coloured circles represent the electron density around each ion due to Debye screening. (a) Two nearby ions are initially well shielded by small Debye clouds of cold electrons. The screened ions interact weakly. (b) If the electrons are suddenly heated, the Debye spheres get larger and the ions are screened less effectively. They suddenly repel more strongly. (c) The ions repel each other and fly apart. Their potential energy is converted into kinetic energy.

This correlation heating mechanism operates on ion-plasma-frequency time scales ${t\sim 1/\omega _{pi}}$, which is much faster than energy exchange due to interspecies Coulomb collisions (More 1980, 1981,  1983).

This curious effect is similar to ionisation or ‘proximity’ heating, in which suddenly changing the charge state of the ions causes them to push away from each other and heat up (Hagelstein 1981; Boercker & More 1986; More 1986; Evans 2008; Baggott, Rose & Mangles 2021). Monte Carlo calculations by Hagelstein (1981) show that when lasers are used to ionise neon gas, the ion temperature could increase by more than an order of magnitude due to ionisation heating. This in turn is reminiscent of the rapid conversion of potential energy into ion kinetic energy by Coulomb explosions (Fleischer, Price & Walker 1965), which has been proposed as a mechanism for creating fast deuterons for nuclear fusion (Last & Jortner 2001; Ron, Last & Jortner 2012).

Another example of correlation heating is disorder-induced heating of ultracold neutral plasmas (e.g. Murillo 2001; Kuzmin & O'Neil 2002; Gericke & Murillo 2003; Simien et al. 2004; Cummings et al. 2005; Killian et al. 2007; Lyon & Rolston 2016). These plasmas are generated when neutral atoms in a magneto-optical trap are cooled to low temperatures, typically $1\,{\rm mK}$ or even lower, and suddenly ionised by lasers. The newly formed ions quickly transition from an uncorrelated, disordered initial state to a state with correlations imposed by their strong Coulomb repulsion. This change in the particle correlations causes a tremendous increase in the ion temperature, which can climb to $1\,{\rm K}$ or more.

However, ultracold neutral plasmas are low-temperature systems in which quantum effects are important. The correlation heating discussed in this paper, caused by changing Debye shielding, is a universal phenomenon, which could occur in any hot, classical plasma in which electrons shield ions. The enormous ion temperature increase that can be achieved using ionisation heating or disorder-induced heating is a tantalising prospect for fusion. This is particularly true for fusion approaches that use more strongly coupled plasmas, in which there is a large reservoir of potential energy.

In this paper, we calculate the ion temperature change due to correlation heating after a sudden input of energy to the electrons in a moderately coupled plasma, by which we mean a plasma satisfying the ordering $1\ll \varLambda /\lambda \ll (m_i/m_e)^{1/2}$. Here, $\varLambda$ is the plasma parameter (the number of particles in a Debye sphere) and $\lambda$ is the Coulomb logarithm. The assumption of moderate coupling allows a novel, rigorous expansion of the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) hierarchy, giving a first-principles derivation of correlation heating in a two-component plasma with electrons and ions. Previous work has investigated the effect of changing shielding in one-component plasmas, both theoretically (Morawetz et al. 2001) and numerically (Gericke et al. 2003; Chen et al. 2004).

Our calculations require formulae for two-particle correlations in a two-temperature plasma, first presented by Salpeter (1963). In § 3.1, we present a streamlined derivation of these correlations that employs only basic plasma physics concepts.

We compute the time-dependence of the ion kinetic energy following electron heating and find that it does not overshoot its final value. We also calculate the modification to the ion velocity distribution and find that correlation heating in a moderately coupled plasma is relatively inefficient at heating suprathermal ions, which are crucial for fusion. However, we argue that correlation heating is more effective at creating fast ions if the plasma is more strongly coupled or if there is a light minority ion species present. This may be particularly relevant for fast ignition approaches to ICF. Even if few fast ions are created directly, the combination of correlation heating and subsequent ion–ion collisions could allow the ions in a fusion device, including suprathermal ions, to be energised more rapidly than would be possible using ion–electron collisions alone.

In complementary work (Fetsch, Foster & Fisch 2023) that considers time scales longer than the $1/\omega _{pi}$ time scales resolved here, the same result for the total change in ion temperature when the electrons are suddenly heated is obtained. Fetsch et al. (2023) generalises thermodynamic relations to two-temperature plasma and obtains a quasi-equation of state. This is used to compute temperature changes under compression and heating of either species; for example, initially equal electron and ion temperatures are found to separate under compression. The thermodynamic method also allows calculation of the entropy of each species in a two-temperature plasma, thus providing a particularly effective way to analyse processes that are adiabatic for at least one species.

2 Ordering assumptions

We are interested in classical, homogeneous, non-magnetised plasmas consisting of electrons and a single ion species of charge $+Ze$. Since the ion charge is fixed, ionisation heating does not occur. This allows us to isolate the effect of correlation heating due to changes in the Debye shielding of ions. If the electron component and ion component have had time to Maxwellianise individually, due to intra-species collisions, then the velocity distribution of species $\alpha$ (with $\alpha = i$ for ions and $\alpha = e$ for electrons) is

(2.1)\begin{equation} f_\alpha(\boldsymbol{v}) = n_\alpha\left( \frac{m_\alpha}{2{\rm \pi} T_\alpha} \right)^{3/2}\exp\left( -\frac{m_\alpha \boldsymbol{v}^2}{2 T_\alpha} \right), \end{equation}

where $n_\alpha$, $m_\alpha$ and $T_\alpha$ are the number density, particle mass and temperature, respectively, of species $\alpha$. This defines what we mean by species-dependent temperature; the temperature of a species is the characteristic width of its Maxwellian velocity distribution.

In this section, we give the ordering assumptions used throughout the paper, which allow the plasma to exist in a two-temperature state. They also define a regime in which correlation heating is analytically tractable and they provide a fair model of real plasmas from experimental contexts involving sudden electron heating. The section concludes with parameters from two example plasmas.

2.1 Masses

Two-temperature plasma can exist when $\tau _{ie}$, the time scale of collisional energy exchange between species, is much greater than the time scales of Maxwellianisation for each species. The Maxwellianisation time scales can be estimated using the ion–ion collision time $\tau _{ii}$ and the electron–electron collision time $\tau _{ee}$. It is a standard result that the considerable mass difference between ions and electrons in a plasma leads to a large separation between these collisional time scales, since $\tau _{ee}/\tau _{ii}\sim \tau _{ii}/\tau _{ie}\sim (m_e/m_i)^{1/2}$ (Spitzer 1956; Montgomery & Tidman 1964). Our first assumption is therefore

(2.2)\begin{equation} \text{Assumption 1:}\quad \frac{m_e}{m_i}\ll 1. \end{equation}

Our formalism treats this mass ratio as the smallest quantity in the problem, unlike conventional plasma kinetic theory in which the weak-coupling limit is taken before the mass ratio is sent to zero. It is worth pointing out that these limits do not commute; this fact is more than just a technical curiosity and will be crucial later for understanding how the system responds when the electrons are suddenly heated.

2.2 Temperatures

A small electron–ion mass ratio is not sufficient to ensure that the species temperatures equilibrate slowly compared with all other time scales in the problem. There must also be a restriction on the ion and electron temperatures because collisional energy exchange can occur quickly if the difference between them is large enough. How far apart can the temperatures be?

In a plasma with large plasma parameter $\varLambda \gg 1$ and the mild restriction ${T_e/m_e\gg T_i/m_i}$, the collisional time scales $\tau _{ee}$, $\tau _{ii}$, $\tau _{ie}$ and $\tau _{ei}$ (the time scale of momentum transfer between electrons and ions) can be determined using the Landau–Spitzer formulae in table 1 (Landau 1936; Spitzer 1956). We constrain the temperatures by assuming that these time scales satisfy

(2.3)\begin{equation} \tau_{ee}, \tau_{ei}\ll \tau_{ii}\ll \tau_{ie}. \end{equation}

This ordering was chosen for two reasons. First, it ensures that the temperatures change slowly compared with the Maxwellianisation rates and compared with the inter-species momentum exchange rate.¹ This justifies the ‘two-temperature’ model in which each species is Maxwellian, and has the same mean flow velocity, but with different temperatures.

Table 1. Landau–Spitzer formulae for collisional time scales, following Hazeltine & Waelbroeck (2018). The $\lambda _{\alpha \beta }$ are Coulomb logarithms.

Second, we have chosen to make the electrons Maxwellianise faster than the ions, which will be justified shortly.

Using the formulae in table 1 with $Z\sim 1$ and $\lambda _{ee}\sim \lambda _{ei}\sim \lambda _{ii}$, condition (2.3) means the temperatures must satisfy (Bobylev et al. 1997)

(2.4)\begin{equation} \text{Assumption 2:}\quad \left(\frac{m_e}{m_i}\right)^{1/3}\ll\frac{T_e}{T_i}\ll \left(\frac{m_i}{m_e}\right)^{1/3}. \end{equation}

In a hydrogen plasma, $(m_i/m_e)^{1/3} \approx 12$, so this constraint on the temperatures is satisfied in most situations of interest. From now on, we order $T_e\sim T_i$.

2.3 Plasma parameter

Our final assumption concerns the plasma parameter $\varLambda$, which is the number of particles in a Debye sphere and which measures the strength of interactions in a plasma.

We order the electron and ion plasma frequencies relative to the collisional time scales as follows:

(2.5)\begin{equation} 1/\omega_{pe}\ll \tau_{ee}, \tau_{ei}\ll 1/\omega_{pi} \ll \tau_{ii}. \end{equation}

This means $\varLambda$ must be large but also includes the unusual condition that the ion plasma frequency $\omega _{pi}$ is slower than the electron Maxwellianisation rate $1/\tau _{ee}$. The electron time scales $\omega _{pe}$ and $\tau _{ee}$ are smaller that the corresponding ion time scales $\omega _{pi}$ and $\tau _{ii}$ by a factor of $(m_e/m_i)^{1/2}$, so they are taken to be the smallest time scales in the problem. This is valid as long as $\varLambda$ is not too large; the resulting regime, while narrow, is relevant for certain plasmas in Z-pinch and ICF experiments. In addition, our correlation heating calculations will be simplified because the electrons Maxwellianise before the ions respond to the heating, which means the ions respond independently of where in velocity space the electrons were energised.

The plasma frequency and Maxwellianisation time scale of a species are related by $1/\omega _{p\alpha }\sim \lambda \tau _{\alpha \alpha }/\varLambda$, where $\lambda$ is the Coulomb logarithm and

(2.6)\begin{equation} \varLambda \sim \left( \frac{T}{n^{1/3}e^2} \right)^{3/2}. \end{equation}

Species indices are not necessary here because we have taken $Z\sim 1$ and $T_e\sim T_i$. Therefore, (2.5) is equivalent to

(2.7)\begin{equation} \text{Assumption 3:}\quad 1\ll \frac{\varLambda}{\lambda}\ll\left( \frac{m_i}{m_e} \right)^{1/2}. \end{equation}

We refer to this as the moderate coupling regime because $\varLambda$ is much greater than one but cannot be arbitrarily large. We caution the reader that this term is often used to describe plasmas with $\varLambda \sim 1$. The large plasma parameter will allow the BBGKY hierarchy to be truncated, giving a closed set of equations which we will solve to find the ion evolution during correlation heating.

2.4 Example parameters

In table 2, we give parameters for two example plasmas.

Table 2. Example parameters from two experimental contexts involving fast electron heating. The coupling strengths $\varGamma _\alpha = (4{\rm \pi} n_\alpha /3)^{1/3}q_\alpha ^2/T_\alpha$ are the ratio of the typical potential energy to the typical kinetic energy of a particle of species $\alpha$. The collision time scales are calculated using the formulae in table 1 together with the conventional definitions of the Coulomb logarithms (Richardson 2019).

The first set of parameters is for a hydrogen plasma (Falcon et al. 2015) created in the Z Pulsed Power Facility at Sandia National Laboratories (Sinars et al. 2020) as part of the Z Astrophysical Plasma Properties (ZAPP) collaboration (Rochau et al. 2014). X-rays from the Z-pinch are used to photoionise hydrogen to study spectral lines of hydrogen plasma under conditions typical of white dwarf photospheres.

The plasma exists for approximately $100\,{\rm ns}$ and, after an initial period where the plasma is overionised, the electron temperature and density eventually stabilise around the given values. The ion temperature is not measured but should be lower than the electron temperature, meaning the ions may be weakly or strongly coupled.

The second set of parameters in table 2 is for an ICF hot spot plasma (deuterium–tritium) in the fast ignition scenario (Tabak et al. 2006). Conventionally, ICF uses stagnation of converging flows to heat a central region up to ignition. In fast ignition, however, it is envisaged that the target would first be imploded with a uniform density, before ignition would be triggered in a central hot spot with a laser pulse (Tabak et al. 1994). This laser pulse would heat the hot spot indirectly by creating fast electrons in the surrounding plasma that rapidly propagate inwards, although other schemes, for example, involving ion beams, have also been proposed (Bychenkov et al. 2001; Logan et al. 2006). Fast ignition offers various advantages over the conventional approach to ICF, including higher gain, lower driver energy and relaxation of symmetry requirements (Badziak, Jabłoński & Wołowski 2007). Fast ignition requirements in two-temperature plasmas have been studied, for example, by Eliezer et al. (2015).

Both example plasmas have coupling strengths that are less than one without being extremely small. Also, both plasmas have electron–electron and electron–ion collision times of the same order as the ion plasma frequency, instead of being much longer as is usually the case in very weakly coupled systems. Since the exact parameters may be varied in the experiment, the plasma in these ICF and Z-pinch contexts may be well described by a moderate-coupling ordering.

3 Prerequisite: pair correlations in a two-temperature plasma

Expressions for the steady-state pair correlations in a two-temperature plasma will be an essential ingredient for the analytical work in this paper.

Pair correlations describe the relative positions of particles in a system. In a thermal equilibrium plasma, assuming large plasma parameter $\varLambda \gg 1$, the pair correlations take a well-known form (Krall & Trivelpeice 1973),

(3.1)\begin{equation} G_{\alpha\beta}(\boldsymbol{r}) ={-}\frac{q_\alpha q_{\beta}}{T}\frac{{\rm e}^{{-}k_D r}}{r}, \end{equation}

where $G_{\alpha \beta }(\boldsymbol {r})$ is the spatial pair correlation between species $\alpha$ and species $\beta$ (see Appendix A for the definitions and notation that we use for pair correlations and distribution functions). Here, $k_D$ is the inverse Debye length, defined by $k_D^2 = k_e^2 + k_i^2$ with

(3.2)\begin{equation} k_\alpha^2 = \frac{4{\rm \pi} n_\alpha q_\alpha^2}{T_\alpha}. \end{equation}

The generalisation of (3.1) to a two-temperature plasma was first found by Salpeter (1963) to study Thomson scattering of radio waves in the ionosphere, which has hotter electrons than ions. Salpeter's results corrected earlier work (Kadomtsev 1958; Renau 1962; Renau, Camnitz & Flood 1962; Renau 1964) and were later rederived using various other methods (Ecker & Kröll 1964; Schram & Kegel 1965; Boercker & More 1986; Seuferling, Vogel & Toepffer 1989; Fetsch et al. 2023). Correlation functions in two-temperature plasmas were thoroughly studied by Seuferling (1987). There has been recent interest in extending Salpeter's formulae to more strongly coupled plasmas (Shaffer, Tiwari & Baalrud 2017; Triola 2022).

In this section, we present two derivations of the steady-state pair correlations in a two-temperature plasma. The first derivation is simpler than other published methods and is based on Salpeter's original approach, using arguments in the style of Debye & Hückel (1923). The second derivation is more rigorous and systematic, and introduces the formalism that will be needed later to obtain analytical results for correlation heating.

3.1 Debye–Hückel method

The key idea is to exploit the fact that the mass ratio $m_e/m_i$ is a small parameter. As long as the ions are not too hot, which is ensured by assumption (2.4), the electron thermal speed is much greater than the ion thermal speed; compared with the electrons, the ions move extremely slowly and are essentially stationary. We therefore begin by finding the steady-state response of the electrons to a set of fixed ions. This will determine the electron–ion and ion–ion correlations.

Assume the electron distribution is a uniform Maxwellian $f_{e0}(\boldsymbol {v})$ plus a small perturbation $f_{e1}(\boldsymbol {v},\boldsymbol {r})$ with vanishing spatial average. The perturbation satisfies the linearised Vlasov equation,

(3.3)\begin{equation} \left( \frac{\partial}{\partial t} + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} \right) f_{e1} + \frac{e}{m_e}(\boldsymbol{\nabla}\varphi)\boldsymbol{\cdot}\frac{\partial f_{e0}}{\partial \boldsymbol{v}} = 0. \end{equation}

If the average potential is taken to be zero, there is a simple steady-state solution

(3.4)\begin{equation} f_{e1} = \frac{e\varphi}{T_e}f_{e0}. \end{equation}

Following Salpeter (1963), we have disregarded any terms that depend on time or that do not involve the potential $\varphi (\boldsymbol {r})$. Equation (3.4) is the usual linearised adiabatic response formula. We could also justify this result by arguing that the electrons reach thermal equilibrium around essentially fixed ions, which means the ion potential can be treated as an external potential acting on the system consisting of the electrons only. Then the Boltzmann distribution for the electrons $f_e(\boldsymbol {r},\boldsymbol {v}) \propto f_{e0}(\boldsymbol {v})\exp (e\varphi (\boldsymbol {r})/T_e)$ leads to (3.4), assuming weak interactions $e\varphi \ll T_e$. It will be convenient to work with Fourier transforms (our convention is given in Appendix A) from now on, so we have

(3.5)\begin{equation} \tilde{f}_{e1}(\boldsymbol{k},\boldsymbol{v}) = \frac{e\tilde{\varphi}(\boldsymbol{k})}{T_e}f_{e0}(\boldsymbol{v}). \end{equation}

To completely specify the electron response to the ions, we need to calculate the potential $\tilde {\varphi }$ for any given ion distribution. It satisfies Poisson's equation

(3.6)\begin{equation} k^2\tilde{\varphi} = 4{\rm \pi}\,(\tilde{\rho}_e + \tilde{\rho}_i), \end{equation}

where $\rho _e$ and $\rho _i$ are the charge densities of the electrons and ions, respectively. The only part of the electron distribution with spatial variation is $f_{e1}$, so $\tilde {\rho }_e(\boldsymbol {k}) = \int f_{e1}(\boldsymbol {k},\boldsymbol {v})\,{\rm d}\boldsymbol {v}$ for non-zero $\boldsymbol {k}$. Then (3.5) gives

(3.7)\begin{equation} \tilde{\rho}_{e} ={-}\frac{n_ee^2\tilde{\varphi}}{T_e}. \end{equation}

Poisson's equation becomes

(3.8)\begin{equation} (k^2+k_e^2)\tilde{\varphi} = 4{\rm \pi} \tilde{\rho}_i, \end{equation}

where $k_e$ is defined by (3.2). Therefore, if the ion charge density is

(3.9)\begin{equation} \rho_i(\boldsymbol{r}) = \sum_{\mathrm{ions\, }i} Ze\,\delta(\boldsymbol{r}-\boldsymbol{r}_i), \end{equation}

corresponding to a collection of fixed ions, then the resulting potential is

(3.10)\begin{equation} \varphi(\boldsymbol{r}) = \sum_{\mathrm{ions }\,i}\frac{Ze}{|\boldsymbol{r}-\boldsymbol{r}_i|}\,{\rm e}^{{-}k_e|\boldsymbol{r}-\boldsymbol{r}_i|} + \mathrm{const.}, \end{equation}

where the constant ensures that the average potential is zero. We can think of the ions as interacting via a screened Coulomb potential (Yukawa potential)

(3.11)\begin{equation} \varphi^{(s)}(\boldsymbol{r}) = \frac{Ze}{r}\,{\rm e}^{{-}k_er}, \end{equation}

where $r=|\boldsymbol {r}|$, with shielding length set by the electron temperature.

For any given ion distribution, we can use (3.5) and (3.8) to calculate the electron response. The resulting electron charge density is

(3.12)\begin{equation} \tilde{\rho}_{e} ={-}\left(\frac{k_e^2}{k^2+k_e^2}\right)\tilde{\rho}_i. \end{equation}

It is worth noting a subtle point which will be important later. In writing (3.4), we assumed that the electron distribution is a steady-state solution of the linearised Vlasov equation. In reality, the electron distribution will fluctuate rapidly as the electrons move around and take part in plasma oscillations; by ignoring these fluctuations, we have essentially averaged the electron distribution over a time period much longer than the fluctuation time scale and much shorter than the time scale of ion motion (these time scales are well separated because of the small mass ratio $m_e/m_i$). The fluctuating part of the electron distribution depends on the fluctuating part of the potential, which is independent of the slowly varying ion distribution. Therefore, the electron fluctuations are independent of the ion distribution and they do not affect the electron–ion or ion–ion correlations. They will, however, be relevant for the electron–electron correlation.

Knowing the response of electrons to fixed ions, we can now work out the ion–ion correlation. The ions are a system of particles with temperature $T_i$ interacting via the screened potential in (3.11), and equilibrium statistical mechanics can be used to calculate the pair correlations for such a system.²

Instead of using statistical mechanics and calculating a partition function, we will follow Debye & Hückel (1923) and assume that the ion–ion correlation is proportional to the excess ion density around a single, fixed ion at the origin. Around the fixed ion, the Boltzmann response formula gives

(3.13)\begin{equation} \rho_i(\boldsymbol{r}) \propto \exp\left( -\frac{Ze\varphi(\boldsymbol{r})}{T_i} \right), \end{equation}

or

(3.14)\begin{equation} \rho_i(\boldsymbol{r}) = Zen_i\left( 1 - \frac{Ze\varphi(\boldsymbol{r})}{T_i} \right), \end{equation}

assuming $Ze\varphi /T_i\ll 1$. So, when there is a fixed ion at the origin, Poisson's equation (3.8) becomes

(3.15)\begin{equation} (k^2+k_e^2+k_i^2)\tilde{\varphi} = 4{\rm \pi} Ze, \end{equation}

which has solution

(3.16)\begin{equation} \tilde{\varphi}(\boldsymbol{k}) = \frac{4{\rm \pi} Ze}{k^2+k_e^2+k_i^2}. \end{equation}

In real space, this is

(3.17)\begin{equation} \varphi(\boldsymbol{r}) = \frac{Ze}{r}\, {\rm e}^{{-}k_Dr}, \end{equation}

where $k_D$ is now the total inverse Debye length because both the ions and the electrons are taking part in the shielding of the test ion at the origin. Then, the ion and electron densities around this fixed ion are

(3.18)\begin{gather} \rho_i(\boldsymbol{r}) = Zen_i\left( 1 - \frac{Z^2e^2}{T_i}\frac{{\rm e}^{{-}k_Dr}}{r} \right), \end{gather}

(3.19)\begin{gather}\rho_e(\boldsymbol{r}) ={-}en_e\left( 1 + \frac{Ze^2}{T_e}\frac{{\rm e}^{{-}k_Dr}}{r} \right). \end{gather}

We can read off the ion–ion correlation, which is

(3.20)\begin{equation} G_{ii}(\boldsymbol{r}) ={-}\frac{Z^2e^2}{T_i}\frac{{\rm e}^{{-}k_D r}}{r}. \end{equation}

The excess electron density around a fixed ion also tells us the electron–ion correlation because the ions move so slowly that they can all be approximated as stationary. So, we can also read off

(3.21)\begin{equation} G_{ei}(\boldsymbol{r}) = \frac{Ze^2}{T_e}\frac{{\rm e}^{- k_D r}}{r}. \end{equation}

Now we just need the electron–electron correlation. The Debye–Hückel method of fixing a particle at the origin and letting the other particles reach thermal equilibrium worked well for finding the electron density around any given ion, since the ions all appear virtually stationary to the electrons. However, the ion density around an electron cannot be found this way.

Nevertheless, the fact that the electrons form identical screening clouds around each ion suggests that we should be able to calculate the electron spatial distribution from the ion spatial distribution, which is encoded in the ion–ion correlations. The electron and ion spatial distributions are related by (3.12); taking the modulus squared, then the thermal average gives

(3.22)\begin{equation} \langle|\tilde{\rho}_e|^2\rangle = \left(\frac{k_e^2}{k^2+k_e^2}\right)^2\langle|\tilde{\rho}_i|^2\rangle. \end{equation}

This equation links fluctuations in the electron charge density with fluctuations in the ion charge density. It is well known that charge density fluctuations are related to pair correlations. We use a formula derived in Appendix B,

(3.23)\begin{equation} \langle|\tilde{\rho}_\alpha|^2\rangle = Vn_\alpha q_\alpha^2 \left( 1 + n_\alpha \tilde{G}_{\alpha\alpha} \right), \end{equation}

where $V$ is the system volume, to find

(3.24)\begin{equation} \tilde{G}_{ee}(\boldsymbol{k}) \stackrel{?}{=} \frac{4{\rm \pi} Z e^2}{T_e}\frac{k_e^2}{(k^2+k_D^2)(k^2+k_e^2)}-\frac{1}{n_e}. \end{equation}

However, this result must be wrong! If the ion charge is taken to zero $Z\to 0$ with $n_e$ fixed, then the ions become unaffected by Coulomb forces and uniformly distributed. Equation (3.24) gives $G_{ee} = -1/n_e$ in this limit, corresponding to ${\langle |\tilde {\rho }_e|^2\rangle = 0}$. That cannot be correct! The electrons screen each other, and have non-zero charge fluctuations $\langle |\tilde {\rho }_e|^2\rangle$, even when $Z\to 0$.

The problem is that, as discussed earlier, we ignored the rapid fluctuations in the electron distribution that correspond to the electrons moving around and taking part in plasma oscillations. It is precisely these fluctuations that are missing from the pair correlation. They are independent of the ions, so the missing term in the electron–electron correlation is the correlation that remains when $Z\to 0$. In this limit, the ions become a uniform neutralising background and we are left with a one-component plasma of electrons at temperature $T_e$. The pair correlation for such a system is, using (3.1),

(3.25)\begin{equation} G_{ee}^{(Z\to 0)}(\boldsymbol{r}) ={-}\frac{e^2}{T_e}\frac{{\rm e}^{{-}k_er}}{r}. \end{equation}

A constant term, independent of the ions, should be added to the part of the electron–electron correlation arising from their interaction with the ions, which is given in (3.24), so that the $Z\to 0$ limit is consistent with (3.25).

In the end, the correct pair correlations in a two-temperature plasma are

(3.26)\begin{gather} G_{ii}(\boldsymbol{r}) ={-}\frac{Z^2e^2}{T_i}\frac{{\rm e}^{{-}k_D r}}{r} , \end{gather}

(3.27)\begin{gather}G_{ei}(\boldsymbol{r}) = \frac{Ze^2}{T_e}\frac{{\rm e}^{- k_D r}}{r} , \end{gather}

(3.28)\begin{gather}G_{ee}(\boldsymbol{r}) ={-}\frac{e^2}{T_e}\left[ \frac{T_i}{T_e}\frac{{\rm e}^{{-}k_D r}}{r} + \left( \frac{T_e - T_i}{T_e}\right)\frac{{\rm e}^{{-}k_e r}}{r} \right]. \end{gather}

When the temperatures $T_e$ and $T_i$ are equal, these results reduce to the usual equilibrium formula (3.1).

This derivation is simple and intuitive. However, just like the thermal equilibrium correlations (Krall & Trivelpeice 1973), Salpeter's two-temperature correlations can be derived using the BBGKY hierarchy. We now present this alternative, more systematic method, which introduces the formalism that we will use to analyse the response of a two-temperature plasma to sudden electron heating.

3.2 Kinetic method

3.2.1 BBGKY equations

The first two equations in the BBGKY hierarchy for a classical, unmagnetised plasma are (Balsecu 1988, p. 102)

(3.29)\begin{align} & \left(\frac{\partial}{\partial t} + \mathcal{L}(1)\right) f_\alpha(1) ={-}\sum_\beta \left( \int \,{\rm d}(2) \mathcal{V}_{\alpha\beta}(12)\,g_{\alpha\beta}(12) \right),\end{align}

(3.30)\begin{align} & \left(\frac{\partial}{\partial t} + \mathcal{L}(12)\right)g_{\alpha\beta}(12) + \sum_\gamma \left(\int \,{\rm d}(3) \mathcal{V}_{\alpha\gamma}(13) f_\alpha(1) g_{\gamma\beta}(32) \right)\nonumber\\ & \qquad + \sum_\gamma \left(\int \,{\rm d}(3) \mathcal{V}_{\gamma\beta}(32)\,f_\beta(2)\,g_{\alpha\gamma}(13) \right)\nonumber\\ & \quad ={-}\mathcal{V}_{\alpha\beta}(12)\, \left[f_\alpha(1)f_\beta(2) + g_{\alpha\beta}(12)\right]\nonumber\\ & \qquad- \sum_\gamma \left(\int \,{\rm d}(3)\, \left[\mathcal{V}_{\alpha\gamma}(13) + \mathcal{V}_{\beta\gamma}(23)\right]\,h_{\alpha\beta\gamma}(123) \right). \end{align}

In Appendix A, we give the definitions of the distribution functions $f_\alpha$, $g_{\alpha \beta }$ and $h_{\alpha \beta \gamma }$, and a detailed explanation of our notation. In summary, if the position and velocity of particle $1$ are $\boldsymbol {r}_1$ and $\boldsymbol {v}_1$, respectively, we write $(1)$ for the phase-space coordinates $(\boldsymbol {r}_1, \boldsymbol {v}_1)$ and $\int \,{\rm d}(1) = \int \,{\rm d}\boldsymbol {r}_1\,{\rm d}\boldsymbol {v}_1$ for integration over these coordinates. When a function $f$ depends on the phase-space coordinates of multiple particles, we write $f(12\ldots )$.

The operators $\mathcal {V}$ and $\mathcal {L}$ are defined as follows. The symmetric operator $\mathcal {V}_{\alpha \beta }(12)$ describes the effect of Coulomb interactions between particles $1$ and $2$:

(3.31)\begin{equation} \mathcal{V}_{\alpha\beta}(12) = q_\alpha q_\beta \frac{\boldsymbol{r}_1 - \boldsymbol{r}_2}{|\boldsymbol{r}_1-\boldsymbol{r}_2|^3}\boldsymbol{\cdot}\left(\frac{1}{m_\alpha}\frac{\partial}{\partial\boldsymbol{v}_1} - \frac{1}{m_\beta}\frac{\partial}{\partial\boldsymbol{v}_2}\right). \end{equation}

The $\mathcal {L}(1)$ operator describes the evolution of a distribution function as particle $1$ moves in the mean electric field:

(3.32)\begin{equation} \mathcal{L}(1)a(1) = \left(\boldsymbol{v}_1\boldsymbol{\cdot}\frac{\partial}{\partial\boldsymbol{r}_1} + \sum_\beta\int \,{\rm d}(2) \mathcal{V}_{\alpha\beta}(12)f_\beta(2)\right)a(1), \end{equation}

for any function $a(1)$. Finally, $\mathcal {L}(12) = \mathcal {L}(1)+\mathcal {L}(2)$.

When a two-temperature plasma has a large plasma parameter $\varLambda \gg 1$, two simplifications to the BBGKY equations (3.29)–(3.30) are possible. First, we can neglect $g_{\alpha \beta }(12)$ compared with $f_\alpha (1)f_\beta (2)$ on the right-hand side of (3.30). This means we ignore initial correlations between particles as they interact with each other and produce new correlations, which is valid as long as the particles are not too close together.³ Second, terms involving the triple correlation $h_{\alpha \beta \gamma }(123)$ in (3.30) may be dropped, which means we neglect the effect of three-particle collisions.⁴

With these two simplifications, we have a closed pair of equations for $f_{\alpha }(1)$ and $g_{\alpha \beta }(12)$:

(3.33)\begin{align}& \left(\frac{\partial}{\partial t} + \mathcal{L}(1)\right) f_\alpha(1) ={-}\sum_\beta \left( \int \,{\rm d}(2) \mathcal{V}_{\alpha\beta}(12)g_{\alpha\beta}(12) \right), \end{align}

(3.34)\begin{align} & \left(\frac{\partial}{\partial t} + \mathcal{L}(12)\right)g_{\alpha\beta}(12) +\sum_\gamma \left(\int \,{\rm d}(3) \mathcal{V}_{\alpha\gamma}(13) f_\alpha(1)\,g_{\gamma\beta}(32) \right)\nonumber\\ & \qquad + \sum_\gamma \left(\int \,{\rm d}(3) \mathcal{V}_{\gamma\beta}(32)\,f_\beta(2)\,g_{\alpha\gamma}(13) \right)\nonumber\\ & \quad ={-}\mathcal{V}_{\alpha\beta}(12)\,f_\alpha(1)f_\beta(2). \end{align}

We now assume the plasma is homogeneous, which means $f_\alpha (1)$ is independent of position and $g_{\alpha \beta }(12)$ depends on $\boldsymbol {r}_1$ and $\boldsymbol {r}_2$ only through the relative separation ${\boldsymbol {r} = \boldsymbol {r}_1 - \boldsymbol {r}_2}$. Then the second term vanishes from the definition (3.32) of $\mathcal {L}(1)$. Taking the Fourier transform of (3.34) with respect to the relative separation $\boldsymbol {r}$ gives an equation for $\tilde {g}_{\alpha \beta }(\boldsymbol {k},\boldsymbol {v}_1, \boldsymbol {v}_2) = \tilde {g}_{\alpha \beta }(12)$,

(3.35)\begin{align} & \left(\frac{\partial}{\partial t} + \tilde{\mathcal{L}}(12)\right)\tilde{g}_{\alpha\beta}(12) + \sum_\gamma \left(\int \,{\rm d}[3] \tilde{\mathcal{V}}_{\alpha\gamma}(13)f_\alpha(1)\tilde{g}_{\gamma\beta}(32) \right)\nonumber\\ & \qquad + \sum_\gamma \left(\int \,{\rm d}[3] \tilde{\mathcal{V}}_{\gamma\beta}(32) f_\beta(2) \tilde{g}_{\alpha\gamma}(13) \right) =-\tilde{\mathcal{V}}_{\alpha\beta}(12) f_\alpha(1)f_\beta(2), \end{align}

where the integral terms are obtained in this form after a change of variables and application of the convolution theorem. Here,

(3.36)\begin{gather} \tilde{\mathcal{V}}_{\alpha\beta}(12) ={-} \frac{4{\rm \pi} q_\alpha q_\beta }{k^2}\,{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\left(\frac{1}{m_\alpha}\frac{\partial}{\partial\boldsymbol{v}_1} - \frac{1}{m_\beta}\frac{\partial}{\partial\boldsymbol{v}_2} \right), \end{gather}

(3.37)\begin{gather}\tilde{\mathcal{L}}(12)= {\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}_1 - \boldsymbol{v}_2) \end{gather}

are the Fourier-space versions of the operators introduced earlier and $\int \,{\rm d}[3] = \int \,{\rm d}\boldsymbol {v}_3$ (see Appendix A).

3.2.2 Electron–ion correlation

As before, we need to find a way to exploit the small mass ratio. The only equation that involves both $m_e$ and $m_i$ is the equation for the mixed pair correlation $\tilde {g}_{ei}$. We therefore begin with this equation by expanding the operators in $m_e/m_i\ll 1$. Estimating the sizes of the various terms in (3.35), we have

(3.38)\begin{equation} \tilde{\mathcal{L}}(12)\,\tilde{g}_{ei}(12) = {\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}_1 - \boldsymbol{v}_2)\,\tilde{g}_{ei}(12) \approx {\rm i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}_1 \tilde{g}_{ei}(12) \end{equation}

and

(3.39)\begin{equation} \int \,{\rm d}[2] \tilde{\mathcal{V}}_{\alpha\beta}(12) f_\alpha(1) \sim \frac{q_\alpha q_\beta}{ k_D\sqrt{m_\alpha T_\alpha}}f_\alpha(1). \end{equation}

The second estimate means

(3.40)\begin{equation} \int \,{\rm d}[3] \tilde{\mathcal{V}}_{e\gamma}(13) f_e(1) \tilde{g}_{i\gamma}(32)\gg \int \,{\rm d}[3] \tilde{\mathcal{V}}_{\gamma i}(32) f_i(2) \tilde{g}_{e\gamma}(13), \end{equation}

because $m_eT_e \ll m_iT_i$ by assumption (2.4). Writing out the operators explicitly and dividing by ${\rm i}\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {v}_1$, we find

(3.41)\begin{equation} \tilde{g}_{ei}(12) + \frac{4{\rm \pi} e^2}{T_e}\frac{1}{k^2}f_e(1)\int \,{\rm d}[3] [\tilde{g}_{ei}(32) - Z \tilde{g}_{ii}(32)] = \frac{4 {\rm \pi}Ze^2}{T_e}\frac{1}{k^2}f_e(1)f_i(2). \end{equation}

Integrating ${\rm d}[1]$ then gives an expression for $\int \, {\rm d}[3] \tilde {g}_{ei}(32)$, which allows us to solve (3.41) for $\tilde {g}_{ei}(12)$ at steady state. We obtain

(3.42)\begin{equation} \tilde{g}_{ei}(12) = \frac{1}{n_e} \frac{Zk_e^2}{k^2+k_e^2}f_e(1)\left(f_i(2) + \int\, {\rm d}[3]\, \tilde{g}_{ii}(32) \right), \end{equation}

which relates the electron–ion correlation to the ion–ion correlation. We will now use this relation to eliminate $\tilde {g}_{ei}$ from the equations for $\tilde {g}_{ee}$ and $\tilde {g}_{ii}$.

3.2.3 Electron–electron and ion–ion correlations

The electron–electron and ion–ion correlations both satisfy an equation of the form

(3.43)\begin{align} & \left(\frac{\partial}{\partial t} + \tilde{\mathcal{L}}(12)\right)\tilde{g}_{\alpha\alpha}(12) + \sum_\gamma \left(\int\, {\rm d}[3] \tilde{\mathcal{V}}_{\alpha\gamma}(13) f_\alpha(1) \tilde{g}_{\gamma \alpha}(32) \right)\nonumber\\ & \qquad + \sum_\gamma \left(\int\, {\rm d}[3] \tilde{\mathcal{V}}_{\gamma \alpha}(32) f_\alpha(2) \tilde{g}_{\alpha\gamma}(13) \right) =-\tilde{\mathcal{V}}_{\alpha\alpha}(12)\,f_\alpha(1)f_\alpha(2). \end{align}

Writing the operators out explicitly, at steady state, we have

(3.44)\begin{align} & {\rm i} \boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}_1-\boldsymbol{v}_2)\tilde{g}_{\alpha\alpha}(12) + \sum_\gamma\frac{4{\rm \pi} q_\alpha q_\gamma}{T_\alpha}\frac{{\rm i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}_1}{k^2}f_\alpha(1)\int\, {\rm d}[3] \tilde{g}_{\gamma\alpha}(32)\nonumber\\ & \qquad - \sum_\gamma\frac{4{\rm \pi} q_\alpha q_\gamma}{T_\alpha}\frac{{\rm i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}_2}{k^2}f_\alpha(2)\int\, {\rm d}[3] \tilde{g}_{\alpha\gamma}(13)\nonumber\\ & \quad ={-}\frac{4 {\rm \pi}q_\alpha^2}{T_\alpha}\frac{1}{k^2}\,{\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}_1 - \boldsymbol{v}_2)f_\alpha(1)f_\alpha(2). \end{align}

This equation can easily be solved if we assume that $\tilde {g}_{\alpha \beta }(\boldsymbol {k},\boldsymbol {v}_1,\boldsymbol {v}_2)$ depends on the wavevector $\boldsymbol {k}$ only through its magnitude $k$. For a method that avoids this assumption, which proves the solution is unique, see Seuferling (1987). Then, for (3.44) to hold for all $\boldsymbol {k}$ and for all $\boldsymbol {v}_1$ and $\boldsymbol {v}_2$, we must have

(3.45)\begin{equation} \tilde{g}_{\alpha\alpha}(12) + \sum_\gamma\frac{4{\rm \pi} q_\alpha q_\gamma}{T_\alpha}\frac{1}{k^2}f_\alpha(1)\int\, {\rm d}[3] \tilde{g}_{\gamma\alpha}(32) ={-}\frac{4 {\rm \pi}q_\alpha^2}{T_\alpha}\frac{1}{k^2}f_\alpha(1)f_\alpha(2). \end{equation}

Equation (3.45) is almost identical to (3.41), which determined $\tilde {g}_{ei}$. The method used to solve that equation can be applied here, leading to

(3.46)\begin{equation} \tilde{g}_{ee}(12) ={-}\frac{1}{n_e} \frac{k_e^2}{k^2+k_e^2}f_e(1)\left(f_e(2) - Z\int\, {\rm d}[3] \tilde{g}_{ie}(32) \right) \end{equation}

for the electrons and

(3.47)\begin{equation} \tilde{g}_{ii}(12) ={-}\frac{1}{n_i} \frac{k_i^2}{k^2+k_i^2}f_i(1)\left(f_i(2) - \frac{1}{Z}\int\, {\rm d}[3] \tilde{g}_{ei}(32) \right) \end{equation}

for the ions. Using our solution (3.42) for $\tilde {g}_{ei}$ in terms of $\tilde {g}_{ii}$, together with ${\tilde {g}_{ie}(\boldsymbol {k},\boldsymbol {v}_1,\boldsymbol {v}_2) = \tilde {g}_{ei}(-\boldsymbol {k},\boldsymbol {v}_2,\boldsymbol {v}_1})$, gives

(3.48)\begin{equation} \tilde{g}_{ee}(12) = \left[\frac{1}{n_e^2}\left( \frac{Zk_e^2}{k^2+k_e^2} \right)^2 \left( n_i + \int\, {\rm d}[34] \tilde{g}_{ii}(34) \right) - \frac{1}{n_e}\frac{k_e^2}{k^2+k_e^2}\right]f_e(1)f_e(2). \end{equation}

Here we can clearly see that the electron–electron correlation consists of two parts, one due to their interactions with the ions and one that is independent of the ions, as explained in § 3.1.

Next we will calculate $\tilde {g}_{ii}$, but it is first worth highlighting what has been achieved so far. Given any model for the ion–ion correlation, (3.42) and (3.48) determine the electron–ion and electron–electron correlations. These relations hold as long as the electrons interact weakly, so they constrain the electron correlations in a plasma with weakly or strongly coupled ions. This is possible because the electrons form completely symmetric, identical screening clouds around the ions, so we expect that knowing how the ions are distributed should tell us how the electrons are distributed. Boercker & More (1986) obtained the same relations by using the electron Helmholtz free energy to construct an effective ion–ion interaction; this procedure is rigorously justified by Rosenfeld (1994).

We now solve (3.47), which assumes the ions interact weakly, for $\tilde {g}_{ii}$. Eliminating $\tilde {g}_{ei}$ using (3.42) gives

(3.49)\begin{equation} \tilde{g}_{ii}(12) ={-}\frac{1}{n_i} \frac{k_i^2}{k^2+k_i^2}f_i(1)\left(\frac{k^2}{k^2+k_e^2}f_i(2) - \frac{k_e^2}{k^2+k_e^2}\int\, {\rm d}[3] \tilde{g}_{ii}(32) \right). \end{equation}

Integrating ${\rm d}[1]$ and substituting into (3.49) leads to

(3.50)\begin{equation} \tilde{g}_{ii}(12) ={-}\frac{1}{n_i}\frac{k_i^2}{k^2+k_D^2}f_i(1)f_i(2). \end{equation}

We have calculated the ion–ion correlation, which now determines both the other correlations via (3.42) and (3.48). The full set of pair correlations is

(3.51)\begin{gather} \tilde{g}_{ii}(12) ={-}\frac{4{\rm \pi} Z^2 e^2}{T_i}\frac{1}{k^2+k_D^2}f_i(1)f_i(2), \end{gather}

(3.52)\begin{gather}\tilde{g}_{ei}(12) = \frac{4{\rm \pi} Z e^2}{T_e}\frac{1}{k^2+k_D^2}f_e(1)f_i(2), \end{gather}

(3.53)\begin{gather}\tilde{g}_{ee}(12) ={-}\frac{4{\rm \pi} e^2}{T_e}\left[\frac{T_i}{T_e}\frac{1}{k^2+k_D^2} + \left( \frac{T_e-T_i}{T_e} \right)\frac{1}{k^2+k_e^2}\right]f_e(1)f_e(2). \end{gather}

Of course, these are consistent with (3.26)–(3.28) upon integrating over velocities and inverse-Fourier transforming. We now have the spatial and velocity parts of the pair correlations; the velocity dependence enters through a product of Maxwellians, as for systems at thermal equilibrium.

4 Response of plasma to sudden electron heating

4.1 Correlation heating mechanism

The main objective of this paper is to answer the following question. Suppose some energy $Q$ per particle is instantaneously supplied to the electrons in a moderately coupled two-temperature plasma. The hotter electrons will now shield the ions less; what effect does this have on the ion temperature?

The answer is there is a small but rapid increase in the ion temperature, on the ion-plasma-frequency time scale $t\sim 1/\omega _{pi}$. To understand why, consider the extreme case of very cold electrons which screen the ions perfectly. The ions will be distributed randomly, as they do not interact. Suppose the ions are also very cold, so they are essentially stationary.

Now imagine the electrons are suddenly heated, which means their Debye length increases. As a result, the ions are not so well screened and will repel each other. This acceleration increases their average kinetic energy. The effect is particularly noticeable for pairs, or larger clusters, of ions that were initially quite close together: the sudden removal of their shielding clouds means they repel each other and fly apart. The ions were initially stationary but have now been given some kinetic energy, so their temperature must have increased. This mechanism is most effective at producing fast ions if the electron shielding is stripped away rapidly.

Disorder-induced heating of ultracold neutral plasmas has exactly the same underlying mechanism, except the starting point is a neutral gas rather than well-shielded ions (see e.g. Acciarri, Moore & Baalrud 2022). In this context, a common explanation is that the ions transition from an initial state of random, uncorrelated positions to a final state with correlations imposed by their Coulomb repulsion. This is depicted in figure 2; their repulsion after the electrons are heated means the ions push away from each other and are generally not found close together any more. This means their spatial correlation has changed and has become more ordered, leading to heating; the time scale for the ion–ion correlation to evolve is $t\sim 1/\omega _{pi}$, which is the time it takes a thermal ion to traverse a Debye length.

Figure 2. Cartoon illustration of ion positions before and after ionisation and disorder-induced heating. (a) Before heating: random, uncorrelated positions. (b) After heating: repulsion leads to more ordered, correlated positions.

This temperature increase can also be understood as a consequence of the Second Law of Thermodynamics; it is required to compensate for the decrease in ion entropy when they become more spatially ordered. For an entropic perspective on correlation heating, see Fetsch et al. (2023).

4.2 Multiple time scales

The response of the plasma occurs in stages over multiple time scales.

Before the heating at $t=0$, the pair correlations in Fourier space are given by (3.51)–(3.53):

(4.1)\begin{gather} \tilde{g}_{ii}(12) ={-}\frac{4{\rm \pi} Z^2 e^2}{T_{i0}}\frac{1}{k^2+k_{D0}^2}f_{i0}(1)f_{i0}(2), \end{gather}

(4.2)\begin{gather}\tilde{g}_{ei}(12) = \frac{4{\rm \pi} Z e^2}{T_{e0}}\frac{1}{k^2+k_{D0}^2}f_{e0}(1)f_{i0}(2), \end{gather}

(4.3)\begin{gather}\tilde{g}_{ee}(12) ={-}\frac{4{\rm \pi} e^2}{T_{e0}}\left[\frac{T_{i0}}{T_{e0}}\frac{1}{k^2+k_{D0}^2} + \left( \frac{T_{e0}-T_{i0}}{T_{e0}} \right)\frac{1}{k^2+k_{e0}^2}\right]f_{e0}(1)f_{e0}(2), \end{gather}

where the ‘$0$’ subscript denotes the initial value of a quantity. This includes the initial electron and ion velocity distributions, which are $f_{e0}$ and $f_{i0}$, respectively.

First, the electrons are instantaneously given a larger speed. Their velocity distribution will change in a way that depends on the exact mechanism of the heating, but the spatial pair correlations cannot change because this instantaneous heating takes place with fixed particle positions.

Next, the electron–electron and electron–ion correlations will adjust on an electron-plasma-frequency time scale $t\sim 1/\omega _{pe}$. The increase in temperature makes the electron Debye length larger, so the electron shielding clouds must rapidly expand.

Then, over the electron–electron collision time scale $t\sim \tau _{ee}$, the electron velocity distribution will become Maxwellian again, with new temperature $T_{e1}$. The electron correlations $g_{ee}$ and $g_{ei}$ will continually adjust as this occurs and will approach their steady-state solutions given constant $f_i$ and $g_{ii}$. At the end of this stage of the process, the electrons are Maxwellian, and the steady-state solutions for the correlations are given in (3.42) and (3.48):⁵ we must have

(4.4)\begin{gather} \tilde{g}_{ei}(12) = \frac{1}{n_e} \frac{Zk_{e1}^2}{k^2+k_{e1}^2}f_{e1}(1)\left(f_{i0}(2) + \int\, {\rm d}[3]\, \tilde{g}_{ii}(32) \right), \end{gather}

(4.5)\begin{gather}\tilde{g}_{ee}(12) = \left[\frac{1}{n_e^2}\left( \frac{Zk_{e1}^2}{k^2+k_{e1}^2} \right)^2 \left( n_i + \int\, {\rm d}[34] \tilde{g}_{ii}(34) \right) - \frac{1}{n_e}\frac{k_{e1}^2}{k^2+k_{e1}^2}\right]f_{e1}(1)f_{e1}(2), \end{gather}

where $k_{e1}$ is the new inverse-Debye length of the electrons at temperature $T_{e1}$ and $f_{e1}$ is their new velocity distribution. Note that these two formulae assumed Maxwellian electrons but did not require Maxwellian ions, so they do not involve an ion temperature.

The ion–ion correlation changes on the ion-plasma-frequency time scale $t\sim 1/\omega _{pi}$, so the ions have not yet had time to adjust. We still have

(4.6)\begin{equation} \tilde{g}_{ii}(12) ={-}\frac{4{\rm \pi} Z^2 e^2}{T_{i0}}\frac{1}{k^2+k_{D0}^2}f_{i0}(1)f_{i0}(2). \end{equation}

The assumption of instantaneous heating is an idealisation that is not actually necessary. Our calculations are valid as long as the electrons are heated fast enough that they can Maxwellianise before the ions respond.

Next, over the ion-plasma-frequency time scale, the ion–ion correlation will change. The ions are screened less by the hotter electrons, so they push each other further apart.

This is the critical time scale on which collisionless heating occurs. At any moment, just as in (4.4) and (4.5), the electron distribution $f_e$ and correlations $g_{ei}$ and $g_{ee}$ may be assumed to have relaxed to their steady-state solutions given constant $f_i$ and $g_{ii}$, because the electron evolution time scales ($1/\omega _{pe}$ and $\tau _{ee}$) are much shorter than the time scale of interest ($1/\omega _{pi}$).

Finally, on the ion–ion collisional time scale $t\sim \tau _{ii}$, the ion distribution will become Maxwellian again. This step is important because correlation heating causes the ion distribution to deviate from a Maxwellian. The correlations all adjust continually as this occurs, but we will show that there is no further change in the ion kinetic energy during this last stage. If the final electron and ion temperatures are $T_{e2}$ and $T_{i2}$, respectively, then at the end, we must have correlations in their steady-state form again:

(4.7)\begin{gather} \tilde{g}_{ii}(12) ={-}\frac{4{\rm \pi} Z^2 e^2}{T_{i2}}\frac{1}{k^2+k_{D2}^2}f_{i2}(1)f_{i2}(2), \end{gather}

(4.8)\begin{gather}\tilde{g}_{ei}(12) = \frac{4{\rm \pi} Z e^2}{T_{e2}}\frac{1}{k^2+k_{D2}^2}f_{e2}(1)f_{i2}(2), \end{gather}

(4.9)\begin{gather}\tilde{g}_{ee}(12) ={-}\frac{4{\rm \pi} e^2}{T_{e2}}\left[\frac{T_{i2}}{T_{e2}}\frac{1}{k^2+k_{D2}^2} + \left( \frac{T_{e2}-T_{i2}}{T_{e2}} \right)\frac{1}{k^2+k_{e2}^2}\right]f_{e2}(1)f_{e2}(2). \end{gather}

There is also a much longer time scale of temperature equilibration. However, we are not interested in processes that occur over such long times. Figure 3 summarises the entire correlation heating process.

Figure 3. Timeline of the system's response to sudden electron heating.

4.3 Ion evolution

To calculate how much the ion temperature increases, we start with (3.33) for the ion distribution, repeated here:

(4.10)\begin{equation} \left(\frac{\partial}{\partial t} + \mathcal{L}(1)\right) f_i(1) ={-}\sum_\beta \left( \int\, {\rm d}(2) \mathcal{V}_{\beta i}(21)\,g_{\beta i}(21) \right). \end{equation}

On the right-hand side, particles $1$ and $2$ and their corresponding species labels were swapped for convenience. In a homogeneous plasma this is equivalent to

(4.11)\begin{equation} \frac{\partial f_i(1)}{\partial t} ={-}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\int\, {\rm d}[2] \left(\tilde{\mathcal{V}}_{ei}(21)^*\,\tilde{g}_{ei}(21) + \tilde{\mathcal{V}}_{ii}(21)^*\,\tilde{g}_{ii}(21) \right). \end{equation}

As explained, we use the steady-state solution (3.42) for $g_{ei}$ given Maxwellian electrons and fixed $f_i$ and $g_{ii}$, together with isotropy in the $\boldsymbol {k}$-space integration, to eliminate the electrons from this equation. We obtain

(4.12)\begin{equation} \frac{\partial f_i(1)}{\partial t} ={-}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\int\, {\rm d}[2] \left( \frac{4{\rm \pi} Z^2e^2}{m_i}\frac{{\rm i}\boldsymbol{k}}{k^2+k_e^2}\boldsymbol{\cdot}\frac{\partial}{\partial\boldsymbol{v}_1}\tilde{g}_{ii}(21) \right). \end{equation}

In real space, this reads

(4.13)\begin{equation} \frac{\partial f_i(1)}{\partial t} ={-} \int\, {\rm d}(2)\, \mathcal{V}^s(12)\,g_{ii}(12), \end{equation}

where we have defined an effective interaction operator

(4.14)\begin{equation} \mathcal{V}^s(12)={-}\frac{1}{m_i}\frac{\partial}{\partial(\boldsymbol{r}_1 - \boldsymbol{r}_2)}\left( \frac{Z^2e^2}{|\boldsymbol{r}_1 - \boldsymbol{r}_2|}\,{\rm e}^{{-}k_e|\boldsymbol{r}_1 - \boldsymbol{r}_2}| \right) \boldsymbol{\cdot}\left(\frac{\partial}{\partial\boldsymbol{v}_1} - \frac{\partial}{\partial\boldsymbol{v}_2}\right). \end{equation}

Ion–ion correlations drive changes in the ion distribution function as if the ions were a gas of particles interacting via the screened Coulomb potential ${\varphi ^{(s)}(\boldsymbol {r}) = Z^2e^2 \exp (-k_er)/r}$.

The ion distribution function changes when the ion correlation has not yet reached its steady-state form. To calculate how much the ion kinetic energy changes, we will need to know how $g_{ii}$ evolves with time.

When the ions are weakly coupled, the BBGKY equation (3.35) for $\tilde {g}_{ii}$ is

(4.15)\begin{align} & \left(\frac{\partial}{\partial t} + \tilde{\mathcal{L}}(12)\right)\tilde{g}_{ii}(12) + \sum_\gamma \left(\int\, {\rm d}[3] \tilde{\mathcal{V}}_{i\gamma}(13)f_i(1)\tilde{g}_{\gamma i}(32) \right)\nonumber\\ & \qquad+ \sum_\gamma \left(\int\, {\rm d}[3] \tilde{\mathcal{V}}_{\gamma i}(32)f_i(2)\tilde{g}_{i\gamma}(13) \right) =-\tilde{\mathcal{V}}_{ii}(12)\,f_i(1)f_i(2). \end{align}

The electron–ion correlation appears in the three-particle interaction terms, for example,

(4.16)\begin{align} \sum_\gamma \left(\int\, {\rm d}[3] \tilde{\mathcal{V}}_{i\gamma}(13)\,f_i(1)\,\tilde{g}_{\gamma i}(32) \right) & ={-} \frac{4{\rm \pi} Z^2e^2}{m_i}\frac{{\rm i}\boldsymbol{k}}{k^2}\boldsymbol{\cdot}\frac{\partial f_i(1)}{\partial\boldsymbol{v}_1} \int\, {\rm d}[3]\,\tilde{g}_{ii}(32)\nonumber\\ & \quad + \frac{4{\rm \pi} Ze^2}{m_i}\frac{{\rm i}\boldsymbol{k}}{k^2}\boldsymbol{\cdot}\frac{\partial f_i(1)}{\partial\boldsymbol{v}_1} \int\, {\rm d}[3]\,\tilde{g}_{ei}(32). \end{align}

Substituting (4.4) in for $\tilde {g}_{ei}$ as before gives

(4.17)\begin{align} \sum_\gamma \left(\int\, {\rm d}[3] \tilde{\mathcal{V}}_{i\gamma}(13)f_i(1)\tilde{g}_{\gamma i}(32) \right) & ={-} \frac{4{\rm \pi} Z^2e^2}{m_i}\frac{{\rm i}\boldsymbol{k}}{k^2+k_e^2}\boldsymbol{\cdot}\frac{\partial f_i(1)}{\partial\boldsymbol{v}_1} \int\, {\rm d}[3]\,\tilde{g}_{ii}(32)\nonumber\\ & \quad + \frac{4{\rm \pi} Z^2e^2}{m_i}\frac{k_e^2}{k^2+k_e^2}\frac{{\rm i}\boldsymbol{k}}{k^2}\boldsymbol{\cdot}\frac{\partial f_i(1)}{\partial\boldsymbol{v}_1} f_i(2). \end{align}

Using the same substitution in the other terms, we therefore find

(4.18)\begin{align} & \left(\frac{\partial}{\partial t} + \tilde{\mathcal{L}}(12)\right)\tilde{g}_{ii}(12) - \frac{4{\rm \pi} Z^2e^2}{m_i}\frac{{\rm i}\boldsymbol{k}}{k^2+k_e^2}\boldsymbol{\cdot}\frac{\partial f_i(1)}{\partial\boldsymbol{v}_1} \int\, {\rm d}[3]\tilde{g}_{ii}(32)\nonumber\\ & \qquad- \frac{4{\rm \pi} Z^2e^2}{m_i}\frac{{\rm i}\boldsymbol{k}}{k^2+k_e^2}\boldsymbol{\cdot}\frac{\partial f_i(2)}{\partial\boldsymbol{v}_2} \int\, {\rm d}[3]\;\tilde{g}_{ii}(13)\nonumber\\ & \quad ={+} \frac{4{\rm \pi} Z^2e^2}{m_i}\frac{{\rm i}\boldsymbol{k}}{k^2+k_e^2}\boldsymbol{\cdot}\left(\frac{\partial }{\partial\boldsymbol{v}_1} - \frac{\partial }{\partial\boldsymbol{v}_2}\right) f_i(1)f_i(2). \end{align}

Returning to real space, we now have a system of equations determining the evolution of $f_i$ and $g_{ii}$:

(4.19)\begin{equation} \frac{\partial f_i(1)}{\partial t} ={-} \int\, {\rm d}(2) \mathcal{V}^s(12)g_{ii}(12), \end{equation}

(4.20)\begin{align} & \left(\frac{\partial}{\partial t} + \mathcal{L}(12)\right)g_{ii}(12) + \int\, {\rm d}(3) \mathcal{V}^s(13)\,f_i(1)\,g_{ii}(32)\nonumber\\ & \qquad + \int\, {\rm d}(3) \mathcal{V}^s(32)f_i(2)g_{ii}(13) ={-}\mathcal{V}^s(12)f_i(1)f_i(2). \end{align}

The ion distribution function and ion–ion pair correlation evolve according to (4.19) and (4.20). The electrons do not enter these equations anywhere, except through the temperature $T_e$ which controls the screening length $k_e$. These equations are precisely the equations that would describe a gas of particles interacting via screened Yukawa potentials, although here the screening length $k_e$ is also a dynamical variable that changes with time. The ions therefore behave like a weakly coupled Yukawa one-component plasma (YOCP). While it is common to model the ions as a YOCP in theoretical and computational work (e.g. Hamaguchi, Farouki & Dubin 1997; Murillo 2001; Gericke et al. 2003; Murillo 2009; Lyon, Bergeson & Murillo 2013; Langin et al. 2016; Sprenkle et al. 2022), we are not aware of any other rigorous demonstration, starting from the BBGKY hierarchy, that this approximation is valid for a two-temperature plasma. This result, which is consistent with the physical picture developed in § 3.1, will allow us to calculate the temperature changes during correlation heating in a very simple way.

5 Temperature changes due to correlation heating

In this section, we calculate how the temperatures of each species change in response to a sudden input of energy to the electrons.

First, we review the sequence of temperatures defined in § 4.2. The electron and ion temperatures start out at $T_{e0}$ and $T_{i0}$, respectively, before each electron's kinetic energy is suddenly increased by $Q$. In general, the electron distribution may not be Maxwellian immediately after this heating; however, their total kinetic energy must change from $3N_eT_{e0}/2$ to $3N_eT_{e0}/2 + N_e Q$. Next, the electron shielding clouds expand and they reach a new Maxwellian distribution at temperature $T_{e1}$, all while the ion temperature remains $T_{i0}$. Then, the ions move further apart (the ion–ion correlation changes) and their kinetic energy increases due to correlation heating. Eventually, the ions re-Maxwellianise, and the final temperatures are $T_{e2}$ and $T_{i2}$.

We expect the changes in particle correlations to cause temperature changes that are small in $1/\varLambda$ because this heating relies on conversion of potential energy into kinetic energy, and the potential energy in a moderately coupled plasma is a factor of $1/\varLambda$ smaller than the kinetic energies. This is consistent with the fact that, for an ideal plasma with $\varLambda \to \infty$, the species temperatures are completely decoupled and there should be no correlation heating. So, $T_{e1}$ and $T_{e2}$ are close to $T_{e0}+2Q/3$, while $T_{i2}$ is close to $T_{i0}$.

5.1 Conservation of energy

The BBGKY equations for a system of particles with interaction potential $\varphi _{\alpha \beta }$ between species $\alpha$ and species $\beta$ conserve energy. The total plasma energy $U$ is the sum of the kinetic energy of each species and the potential energy, which depends on the relative positions of the particles and must therefore be calculated from the pair correlations. It is given by (Landau & Lifshitz 1980)

(5.1)\begin{equation} U = \sum_\alpha V\int\, {\rm d}\boldsymbol{v}\frac{1}{2}m_\alpha\boldsymbol{v}^2 f_\alpha(\boldsymbol{v}) + \frac{1}{2}\sum_{\alpha}\sum_{\beta} Vn_\alpha n_\beta\int\, {\rm d}\boldsymbol{r}\, \varphi_{\alpha\beta}(\boldsymbol{r})\, G_{\alpha\beta}(\boldsymbol{r}), \end{equation}

where $V$ is the system volume. For example, the energy of a steady-state two-temperature plasma, which is found using the pair correlations (3.26)–(3.28) and Coulomb potential $\varphi _{\alpha \beta }(\boldsymbol {r}) = q_\alpha q_\beta /r$, is

(5.2)\begin{equation} U = \frac{3}{2}N_eT_e + \frac{3}{2}N_iT_i -\frac{V}{8{\rm \pi}}\left[T_i(k_D^3 - k_e^3) + T_ek_e^3\right]. \end{equation}

This calculation was carried out by More (1980), although it seems the correct formula (5.2) was first published by Triola (2022). An alternative derivation was provided by Fetsch et al. (2023). Note that when $Z\sim 1$ and $T_e\sim T_i$, the potential energy correction is ${\sim NT/\varLambda}$, as expected.

5.2 Ion temperature $T_{i2}$

A central aim of this paper is to calculate the ion temperature after correlation heating. It is not possible to find $T_{i2}$ just using conservation of the total plasma energy, because the final energy depends on two unknowns, $T_{e2}$ and $T_{i2}$; one conservation equation is not enough to determine both temperatures.

However, at the end of § 4, we found that the ion distribution function and ion–ion pair correlation satisfy the same equations as a weakly coupled YOCP with screening length $k_e$. Over the course of the ion correlation adjustment and subsequent re-Maxwellianisation, the electron temperature changes from $T_{e1}$ to $T_{e2}$, so this screening length is not constant. However, the variation in $k_e$ is small in $1/\varLambda$ and is therefore a higher-order effect which can be neglected for the purposes of obtaining a leading-order expression for the ion heating.

Thus, the screening length can be treated as fixed, $k_e=k_{e1}$. The ion equations (4.19) and (4.20) then have a conserved energy,

(5.3)\begin{equation} E = \int\, {\rm d}\boldsymbol{v} \frac{1}{2}\,m_i\boldsymbol{v}^2 f_i(\boldsymbol{v}) + \frac{1}{2}Vn_i^2\int\, {\rm d}\boldsymbol{r} \left(\frac{Z^2e^2}{r}\,{\rm e}^{{-}k_{e1}r}\right) G_{ii}(\boldsymbol{r}), \end{equation}

which is the energy the ions would have if they were a YOCP. Initially, before the ion correlation adjusts, we have $G_{ii}(\boldsymbol {r},t=0)=-(Z^2e^2/T_{i0})({\rm e}^{-k_{D0}r}/r)$ and the energy is

(5.4)\begin{equation} E_1 = \frac{3}{2}N_iT_{i0} - \frac{1}{8{\rm \pi}}VT_{i0}\frac{k_{i0}^4}{k_{D0}+k_{e1}}. \end{equation}

After correlation heating, the ion correlation is $G_{ii}(\boldsymbol {r},t=\infty ) = -(Z^2e^2/T_{i2})({\rm e}^{-k_{D2}r}/r)$, which is equivalent to $G_{ii}(\boldsymbol {r},t=\infty ) = -(Z^2e^2/T_{i0})({\rm e}^{-k_{D1}r}/r)$ up to corrections small in $1/\varLambda$. The final energy is therefore

(5.5)\begin{equation} E_2 = \frac{3}{2}N_iT_{i2} - \frac{1}{8{\rm \pi}}VT_{i0}\frac{k_{i0}^4}{k_{D1}+k_{e1}}. \end{equation}

Equating $E_1$ and $E_2$ allows us to determine the change in ion temperature, ${T_{i2} - T_{i1} = \Delta T_i}$:

(5.6)\begin{equation} \Delta T_i = \frac{Z^2e^2}{3}k_{i0}^2\left(\frac{1}{k_{D1}+k_{e1}} - \frac{1}{k_{D0}+k_{e1}} \right). \end{equation}

This is one of our main results.

1. (i) When the electrons are not heated at all, $k_{D1}=k_{D0}$ and so $\Delta T_i = 0$. Nothing happens.

2. (ii) If the electrons are heated, then $k_{e1}< k_{e0}$. In this case, $\Delta T_i > 0$ and the ions also heat up. However, if the electrons are suddenly cooled, then the ions also cool.

3. (iii) As expected, the change in the ion temperature is small in a moderately coupled plasma because $\Delta T_i \sim e^2 k_D \sim T/\varLambda$. However, the temperature change is not small in the mass ratio $m_e/m_i$; it is non-zero, even though we took $m_e/m_i \to 0$ in the derivation.

4. (iv) If the electrons become much hotter than the ions, meaning $T_e \gg T_i$, then the ion heating approaches a finite limit that depends on the initial conditions:

   (5.7)\begin{equation} \Delta T_i = \frac{Z^2e^2}{3}k_{i0}^2\left( \frac{1}{k_{i0}} - \frac{1}{k_{D0}} \right). \end{equation}

5. (v) However, if the electrons are suddenly cooled so that $T_e\ll T_i$, then ${\Delta T_i \to 0}$. Physically, the reason is that cooling the electrons causes them to form small shielding clouds around each of the ions, which become perfectly screened. Then, the ions no longer interact and so must continue moving in a straight line at a constant speed; there is no heating.

Expression (5.6) is plotted in figure 4 for different ratios of initial electron temperature to initial ion temperature. The figure shows that $\Delta T_i$ approaches a constant, maximum achievable value when the electron heating is large, and that it vanishes when the electrons are cooled to zero temperature. It shows that correlation heating has the greatest effect on the ions when the electrons are started at a low temperature.

Figure 4. Plot of fractional change in ion temperature $\Delta T_i/T_{i0}$ against the change in electron temperature, for three different initial electron temperatures. Here, $Z=1$ and we normalise using ${\varLambda _i = n_i / k_{i0}^3}$, which means the actual temperature change is roughly a factor of the plasma parameter smaller than the values on these curves. Note that the curves do not continue to arbitrarily negative abscissae because $T_{e1}$ cannot be reduced below zero.

5.3 Electron temperature $T_{e1}$

We now work out how the electron temperature changes. This occurs in two stages: first, they are suddenly heated and they re-Maxwellianise before the ions respond. At this point, the electron temperature is $T_{e1}$. Then, once the ion correlations have adjusted, the electrons will reach a new temperature $T_{e2}$.

The potential energy is easiest to calculate using Fourier-space pair correlations. The energy formula (5.1) is equivalent to

(5.8)\begin{equation} U = \sum_\alpha V\int\, {\rm d}\boldsymbol{v}\frac{1}{2}m_\alpha\boldsymbol{v}^2 f_\alpha(\boldsymbol{v}) + \frac{1}{2}\sum_{\alpha}\sum_{\beta} Vn_\alpha n_\beta\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3} \tilde{\varphi}_{\alpha\beta}(\boldsymbol{k}) \tilde{G}_{\alpha\beta}(\boldsymbol{k}), \end{equation}

with $\tilde {\varphi }_{\alpha \beta }(\boldsymbol {k}) = 4{\rm \pi} q_\alpha q_\beta / k^2$ for the Coulomb potential. This expression involves all three correlations $G_{ee}$, $G_{ei}$ and $G_{ii}$ in the potential energy. However, we have already seen that the electron–electron and electron–ion correlations are determined by the slowly varying ion correlation. Assuming each species is Maxwellian, we now use (4.4) and (4.5) integrated over velocities,

(5.9)\begin{gather} \tilde{G}_{ei}(\boldsymbol{k}) = \frac{1}{n_i}\frac{k_e^2}{k^2+k_e^2}\left( 1 + n_i\tilde{G}_{ii}(\boldsymbol{k}) \right), \end{gather}

(5.10)\begin{gather}\tilde{G}_{ee}(\boldsymbol{k}) = \frac{1}{n_i} \left(\frac{k_e^2}{k^2+k_e^2}\right)^2 \left( 1 + n_i\tilde{G}_{ii}(\boldsymbol{k}) \right) - \frac{1}{n_e}\frac{k_e^2}{k^2+k_e^2}, \end{gather}

to express the potential energy in terms of the ion correlation only. We find

(5.11)\begin{equation} U = \frac{3}{2}N_eT_e + \frac{3}{2}N_iT_i - \frac{1}{2}N_e e^2 \left(\frac{3Z}{2}+1\right)k_e+ 2{\rm \pi} n_i^2 VZ^2e^2\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{k^2\tilde{G}_{ii}(\boldsymbol{k})}{(k^2+k_e^2)^2}. \end{equation}

This formula for the energy of a two-temperature plasma with arbitrary ion correlation is interesting in its own right. It uses the fact that the electrons interact weakly but does not assume the same about the ions, as explained in § 3.2.3.

As previously discussed, the ions in a two-temperature plasma are often approximated using a YOCP model (e.g. Hamaguchi et al. 1997; Murillo 2001; Gericke et al. 2003; Murillo 2009; Lyon et al. 2013; Langin et al. 2016; Sprenkle et al. 2022). When such a plasma is not at steady state, the electron temperature does not necessarily remain constant, so the screening length is a dynamical variable which should be evolved along with the ion distribution function. Conservation of energy, using the potential energy in (5.11), allows the electron temperature to be calculated at any given time from the ion–ion correlation, without reference to the fast electron dynamics. The fact that the effective interaction range of the ions is a dynamical variable that changes with time could be important when studying strongly coupled plasmas.

Now using the initial ion correlation $G_{ii}(\boldsymbol {r},t=0) = -(Z^2e^2/T_{i0})({\rm e}^{-k_{D0}r}/r)$ in (5.11), we find the initial energy is

(5.12)\begin{equation} U_0 = \frac{3}{2}N_e T_{e0} + \frac{3}{2}N_i T_{i0} - \frac{1}{2}N_ee^2\left( \frac{3Z}{2} + 1 \right)k_{e0} - \frac{1}{4}N_iZ^2e^2k_{i0}^2\frac{2k_{D0}+k_{e0}}{(k_{D0}+k_{e0})^2}, \end{equation}

and the energy after electron heating, but before the ion correlations adjust, is

(5.13)\begin{equation} U_1 = \frac{3}{2}N_e T_{e1} + \frac{3}{2}N_i T_{i0} - \frac{1}{2}N_ee^2\left( \frac{3Z}{2} + 1 \right)k_{e1} - \frac{1}{4}N_iZ^2e^2k_{i0}^2\frac{2k_{D0}+k_{e1}}{(k_{D0}+k_{e1})^2}. \end{equation}

In the potential energy term, $k_{e1}$ depends on $T_{e1}$, which is the unknown that we are trying to find. However, the potential energy term is small in $1/\varLambda$, so we can use the leading order formula

(5.14)\begin{equation} k_{e1} = \left(\frac{4{\rm \pi} n_e e^2}{T_{e0}+2Q/3}\right)^{1/2}, \end{equation}

since the difference between $T_{e1}$ and $T_{e0}+2Q/3$ is higher-order in $1/\varLambda$.

We find the unknown temperature $T_{e1}$ using conservation of energy. Equating the final energy after the ion adjustment to the energy just after heating, $U_1 = U_0 + N_e Q$, gives

(5.15)\begin{align} T_{e1} - T_{e0} & = \frac{2}{3} Q - \frac{1}{2}N_ee^2\left( \frac{3Z}{2} + 1 \right)(k_{e0}-k_{e1})\nonumber\\ & \quad+ \frac{1}{4}N_iZ^2e^2k_{i0}^2\left(\frac{2k_{D0}+k_{e1}}{(k_{D0}+k_{e1})^2} - \frac{2k_{D0}+k_{e0}}{(k_{D0}+k_{e0})^2} \right). \end{align}

Expression (5.15) is not particularly enlightening, but it can be used to prove that ${T_{e1}-T_{e0} < 2Q/3}$, so the electrons do not heat up quite as much as they would if there were an ideal gas (the opposite is true if they are cooled).

5.4 Electron temperature $T_{e2}$

The last unknown temperature to calculate is the final electron temperature $T_{e2}$. Now that $T_{i2}$ is known, conservation of the total plasma energy can be used to calculate $T_{e2}$.

Using (5.11) again, the final energy after the ions have adjusted is

(5.16)\begin{equation} U_2 = \frac{3}{2}N_e T_{e2} + \frac{3}{2}N_i T_{i2} - \frac{1}{2}N_ee^2\left( \frac{3Z}{2} + 1 \right)k_{e2} - \frac{1}{4}N_iZ^2e^2k_{i2}^2\frac{2k_{D2}+k_{e2}}{(k_{D2}+k_{e2})^2}, \end{equation}

where $k_{D2}^2 = k_{i2}^2 + k_{e2}^2$. Approximating $k_{i2}=k_{i0}$, $k_{e2} = k_{e1}$ and $k_{D2} = k_{D1}$ to leading order in the potential energy term, we equate $U_2 = U_1$ and obtain

(5.17)\begin{equation} T_{e2}-T_{e1} ={-}\frac{Z^2e^2}{6}k_{i0}^2k_{e1}\left(\frac{1}{(k_{D1}+k_{e1})^2} - \frac{1}{(k_{D0}+k_{e1})^2} \right). \end{equation}

If $k_{e1}< k_{e0}$, so the electrons are heated initially, then this expression is negative and the electron temperature drops back down slightly (that is, by an amount small in $1/\varLambda$).

In figure 5, we summarise the temperature changes that occur during correlation heating. If the electrons are initially heated, then the ions heat up due to correlation heating, while the electrons cool during the process. The opposite is true if the electrons are initially cooled; in either case, correlation heating causes the two temperatures to move towards each other.

Figure 5. Temperature changes for each species during correlation heating.

6 Time evolution of the ion kinetic energy

In the previous section, we calculated the change in ion temperature due to correlation heating, assuming the final state has pair correlations given by the Salpeter formulae for a steady-state, two-temperature plasma. In this section, we find the time dependence of the ion kinetic energy on the ion-plasma-frequency time scale $t\sim 1/\omega _{pi}$.

This calculation is interesting for two reasons. First, it shows that the ion heating occurs entirely on this fast time scale. This means that, when the ions re-Maxwellianise on the ion–ion collision time scale $t\sim \tau _{ii}$, they do not undergo any further heating.

Second, it allows us to investigate whether the kinetic energy oscillates or overshoots its final value. Such an overshoot could facilitate fusion reactions within a short window of time. Indeed, in strongly coupled plasmas, it has been observed in simulations (Morozov & Norman 2003; Pohl, Pattard & Rost 2004; Niffenegger, Gilmore & Robicheaux 2011; Sprenkle et al. 2022) and experiments (Chen et al. 2004; Laha et al. 2006; Lyon & Bergeson 2011; Langin et al. 2016) that the kinetic energy does oscillate about its final value if the correlations are suddenly modified. However, we are able to prove that, in a moderately coupled plasma, the kinetic energy may not be a monotonic function of time but it does not overshoot. A similar problem was solved for a YOCP by Morawetz et al. (2001) using the Kadanoff–Baym equations of quantum statistical theory.

To evolve the ion–ion correlation on this time scale, we need to solve (4.20), repeated here:

(6.1)\begin{align} & \left(\frac{\partial}{\partial t} + \mathcal{L}(12)\right)g_{ii}(12) + \int\, {\rm d}(3) \mathcal{V}^s(13)f_i(1)g_{ii}(32)\nonumber\\ & \qquad+ \int\, {\rm d}(3) \mathcal{V}^s(32) f_i(2)\,g_{ii}(13) ={-}\mathcal{V}^s(12)\,f_i(1)f_i(2), \end{align}

with fixed ion distribution $f_i = f_{i0}$ and screening length $k_e = k_{e1}$.

The reason the ion distribution can be fixed is that it changes on the much slower ion–ion collision time scale; we invoke the usual Bogoliubov time scale hierarchy.⁶ Of course, we are claiming that the ions heat up on the plasma-frequency time scale $t\sim 1/\omega _{pi}$, so the distribution function $f_i$ cannot actually be constant! However, we will show that this change is small in $1/\varLambda$ and is therefore a higher-order effect that can be neglected (remember that in (6.1), terms which are small in $1/\varLambda$ have already been dropped).

The reason the electron temperature, and therefore the screening length $k_e$, can be fixed is similar and has already been discussed: the variation is small in $1/\varLambda$ so can be neglected.

To proceed, we will use a Green's function approach to solve

(6.2)\begin{align} & \left(\frac{\partial}{\partial t} + \mathcal{L}(12)\right)g_{ii}(12) + \int\, {\rm d}(3) \mathcal{V}^s(13)f_{i0}(1)\,g_{ii}(32)\nonumber\\ & \qquad+ \int\, {\rm d}(3) \mathcal{V}^s(32)f_{i0}(2)g_{ii}(13) = S(12), \end{align}

where $\mathcal {V}^s$ is assumed to have fixed screening length $k_e=k_{e1}$ from now on, and the source term is $S(12) = -\mathcal {V}^s(12)f_{i0}(1)f_{i0}(2)$.

Since $\mathcal {L}(12) = \mathcal {L}(1) + \mathcal {L}(2)$, we can write (6.2) as

(6.3)\begin{equation} \left(\frac{\partial}{\partial t} + \mathcal{O}(1) + \mathcal{O}(2)\right)g_{ii}(12) = S(12), \end{equation}

where the operator $\mathcal {O}(1)$ is defined by

(6.4)\begin{equation} \mathcal{O}(1)a(1) = \mathcal{L}(1)a(1) + \int\, {\rm d}(3) \mathcal{V}^s(13)f_{i0}(1)a(3) \end{equation}

for any function $a(1)$.

Let $U(12,t)$ be the solution of

(6.5)\begin{equation} \left( \frac{\partial}{\partial t} + \mathcal{O}(1) \right)U(12,t) = \delta(t)\delta(12), \end{equation}

where $\delta (12) = \delta (\boldsymbol {v}_1-\boldsymbol {v}_2)\delta (\boldsymbol {r}_1-\boldsymbol {r}_2)$ and $U(12,t)=0$ for $t<0$. This means $U(12,t)$ is the Green's function for the linearised Vlasov equation with Yukawa interactions.

Equation (6.5) is first-order in time, so $U(12,t)$ can also be characterised as the function that satisfies

(6.6)\begin{equation} \left( \frac{\partial}{\partial t} + \mathcal{O}(1) \right)U(12,t) = 0 \end{equation}

for $t>0$, with initial condition $\lim _{t\to 0^+}U(12,t) = \delta (12)$.

Then, the Green's function for (6.2) or (6.3) is $U(13,t)U(24,t)$. This means

(6.7)\begin{equation} \left( \frac{\partial}{\partial t} + \mathcal{O}(1) + \mathcal{O}(2)\right)U(13,t)U(24,t) = \delta(t)\delta(13)\delta(24). \end{equation}

To prove this, we simply observe that

(6.8)\begin{equation} \left( \frac{\partial}{\partial t} + \mathcal{O}(1) + \mathcal{O}(2)\right)U(13,t)U(24,t) = 0, \end{equation}

for $t>0$ and $\lim _{t\to 0^+}U(13,t)U(24,t) = \delta (13)\delta (24)$.

Therefore, the general solution to (6.2) is

(6.9)\begin{align} g_{ii}(12,t) & = \int\, {\rm d}(34)g_{ii}(34,t=0)U(13,t)U(24,t)\nonumber\\ & \quad + \int_0^t {\rm d} t' \int\, {\rm d}(34)\, S(34)U(13,t-t')U(24,t-t'). \end{align}

The source term is time independent, so we can also write

(6.10)\begin{align} g_{ii}(12,t) & = \int\, {\rm d}(34)g_{ii}(34,t=0)U(13,t)U(24,t)\nonumber\\ & \quad+ \int_0^t {\rm d} t' \int\, {\rm d}(34) S(34)U(13,t')U(24,t'). \end{align}

We still need to determine the Green's function $U(12,t)$. It can found using Fourier and Laplace transforms, following the usual solution for the Vlasov equation with one species (see pp. 44–46 of Ichimaru 1992) but with the Coulomb potential $\tilde {\varphi }=4{\rm \pi} Z^2e^2/k^2$ replaced with $\tilde {\varphi }^{(s)} = 4{\rm \pi} Z^2e^2/(k^2+k_{e1}^2)$. Its Fourier transform is

(6.11)\begin{align} \tilde{U}(\boldsymbol{k}, \boldsymbol{v}_1, \boldsymbol{v}_2, t) & = \int_{\mathcal{B}}\frac{{\rm d}\omega}{2{\rm \pi}}{\rm e}^{-{\rm i}\omega t}\left(\frac{{\rm i}\,\delta(\boldsymbol{v}_1-\boldsymbol{v}_2)}{\omega - \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}_1} + \frac{1}{n_i}\frac{k_{i0}^2}{k^2+k_{e1}^2}\frac{1}{\epsilon(\omega,\boldsymbol{k})}\right.\nonumber\\ & \quad \times \left.\frac{{\rm i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}_1}{\omega - \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}_1}\frac{1}{\omega - \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}_2}f_{i0}(\boldsymbol{v}_1)\right).\end{align}

Here, the contour $\mathcal {B}$ is the usual Bromwich contour for inverse Laplace transforms, meaning it is a straight line in the complex plane from $i\sigma -\infty$ to $i\sigma +\infty$, where $\sigma$ is chosen so that all the poles of the integrand lie below the contour. The dielectric function $\epsilon (\omega,\boldsymbol {k})$ is given by

(6.12)\begin{equation} \epsilon(\omega ,\boldsymbol{k}) = 1 - \frac{1}{n_{{\rm i}}}\frac{k_{i0}^2}{k^2+k_{e1}^2} \int_{\mathcal{L}} \,\mathrm{d}\boldsymbol{v} \frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}}{\omega -\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}} f_{i0}(\boldsymbol{v}), \end{equation}

where $\mathcal {L}$ is the usual Landau contour (Landau 1965). In Fourier space, (6.10) reads

(6.13)\begin{align} \tilde{g}_{ii}(\boldsymbol{k},\boldsymbol{v}_1,\boldsymbol{v}_2,t) & = \int\, {\rm d}\boldsymbol{v}_3{\rm d}\boldsymbol{v}_4\,\tilde{g}_{ii}(\boldsymbol{k},\boldsymbol{v}_3,\boldsymbol{v}_4,t=0)\tilde{U}(\boldsymbol{k},\boldsymbol{v}_1,\boldsymbol{v}_3,t)\tilde{U}(-\boldsymbol{k},\boldsymbol{v}_2,\boldsymbol{v}_4,t)\nonumber\\ & \quad + \int_0^t \,{\rm d} t' \int\, {\rm d}\boldsymbol{v}_3{\rm d}\boldsymbol{v}_4 \tilde{S}(\boldsymbol{k},\boldsymbol{v}_3,\boldsymbol{v}_4)\tilde{U}(\boldsymbol{k},\boldsymbol{v}_1,\boldsymbol{v}_3,t')\tilde{U}(-\boldsymbol{k},\boldsymbol{v}_2,\boldsymbol{v}_4,t'). \end{align}

To find the heating that results from this change in correlations, we should use this solution for $g_{ii}$ to evolve $f_i$ according to

(6.14)\begin{equation} \frac{\partial f_i(\boldsymbol{v}_1,t)}{\partial t} ={-}\frac{T_{i0}k_{i0}^2}{n_im_i}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3} \frac{1}{k^2+k_{e1}^2}{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial}{\partial\boldsymbol{v}_1}\int\, {\rm d}\boldsymbol{v}_2\tilde{g}_{ii}(\boldsymbol{k},\boldsymbol{v}_1,\boldsymbol{v}_2,t), \end{equation}

which is just (4.12) with fixed screening length $k_e = k_{e1}$. In fact, we are really interested in the change in the ion kinetic energy density $K_i = \int \, {\rm d}\boldsymbol {v}_1 \tfrac {1}{2}m_i\boldsymbol {v}_1^2 f_i(\boldsymbol {v}_1,t)$, given by

(6.15)\begin{equation} \frac{{\rm d} K_i}{{\rm d} t} = \frac{T_{i0}k_{i0}^2}{n_i}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{1}{k^2+k_{e1}^2}\int\, {\rm d}\boldsymbol{v}_1 {\rm i} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}_1\int\, {\rm d}\boldsymbol{v}_2 \tilde{g}_{ii}(\boldsymbol{k},\boldsymbol{v}_1,\boldsymbol{v}_2,t). \end{equation}

Using (6.13) and (6.11), this becomes

(6.16)\begin{equation} \frac{{\rm d} K_i}{{\rm d} t} = \frac{{\rm d} K_i^{(1)}}{{\rm d} t} + \frac{{\rm d} K_i^{(2)}}{{\rm d} t}, \end{equation}

where

(6.17)\begin{align} \frac{{\rm d} K_i^{(1)}}{{\rm d} t} & = \frac{T_{i0}k_{i0}^4}{n_i^2}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{1}{k^2+k_{D0}^2}\frac{1}{k^2+k_{e1}^2}\int\, {\rm d}\boldsymbol{v}_3{\rm d}\boldsymbol{v}_4 f_i(\boldsymbol{v}_3)f_i(\boldsymbol{v}_4)\nonumber\\ & \quad \times\int\, {\rm d}\boldsymbol{v}_1{\rm d}\boldsymbol{v}_2({\rm i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{v}_1) \int_\mathcal{B}\frac{{\rm d}\omega}{2{\rm \pi}}\int_{\mathcal{B}'}\frac{{\rm d}\omega'}{2{\rm \pi}}\, {\rm e}^{-{\rm i}(\omega+\omega')t}\nonumber\\ & \quad \times\left( \frac{\delta(\boldsymbol{v}_1-\boldsymbol{v}_3)}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1} + \frac{1}{n_i}\frac{k_{i0}^2}{k^2+k_{e1}^2}\frac{1}{\epsilon(\omega,\boldsymbol{k})}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1}\frac{1}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_3} f_{i0}(\boldsymbol{v}_1)\right)\nonumber\\ & \quad \times\left( \frac{\delta(\boldsymbol{v}_2-\boldsymbol{v}_4)}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2} - \frac{1}{n_i}\frac{k_{i0}^2}{k^2+k_{e1}^2}\frac{1}{\epsilon(\omega',\boldsymbol{k})}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2}\frac{1}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_4} f_{i0}(\boldsymbol{v}_2)\right), \end{align}

(6.18)\begin{align} \frac{{\rm d} K_i^{(2)}}{{\rm d} t} & = \frac{T_{i0}k_{i0}^4}{n_i^2}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{1}{(k^2+k_{e1}^2)^2}\int\, {\rm d}\boldsymbol{v}_3{\rm d}\boldsymbol{v}_4 \,f_i(\boldsymbol{v}_3)f_i(\boldsymbol{v}_4)({\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}_3-\boldsymbol{v}_4))\nonumber\\ & \quad \times\int\, {\rm d}\boldsymbol{v}_1{\rm d}\boldsymbol{v}_2\,({\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1) \int_0^t\, {\rm d} t'\int_\mathcal{B}\frac{{\rm d}\omega}{2{\rm \pi}}\int_{\mathcal{B}'}\frac{{\rm d}\omega'}{2{\rm \pi}}\, {\rm e}^{-{\rm i}(\omega+\omega')t'}\nonumber\\ & \quad \times\left( \frac{\delta(\boldsymbol{v}_1-\boldsymbol{v}_3)}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1} + \frac{1}{n_i}\frac{k_{i0}^2}{k^2+k_{e1}^2}\frac{1}{\epsilon(\omega,\boldsymbol{k})}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1}\frac{1}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_3} f_{i0}(\boldsymbol{v}_1)\right)\nonumber\\ & \quad \times\left( \frac{\delta(\boldsymbol{v}_2-\boldsymbol{v}_4)}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2} - \frac{1}{n_i}\frac{k_{i0}^2}{k^2+k_{e1}^2}\frac{1}{\epsilon(\omega',\boldsymbol{k})}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2}\frac{1}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_4} f_{i0}(\boldsymbol{v}_2)\right). \end{align}

The first term ${\rm d} K_i^{(1)}/{\rm d} t$ involves the initial correlation, and therefore the parameter $k_{D0}$, whereas the second does not. While it is common in plasma kinetic theory to ignore the initial condition when evolving pair correlations, because any memory of the initial state is quickly lost, in situations involving ultrafast relaxation, this Markovian approximation is not valid (Semkat, Kremp & Bonitz 1999; Morawetz et al. 2001). Correlation heating arises because the initial ion–ion correlation is not correct for the new electron temperature and it occurs in the brief time before memory of the initial state is lost. If there is no electron heating, so that $k_{D0} = k_{D1}$, then ${\rm d} K_i^{(1)}/{\rm d} t$ and ${\rm d} K_i^{(2)}/{\rm d} t$ must cancel.

In Appendix C, the integrals in (6.17) and (6.18) are simplified separately and the result is

(6.19)\begin{equation} \frac{{\rm d} K_i}{{\rm d} t} ={-}\frac{{\rm d}}{{\rm d} t}\left[\frac{1}{2}T_{i0}(k_{D0}^2-k_{D1}^2)\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{(k^2 + k_{D1}^2)|\mathcal{F}(t,\boldsymbol{k})|^2}{(k^2+k_{e1}^2)(k^2+k_{D0}^2)}\right]. \end{equation}

Indeed, when $k_{D0}=k_{D1}$, this vanishes. Here, the function $\mathcal {F}(t,\boldsymbol {k})$ is defined by

(6.20)\begin{equation} \mathcal{F}(t,\boldsymbol{k}) = \int_\mathcal{C}\frac{{\rm d}\omega}{2{\rm \pi}}\, {\rm e}^{-{\rm i}\omega t} \frac{1}{\omega} \frac{\epsilon(\omega,\boldsymbol{k})-1}{\epsilon(\omega,\boldsymbol{k})}, \end{equation}

where the contour $\mathcal {C}$ passes below the pole at the origin $\omega = 0$ but above all the poles of $1/\epsilon$, as shown in figure 6.

Figure 6. Contour $\mathcal {C}$ used to define $\mathcal {F}(t,\boldsymbol {k})$.

We can check that this is consistent with our ion temperature change in § 5.2. Integrating (6.19) gives the change in ion kinetic energy density,

(6.21)\begin{equation} \Delta K_i(t) ={-}\frac{1}{2}T_{i0}(k_{D0}^2-k_{D1}^2)\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{(k^2 + k_{D1}^2)\left(|\mathcal{F}(t,\boldsymbol{k})|^2 - |\mathcal{F}(0,\boldsymbol{k})|^2\right)}{(k^2+k_{e1}^2)(k^2+k_{D0}^2)}. \end{equation}

Since the contour $\mathcal {C}$ can run entirely below the real axis, the integrand in (6.20) is exponentially suppressed at large times. Therefore, $\mathcal {F}(\infty, \boldsymbol {k})=0$ and the total change over this time scale is

(6.22)\begin{equation} \Delta K_i(\infty) = \frac{1}{2}T_{i0}(k_{D0}^2-k_{D1}^2)\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{(k^2 + k_{D1}^2)|\mathcal{F}(0,\boldsymbol{k})|^2}{(k^2+k_{e1}^2)(k^2+k_{D0}^2)}. \end{equation}

To evaluate this integral, we need

(6.23)\begin{equation} \mathcal{F}(0,\boldsymbol{k}) = \int_\mathcal{C}\frac{{\rm d}\omega}{2{\rm \pi}} \frac{1}{\omega} \frac{\epsilon(\omega,\boldsymbol{k})-1}{\epsilon(\omega,\boldsymbol{k})}. \end{equation}

Now that there is no longer an exponential factor, the integrand decays as $|\omega |\to \infty$. So, completing the contour in the upper half-plane as shown in figure 7 encloses a single pole at the origin $\omega = 0$:

(6.24)\begin{equation} \mathcal{F}(0,\boldsymbol{k}) = {\rm i} \left( \frac{\epsilon(0,\boldsymbol{k})-1}{\epsilon(0,\boldsymbol{k})}\right) = {\rm i}\left(\frac{k_{i0}^2}{k^2 + k_{D1}^2}\right). \end{equation}

Using this in (6.22), we find

(6.25)\begin{equation} \Delta K_i(\infty) = \frac{1}{8{\rm \pi}}T_{i0}k_{i0}^4\left( \frac{1}{k_{D1}+k_{e1}} - \frac{1}{k_{D0} + k_{e1}} \right).\end{equation}

This change in ion kinetic energy density is exactly equivalent to the ion temperature change calculated in (5.6). Therefore, the ions only change their kinetic energy on the extremely fast ion-plasma-frequency time scale.

Figure 7. Completing the contour in the upper half-plane.

Equation (6.21) is equivalent to

(6.26)\begin{equation} \Delta K_i(t) = \Delta K_i(\infty) - \frac{1}{2}T_{i0}(k_{D0}^2-k_{D1}^2)\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{(k^2 + k_{D1}^2)|\mathcal{F}(t,\boldsymbol{k})|^2}{(k^2+k_{e1}^2)(k^2+k_{D0}^2)} . \end{equation}

Although $\mathcal {F}(t,\boldsymbol {k})$ is still an unknown function, the fact that the integrand is positive means that when the electrons are heated, so $k_{D0}>k_{D1}$, we must have ${\Delta K_i(t)<\Delta K_i(\infty )}$. Therefore, the ion temperature cannot overshoot its final value. A similar result for the ultrafast relaxation of a weakly coupled one-component plasma with zero initial correlation has been found using molecular dynamics simulations (Zwicknagel 1999).

Equation (6.26) is plotted in figure 8 for different combinations of electron temperature before and after heating. In figure 8(a), the two curves with $T_{e1} = 10T_{i0}$ are not monotonically increasing; at later times, $\Delta K_i$ contains a slight wobble.

Figure 8. Change in ion kinetic energy for different combinations of $T_{e0}/T_{i0}$ and $T_{e1}/T_{i0}$, normalised using $\varLambda _i = n_i/k_{i0}^3$. In panel (a), the electrons are heated, while in panel (b) they are cooled.

7 Modification to the ion distribution

Fusion reactions involve suprathermal ions that have velocities well above the thermal speed. So, to increase the fusion rate, it is desirable to increase the number of fast ions in the tail of the distribution. In this section, we investigate how efficient correlation heating is at creating suprathermal ions, by calculating the change in the ion velocity distribution resulting from correlation heating on the ion-plasma-frequency time scale. There is no reason to expect the ions to remain Maxwellian; molecular dynamics simulations of an initially uncorrelated YOCP by Murillo (2006) showed strong deviations from a Maxwellian in the early stages of the ultrafast relaxation of the pair correlation.

Using the same approach as in the previous section, we can integrate (6.14) from $t=0$ to $t=\infty$ to find the change in the ion distribution ${\delta f(\boldsymbol {v}) = f_i(\boldsymbol {v},t=\infty ) - f_i(\boldsymbol {v},t=0)},$ which is

(7.1)\begin{equation} \delta f(\boldsymbol{v}) = \delta f^{(1)}(\boldsymbol{v}) + \delta f^{(2)}(\boldsymbol{v}), \end{equation}

where

(7.2)\begin{align} \delta f^{(1)}(\boldsymbol{v}_1) & ={-}\frac{T_{i0}k_{i0}^4}{n_i^2m_i}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{1}{k^2+k_{D0}^2}\frac{1}{k^2+k_{e1}^2}{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial}{\partial \boldsymbol{v}_1}\int\, {\rm d}\boldsymbol{v}_3{\rm d}\boldsymbol{v}_4 f_i(\boldsymbol{v}_3)f_i(\boldsymbol{v}_4)\nonumber\\ & \quad \times\int\, {\rm d}\boldsymbol{v}_2\, \int_0^\infty \,{\rm d} t\int_\mathcal{B}\frac{{\rm d}\omega}{2{\rm \pi}}\int_{\mathcal{B}'}\frac{{\rm d}\omega'}{2{\rm \pi}}\, {\rm e}^{-{\rm i}(\omega+\omega')t} \nonumber\\ & \quad \times\left( \frac{\delta(\boldsymbol{v}_1-\boldsymbol{v}_3)}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1} + \frac{1}{n_i}\frac{k_{i0}^2}{k^2+k_{e1}^2}\frac{1}{\epsilon(\omega,\boldsymbol{k})}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1}\frac{1}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_3} f_{i0}(\boldsymbol{v}_1)\right)\nonumber\\ & \quad \times\left( \frac{\delta(\boldsymbol{v}_2-\boldsymbol{v}_4)}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2} - \frac{1}{n_i}\frac{k_{i0}^2}{k^2+k_{e1}^2}\frac{1}{\epsilon(\omega',\boldsymbol{k})}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2}\frac{1}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_4} f_{i0}(\boldsymbol{v}_2)\right), \end{align}

(7.3)\begin{align} \delta f^{(2)}(\boldsymbol{v}_1) & ={-}\frac{T_{i0}k_{i0}^4}{n_i^2m_i}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{1}{(k^2+k_{e1}^2)^2}{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\frac{\partial}{\partial \boldsymbol{v}_1}\int\, {\rm d}\boldsymbol{v}_3{\rm d}\boldsymbol{v}_4 f_i(\boldsymbol{v}_3)f_i(\boldsymbol{v}_4)\nonumber\\ & \quad \times\left({\rm i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{v}_3-\boldsymbol{v}_4)\right)\int\, {\rm d}\boldsymbol{v}_2\,\int_0^\infty {\rm d} t \int_0^t {\rm d} t'\int_\mathcal{B}\frac{{\rm d}\omega}{2{\rm \pi}}\int_{\mathcal{B}'}\frac{{\rm d}\omega'}{2{\rm \pi}}\, {\rm e}^{-{\rm i}(\omega+\omega')t'}\nonumber\\ & \quad \times\left( \frac{\delta(\boldsymbol{v}_1-\boldsymbol{v}_3)}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1} + \frac{1}{n_i}\frac{k_{i0}^2}{k^2+k_{e1}^2}\frac{1}{\epsilon(\omega,\boldsymbol{k})}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_1}\frac{1}{\omega - \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_3} f_{i0}(\boldsymbol{v}_1)\right)\nonumber\\ & \quad \times\left( \frac{\delta(\boldsymbol{v}_2-\boldsymbol{v}_4)}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2} - \frac{1}{n_i}\frac{k_{i0}^2}{k^2+k_{e1}^2}\frac{1}{\epsilon(\omega',\boldsymbol{k})}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_2}\frac{1}{\omega' + \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}_4} f_{i0}(\boldsymbol{v}_2)\right). \end{align}

These integrals are simplified in Appendix D and the result is

(7.4)\begin{equation} \delta f(\boldsymbol{v}) = \frac{T_{i0}k_{i0}^2}{2n_im_i}(k_{D0}^2-k_{D1}^2)\,\frac{\partial}{\partial \boldsymbol{v}}\boldsymbol{\cdot}\left[\frac{\boldsymbol{v}f_{i0}(\boldsymbol{v})}{\boldsymbol{v}^2}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{1-|\epsilon(0,\boldsymbol{k})/\epsilon(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}, \boldsymbol{k})|^2}{(k^2+k_{D0}^2)(k^2+k_{D1}^2)}\right]. \end{equation}

Since this takes the form of a velocity-space divergence, the perturbation expansion adopted here conserves the total number of ions. It is also reassuring that, by multiplying by $\tfrac {1}{2} m_i\boldsymbol {v}^2$, integrating over $\boldsymbol {v}$ and making use of the integral

(7.5)\begin{equation} \int\, {\rm d}\boldsymbol{v}\, \frac{f_{i0}(\boldsymbol{v})}{|\epsilon(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}, \boldsymbol{k})|^2} = n_i\left(\frac{k^2+k_{e1}^2}{k^2+k_{D1}^2}\right), \end{equation}

the same result is obtained for the total ion kinetic energy change that we found in (5.6) and (6.25). Equation (7.4) is plotted in figure 9 for different combinations of the initial and final electron temperatures $T_{e0}$ and $T_{e1}$. We see that heating the electrons has a larger effect on the ion distribution than cooling them. In addition, the effect is larger when the electrons are initially colder.

Figure 9. (a) Change in ion speed distribution $4{\rm \pi} v^2 \delta f(v)$, normalised using ${\varLambda _i = n_i/k_{i0}^3}$, for three different combinations of $T_{e0}/T_{i0}$ and $T_{e1}/T_{i0}$. In each case, the electrons are heated. (b) Cross-section, in the plane $v_y=v_z=0$, of $\delta f(\boldsymbol {v})/f_{i0}(\boldsymbol {v})$ for the same three cases. (c) Similar to panel (a), except now three cases are presented in which the electrons are cooled. (d) Similar to panel (b), except the electrons are cooled.

An interesting feature is that if the electrons are made hotter than the ions, then $\delta f(\boldsymbol {v})$ develops a small bump. This bump is located around the ion sound speed, which is of the order of $v=v_{\mathrm {th}i}\,k_{i0}/k_{e1}$, where $v_{\mathrm {th}i}=(2T_{i0}/m_i)^{1/2}$. The reason is that when the electrons are made very hot, correlation heating leads to the excitation of long-lived ion-acoustic modes, which damp partially on the ions. The ion distribution therefore flattens slightly at speeds in resonance with these modes, which appears as a bump in $\delta f(\boldsymbol {v})$.

We also see that at large velocities, $\delta f(\boldsymbol {v})/f_{i0}(\boldsymbol {v})$ has a constant limit. It is simple to justify this behaviour analytically. If $|\boldsymbol {v}|$ is large enough, we can replace the dielectric function $\epsilon (\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {v},\boldsymbol {k})\to 1$.⁷ Then,

(7.6)\begin{equation} \delta f(\boldsymbol{v}) = \frac{T_{i0}k_{i0}^2}{2n_im_i}(k_{D0}^2-k_{D1}^2)\,\frac{\partial}{\partial \boldsymbol{v}}\boldsymbol{\cdot}\left[\frac{\boldsymbol{v}f_{i0}(\boldsymbol{v})}{\boldsymbol{v}^2}\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{1-|\epsilon(0,\boldsymbol{k})|^2}{(k^2+k_{D0}^2)(k^2+k_{D1}^2)}\right]. \end{equation}

At large velocities, the dominant contribution comes from the derivative acting on $f_i(\boldsymbol {v}):$

(7.7)\begin{equation} \frac{\partial}{\partial \boldsymbol{v}}\boldsymbol{\cdot}\left[\frac{\boldsymbol{v}f_{i0}(\boldsymbol{v})}{\boldsymbol{v}^2}\right] \sim{-}\frac{m_i}{T_{i0}}f_{i0}(\boldsymbol{v}), \end{equation}

so

(7.8)\begin{equation} \delta f(\boldsymbol{v}) = \frac{1}{2}\frac{k_{i0}^2}{n_i}(k_{D0}^2-k_{D1}^2)f_{i0}(\boldsymbol{v})\int \frac{{\rm d}\boldsymbol{k}}{(2{\rm \pi})^3}\frac{|\epsilon(0,\boldsymbol{k})|^2-1}{(k^2+k_{D0}^2)(k^2+k_{D1}^2)}. \end{equation}

This $\boldsymbol {k}$-integral may be evaluated using $\epsilon (0,\boldsymbol {k}) = (k^2+k_{D1}^2)/(k^2+k_{e1}^2)$ but the result is not very enlightening. The important conclusion is that $\delta f(\boldsymbol {v})\propto f_{i0}(\boldsymbol {v})$, which means the number of ions at large velocities is increased by a constant factor.

This result has the following implication: using correlation heating to quickly supply a given amount of energy to the ions in a moderately coupled plasma leads to an ion distribution with a depleted tail relative to a Maxwellian at the same energy. Over a longer time scale, the tail will be filled out by collisions.

To see this, suppose the initial ion distribution is

(7.9)\begin{equation} f_{i0}(\boldsymbol{v}) = n_i\left( \frac{m_i}{2{\rm \pi} T_{i0}} \right)^{3/2}\exp\left( -\frac{m_i\boldsymbol{v}^2}{2 T_{i0}} \right). \end{equation}

If the distribution remains Maxwellian but the temperature is increased by a small amount $\delta T_i \sim T_i/\varLambda$, then the modification to the distribution at large velocities is

(7.10)\begin{equation} \delta f(\boldsymbol{v}) \sim \left( \frac{m_i\boldsymbol{v}^2}{2T_{i0}^2}\delta T_{{\rm i}}\right)f_{i0}(\boldsymbol{v}) \end{equation}

for $|\boldsymbol {v}|\gg (T_{i0}/m_i)^{1/2}$. For sufficiently large velocities, this modification is always greater than the change in (7.8) that can be achieved using correlation heating. This is not surprising: widening a Gaussian distribution increases the area under the tail more effectively than multiplication of the distribution by a constant factor. Therefore, correlation heating does not efficiently energise suprathermal ions on the ion-plasma-frequency time scale in a moderately coupled plasma. In the next section, we argue that this conclusion may not hold in strongly coupled systems.

8 Discussion

8.1 Key formulae

If some energy $Q$ is suddenly supplied to the electrons, how much energy and power is transferred to the ions as a result of the change in shielding? The calculations in § 5 allow us to answer this question quantitatively for moderately coupled plasmas. If the electrons become sufficiently hot, the amount of energy transferred is independent of $Q$. In this limit, the average energy supplied per ion is

(8.1)\begin{equation} \text{Energy in eV per ion} \approx 30 \frac{ Z^3 n_{26}^{1/2}}{T_{1}^{1/2}}\left( 1 - \frac{1}{\sqrt{1+(1/Z)(1/\tau)}} \right). \end{equation}

The time taken for the ion kinetic energy to increase is roughly $t\approx 1.5/\omega _{pi}$ (see figure 8). So, the power delivered to the ions is

(8.2)\begin{equation} \text{Power in eV per fs per ion} \approx 76 \frac{ Z^4 n_{26}}{A^{1/2} T_{1}^{1/2}}\left( 1 - \frac{1}{\sqrt{1+(1/Z)(1/\tau)}} \right). \end{equation}

Here, $Z$ is the ion charge state, $n_{26}$ is the ion number density normalised to $10^{26}\,{\rm cm}^{-3}$, $A$ is the ion mass normalised to the proton mass, $T_1$ is the initial ion temperature in ${\rm keV}$, and $\tau = T_{e0}/T_{i0}$ is the ratio between the initial electron and ion temperatures. These formulae give a practical way to estimate the maximum amount of correlation heating possible in any given system.

8.2 How large can correlation heating be?

The moderate coupling ordering $1\ll \varLambda /\lambda \ll (m_i/m_e)^{1/2}$ is required for the rigorous analytical description of correlation heating presented in this paper. However, since the ion heating is $\Delta T \sim T/\varLambda$, this ordering precludes substantial heating, which might be possible in more strongly coupled plasmas. It is therefore important to consider the maximum heating attainable when the moderate coupling restriction is removed.

The two key energy scales in a plasma are the temperature $T$ (assuming $T_e\sim T_i$) and the typical potential energy between nearby particles, $\varphi \sim e^2n^{1/3}$. The plasma parameter scales as $\varLambda \sim (T/\varphi )^{3/2}$, so in a moderately coupled plasma, $\Delta T \sim (\varphi /T)^{3/2}\, T$. As the coupling strength $\varphi /T$ increases, correlation heating gets larger. A natural question to ask is whether this scaling continues to hold in strongly coupled plasmas in which $\varphi \gg T$ initially. Is it possible to get significant heating?

Consider the simplest case of electrons heated to a large, essentially infinite, temperature. They become a uniform neutralising background and have no influence over the subsequent ion relaxation. As the ion–ion correlation adjusts, the ion energy is conserved,

(8.3)\begin{equation} U = \frac{3}{2}N_iT_i + \frac{1}{2}Vn_i^2\int\, {\rm d}\boldsymbol{r} \left(\frac{Z^2e^2}{r}\right)G_{ii}(\boldsymbol{r}) = \mathrm{const.}, \end{equation}

as in (5.3). This means (Gericke & Murillo 2003)

(8.4)\begin{equation} \Delta T_i ={-} \frac{1}{3}n_i\int\, {\rm d}\boldsymbol{r} \left(\frac{Z^2e^2}{r}\right)\Delta G_{ii}(\boldsymbol{r}), \end{equation}

where $\Delta G_{ii}(\boldsymbol {r})$ is the change in the ion–ion correlation. It is clear that to obtain heating, we want $\Delta G_{ii}<0$ (for most values of $\boldsymbol {r}$, at least). Since the ion–ion correlation $G_{ii}(\boldsymbol {r})$ is negative for most $\boldsymbol {r}$, this means its absolute value must increase. Therefore, the ions have to become more ordered, or more correlated, to heat up.

This means it is impossible for ions that were initially strongly coupled, $\varphi \gg T$, to heat up so much that they become weakly coupled, $\Delta T \gg \varphi$. If this were the case, the ions would have a smaller $|G_{ii}|$ in their weakly coupled final state than in their initial state, so they would become less ordered and (8.4) implies $\Delta T_i<0$, which is a contradiction. Correlation heating occurs as the ions become more ordered; however, an increase in temperature encourages spatial disorder and opposes further heating. This negative feedback places an upper limit on how hot the ions can become. Therefore, the scaling $\Delta T\sim (\varphi /T)^{3/2}T$ must be modified in strongly coupled plasmas. For comparison, in disorder-induced heating, the ion temperature generally rises until $\varphi /T\sim$ $1$–$4$ (Murillo 2001; Chen et al. 2004; Simien et al. 2004; Cummings et al. 2005).

Strongly coupled ions have an equilibrium pair correlation that approximately represents a sphere that perfectly excludes other ions, with radius of the order of the inter-particle spacing (Ichimaru 1982). In this picture, the electron temperature does not affect the ion correlation, which implies $|\Delta G_{ii}|\ll |G_{ii}|$ so $\Delta T \ll \varphi$. However, significant heating ${\Delta T \gg T}$ could still be achievable.

Equilibrium correlations of strongly coupled YOCPs and two-temperature plasmas have been found numerically. The change in ion correlation when the electron temperature is modified, $\Delta G_{ii}$, is visible in plots of the pair correlation (radial distribution function) by Daughton, Murillo & Thode (2000) and Shaffer et al. (2017) for different electron temperatures. These results suggest significant heating could occur. For example, if the correlation changes by roughly $1\,\%$ in a plasma with coupling strength $\varphi /T\sim 100$ (as in Daughton et al. (2000) for an electron temperature change by a factor of $9$), then $\Delta T \sim T$ may be possible. We conclude that correlation heating could be large if the electrons and ions are initially strongly coupled.

8.3 Application of correlation heating

In magnetic confinement fusion, it is common for plasma to exist in a two-temperature state because of heating or cooling mechanisms that target one species preferentially and drive the electron and ion temperatures apart.

For example, devices may be operated in a hot-ion mode in which the ion temperature is maintained significantly above the electron temperature (Clarke 1980; Jin, Reiman & Fisch 2021; Kolmes et al. 2021) by mechanisms such as ion cyclotron heating, neutral beam heating or alpha channelling (Fisch & Rax 1992; Fisch & Herrmann 1994, 1999). The supershot plasma regime in the Tokamak Fusion Test Reactor (TFTR) (Strachan et al. 1987) and the Hot Ion H-mode regime in the Joint European Torus (JET) (Keilhacker et al. 1999) were both characterised by high fusion rates and much hotter ions than electrons.

As another example, electrons are cooled by Bremsstrahlung, which is particularly significant in reactors operating at the high temperatures required for advanced fuels such as proton-${}^{11}$Boron (Rider 1995, 1997). In such reactors, Bremsstrahlung makes hotter ions than electrons a requirement (Cai et al. 2022).

In inertial confinement fusion, a two-temperature state can arise when the electron temperature increases so rapidly that collisional energy exchange is too slow to enforce equal species temperatures. ICF plasmas receive a sudden input of energy from powerful lasers and may have high coupling strengths, so correlation heating is most relevant in this context.

The fact that increasing the electron temperature causes the ion temperature to also rise can be captured by equation of state models that allow different electron and ion temperatures while also incorporating strong-coupling effects (e.g. Fetsch et al. 2023). The novel possibility unlocked by correlation heating is fast, collisionless ion heating, which can be used (once, unless the electrons cool back down) to step the ion temperature upwards faster than would be possible just using ion–electron Coulomb collisions.

This could be important in fast ignition approaches to ICF, in which the hotspot must usually be heated to ignition in approximately $20\,{\rm ps}$ (Tabak et al. 1994; Atzeni 1999). The heating needs to be rapid so that there is no time for pressure equilibration within the fuel, meaning it needs to be faster than the hotspot confinement time (Tabak et al. 2006). It is usually assumed that the ions heat up on the ion–electron collision time scale.

Correlation heating could be leveraged to allow faster heating of the ions, which would be most effective if laser heating of the electrons occurred on an ion-plasma-frequency time scale. For this to be a more efficient way to achieve ignition, it is important that energy is supplied to the suprathermal tail of the ion distribution. The calculations in § 7 for a moderately coupled plasma suggest that this does not occur. However, this result may not hold in all cases.

The simple picture given in § 1 of ions flying apart after they are suddenly stripped of their shielding clouds suggests that suprathermal ions could be created if a group of ions all work together to push another ion in one direction. This would occur when there are significant fluctuations in the number of ions in any given region, before the heating. Large fluctuations, corresponding to a large clumps of well-shielded ions in close proximity, would have large potential energy immediately after the electron screening gets removed. This potential energy would get evenly shared between the ions when they fly apart. In this picture, correlation heating is similar to localised Coulomb explosions wherever there are small-scale density fluctuations in the plasma.

This picture suggests that fast ion creation could be significantly enhanced if the plasma is initially strongly coupled, so that it is dense and so that the electrons are cold enough to provide efficacious shielding before they are heated. However, if the plasma is too strongly coupled, then the ions can begin to form a regular spatial structure, which makes spatial inhomogeneities that create fast ions less likely.

Furthermore, correlation heating could be very effective at accelerating light minority ions to suprathermal speeds. If the repulsion between an initially stationary pair of ions is suddenly switched on, which is a simple model for the removal of strong shielding, then the particles fly apart but the majority of the interaction energy is received by the light ion. Note that the arguments in § 8.2, leading to the conclusion that correlation heating cannot heat strongly coupled ions enough that they become weakly coupled, do not apply to a minority ion species.

It is also worth pointing out that in systems with sufficiently hot electrons, the ion–ion collision time can be shorter than the ion–electron collision time, even for a fast ion. In such cases, after correlation heating supplies additional energy to the ions, the ion distribution re-Maxwellianises and this may fill out the tail of the distribution faster than ion–electron collisions would have.

For all of these reasons, it may be that there is an advantage to heating the plasma even faster than the typical laser duration of ${\sim }20\,{\rm ps}$ envisioned for fast ignition. This is within reach of modern ultrafast lasers, which can operate at the femtosecond level. A short pulse to leverage correlation heating could also be combined with a longer pulse that heats the plasma on the usual collisional time scale. The faster heating that results may allow a smaller hotspot, which would require less energy input to reach ignition – such intriguing possibilities merit further investigation.

9 Conclusion

In this paper, we used a novel expansion of the BBGKY hierarchy to rigorously derive the ion temperature increase that can be achieved using correlation heating in moderately coupled plasmas. The resulting formulae are compact and comprehensible, establishing a framework that makes it possible to quickly estimate how much potential energy could be extracted by suddenly stripping ions of their Debye clouds.

We began in § 3 with a detailed review of the pair correlations in a two-temperature plasma, including a new, simple derivation that only uses elementary plasma physics concepts.

In a departure from conventional plasma kinetic theory, we focused on moderately coupled plasmas by taking the mass ratio $m_e/m_i$ to be the smallest parameter in the problem, rather than $1/\varLambda$. Interchanging these two limits means the electron collision frequencies become faster than the ion plasma frequency, so the usual Bogoliubov time scale hierarchy is modified. This is crucial for understanding how the system relaxes on short time scales after rapid heating. We prove that the ions behave like a Yukawa one-component plasma, with electrons only entering via an effective screening length which can be calculated using conservation of energy.

The formalism is applied in §§ 5–7 to find the temperature changes that occur for both ions and electrons, the time-dependence of the ion kinetic energy, and the modification to the ion distribution function. The ion kinetic energy does not overshoot its final value, unlike in strongly coupled plasmas. Correlation heating is shown to be inefficient at heating fast ions in a moderately coupled plasma; however, we argue that this conclusion may not hold in strongly coupled systems, particularly if there if a light minority ion species present. This may be relevant for the fast-ignition approach to ICF, where a femtosecond laser pulse would likely be necessary to take advantage of correlation heating.

This work also uncovered potential advantages of integrating concepts from the fields of strongly coupled and ultracold plasmas into fusion research (Bergeson et al. 2019). We studied a process in which energy is stored in correlations and quickly converted into particle kinetic energy. An alternative possibility is to release the correlation energy by compressing the plasma (Avinash 2010; Avinash & Kaw 2014; Fetsch et al. 2023). Similar mechanisms have been proposed for different types of correlation; for example, compressing a plasma with correlations arising from Langmuir waves can energise hot electrons (Schmit, Dodin & Fisch 2010; Schmit & Fisch 2013; Schmit et al. 2013), while compressing a plasma with correlations arising from turbulence can cause sudden dissipation of the turbulent energy as heat (Davidovits & Fisch 2016a,b, 2017, 2019). The heating discussed in this paper depends on intrinsic correlations that are present in any plasma in which electrons shield ions, and is therefore of universal interest. Other applications could be found as strongly coupled plasmas are further studied theoretically, numerically and in future experiments on ultracold plasmas.