Self-consistent strong screening applied to thermonuclear reactions

Abstract

Self-consistent strong plasma screening around light nuclei is implemented in the Big Bang nucleosynthesis (BBN) epoch to determine the short-range screening potential, $e\phi(r)/T\geq 1$, relevant for thermonuclear reactions. We numerically solve the non-linear Poisson-Boltzmann equation incorporating Fermi-Dirac statistics, adopting a generalized screening mass to find the electric potential in the cosmic BBN electron-positron plasma for finite-sized ⁴He nuclei as an example. Although the plasma follows Boltzmann statistics at large distances, Fermi-Dirac statistics is necessary when work performed by ions on electrons is comparable to their rest mass energy. While strong screening effects are generally minor due to the high BBN temperatures, they can enhance the fusion rates of high-$Z>2$ elements while leaving fusion rates of lower-$Z\leq 2$ elements relatively unaffected. Our results also reveal a pronounced spatial dependence of the strong screening potential near the nuclear surface. These findings about the electron-positron plasma’s role refine BBN theory predictions and offer broader applications for studying weakly coupled plasmas in diverse cosmic and laboratory settings.

1 Introduction

1.1 The Big Bang nucleosynthesis epoch

We address the formation of light elements during the Big Bang nucleosynthesis (BBN) epoch in the temperature range $86\,\mathrm{keV}>T>20\,\mathrm{keV}$ (Pitrou et al., 2018). In this temperature regime, electron-positron pairs ($e^{-}e^{+}$ plasma) are abundant, as highlighted in various recent studies (Wang et al., 2011; Hwang et al., 2021; Rafelski et al., 2023). We determine the magnitude and effect of the self-consistent strong field screening for finite-sized stationary nuclei in the early Universe induced by this exotic $e^{-}e^{+}$ plasma. These effects are of interest for precision BBN calculations at the 0.1% level and are of general interest in exploring the properties of plasma theory.

Here, we explore various static and nonlinear models of charge screening, a collective effect within the plasma that alters the potential between nuclei. Plasma screening involves electrons surrounding an ion’s charge $Ze$ (elementary charge $e>0$), which effectively ’screens’ or diminishes the influence of other nuclear charges beyond their radius, lowering the Coulomb barrier. In the context of nuclear reactions, this reduction in the Coulomb barrier facilitates increased penetration probability. This, in turn, boosts the rates of thermonuclear reactions and consequently alters the abundance of lighter elements formed in the early universe.

Traditionally, most BBN plasma studies assume the “weak-field” limit where the electromagnetic potential energy $\phi(r)$ is small compared to the thermal energy $T$

$$\frac{q\phi(r)}{T}\ll 1\,.$$ (1)

In this case, the movement of plasma particles is dominated by thermal fluctuations rather than the Coulomb force. Plasmas satisfying this condition are weakly coupled, indicating plasma effects will lead to linear corrections to the potential, such as in Debye-Hückel theory (Debye & Hückel, 1923).

Due to the relatively low baryon number density $n_{B}$ in the early Universe (Grayson et al., 2023), the internuclear distance $a=n_{B}^{-1/3}$ is large and the macroscopic properties of the early Universe plasma satisfy the weak field limit for the inter-nuclear spacing $a$ because, at this distance, the potential is much weaker than the thermal energy, which is below $86\,$keV

$$\frac{Ze\phi(a)}{T}\ll 1\,,\quad\text{with}\quad a=n_{B}^{-1/3}\,,$$ (2)

for light nuclei with charge $Z$. The weak field limit can accurately describe the electromagnetic fields in the plasma at large distances on the order of the inter-nuclear spacing $a$ but not at short distances where $Ze\phi/T$ could be larger than one.

Although the BBN plasma is weakly coupled, nonlinear corrections to the short-distance electromagnetic potential may be relevant to quantum tunneling in thermonuclear reactions. This is because the Gamow energy $E_{G}$ at which nuclei are most likely to tunnel is higher than the BBN thermal energy (Shaviv & Shaviv, 1996). The internuclear distance corresponding to Gamow energy is on the order of femtometers, such that for the short-range potential

$$\frac{Ze\phi(r_{E_{G}})}{T}>1\,,$$ (3)

where $r_{E_{G}}$ is the classical turning point at the Gamow energy.

$$E_{G}=\left(\frac{(\pi TZ_{1}Z_{2}\alpha)^{2}\mu_{r}}{2}\right)^{1/3}\,.$$ (4)

The reduced mass of the colliding light nuclei is $\mu_{r}$, and $Z_{1}$ and $Z_{2}$ are the respective charges of the nuclei. The condition Eq. (3) indicates that although the BBN plasma can be treated globally as weakly coupled, locally, the short-range potential can have nonlinear corrections simply because the electrostatic energy close to a nucleus can be much higher than the thermal energy. Naively using Boltzmann statistics, one expects exponential enhancement of the charge density (Grayson et al., 2023). We anticipate that these nonlinear corrections are relevant because they are on the order of the weak-field corrections due to dynamic motion and damping.

1.2 Methods of evaluation

Plasma screening effects were first considered in 1954 by (Salpeter, 1954), who proposed evaluating the enhancement of nuclear reactions by employing the static Debye-Hückel potential (Debye & Hückel, 1923; Salpeter & van Horn, 1969; Famiano et al., 2016). These applications focus on collision-less plasma only. Later, this approach was generalized for nuclei moving in the plasma (Hwang et al., 2021; Carraro et al., 1988; Gruzinov, 1998; Opher & Opher, 2000; Yao et al., 2017) i.e., ‘dynamic’ screening. In our previous work, we addressed scattering damping (Formanek et al., 2021) in the quantum electrodynamic (QED) $e^{-}e^{+}\gamma$ plasma (Grayson et al., 2023) where the BBN reaction network occurs. Similar numerical simulations were performed to study damping in the $e^{-}e^{+}\gamma$ plasma (Sasankan et al., 2020; Kedia et al., 2021).

