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 Reference Spitzer1956; Montgomery & Tidman Reference Montgomery and Tidman1964). 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 Reference Landau1936; Spitzer Reference Spitzer1956). 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.Footnote ¹ 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 (Reference Hazeltine and Waelbroeck2018). 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. Reference Bobylev, Potapenko and Sakanaka1997)

(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 Reference Richardson2019).

The first set of parameters is for a hydrogen plasma (Falcon et al. Reference Falcon, Rochau, Bailey, Gomez, Montgomery, Winget and Nagayama2015) created in the Z Pulsed Power Facility at Sandia National Laboratories (Sinars et al. Reference Sinars, Sweeney, Alexander, Ampleford, Ao, Apruzese, Aragon, Armstrong, Austin and Awe2020) as part of the Z Astrophysical Plasma Properties (ZAPP) collaboration (Rochau et al. Reference Rochau, Bailey, Falcon, Loisel, Nagayama, Mancini, Hall, Winget, Montgomery and Liedahl2014). 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. Reference Tabak, Hinkel, Atzeni, Campbell and Tanaka2006). 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. Reference Tabak, Hammer, Glinsky, Kruer, Wilks, Woodworth, Campbell, Perry and Mason1994). 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. Reference Bychenkov, Rozmus, Maksimchuk, Umstadter and Capjack2001; Logan et al. Reference Logan, Bangerter, Callahan, Tabak, Roth, Perkins and Caporaso2006). 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 Reference Badziak, Jabłoński and Wołowski2007). Fast ignition requirements in two-temperature plasmas have been studied, for example, by Eliezer et al. (Reference Eliezer, Henis, Nissim, Pinhasi and Val2015).

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.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 (Reference Salpeter1963), 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.Footnote ²

Instead of using statistical mechanics and calculating a partition function, we will follow Debye & Hückel (Reference Debye and Hückel1923) 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 Reference Krall and Trivelpeice1973), 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 Reference Balsecu1988, 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.Footnote ³ 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.Footnote ⁴

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 (Reference Seuferling1987). 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 (Reference Boercker and More1986) 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 (Reference Rosenfeld1994).

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.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 Reference Acciarri, Moore and Baalrud2022). 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. (Reference Fetsch, Foster and Fisch2023).

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):Footnote ⁵ 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 Reference Hamaguchi, Farouki and Dubin1997; Murillo Reference Murillo2001; Gericke et al. Reference Gericke, Murillo, Semkat, Bonitz and Kremp2003; Murillo Reference Murillo2009; Lyon, Bergeson & Murillo Reference Lyon, Bergeson and Murillo2013; Langin et al. Reference Langin, Strickler, Maksimovic, McQuillen, Pohl, Vrinceanu and Killian2016; Sprenkle et al. Reference Sprenkle, Silvestri, Murillo and Bergeson2022), 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.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 Reference Landau and Lifshitz1980)

(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 (Reference More1980), although it seems the correct formula (5.2) was first published by Triola (Reference Triola2022). An alternative derivation was provided by Fetsch et al. (Reference Fetsch, Foster and Fisch2023). 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. Reference Hamaguchi, Farouki and Dubin1997; Murillo Reference Murillo2001; Gericke et al. Reference Gericke, Murillo, Semkat, Bonitz and Kremp2003; Murillo Reference Murillo2009; Lyon et al. Reference Lyon, Bergeson and Murillo2013; Langin et al. Reference Langin, Strickler, Maksimovic, McQuillen, Pohl, Vrinceanu and Killian2016; Sprenkle et al. Reference Sprenkle, Silvestri, Murillo and Bergeson2022). 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.