To estimate the damped-dynamic enhancement of thermonuclear reactions during the BBN epoch, we can use the approximate analytic damped-dynamic screening potential denoted as $\phi_{\text{DD}}(r)$ in (Grayson et al., 2023). For weak screening with potential $e\phi/T\ll 1$, the reaction rate enhancement is related to the usual Salpeter enhancement factor (Salpeter, 1954), which only depends on the value of $\phi_{\text{DD}}(r)$ at the origin $r=0$

$$e\phi_{\text{DD}}(0)=Z\alpha\left(m_{D}c^{2}-\frac{3}{2}\beta\sigma_{0}\right)\,,$$ (5)

where $\alpha\approx 1/137$ is the fine structure constant, $\beta$ is the thermal velocity of ions in the plasma

$$\beta=\sqrt{\frac{2T}{M}}\,$$ (6)

with $M$ being the mass of nuclei and $T$ is the temperature of the BBN electron-positron plasma, with $k_{B}=1$. Then $m_{D}$ is the Debye screening mass, and $\sigma_{0}$ is the static conductivity of the early universe during BBN. The Debye screening mass is related to the Debye length $\lambda_{D}$ as

$$m_{D}c^{2}=\frac{\hbar c}{\lambda_{D}}\,,$$ (7)

where $\hbar c=197.3\ \mathrm{MeV}\ \mathrm{fm}$. From now on, we will absorb $c$ into $ct\rightarrow t$, such that $m_{D}c^{2}\rightarrow m_{D}$ has units of energy. Eq. (5) is only valid in the weak damping limit $\omega_{p}<\kappa$, with $\kappa$ being the average rate in of collisions and $\omega_{p}$ the plasma frequency. The thermal velocity of ions during BBN is very small due to their large mass $\beta<10^{-2}$. The conductivity of the early universe $\sigma_{0}$ is estimated in (Grayson et al., 2023) to be $0.23\,$keV. Therefore, the dynamic correction in the second term of Eq. (5) is small compared to the static result in the first term, $m_{D}\approx 2.7\,$keV. For this reason, the dynamic contribution to screening will be neglected in this work.

Approaches to strong screening have been previously developed in astrophysical plasmas such as solar plasmas and in ¹²C plasmas within white dwarfs. This research began with Salpeter in 1954, modeling a strongly coupled plasma with a uniform density of electrons that can completely screen nuclei. The following papers continue this research (Salpeter, 1954; Dewitt et al., 1973; Itoh et al., 1977, 1979; Ichimaru, 1982; Kravchuk & Yakovlev, 2014). Strong screening models are applicable for plasmas where the Coulomb energy is much larger than the thermal energy, even at distances larger than the ion separation. In this work, we study ‘intermediate screening,’ which we call self-consistent strong screening since it captures both the strong and weak field regimes and involves a self-consistent determination of the potential.

The strong screening potential is described by the Poisson-Boltzmann equation (Dzitko et al., 1995; Gruzinov & Bahcall, 1998; Brüggen & Gough, 1997, 2000; Bi et al., 2000; Liolios, 2004; Luo et al., 2020). The Poisson-Boltzmann equation is also used in physical chemistry to calculate the electromagnetic potential in electrolytic solutions (Fogolari et al., 2002). Traditionally, using the Boltzmann approximation assumes point-like charges, which cause the screening potential to diverge near their locations (Gruzinov & Bahcall, 1998). In the past, this issue was circumvented by calculating strong screening using a cluster expansion of weak field screening (Graboske et al., 1973), by considering ion correlations (Dewitt et al., 1973; Itoh et al., 1977), and employing the density matrix equation in quantum statistical mechanics (Gruzinov & Bahcall, 1998; Elsing et al., 2022). We solve the Poisson-Boltzmann equation directly by implementing a finite-sized source charge density and using Fermi-Dirac statistics. Other solution methods to the Poisson-Boltzmann equation using Fermi-Dirac statistics are explored in the context of stellar fusion by (Brüggen & Gough, 2000), where the approximate screening potential is found analytically using matching conditions. Numerical studies of the Poisson-Boltzmann equation were performed in (Cowan & Kirkwood, 1958) to compare with quantum mechanical calculations.

1.3 Novel approaches and results

Our study introduces several advancements and findings in plasma screening during the Big Bang nucleosynthesis (BBN) epoch. Our main result is the short-range potential around light nuclei. We introduce finite-sized nuclei into the Poisson-Boltzmann equation, eliminating singular behavior mentioned in (Gruzinov & Bahcall, 1998) that typically arises with point-like charges. Additionally, we formulate a generalized screening mass that characterizes the strength of polarization within the plasma, providing a framework for comparing various screening models. Our results indicate that Fermi-Dirac statistics is necessary when the work performed by the external potential of ions on electrons is comparable to their rest mass energy, even though the plasma adheres to Boltzmann statistics at larger distances. We establish an analytic strong screening mass in the ultrarelativistic limit for Fermi-Dirac statistics, which could be instrumental in developing a dynamic theory of strong screening in future studies.

While strong screening is generally small due to the high temperatures prevalent during BBN, it is interesting since it enhances the fusion rates of high-$Z$ elements while leaving elements with lower-$Z$ relatively unaffected. Our findings reveal that although the overall screening effect diminishes with decreasing temperature, the relative enhancement of strong screening over weak screening grows as the temperature decreases. We also find that the strong screening potential exhibits a pronounced spatial dependence near the nuclear surface. The step nature of this potential shown in Fig. 2 implies that relying solely on the screening energy near the origin can lead to an overestimate of the reaction rate enhancement due to strong screening. This overestimate was previously noted in (Itoh et al., 1977).

In Section 2, we review the theoretical background of self-consistent strong screening. Section 3 presents numerical Solutions to the Poisson-Boltzmann equation for ⁴He ions in the BBN plasma. In Section 4, we estimate the strong screening enhancement factor for ⁴He-⁴He and ¹²C-¹²C Collisions by evaluating the WKB tunneling probability of the strong screening potential. Section 5 reviews our results and discusses their implications for BBN and fusion.

2 Self-Consistent Strong Screening of Light Nuclei

We find the screened potential $\phi(r)$ in plasma by solving the Poisson equation for the induced polarization charge density $\rho_{\mathrm{ind}}(r)$ and the external charge density of the light nuclei $\rho_{\mathrm{ext}}(r)$

$$-\nabla^{2}\phi(r)=\rho_{\mathrm{tot}}(r)/\varepsilon_{0}=[\rho_{\mathrm{ext}}(r)+\rho_{\mathrm{ind}}(r)]/\varepsilon_{0}\,,$$ (8)

where $\varepsilon_{0}$ is the vacuum permittivity. The equilibrium-induced charge density is the difference between the charge density of electrons and positrons

$$\rho_{\mathrm{ind}}(r)=en_{+}(r)-en_{-}(r)\,,$$ (9)

where $n_{\pm}(r)$ represents the number density of electrons and positrons. The induced charge density can be calculated by assuming some statistical distribution for charges in the plasma. For the choice of Boltzmann statistics, one recovers the usual Poisson-Boltzmann equation.

We introduce the re-scaled potential

$$\Phi(r)\equiv\frac{e\phi(r)}{T}\,,$$ (10)

to create a dimensionless variable $\Phi$ which compares the Coulomb energy $e\phi$ to the plasma temperature $T$. The re-scaled Poisson equation reads

$$-\nabla^{2}\Phi(r)-e\rho_{\mathrm{ind}}(r)/(\varepsilon_{0}T)=e\rho_{\mathrm{ext}}(r)/(\varepsilon_{0}T)\,.$$ (11)

We introduce the re-scaled external charge distribution $P_{\mathrm{ext}}$ to simplify the right-hand side

$$P_{\mathrm{ext}}(r)\equiv e\frac{\rho_{\mathrm{ext}}(r)}{\varepsilon_{0}T}=\frac{4\pi Z\alpha\hbar c}{\pi^{3/2}R^{3}T}e^{-\frac{r^{2}}{R^{2}}}\,,$$ (12)

where we chose to model the charge distribution of a nucleus as a Gaussian and $R$ is the root-mean-squared charge radius

$$R=\sqrt{\frac{2}{3}}R_{\alpha}\,.$$ (13)

For the radius of ⁴He we use the charge radius $R_{\alpha}=1.67824\,$fm (Krauth et al., 2021). For the radius of ¹²C, we used an equilateral cluster of alphas (Smith et al., 2020). We rewrite the Poisson equation Eq. (8) in terms of an effective screening mass to make the screening strength more explicit. The screening mass is defined as

$$\frac{m_{s}^{2}(\Phi)}{(\hbar c)^{2}}\equiv-\frac{e\rho_{\mathrm{ind}}(r)/(\varepsilon_{0}T)}{\Phi(r)}\,.$$ (14)

Using this screening mass Eq. (11) takes the form

$$-\nabla^{2}\Phi(r)+\frac{m_{s}^{2}(\Phi)}{(\hbar c)^{2}}\Phi(r)=P_{\mathrm{ext}}(r)\,.$$ (15)

The screening mass $m_{s}^{2}$ is related to the usual Deybe mass but exhibits additional nonlinear behavior due to its dependence on the potential. Using this generalized screening mass, we will compare Boltzmann and Fermi-Dirac equilibrium distributions for the plasma.

2.1 Boltzmann Statistics

The full equilibrium Boltzmann distribution necessary to describe the short-range potential is (Hakim, 1967; Groot et al., 1980)

$$f_{B}^{\pm}(x,p)=e^{-\left\{u_{\mu}[p^{\mu}+qA^{\mu}(x)]\right\}/T}\,.$$ (16)

The charge of a plasma particle is $q$, $A^{\mu}(x)$ is the electromagnetic 4-potential, and $T$ is the plasma temperature. During the BBN epoch $86\,\mathrm{keV}>T>20\,\mathrm{keV}$, the number of electrons and positrons are almost equal to each other due to the charge neutrality (Grayson et al., 2023). In this temperature range, we set the chemical potential to zero for our study. We assume the plasma is at rest such that its 4-velocity reads $u^{\mu}=(1,0,0,0)$, then

$$u_{\mu}A^{\mu}(x)=\phi(x)\,.$$ (17)

The potential term in Eq. (16) that accounts for the energy change due to the rest frame potential (or equivalently uses the kinetic vs canonical momentum) to describe the rest energy of electrons and positrons in the plasma (Hakim, 1967; Groot et al., 1980). This phase space density leads to a number density that is the normal expression but enhanced by the exponential factor in potential

$$n_{\pm}(x)=2\int\frac{d^{3}\boldsymbol{p}}{(\hbar c)^{3}(2\pi)^{3}p^{0}}p^{0}f_{B}^{\pm}(x,p)=n_{\mathrm{eq}}e^{-qu\cdot A(x)/T}\,,$$ (18)

where the equilibrium density is defined as

$$n_{eq}\equiv 2\int\frac{d^{3}\boldsymbol{p}}{(\hbar c)^{3}(2\pi)^{3}}e^{-u\cdot p/T}\,.$$ (19)

The induced equilibrium charge density is used to calculate the screening mass Eq. (14).

Calculating the induced charge density for the Boltzmann case using Eq. (9) and Eq. (18) and then inserting them into the definition for the screening mass Eq. (14), one finds

$$m_{s,\mathrm{Boltz}}^{2}(\Phi)=m_{D}^{2}\frac{\sinh\left[\Phi(r)\right]}{\Phi(r)}\,,$$ (20)

since the exponential factors in Eq. (18) combine to give hyperbolic sine. In this expression, the Debye screening mas $m_{D}$ has its standard definition

$$\frac{1}{\lambda_{D}^{2}}=\frac{m_{D}^{2}}{(\hbar c)^{2}}=\frac{8\pi\alpha\hbar c}{T}\,n_{\mathrm{eq}}\,.$$ (21)

$\lambda_{D}$ is the characteristic distance over which the linearized field is screened.

2.2 Fermi-Dirac Statistics

We anticipate that the induced charge density near the nuclear surface where $e\phi\geq m_{e}$ may become large enough for degeneracy effects to become important. To model this correctly, we use the relativistic Fermi-Dirac distribution

$$f_{F}^{\pm}(x,p)=\frac{1}{e^{\left[u_{\mu}(p^{\mu}+qA^{\mu}(x))\right]/T}+1}\,.$$ (22)

We then sum the induced charge densities of electrons and positrons as in Eq. (9). After simplification, one finds

$$m^{2}_{s,\mathrm{FD}}(\Phi)=\frac{8\pi\alpha}{T\Phi(r)}\int\frac{d^{3}\boldsymbol{p}}{(2\pi)^{3}}\frac{\sinh\left[\Phi(r)\right]}{\cosh\left(p_{0}/T\right)+\cosh\left[\Phi(r)\right]}\,,$$ (23)

which does not have a simple analytic form in the full relativistic limit due to the dependence on the re-scaled potential $\Phi$ (Frolov, 2023). An analytic solution to Eq. (23) exists in the ultrarelativistic limit where one can approximate $p_{0}=\sqrt{m^{2}+\boldsymbol{p}^{2}}\approx|\boldsymbol{p}|$ or equivalently $m/T\ll 1$. Eq. (23) is then integrated analytically into the known form (Elze et al., 1980) valid to all orders in $\Phi(r)$

$$\begin{split}m^{2}_{s,\mathrm{FD}}(\Phi)&\approx\frac{4\alpha T^{2}}{3\pi}\left[\pi^{2}+\left(\Phi(r)+\frac{\mu}{T}\right)^{2}\right]\,.\end{split}$$ (24)

This is the usual ultra-relativistic Debye mass plus a quadratic dependence in $\phi$. This screening mass has much slower quadratic growth than the exponential dependence in $\Phi$ predicted by the Boltzmann distribution Eq. (20).

Keeping higher order terms in a low density expansion $z={m}/({\Phi T})\ll 1$ [see (Birrell et al., 2024)] and for $\mu=0$ one finds (Kodama, 2002)

$$m^{2}_{s,\mathrm{FD}}(\Phi)\approx\frac{4\alpha}{3\pi}\left(\frac{\pi^{2}}{2}T^{2}\frac{2-z^{2}}{\sqrt{1-z^{2}}}+(\Phi T)^{2}(1-z^{2})^{3/2}+\frac{7\pi^{4}T^{4}}{40(\Phi T)^{2}}\frac{z^{4}}{(1-z^{2})^{5/2}}\right)\,.$$ (25)

Figure 1: The effective screening mass for both Fermi Eq. (22) and Boltzmann Eq. (16) statistics are shown as a blue-dashed line and a red-dotted line, respectively. We also show the ultrarelativistic expansion as an orange dashed-dotted line and the series expansion derived in (Kodama, 2002) shown as a purple long-dashed line. Strong screening predicts a much larger screening effect at large $\Phi=e\phi/T$ values than weak screening. The Gamow energy $E_{G}\approx 390$ keV, the most probable tunneling energy, is shown as a gray vertical line, and the value of the potential at the origin $e\phi(0)\approx 2300\,$keV is shown as a dashed gray line.

In Fig. 1, we compare the effective screening mass Eq. (14) for a Fermi-Dirac distribution and a Boltzmann distribution for varying values of the potential $\phi$. In general, a larger screening mass indicates more screening charge will be present for a particular external potential $\phi$, leading to a larger enhancement of nuclear reaction rates. In Fig. 1, the Boltzmann distribution vastly over-predicts the screening mass and, thus, the screening effect. This is because it does not include the stacking of fermion states of electrons and positrons in the polarization cloud. The domain of the screening mass relevant for quantum tunneling lies between the gamow energy and the potential value near the origin $E_{G}<e\phi<e\phi(0)$. In Fig. 1, we see that in this region, strong screening predicts around a $10>m^{2}_{s,\mathrm{FD}}/m_{D}>100$ times larger screening mass as compared to weak screening. This is a much larger change in the screening effect than predicted by dynamic collision-less screening (Hwang et al., 2021) and damped dynamic screening (Grayson et al., 2023). For large values of potential $\phi$, the Fermi-Dirac screening mass $m_{s,\mathrm{FD}}$ approaches the ultrarelativistic limit Eq. (24). This indicates that the analytic expression for the ultra-relativistic Fermi-Dirac screening mass fully captures the short-distance potential and the behavior of a high-temperature strong field plasma in general. This could potentially lead to huge simplification of analytic methods. To numerically solve Eq. (15) we must use the relativistic Fermi-Dirac screening mass Eq. (23) instead of the analytic result Eq. (24), since our boundary conditions are given at large distances where Eq. (24) is invalid.

3 Solving the strong field Poisson equation self consistently

To find the numerical solution to Eq. (15) we start the integration of Eq. (15) at large values of $r$ where the weak field condition Eq. (2) applies since there we know the analytic solution [see Appendix A, Eq. (A3)] up to an overall normalization due to the additional screening that will occur at short ranges. For this calculation, we typically began integration at 1 Å, using the weak field analytic solution,

$$\Phi_{r\gg\hbar c/m_{s}}(r)\approx\frac{Q_{\mathrm{eff}}\,e^{-m_{D}r}}{r}\,.$$ (26)

as a boundary condition at large distances [for better accuracy we use the Debye-Hückel field for Gaussian chargers Eq. (A2)]. One can integrate inwards with any standard solver to small values of $r$ varying a shooting parameter $Q_{\mathrm{eff}}$, which represents the effective charge at large distances. The solution will either diverge to positive or negative infinity based on whether the effective charge has been chosen to be too large or too small. One can then iterate the shooting algorithm, increasing or decreasing the effective charge seen at large distances until the divergence point converges sufficiently close to $r=0$.

The numerical solution can be checked for consistency as discussed in Appendix B. However, the precision of the solution was mainly determined by ensuring the stability of the solution when varying parameters such as step size, boundary conditions, and initial guesses for the shooting algorithm. The overall accuracy of the solution is determined by the choice of boundary conditions, which we assume is the analytic solution at large distances Eq. (A2).

For the Fermi-Dirac case, there is additional complexity since the density cannot be easily integrated analytically over momentum $p$. However, since the dependence on the potential is trivial, one can parameterize a numerical solution for various values of potential $\Phi$ and then use the resulting parameterization in the Poisson equation Eq. (15). This parameterization is plotted in Fig. 1.

Figure 2: The potential $e\phi_{\text{ind}}$ due to the induced screening charge density for both Fermi Eq. (22) and Boltzmann Eq. (16) statistics, including strong screening shown as a blue dashed line and a red dotted line respectively. The black solid line is the screening potential for the weak-field limit Eq. (A2). This calculation was done at temperature $T=86\,$keV at the beginning of BBN when the screening effect is largest. The gray area shows the nuclear interior at a radius of $R=\sqrt{2/3}R_{\alpha}$. The potential calculated using Boltzmann statistics levels off near this radius at $e\phi_{\text{ind}}=170\,$keV, outside the scope of this plot.

Figure 3: The potential $e\phi_{\text{ind}}$ due to the induced screening charge density for Fermi-Dirac strong screening Eq. (22) at various temperatures as solid lines ranging from blue to red. The weak screening potentials are shown as dashed lines ranging from blue to red. Overall screening decreases with temperature $T$, but the difference between weak and strong becomes larger for small $T$.

We are interested in the potential due to screening charges $\phi_{\text{ind}}$ where

$$\phi_{\text{ind}}(r)\equiv\phi(r)-\phi_{\text{vac}}(r)\,,$$ (27)

with $\phi$ being the total potential and $\phi_{\text{vac}}$ the potential in vacuum. This quantity is plotted for an alpha particle (⁴He) in Fig. 2. Strong screening predicts a larger induced screening potential due to the enhancement of screening charge density Eq. (18). This effect is more pronounced near the origin where $e\phi/T$ is large. For an alpha particle and $T=86\,$keV, this value is around $e\phi(0)/T\approx 27$. The Boltzmann screening potential, shown as a red dotted line, calculated using Eq. (20) over-predicts the screening effects by a factor of up to $10^{3}$ since the density quickly becomes large enough that a Fermi-Dirac distribution is required to describe it correctly. The Fermi-Dirac screening potential, shown as a blue dashed line, displays the usual step behavior. The black solid line shows the weak screening potential Eq. (A2), which is nearly constant up to large distances and equal to its value at the origin

$$e\phi_{\text{weak}}(0)=Z\alpha m_{D}\,.$$ (28)

The overall screening effect in the potential Fig. 2 is much less than what is naively expected by Eq. (14) since the potential at the origin is related linearly to the screening mass in the weak screening limit. In strong screening, the screening mass in Eq. (14) and plotted in Fig. 1 can be used to accurately predict the induced charge density $\rho_{\text{ind}}$ but not the potential, which can only be found after solving Eq. (15).

In Fig. 3, the screening potential $e\phi_{\text{ind}}$ is plotted at decreasing temperatures, the BBN temperature range being $T=86.2-50\,$keV (Pitrou et al., 2018). Temperatures lower than $20\,$keV are not considered because, at this temperature in the early Universe, positrons begin to disappear (Grayson et al., 2023). While the overall strength of screening decreases with temperature, the difference between weak and strong screening is much more pronounced at low temperatures, shown in blue, due to the large value of $e\phi/T$ near the origin. Thus, it is much more important to consider the effect of strong screening at low temperatures.

4 Enhancement factor for nuclear reaction rate

The thermonuclear reaction rate for a reaction process involving two nuclei is given by (Rose, 1998)

$$R_{12}=\chi n_{eq}^{1}n_{eq}^{2}I$$ (29)

where $\chi$ is a symmetry factor corresponding to $1/2$ for identical particles and $1$ for distinguishable, $n_{eq}^{1}$ and $n_{eq}^{2}$ are the equilibrium densities of the reacting light nuclei, which in the early universe we can assume to follow Boltzmann distribution since their mass is much larger than the temperature which is in the range of $86-50\,$keV (Pitrou et al., 2018). The probability of the reaction occurring is determined by the integral

$$I=\sqrt{\frac{8}{\mu_{r}\pi T^{3}}}\int_{E_{\mathrm{th}}}^{\infty}dE\,\sigma(E)E\exp{\left(-E/T\right)}\,,$$ (30)

where $\sigma(E)$ is the reaction cross-section, $\mu_{r}$ is the reduced mass of the two colliding ions, $E_{th}$ is the threshold energy for the reaction, and we assume a Boltzmann distribution for the relative energy $E$ of the collision. In a plasma environment, the screened cross-section is

$$\begin{split}\sigma_{\text{sc}}(E)&=\!\frac{S(E)}{E}\exp{\left(\!-\frac{2\sqrt{2\mu_{r}}}{\hbar c}\!\!\int_{R}^{r_{c}}\!\!\!dr\sqrt{U_{\text{sc}}(r)\!-\!E}\right)}\,.\end{split}$$ (31)

where $S(E)$ is the usual tabulated S-factor, which removes the $1/E$ dependence of the cross-section and the Coulomb penetration probability from the cross-section. The tunneling probability through the Coulomb barrier for the s-wave potential from the classical turning point $r_{c}$ to the nuclear radius $R$ is given by the WKB approximation. $U_{\text{sc}}(r)$ is the inter-nuclear potential energy of the screened Coulomb potential.

We are interested in finding the effect of the screened reaction rate involving the integral

$$I_{\text{sc}}=\sqrt{\frac{8}{\mu_{r}\pi T^{3}}}\int_{E_{\mathrm{th}}}^{\infty}dES(E)e^{-g_{\text{sc}}(E)}\,,$$ (32)

where

$$g_{\text{sc}}(E)\equiv E/T+\frac{2\sqrt{2\mu_{r}}}{\hbar c}\int_{R}^{r_{c}}dr\sqrt{U_{\text{sc}}(r)-E}\,.$$ (33)

To calculate the integral in Eq. (32), we use the usual saddle point approximation to isolate the contribution of Coulomb penetration probability.

$$I_{\text{sc}}\approx\sqrt{\frac{16}{\mu_{r}T^{3}g_{\text{sc}}^{\prime\prime}(E_{G})}}S(E_{G})e^{-g_{\text{sc}}(E_{G})}\,,$$ (34)

with Gamow energy $E_{G}$ defined in Eq. (4). We numerically calculate Eq. (32) using the saddle point approximation and compare it to various screening potential energy models, specifically weak screening to strong screening. We do not compare to Boltzmann statistics since it vastly over-predicts the screening effect.

The potential energy between the two nuclei is related to the potential $\phi$ calculated in Sect. 3 by the self energy

$$U_{\text{sc}}(r)=\frac{1}{2}\int d^{3}r^{\prime}\rho(r^{\prime})\phi(r-r^{\prime})-U(r\rightarrow\infty)=\frac{1}{2}\int d^{3}r^{\prime}\left[\rho_{1}(r^{\prime})\phi_{2}(r-r^{\prime})+\rho_{2}(r^{\prime})\phi_{1}(r-r^{\prime})\right]\,.$$ (35)

Where we have subtracted the in-plasma self-energies of the two nuclei separated at $r\rightarrow\infty$

$$U(r\rightarrow\infty)=\frac{1}{2}\int d^{3}r^{\prime}\left[\rho_{1}(r^{\prime})\phi_{1}(r-r^{\prime})+\rho_{2}(r^{\prime})\phi_{2}(r-r^{\prime})\right]\,.$$ (36)

The total charge density is the sum of the induced plus the external charge

$$U_{\text{sc}}(r)=\frac{1}{2}\int d^{3}r^{\prime}\left\{\left[{\rho_{\text{ext}}}_{1}(r^{\prime})+{\rho_{\text{ind}}}_{1}(r^{\prime})\right]\phi_{2}(r-r^{\prime})+\left(1\leftrightarrow 2\right)\right\}\,.$$ (37)

Then, we use Eq. (14) to rewrite the induced charge density in terms of the total potential

$$U_{\text{sc}}(r>R)=\frac{1}{2}Z_{1}\phi_{2}(r)+\frac{1}{2}\int d^{3}r^{\prime}\left[m_{s}^{2}(\phi_{1})\phi_{1}(r^{\prime})\phi_{2}(r-r^{\prime})\right]+\left(1\leftrightarrow 2\right)\,.$$ (38)

We approximate $\rho_{\text{ext}}$ as a point charge at distances larger than $R$, relevant for the tunneling integral in Eq. (31). Often, in the weak screening limit, one only considers the charge of colliding nuclei with the total potential of the other. If the charges are approximately symmetric, the screening potential is then

$$U_{\text{sc}}(r)\approx Z_{2}\phi_{1}(r)\,.$$ (39)

Since we only want to compare weak screening to strong screening, for simplicity, we neglect the second term in Eq. (38). This approximation is only good for small $Z$ when the value of $m_{s}$ is much smaller than the vacuum potential. To calculate the exact potential energy in this model one would need to calculate the screening of two nuclei as they approach since the the principle of superposition does not apply when nonlinear screening is present. In the weak screening limit, the second term in Eq. (38) is known to lead to a factor $3/2$ in the Salpeter’s usual enhancement factor Eq. (44) (Brüggen & Gough, 1997).

We will compare Eq. (32) to the same integral in vacuum, assuming the electric potential used to find the potential energy is given by Eq. (A3)

$$I_{\text{vac}}=\sqrt{\frac{8}{\mu_{r}\pi T^{3}}}\int_{E_{\mathrm{th}}}^{\infty}dES(E)e^{-g_{\text{vac}}(E)}\,,$$ (40)

where

$$g_{\text{vac}}(E)\equiv E/T+\frac{2\sqrt{2\mu_{r}}}{\hbar c}\int_{R}^{r_{c}}dr\sqrt{U_{\text{vac}}(r)-E}\,.$$ (41)

In this expression, $U_{\text{vac}}$ is the vacuum internuclear potential. We represent the enhancement of the thermonuclear reaction rate as

$$\mathcal{F}_{\text{sc}}\equiv\frac{R_{\text{sc}}}{R_{\text{vac}}}=\frac{I_{\text{sc}}}{I_{\text{vac}}}\,.$$ (42)

In Appendix C, we mention the effect of finite size on reaction rates, which is much larger than screening and essential to include.

When the Debye length $\lambda_{D}$ is very large compared to classical turning point $r_{c}$, the screening charge density can be approximated to be constant near the origin. In that case, the screening effect can be described as a shift in the energy of the reaction rate leading to a simple exponential enhancement referred to as the Salpeter enhancement factor given by the change in the potential energy at the origin.

$$\mathcal{F}_{\text{sc}}(0)=\exp{\left(\frac{U_{vac}(0)-U_{\text{sc}}(0)}{T}\right)}=\exp{\left(\frac{H(0)}{T}\right)}\,.$$ (43)

$H(0)$ is often used in literature to represent the electric potential energy of the induced charge density. Alternative derivations of Eq. (44) lead to different forms of the Slapeter enhancement factor, summarized in (Bahcall et al., 2002). Since we compare the enhancement factor for strong screening to weak screening, the exact form of $\mathcal{F}_{\text{sc}}$ is not important since multiplicative factors will cancel. For the weak screening limit using only the leading order potential energy Eq. (39), one finds (Salpeter, 1954)

$$\begin{split}\mathcal{F}_{\text{weak}}(0)=\exp\left(\frac{Z_{1}Z_{2}\alpha m_{D}}{T}\right)\,.\end{split}$$ (44)

Figure 4: Reaction rate enhancement $\mathcal{F}_{\text{sc}}$ [see Eq. (42)] for ⁴He-⁴He scattering as a function of temperature. The weak screening model is shown as a black solid line, and the strong screening model is plotted as a blue dashed line. The approximate enhancement factors using the screening potential energy at the origin Eq. (44) are shown as a red dotted line for weak screening and an orange dashed line for strong screening. The orange shaded region marks the BBN temperature range $T=50-86.2\,$keV.

In Fig. 4, we calculate the enhancement due to screening in the BBN plasma environment by numerically calculating the ratio $\mathcal{F}_{\text{sc}}$ from Eq. (42) using the weak the screening potential Eq. (A2) shown as a black solid line and the strong screening potential, found by solving Eq. (15) with Fermi-Dirac statistics, shown as a blue dashed line. This is done by calculating the WKB tunneling integral in Eq. (33) using the numerical solutions for the potential plotted in Fig. 3 and then estimating the integral Eq. (32) using the saddle point approximation Eq. (34). The blue dashed line and the black solid line include the finite size of the colliding nuclei and the screening potential’s spatial dependence. In the calculation of Eq. (34), the change in the Gamow energy due to screening is calculated for completeness, but, in general, this change is negligible since screening does not change the distribution of ions or the penetration probability significantly.

Strong screening generally predicts a larger enhancement of reaction rates. However, for ⁴He-⁴He collisions, this enhancement is very small, as indicated by the blue dashed line almost tracing the black solid line. This is expected since predicted ⁴He abundances match theoretical measurements (Pitrou et al., 2018). The discrepancy between weak and strong screening becomes a constant at high temperatures $\sim 10^{-6}$, indicating the polarization density is approaching the ultrarelativistic limit of Fermi-Dirac statistics Eq. (24).

The Salpeter enhancement factor $\mathcal{F}_{\text{weak}}(0)$ reasonably approximates the enhancement due to weak screening, as seen by the red dotted line tracing the black solid line. The Salpeter enhancement factor for strong screening is shown as an orange dashed line in Fig. 4 using the value of the strong screening potential in Fig. 3 at the origin. One can see in Fig. 4 that the Salpeter approximation for screening $\mathcal{F}_{\text{strong}}(0)$ overestimates the numerical solution. This is because the induced charge density is not constant near the origin but has a step in the potential created by the Fermi function in the calculation of the strong screening potential. This overestimation of the Salpeter enhancement factor is also found in (Itoh et al., 1977). The Salpeter enhancement factor is a less accurate approximation at low temperatures where the jump in potential is much more pronounced.

Figure 5: Reaction rate enhancement $\mathcal{F}_{\text{sc}}$ [see Eq. (42)] for ¹²C-¹²C scattering as a function of temperature. The weak screening model is shown as a black solid line, and the strong screening model is plotted as a blue dashed line. The approximate enhancement factors using the screening potential energy at the origin Eq. (44) are shown as a red dotted line for weak screening and an orange dashed line for strong screening. The orange-shaded region marks the BBN temperature range $T=50-86.2\,$keV.

As an informative exercise, we show the enhancement factor for ¹²C-¹²C fusion in Fig. 5. This choice of reaction is to connect to stellar fusion in white dwarfs. We can see that at higher $Z$, the strong screening is becoming more important due to its nonlinear dependence on $Z$. Here, the enhancement at large temperatures is compared to weak screening is $\sim 5\times 10^{-4}$. Compared to ⁴He-⁴He, This is a $\sim 500$ times the stronger effect for an increase in $Z$ of only 4. Although this effect is too small to account for issues in the abundance of light elements, such as the Lithium problem (Pitrou et al., 2018), strong screening is very interesting because it can affect reactions with higher Z while leaving low Z reactions unchanged.

The enhancement predicted by strong screening at the origin, shown as the orange dashed line in Fig. 4, deviates more since the step in the screening potential is more pronounced at higher $Z$. In Fig. 5, the strong screening enhancement at low temperatures is offset from the weak field limit because the increase in $\phi/T$ compensates for the decrease in the Debye mass.

5 Conclusion and discussion

5.1 Summary

In this study, we implemented self-consistent strong electron-positron plasma screening to determine the short-range screening potential relevant to thermonuclear reactions during the Big Bang nucleosynthesis (BBN) epoch. This research direction was originally proposed in our previous work on BBN screening (Grayson et al., 2023). Our approach included incorporating finite-sized nuclei into the Poisson-Boltzmann equation and introducing a generalized screening mass to facilitate the comparison of various screening models.

Our results in Fig. 4 and Fig. 5 indicate that strong screening, while generally small due to the high temperatures prevalent during BBN, can enhance the fusion rates of high-$Z$ elements while leaving lower-$Z$ elements relatively unaffected.

One of our key findings is the necessity of using relativistic Fermi-Dirac statistics to describe the polarization cloud accurately at short distances, particularly when the electromagnetic work performed by ions on electrons is comparable to their rest mass energy creating a large polarization density around the nucleus. This is true even when the plasma adheres to Boltzmann statistics at larger distances.

Additionally, we show in Fig. 3 that the spatial dependence of the strong screening potential near the nuclear surface is very pronounced, having a Fermi-function “step shape.” Relying solely on the screening energy near the origin can lead to overestimating reaction rate enhancement due to strong screening.

In this work we accounted for strong screening, the finite size of nuclei, and finite tunneling distances to calculate reaction rate enhancement in the BBN epoch for two simple cases of ⁴He-⁴He and ¹²C-¹²C collisions. These improvements provide a more precise framework for understanding plasma screening effects during BBN and have broader applications in studying weakly coupled plasmas in various cosmic and laboratory settings.

5.2 Outlook

Usually, a formulation of strong field dynamics in the semi-classical regime using the Vlasov-Boltzmann equation is intractable due to the nonlinear exponential behavior of the induced charge density in the plasma. We think that the novel analytic form of the Fermi-Dirac screening mass in the ultrarelativistic limit Eq. (24) could lead to a solvable strong field theory of plasma dynamics in the high-temperature regime. Our strong screening approach could be relevant to the dynamics of heavy quarks in quark-gluon plasma, where light quarks are ultrarelativistic and strongly screen slow-moving heavier quarks (Mrówczyński, 2018; Carrington et al., 2020).

Strong screening has recently gained interest in the context of laser plasma fusion, as demonstrated by studies predicting in the context of the Boltzmann limit large enhancements of fusion reactions as compared to weak screening at low temperatures and high densities (Elsing et al., 2022).

Additionally, we are interested in exploring nuclear fusion catalysis by heavy ions. Strong screening predicts a significant screening cloud in proximity to large $Z$ nuclei. Light elements could collide in an environment near the heavy nuclei where their repulsion is mostly neutralized, leading to a consequential enhancement of fusion reactions.

In future studies, we are interested in exploring how the pion exchange impacts the nuclear attraction to refine our understanding of the nuclear charge distribution and the effective tunneling distance. The effective nuclear radius $R$ used in tunneling is an important factor affecting strong screening. Considering a more accurate model of nuclear size would alter the tunneling distance through potential wall.

Historically, experimental determinations of astrophysical S-factors have found anomalous screening at low temperatures (Shoppa et al., 1993). This anomalous screening was measured again recently (Zhang et al., 2020). We speculate that this anomalous screening at low collision energy could be explained by strong screening effects obtained in this work.

The highly nonlinear response of the strong screening effect to $Z$ suggests that strong screening could help address the observed lithium over-abundance problem by enhancing the fusion of higher-$Z$ elements (Pitrou et al., 2018). However, based on the screening enhancements obtained in this study, an explanation of the observed Lithium abundance would require an additional screening mechanism, such as a combined strong-dynamic screening theory.

6 Acknowledgements

This research did not receive specific grants from funding agencies in the public, commercial, or not-for-profit sectors. However, we thank Tamás Biró for his hospitality at the PP2024 Particles & Plasmas Symposium, where this work was presented.

Appendix A Solving the Poisson-Boltzmann equation in the weak-field limit

In the weak field approximation $Ze\phi/T\ll 1$, the hyperbolic sine term in Eq. (20) can be approximated to first order

$$-\frac{d^{2}}{dr^{2}}\Phi(r)-\frac{2}{r}\frac{d}{dr}\Phi(r)+\frac{1}{\lambda_{D}^{2}}\Phi(r)=P_{\mathrm{ext}}(r)\,,$$ (A1)

for the external charge distribution given in Eq. (12) this has an analytic solution that can be found via Fourier-transform-methods using the boundary condition of vanishing potential for large distances (Grayson et al., 2023)

$$\phi_{\text{weak}}(r)=\frac{Ze}{4\pi\varepsilon_{0}}\frac{e^{R^{2}/(4\lambda_{D})}}{r}\left(\frac{\text{Erf}\left(\frac{r}{R}-\frac{R}{2\lambda_{D}}\right)}{2}e^{-r/\lambda_{D}}+\frac{\text{Erf}\left(\frac{r}{R}+\frac{R}{2\lambda_{D}}\right)}{2}e^{r/\lambda_{D}}-\sinh(r/\lambda_{D})\right)\,.$$ (A2)

This solution is valid for large distances where the potential is small enough compared to temperature that the exponential term in Eq. (18) is close to one. The vacuum result for $m_{D}\rightarrow 0$ or $\lambda_{D}\rightarrow\infty$

$$\phi_{\text{vac}}(r)=\frac{Ze}{4\pi\varepsilon_{0}}\frac{\text{Erf}\left(\frac{r}{R}\right)}{r}\,,$$ (A3)

is often used to model the Coulomb field of ⁴He (Kumar et al., 2022).

Appendix B Checking the numerical solution

One can check for consistency of the solution by integrating the solution to find the total charge and comparing it to the effective charge at large distances represented by the large distance behavior of the solution. To find the total charge predicted by the numerical solution, one rearranges the Poisson equation

$$-\nabla^{2}\Phi(r)+\frac{m_{D}^{2}}{(\hbar c)^{2}}\Phi(r)=\left[P_{\mathrm{ext}}(r)+\left(\frac{m_{D}^{2}}{(\hbar c)^{2}}-\frac{m_{s}^{2}(r)}{(\hbar c)^{2}}\right)\Phi(r)\right]\,.$$ (B1)

We then identify the right-hand side as the source for the left-hand side

$$\begin{split}Q_{\mathrm{eff}}=\int dV\Big{[}P_{\mathrm{ext}}(r)+\left(\frac{m_{D}^{2}}{(\hbar c)^{2}}-\frac{m_{s}^{2}(r)}{(\hbar c)^{2}}\right)\Phi(r)\Big{]}\,.\end{split}$$ (B2)

The large distance solution to the left-hand side is of the Debye-Hückel form

$$\Phi_{r\gg\hbar c/m_{s}}(r)=\frac{Q_{\mathrm{eff}}\,e^{-m_{D}r}}{r}\,.$$ (B3)

Then the consistency of the numerical solution can be checked by comparing the value of the effective charge found from Eq. (B2) and Eq. (B3).

Appendix C Enhancement due to finite size

The textbook formula for the penetration probability of tunneling to some finite distance $R$ at the surface of the nuclei is (Griffiths & Schroeter, 2018)

$$\mathcal{P}(R)=\exp{\left[\frac{2\sqrt{2\mu_{r}E_{G}}}{\hbar c}\left(\frac{\pi}{2}r_{c}-2\sqrt{Rr_{c}}\right)\right]}\,.$$ (C1)

Comparison can also be made with point nuclei charge, where the finite size of the nuclei is neglected, and the tunneling distance goes from the classical turning point to the origin

$$\sigma_{\text{Som}}(E)=\frac{S(E)}{E}\exp{\left[-2\pi\eta(E)\right]}\,,$$ (C2)

where $\eta(E)$ in the exponential is often referred to as the Sommerfeld parameter

$$\eta(E)=Z_{1}Z_{2}\alpha\sqrt{\frac{\mu_{r}}{2E}}\,.$$ (C3)

For the integral Eq. (30) using Eq. (C2) one finds the usual

$$I_{\text{Som}}=4\sqrt{\frac{2E_{G}}{3\mu}}\frac{S(E_{G})}{T}e^{-3E_{G}/T}\,.$$ (C4)

The ratio of the finite size expression to the Sommerfeld parameter is

$$\frac{\mathcal{P}(R)}{e^{-2\pi\eta}}\approx\exp\left(4\sqrt{2Z_{1}Z_{2}\alpha\frac{R\mu_{r}}{\hbar c}}\right)\,,$$ (C5)

which is around 32.8 for alpha-alpha collisions during BBN, a substantial enhancement compared to the effect produced by screening.