Department of Physics, The University of Arizona, Tucson, Arizona 85721, USA
ELI Beamlines Facility, The Extreme Light Infrastructure ERIC, 252 41 Dolni
Brezany, Czech Republic

# Electron-positron plasma in BBN:
damped-dynamic screening

Christopher Graysona, Cheng Tao Yanga,
Martin Formanekb, Johann Rafelskia

(August 17, 2023)

###### Abstract

We characterize in detail the dense $e^{-}e^{+}\gamma$ plasma present during
the Big-Bang Nucleosynthesis (BBN) and explore how it is perturbed
electromagnetically by “impurities, i.e., spatially dispersed protons and
light nuclei undergoing thermal motion. The internuclear electromagnetic
screened potential is obtained (analytically) using the linear response
approach, allowing for the dynamic motion of the electromagnetic field sources
and the damping effects due to plasma component scattering. We discuss the
limits of the linear response method and suggest additional work needed to
improve BBN reaction rates in the primordial Universe. Our theoretical methods
describing the potential between charged dust particles align with previous
studies on planetary and space dusty plasma and could have a significant
impact on the interpretation of standard cosmological model results.

##  1 Introduction

Electron-positron $e^{-}e^{+}$-pair-plasma has a significant presence in the
early Universe during Big-Bang nucleosynthesis (BBN) [1, 2, 3]. Charge
screening is one of the important collective plasma effects modifying the
internuclear potential $\phi(r)$. The electron cloud surrounding the charge of
an ion screens other nuclear charges far from their own radius and reduces the
Coulomb barrier. In nuclear reactions, reducing the Coulomb barrier increases
the penetration probability, enhances thermonuclear reaction rates, and thus
modifies light element abundances in the early universe.

The screening effect in nuclear reactions has a long and distinguished
history. Plasma screening effects were first considered in 1954 by Salpeter
[4], who proposed evaluating the enhancement of nuclear reactions by employing
the static Debye-Hückel potential [5, 6, 7]. These applications focus on
collision-less electron-positron plasma only. Later, this approach was
generalized for nuclei moving in the plasma [2, 8, 9, 10, 11] i.e., ‘dynamic’
screening. In this work, we address the ultra dense nature of the quantum
electrodynamics (QED) $e^{-}e^{+}\gamma$ plasma in which the BBN reaction
network occurs. Only recently has the ultra dense nature of the BBN plasma
been taken into account by the inclusion of collisions in numerical
simulations [12, 13]. The dense nature of the the BBN plasma is important
since the BBN reactions typically take place within the temperature interval
$86\,\mathrm{keV}>\mathrm{T_{BBN}}>50\,\mathrm{keV}$ [14] which is a non-
negligible fraction of the electron mass energy equivalent, $m_{e}c^{2}=511$
keV. The individual densities of electron-positron pairs are higher than the
electron density in the solar core and surpass densities found in most stellar
burning environments. This dense plasma environment indicates the need to
understand the effect of damping due to scattering of plasma components
$e^{-}e^{+}\gamma$ during the BBN epoch.

To improve the understanding of the internuclear screened Coulomb potential
$\phi$ governing BBN reaction rates, this work expands the prior studies [2,
8, 10, 11] carried out in the linear response approximation applicable in the
regime where the modification of the Coulomb potential is a small effect. We
extend the ’dynamic screening’ effort (and in particular that of Hwang et al
[2] which introduced dynamic screening during BBN) to study ‘damped-dynamic’
screening as follows:

  1. 1.

We detail the scattering damping effect of plasma components in the
temperature range beginning near $T\simeq m_{e}c^{2}=511$ keV extending down
to a temperature near $T_{\mathrm{split}}\simeq 20.3$ keV where practically
all $e^{-}e^{+}$-pairs disappeared.

  2. 2.

We obtain analytical results for damped-dynamic screening applicable to the
BBN epoch temperature range which demonstrates a connection between the BBN
plasma and dusty plasma theory.

  3. 3.

We recognize and describe the limits of applicability of the linear response
method to BBN by considering the full equilibrium distribution necessary in
the presence of strong fields. We explore a one-dimensional toy model at
distances a fraction of an Ångstrom in the quantum tunneling regime where the
potential $\phi$ is large compared to the thermal collision energy
$|\phi|/3T>1$, which predicts enhanced screening at small distances.

Such effort to improve the BBN model through screening is well justified for
the following reasons:

  1. 1.

Currently, the most accessible path to observe the early Universe before
recombination is through light element abundances.

  2. 2.

Despite the great initial success of the BBN model regarding the prediction of
the abundances of D, 3 He, 4 He, there remain well-documented discrepancies.
Currently, substantial efforts are being directed toward reconciling the
theoretical BBN model with present-day observations. For a more detailed
discussion, see Refs. [15, 16].

  3. 3.

Enhanced production of light elements during the BBN era could help resolve
two problems: a) Elements Beryllium and Boron cannot be produced in stellar
burning processes while their high abundance suggests additional mechanisms of
production beyond secondary processes after stellar explosions. b)
Understanding the age of stars obtained using the abundance of light elements
(metallicity) may be improved by adaptation of the BBN network to account for
nuclear reactions in ultra-dense QED plasma.

  4. 4.

We note that screening effects are more pronounced for nuclei with greater
charge $Z$ leading to stronger modification in predicted abundances, just
where the most significant current disagreements exist (e.g. the abundance of
7Li). Inclusion of screening effects in our current models of the BBN reaction
network could explain this discrepancy without the need to invoke more exotic
modifications of the Coulomb internuclear potential. An example of such a
modification would be the variation of $\alpha$, the fine-structure constant
[17].

  5. 5.

Recent James Webb space telescope-Hubble space telescope (JWST-HST)
observations implying early onset of galactic structure formation and/or
presence of luminous objects as early as 300 million years after BBN [18, 19,
20] could be seen as evidence supporting nuclear dust matter clustering in
dense QED plasma.

In the standard BBN model, thermonuclear reaction rates are evaluated using
nuclear reactions specific to the vacuum state. In a high-precision BBN model,
the effect of $e^{-}e^{+}$-pair-plasma screening of nuclear charges must be
allowed. We describe in Section 2 the presence of up to several millions of
electron-positron ($e^{-}e^{+}$) pairs per every charged nucleon. This
situation prompts an effort to reconsider BBN in an ultra dense plasma
environment allowing for a dense $e^{-}e^{+}\gamma$ plasma medium.

We account for collisions between plasma components using the current-
conserving Bhatnagar, Gross, and Krook (BGK) collision term [21], which allows
us to study damping in the dynamic plasma. For an in-depth discussion
presented in contemporary covariant notation and preserving current
conservation, see Formanek et al. [22]. Our approach is different from Monte-
Carlo simulation of two particle collisions [12, 13] as in our case particles
participating in scattering also experience simplified static screening which
in terms of collision description would require more than two particle
collisions. In addition, the simplified BGK collision term we employ allows us
to obtain an analytic ‘damped-dynamic’ internuclear potential relevant to the
BBN reaction network in the linear response approximation.

Several well known processes, including Møller, Bhabha, and inverse Compton
scattering, characterize the damping in plasma. These textbook results will be
adapted to the plasma environment: When studied in vacuum, electrons and
positrons interact with each other by exchanging a massless photon. However,
when a photon propagates through a plasma of electrons and positrons, its
properties are modified by interactions with the medium. When evaluating
Moller and Bhabha scattering, we include the temperature-dependent mass of the
photon obtained in plasma theory without damping. We find that the total
relaxation rate during BBN is much larger than the screening mass, which is
the inverse of the Debye screening length $m_{D}=1/\lambda_{D}$. This suggests
that electromagnetic perturbations of the plasma will be overdamped.

Since, as noted above, the computed damping strength is the dominant scale, it
is also the main parameter determining the photon mass. However, introducing
the damping strength in the photon mass would make the calculation of the
damping strength self-referential. The complexity of such a self-consistent
evaluation of damping is beyond the scope of this work. This underscores the
need to develop a self-consistent approach where both damping and photon
properties in plasma are determined in a mutually consistent manner.

The QED $e^{+}e^{-}\gamma$-plasma present during BBN is perturbed by a variety
of comparatively very heavy charged protons $p$ and even heavier isotopes and
light elements. In calculating the electromagnetic potential $\phi$, we
naturally view BBN active participants as heavier, higher charged impurities
i.e. “charged dust”in the QED-plasma. Since theories of such a system have
been developed by others working on charged dust grains in planetary and space
plasma [23, 24, 25, 26] we can validate our results and adapt previously
developed theoretical tools necessary for describing heavy charged impurities
within a plasma of lighter particles.

We provide the first step in merging these two fields by including damping to
analyze the electromagnetic potential during BBN, similar to work done by
Stenflo in 1973 for dusty plasma [24]. We advance this work by applying it to
the BBN epoch, allowing for the finite size of the ”dust grains,” and by
finding the short-distance behavior of the potential in the weak field limit.
Using a general model of simplified collisional damping in linear response, we
calculate the screened electromagnetic potential of light nuclei undergoing
thermal motion subject to damping relying on previous work [22, 27].

Our manuscript is organized as follows: In Section 2 we derive the electron-
positron density and chemical potential to show that during the normal BBN
temperature range $86\,\mathrm{keV}>\mathrm{T_{BBN}}>50\,\mathrm{keV}$ [14]
the Universe was filled with a dense electron-positron pair-plasma dotted with
dispersed baryonic matter dust. In Section 3, we calculate the relaxation rate
$\kappa=1/\tau$, where $\tau$ is the mean time between collisions in the
plasma. This is done by finding the average of the most relevant reaction
rates for $2\leftrightarrow 2$ scattering: Møller, Bhabha, and inverse Compton
scattering.

We proceed in Section 4 to solve the Vlasov-Boltzmann Equation (VBE) in the
weak field limit with damping. Our approach considers small plasma
perturbations in linear response to find the polarization tensor present
during the BBN era. In Section 5, we find an analytic expression for the
damped-dynamic potential applicable to BBN and study the limitation of the
linear response method. We show that additional strong field phenomena beyond
linear response alter the behavior of the screened potential $\phi$. We
discuss our results and describe their interdisciplinary importance in Section
6.

##  2 Electron chemical potential and number density

In this section, we will derive the dependence of electron chemical potential,
and hence $e^{+}e^{-}$ density, as a function of the photon background
temperature $T$ by employing the following physical principles:

  1. 1.

Charge neutrality of the Universe:

$\displaystyle n_{e^{-}}-n_{{e^{+}}}=n_{p}-n_{\overline{p}}\approx\,n_{p},$
(1)

where $n_{\ell}$ denotes the number density of particle type $\ell$.

  2. 2.

Neutrinos decouple (freeze out) at a temperature $T_{f}\simeq 2$ MeV, after
which they free stream through the Universe with an effective temperature [28]

$\displaystyle T_{\nu}(t)=T_{f}\,\frac{a(t_{f})}{a(t)},$ (2)

where $a(t)$ is the–Friedmann–Lemaître–Robertson–Walker (FLRW) Universe scale
factor which is a function of cosmic time $t$, and $t_{f}$ represents the
cosmic time when neutrino freezes out.

  3. 3.

Total comoving entropy is conserved. At $T\leq T_{f}$, the dominant
contributors to entropy are photons, $e^{+}e^{-}$, and neutrinos. In addition,
after neutrino freeze out, neutrino comoving entropy is independently
conserved [28]. This implies that the combined comoving entropy in
$e^{+}e^{-}\gamma$ is also conserved for $T\leq T_{f}$.

Motivated by the fact that comoving entropy in $\gamma$, $e^{+}e^{-}$ is
conserved after neutrino freezeout, we rewrite the charge neutrality
condition, Eq. (1), in the form

$\displaystyle
n_{e^{-}}-n_{{e^{+}}}=X_{p}\frac{n_{B}}{s_{\gamma,e^{\pm}}}s_{\gamma,e^{\pm}},\qquad
X_{p}\equiv\frac{n_{p}}{n_{B}},$ (3)

where $n_{B}$ is the number density of baryons, $s_{\gamma,e^{\pm}}$ is the
combined entropy density in photons, electrons, and positrons. During the
Universe expansion, the comoving entropy and baryon number are conserved
quantities; hence the ratio $n_{B}/s_{\gamma,e^{\pm}}$ is conserved. We have

$\displaystyle\frac{n_{B}}{s_{\gamma,e^{\pm},}}=\left(\frac{n_{B}}{s_{\gamma,e^{\pm}}}\right)_{t_{0}}\\!\\!\\!\\!=\left(\frac{n_{B}}{s_{\gamma}}\right)_{t_{0}}\\!\\!\\!\\!=\left(\frac{n_{B}}{n_{\gamma}}\right)_{t_{0}}\left(\frac{n_{\gamma}}{s_{\gamma}}\right)_{t_{0}},$
(4)

where the subscript $t_{0}$ denotes the present day value, and the second
equality is obtained by observing that the present day $e^{+}e^{-}$-entropy
density is negligible compared to the photon entropy density. We can evaluate
the ratio by giving the present day baryon-to-photon ratio:
$B/N_{\gamma}=n_{B}/n_{\gamma}=0.605\times 10^{-9}$ from Cosmic Microwave
Background (CMB) [29] and the entropy per particle for a massless boson:
$(s/n)_{\mathrm{boson}}\approx 3.602$.

The total entropy density of photons, electrons, and positrons can be written
as

$\displaystyle
s_{\gamma,e^{\pm}}=\frac{2\pi^{2}}{45}g_{\gamma}\,T^{3}+\frac{\rho_{e^{\pm}}+P_{e^{\pm}}}{T}-\frac{\mu_{e}}{T}(n_{e^{-}}-n_{{e^{+}}}),$
(5)

where $\rho_{e^{\pm}}=\rho_{e^{-}}+\rho_{e^{+}}$ and
$P_{e^{\pm}}=P_{e^{-}}+P_{{e^{+}}}$ are the total energy density and pressure
of electrons and positron respectively.

By incorporating Eq. (3) and Eq. (5), the charge neutrality condition can be
expressed as

$\displaystyle\left[1+X_{p}\left(\frac{n_{B}}{n_{\gamma}}\right)_{t_{0}}\left(\frac{n_{\gamma}}{s_{\gamma}}\right)_{t_{0}}\frac{\mu_{e}}{T}\right]\frac{n_{e^{-}}-n_{{e^{+}}}}{T^{3}}$
$\displaystyle\qquad\qquad\qquad=X_{p}\left(\frac{n_{B}}{n_{\gamma}}\right)_{t_{0}}\left(\frac{n_{\gamma}}{s_{\gamma}}\right)_{t_{0}}\left(\frac{2\pi^{2}}{45}g_{\gamma}+\frac{\rho_{e^{\pm}}+P_{e^{\pm}}}{T^{4}}\right).$
(6)

Using Fermi distribution, the number density of electrons over positrons in
the early Universe is given by

$\displaystyle n_{e^{-}}-n_{{e^{+}}}$
$\displaystyle=\frac{g_{e}}{2\pi^{2}}\left[\int_{0}^{\infty}\frac{p^{2}dp}{\exp{\left((E-\mu_{e})\right)/T}+1}\right.\left.-\int_{0}^{\infty}\frac{p^{2}dp}{\exp{\left((E+\mu_{e})/T\right)}+1}\right]$
$\displaystyle=\frac{g_{e}}{2\pi^{2}}\,{T^{3}}\,\tanh(b_{e})M_{e}^{3}\int_{1}^{\infty}\\!\\!\\!\\!\frac{\eta\sqrt{\eta^{2}-1}d\eta}{1+\cosh(M_{e}\eta)/\cosh(b_{e})},$
(7)

where we have introduced the dimensionless variables as follows:

$\displaystyle\eta=\frac{E}{m_{e}},\qquad M_{e}=\frac{m_{e}}{T},\qquad
b_{e}=\frac{\mu_{e}}{T}.$ (8)

Substituting Eq. (2) into Eq. (2) and giving the value of $X_{p}$, then the
charge neutrality condition can be solved to determine $\mu_{e}/T$ as a
function of $M_{e}$ and $T$.

Figure 1: Left axis: The chemical potential of electrons as a function of
temperature by numerically solving Eq. (2) with $n_{p}/n_{B}=0.878$ and
$n_{B}/n_{\gamma}=6.05\times 10^{-10}$. Right axis: the ratio of electron
(positron) number density to baryon density as a function of temperature. The
solid blue line is the electron density, the red dashed line is the positron
density, and the green dotted line is the number density with $\mu_{e}=0$.
When $T=20.3\,\mathrm{keV}$ (the purple vertical line) positron density
decreases rapidly because of the annihilation. The vertical black dotted lines
represent BBN temperature range
$86\,\mathrm{keV}>\mathrm{T_{BBN}}>50\,\mathrm{keV}$.

In Fig. 1 (left axis), we solve Eq. (2) numerically and plot the electron
chemical potential as a function of temperature with the following parameters:
proton concentration $X_{p}=0.878$ from observation [29] and
$n_{B}/n_{\gamma}=6.05\times 10^{-10}$ from CMB. We can see the value of
chemical potential is comparatively small $\mu_{e}/T\approx 10^{-6}\sim
10^{-7}$ during the BBN temperature range, implying an equal number of
electrons and positrons in plasma. The ratio of electron (positron) number
density to baryon density shows that the Universe was filled with an electron-
positron rich plasma during the accepted BBN temperature range. For example,
when the temperature is around $T=70\,\mathrm{keV}$, the density of electrons
and positrons is comparatively large $n_{e^{\pm}}\approx 10^{7}\,n_{B}$. At
$90$ keV, the electron and positron density is near the solar core density,
see Fig. 19 in Ref. [3]. Later when the temperature is around
$T=20.3\,\mathrm{keV}$, the positron density decreases, transforming the pair-
plasma to an electron-baryon plasma.

##  3 The relaxation rate in electron-positron plasma

In electron-positron plasma, the major reactions between photons and
$e^{+}e^{-}$ pairs are inverse Compton scattering, Møller scattering, and
Bhabha scattering:

$\displaystyle e^{\pm}+\gamma\longrightarrow e^{\pm}+\gamma,\qquad
e^{\pm}+e^{\pm}\longrightarrow e^{\pm}+e^{\pm},\qquad
e^{\pm}+e^{\mp}\longrightarrow e^{\pm}+e^{\mp}.$ (9)

In general, to evaluate the reaction rate in two-body reaction $1+2\rightarrow
3+4$ in the Boltzmann approximation we can use reaction cross-section
$\sigma(s)$ and the relation from Ref. [30]:

$\displaystyle R_{12\rightarrow
34}=\frac{g_{1}g_{2}}{32\pi^{4}}\frac{T}{1+I_{12}}\\!\\!\int^{\infty}_{s_{th}}\\!\\!\\!\\!ds\,\sigma(s)\frac{\lambda_{2}(s)}{\sqrt{s}}\\!K_{1}\\!\\!\left({\sqrt{s}}/{T}\right),$
(10)

where $K_{1}$ is the first-order Bessel function and the Källén function
$\lambda_{2}(s)$ is defined as

$\displaystyle\lambda_{2}(s)=\left[s-(m_{1}+m_{2})^{2}\right]\left[s-(m_{1}-m_{2})^{2}\right],$
(11)

with $m_{1/2}$ and $g_{1/2}$ are the masses and degeneracies of initial
interacting particles. The factor $1/(1+I_{12})$ is introduced to avoid double
counting of indistinguishable pairs of particles, we have $I_{12}=1$ for
identical particles and $I_{12}=0$ otherwise.

The two-body cross-section can be obtained by averaging the matrix element
over the Mandelstam variable $t$. We have

$\displaystyle\sigma_{e^{\pm}\gamma}=\frac{1}{16\pi(s-m_{e}^{2})^{2}}\int^{0}_{-(s-m_{e}^{2})^{2}/s}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dt\,|M_{e^{\pm}\gamma}|^{2},$
(12) $\displaystyle\sigma_{e^{-}e^{-}}=\frac{1}{16\pi
s(s-4m_{e}^{2})}\int^{0}_{-(s-4m_{e}^{2})}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dt\,|M_{e^{\pm}e^{\pm}}|^{2},$
(13) $\displaystyle\sigma_{e^{-}e^{+}}=\frac{1}{16\pi
s(s-4m_{e}^{2})}\int^{0}_{-(s-4m_{e}^{2})}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dt\,|M_{e^{\pm}e^{\mp}}|^{2}.$
(14)

The matrix element associated with inverse Compton scattering is given by [31,
32]

$\displaystyle|M_{e^{\pm}\gamma}|^{2}\\!$
$\displaystyle=32\pi^{2}\alpha^{2}\bigg{[}4\left(\frac{m_{e}^{2}}{m_{e}^{2}-s}+\frac{m_{e}^{2}}{m_{e}^{2}-u}\right)^{2}$
$\displaystyle-\frac{4m_{e}^{2}}{m_{e}^{2}-s}-\frac{4m_{e}^{2}}{m_{e}^{2}-u}-\frac{m_{e}^{2}-u}{m_{e}^{2}-s}-\frac{m_{e}^{2}-s}{m_{e}^{2}-u}\bigg{]},$
(15)

and for Møller and Bhabha scattering we have respectively

$\displaystyle|M_{e^{-}e^{-}}|^{2}\\!=64\pi^{2}\alpha^{2}\bigg{[}\frac{s^{2}+u^{2}+8m_{e}^{2}(t-m_{e}^{2})}{2(t-m^{2}_{\gamma})^{2}}$
$\displaystyle+\frac{{s^{2}+t^{2}}+8m_{e}^{2}(u-m_{e}^{2})}{2(u-m_{\gamma}^{2})^{2}}+\frac{\left({s}-2m_{e}^{2}\right)\left({s}-6m_{e}^{2}\right)}{(t-m_{\gamma}^{2})(u-m_{\gamma}^{2})}\bigg{]},$
(16)

and

$\displaystyle|M_{e^{-}e^{+}}|^{2}=64\pi^{2}\alpha^{2}\bigg{[}\frac{s^{2}+u^{2}+8m_{e}^{2}(t-m_{e}^{2})}{2(t-m^{2}_{\gamma})^{2}}$
$\displaystyle+\frac{u^{2}+t^{2}+8m_{e}^{2}(s-m_{e}^{2})}{2(s-m^{2}_{\gamma})^{2}}+\frac{\left({u}-2m_{e}^{2}\right)\left({u}-6m_{e}^{2}\right)}{(t-m^{2}_{\gamma})(s-m^{2}_{\gamma})}\bigg{]},$
(17)

where $s,t,u$ are the Mandelstam variables and we included the photon mass
$m_{\gamma}$ to avoid singularity in reaction matrix elements. The photon mass
in plasma is equal to the plasma frequency, in our case we have [33]

$m^{2}_{\gamma}=\omega^{2}_{p}=8\pi\alpha\int\frac{d^{3}p_{e}}{(2\pi)^{3}}\left(1-\frac{p_{e}^{2}}{3E_{e}^{2}}\right)\frac{f_{-}+f_{+}}{E_{e}},$
(18)

where $f_{-}$ and $f_{+}$ are the single-particle distribution functions for
electrons and positrons respectively, and $E_{e}=\sqrt{p_{e}^{2}+m^{2}_{e}}$
is the electron energy. In the BBN temperature range $m_{e}\gg T$ and if we
consider the non-relativistic electron and positron plasma

$\displaystyle
m^{2}_{\gamma}=\frac{4\pi\alpha}{2m_{e}}\left(\frac{2m_{e}T}{\pi}\right)^{3/2}e^{-m_{e}/T}\cosh\left(\frac{\mu_{e}}{T}\right).$
(19)

Substituting the cross-sections into Eq. (10), the thermal reaction rate per
volume for inverse Compton scattering can be written as

$\displaystyle
R_{e^{\pm}\gamma}=\frac{g_{e}g_{\gamma}}{16\left(2\pi\right)^{5}}T\int_{m_{e}^{2}}^{\infty}\\!\\!\\!\\!ds\frac{K_{1}(\sqrt{s}/T)}{\sqrt{s}}\int^{0}_{-(s-m_{e}^{2})^{2}/s}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dt\,|M_{e^{\pm}\gamma}|^{2},$
(20)

for Møller and Bhabha scattering we have

$\displaystyle
R_{e^{\pm}e^{\pm}}=\frac{g_{e}g_{e}}{16\left(2\pi\right)^{5}}T\\!\\!\int_{4m_{e}^{2}}^{\infty}\\!\\!\\!\\!ds\frac{K_{1}(\sqrt{s}/T)}{\sqrt{s}}\int^{0}_{-(s-4m_{e}^{2})}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dt\,|M_{e^{\pm}e^{\pm}}|^{2},$
(21) $\displaystyle
R_{e^{\pm}e^{\mp}}=\frac{g_{e}g_{e}}{16\left(2\pi\right)^{5}}T\\!\\!\int_{4m_{e}^{2}}^{\infty}\\!\\!\\!\\!ds\frac{K_{1}(\sqrt{s}/T)}{\sqrt{s}}\int^{0}_{-(s-4m_{e}^{2})}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dt\,|M_{e^{\pm}e^{\mp}}|^{2}.$
(22)

In Fig. 2, we show the reaction rates for Møller, Bhabha, and inverse Compton
scattering as a function of temperature. For temperatures $T>12.0$ keV, the
dominant reactions in plasma are Møller and Bhabha scatterings between
electrons and positrons. Thus, we can neglect the inverse Compton scattering
in the BBN temperature range.

It is convenient to define the average relaxation rate for the electron-
positron in the plasma as follows:

$\displaystyle\kappa=\frac{R_{e^{\pm}e^{\pm}}+R_{e^{\pm}e^{\mp}}+R_{e^{\pm}\gamma}}{\sqrt{n_{e^{-}}n_{e^{+}}}}\approx\frac{R_{e^{\pm}e^{\pm}}+R_{e^{\pm}e^{\mp}}}{\sqrt{n_{e^{-}}n_{e^{+}}}},$
(23)

where we neglect the inverse Compton scattering during BBN as discussed. The
density function ${\sqrt{n_{e^{-}}n_{e^{+}}}}$ in the Boltzmann limit is given
by

$\displaystyle{\sqrt{n_{e^{-}}n_{e^{+}}}}=\frac{g_{e}}{2\pi^{3}}T^{3}\left(\frac{m_{e}}{T}\right)^{2}K_{2}(m_{e}/T).$
(24)

Figure 2: The relaxation rate $\kappa$ as a function of temperature in non-
relativistic electron-positron plasma. We present reaction rates for Møller
scattering $R_{e^{\pm}e^{\pm}}$ (blue dashed line), Bhabha scattering
$R_{e^{\pm}e^{\mp}}$ (red dashed line), and inverse Compton scattering
$R_{e^{\pm}\gamma}$ (green dashed line). The average relaxation rate from Eq.
(23) is shown as a black solid line. The vertical black dotted lines represent
BBN temperature range $86\,\mathrm{keV}>\mathrm{T_{BBN}}>50\,\mathrm{keV}$
during which the average relaxation rate is $\kappa=10\sim 12$ keV. The
dominant reactions during BBN are the Møller and Bhabha scatterings. The
purple solid line represents the Debye mass given by Eq. (37)

.

In Fig. 2 that the total relaxation rate $\kappa$ given by Eq. (23) is
approximately constant $\kappa=10-12$ keV during the BBN. This value is much
larger than the Debye mass [see definition below in Eq. (37)]; we expect the
actual value of $\kappa$ to be smaller for a self-consistent calculation in
which damping is included in the reaction rates. For $T<20.3$ keV, the
relaxation rate $\kappa$ decreases rapidly because the positrons disappear. At
$T=12$ keV, the inverse Compton scattering of remaining electrons becomes the
dominant scattering process.

##  4 Early Universe plasma from kinetic theory

We evaluate the properties of the BBN plasma by finding the polarization
tensor using the linear response relation [34]

$\widetilde{j}_{\mathrm{ind}}^{\mu}(k)={\Pi^{\mu}}_{\nu}(k)\widetilde{A}^{\nu}(k)\,.$
(25)

here ‘$k$’ is the conjugate variable to the configuration space ‘$x$’ seen
below. The polarization tensor ${\Pi^{\mu}}_{\nu}(k)$ can be derived from the
induced current $\widetilde{j}_{\mathrm{ind}}^{\mu}(k)$ due to perturbations
$\delta f$ away from equilibrium in the relativistic Vlasov-Boltzmann
transport equations

$\displaystyle(p\cdot\partial)f_{\pm}(x,p)+$ $\displaystyle
qF^{\mu\nu}p_{\nu}\frac{\partial f_{\pm}(x,p)}{\partial
p^{\mu}}=C_{\pm}(x,p)\,,$ (26) $\displaystyle(p\cdot\partial)f_{\gamma}(x,p)$
$\displaystyle=C_{\gamma}(x,p)\,.$ (27)

The subscripts $-$, $+$, and $\gamma$ indicate the transport equation for
electrons, positrons, and photons. These form a system of differential
equations for each distribution function $f(x,p)$. We take the collision term
$C(x,p)$ in the BGK form as given by [21, 22]

$C_{i}(x,p)=\kappa_{i}(p\cdot
u)\left(f_{i}^{\text{eq}}(p)\frac{n_{i}(x)}{{n_{i}^{\text{eq}}}}-f_{i}(x,p)\right)\,,$
(28)

where $u$ is the 4-velocity of the plasma, $\kappa$ is the damping
coefficient, and $n$ is plasma particle the density. The subscripts ”i”
indicate the distributions of fermions and photons, respectively. We
suppressed the 4-momentum subscript for each species
$f_{i}(x,p)=f_{i}(x,p_{i})$ to simplify notation. Plasma constituents collide
on a momentum-averaged time scale $\tau_{\text{rel}}=\kappa^{-1}$, or the
”mean free time” in the plasma. The effect of the BGK collision term Eq. (28)
is to model the sum of all scattering effects on a species as a dissipative
medium effect which returns the system to equilibrium. The BGK collision term
is constructed such that Eq. (26) retains current conservation [21].

Since photons cannot couple directly to the electromagnetic field, they do not
contribute to the polarization tensor at first order in $\delta f$ as
indicated in Eq. (27). We neglect photon influence on the electron-positron
distribution through the scattering term since the rate of inverse Compton
scattering $R_{e^{\pm}\gamma}$ shown in green in Fig. (2) is much smaller than
the total rate $\kappa$ shown black in the BBN temperature range and in
general when noticeable positron abundance is present. Each fermion Boltzmann
equation Eq. (26) can be solved independently. Since the equations for
electrons and positrons are equivalent, except for the sign of the charge,
only one needs to be solved to understand the dynamics.

We do not consider the influence of light nuclei on the polarization tensor
since their density during the BBN epoch is much smaller than that of
electrons and positrons Fig. (1). One can see in Fig. (3) that the separation
of baryons $n_{B}^{-1/3}$ in black is much larger than the size of the
polarizing Debye sphere, so baryons do not participate significantly in
screening.

We take the equilibrium one particle distribution function $f_{\pm}^{(eq)}$ of
electrons and positrons to be the relativistic Fermi-distribution

$f^{(eq)}_{\pm}(p)=\frac{1}{\exp{\left(\frac{\sqrt{\boldsymbol{p}^{2}+m^{2}}}{T}\right)}+1}\,,$
(29)

with chemical potential $\mu=0$. The electron and positron mass will be
indicated by $m$ unless otherwise stated. At temperatures interesting for
nucleosynthesis $T=50-86$ keV, we expect the plasma temperature to be much
less than the mass of the plasma constituents. Only the non-relativistic form
of Eq. (29) will be relevant at these temperature scales

$f^{(eq)}_{\pm}(p)\approx\exp\left(-\frac{m}{T}\left(1+\frac{|\boldsymbol{p}|^{2}}{2m^{2}}\right)\right)\,.$
(30)

Keeping terms up to quadratic order in $|\boldsymbol{p}|/m$ we solve the
Vlasov-Boltzmann equation [Eq. (26)] for the induced current and identify the
polarization tensor. This is done in detail in our previous work in [22], but
we will reiterate a few key steps.

First, we expand Eq. (26) around small perturbations from equilibrium

$f_{\pm}(x,p)={f^{(eq)}_{\pm}}(p)+\delta f_{\pm}(x,p)\,,$ (31)

and solve Eq. (26) for $\delta f_{\pm}(x,p)$ in Fourier space. The induced
current in Fourier space is given by

$\tilde{j}_{\mathrm{ind}}^{\mu}(k)=2\int\frac{d^{4}p}{(2\pi)^{4}}p^{\mu}4\pi\delta_{+}(p^{2}-m^{2})\sum_{i=\,+,\,-}q_{i}\tilde{f}_{i}(k,p)\,,$
(32)

with the factor of two accounting for spin. After inserting Eq. (31), and
specifying $q_{\pm}=\pm e$ the induced current is a function of the
perturbation

$\tilde{j}_{\mathrm{ind}}^{\mu}(k)=2\int\frac{d^{3}p}{(2\pi)^{3}p^{0}}p^{\mu}\Big{(}e\left[\tilde{f}^{(eq)}_{+}(k,p)-\tilde{f}^{(eq)}_{-}(k,p)\right]\\\
+e\left[\delta\tilde{f}_{+}(k,p)-\delta\tilde{f}_{-}(k,p)\right]\Big{)}\\\
=4e\int\frac{d^{3}p}{(2\pi)^{3}p^{0}}p^{\mu}\delta\tilde{f}(k,p)\,,$ (33)

because the equilibrium currents cancel in the weak field limit and the
perturbations add since they differ by the charge $\delta f_{\pm}=\pm e\delta
f^{\prime}$. This term is studied for finite chemical potential in [1]. We
focus on the second term related to the polarization response of the plasma.

The polarization tensor can be obtained from Eq. (33) by computing the
4-momentum integrals over $p$ in the rest frame of the plasma. Once
integrated, the 4-potential $\widetilde{A}^{\nu}(k)$, coming from the
electromagnetic field strength tensor $F^{\mu\nu}$ in Eq. (26), factors out
since it only depends on the 4-wavevector $k$ and Eq. (33) attains the form of
Eq. (25). For details, see Ref. [22].

In the infinite homogeneous plasma filling the early Universe, the
polarization tensor only has two independent components: the longitudinal
polarization function $\Pi_{\parallel}$ parallel to field wave-vector
$\boldsymbol{k}$ in the rest frame of the plasma and the transverse
polarization function $\Pi_{\perp}$ perpendicular to $\boldsymbol{k}$ [35]. In
the non-relativistic limit, these functions are [22]

$\displaystyle\Pi_{\parallel}(\omega,\boldsymbol{k})$
$\displaystyle=-\omega_{p}^{2}\frac{\omega^{2}}{(\omega+i\kappa)^{2}}\frac{1}{1-\frac{i\kappa}{\omega+i\kappa}\left(1+\frac{T|\boldsymbol{k}|^{2}}{m(\omega+i\kappa)^{2}}\right)}\,,$
(34) $\displaystyle\Pi_{\perp}(\omega)$
$\displaystyle=-\omega_{p}^{2}\frac{\omega}{\omega+i\kappa}\,.$ (35)

In these expressions, the plasma frequency $\omega_{p}$ (defined as $m_{L}$ in
Ref. [22]) is related to the Debye screening mass in the non-relativistic
limit as

$\omega_{p}^{2}=m_{D}^{2}\frac{T}{m}\,.$ (36)

Figure 3: The average distance between baryons $n_{B}^{-1/3}$ and the Debye
length $\lambda_{D}$ ($\mu_{e}\neq 0$) as a function of temperature (red solid
line). During the BBN epoch (vertical dotted lines)
$n_{B}^{-1/3}>\lambda_{D}$. For temperature below $T<32.76$ keV we have
$n_{B}^{-1/3}<\lambda_{D}$. For comparison, the Debye length for zero chemical
potential $\mu_{e}=0$ is also plotted as a blue dashed line.

The transverse response $\Pi_{\perp}$ relates to the dispersion of photons in
the plasma. Here we need only consider $\Pi_{\parallel}$ since the vector
potential $\boldsymbol{A}(t,\boldsymbol{x})$ of the traveling ion will be
small in the non-relativistic limit. In this work, we do not consider the
effect of a primordial magnetic field discussed in Ref. [36]. We note that
Debye mass $m_{D}$ is related to the usual Debye screening length of the field
in the plasma as

$1/\lambda_{D}^{2}=m_{D}^{2}=4\pi\alpha\left(\frac{2mT}{\pi}\right)^{3/2}\frac{e^{-m/T}}{2T}\,.$
(37)

This formula describes the characteristic length scale of screening in the
plasma.

##  5 Effective internuclear potential

We calculate the potential of light nuclei in the early Universe electron-
positron plasma by finding the screened scalar potential $\phi$ of a single
traveling nuclei

$\phi(t,\boldsymbol{x})=\int\frac{d^{4}k}{(2\pi)^{4}}e^{-i\omega
t+i\boldsymbol{k}\cdot\boldsymbol{x}}\frac{\widetilde{\rho}_{\text{ext}}(\omega,\boldsymbol{k})}{\varepsilon_{\parallel}(\omega,\boldsymbol{k})(\boldsymbol{k}^{2}-\omega^{2})}\,,$
(38)

where $\widetilde{\rho}_{\text{ext}}(\omega,\boldsymbol{k})$ is the Fourier-
transformed charge distribution of nuclei traveling at a constant velocity, as
found in A, and $\varepsilon_{\parallel}(\omega,\boldsymbol{k})$ is the
longitudinal relative permittivity.111We omit $\varepsilon_{0}$ until the
final result, to avoid confusion with
$\varepsilon_{\parallel}(\omega,\boldsymbol{k})$. The screened potential is
obtained by solving Maxwell equations algebraically in Fourier space in the
presence of a polarizing medium. The relative permittivity can be written in
terms of the polarization tensor as

$\varepsilon_{\parallel}(\omega,\boldsymbol{k})=\left(\frac{\Pi_{\parallel}(\omega,\boldsymbol{k})}{\omega^{2}}+1\right)\,.$
(39)

The potential between two nuclei can be inferred simply from the potential of
a single nucleus, since in the linear response framework Eq. (25) the
electromagnetic field still obeys the principle of superposition.

We can perform the $\omega$ integration in Eq. (38) using the delta function
in the definition of the external charge distribution Eq. (65), whose effect
is to set $\omega=\boldsymbol{\beta_{\text{N}}}\cdot\boldsymbol{k}$ where
$\boldsymbol{\beta}_{N}=\boldsymbol{v}_{N}/c$ is the nuclei velocity. Then we
have

$\phi(t,\boldsymbol{x})=Ze\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}e^{i\boldsymbol{k}\cdot(\boldsymbol{x}-\boldsymbol{\beta_{\text{N}}}t)}\frac{e^{-\boldsymbol{k}^{2}\frac{R^{2}}{4}}}{\boldsymbol{k}^{2}\varepsilon_{\parallel}(-\boldsymbol{\beta_{\text{N}}}\cdot\boldsymbol{k},\boldsymbol{k})}\,,$
(40)

where $R$ is the Gaussian radius parameter. In non-relativistic approximation
the Lorentz factor $\gamma\approx 1$ and we use the convention
$\varepsilon_{\parallel}(-\boldsymbol{\beta_{\text{N}}}\cdot\boldsymbol{k},\boldsymbol{k})$
[23, 24, 25, 37] which gives the correct causality for the potential.
Specifically, without damping, the wakefield occurs behind the moving nucleus.

In the limit of a single nucleus at rest in the plasma, the screened potential
has the usual Debye-Hückel form [5]

$\phi_{\text{stat}}(t,\boldsymbol{x})=\frac{Ze}{4\pi\varepsilon_{0}}\frac{e^{-m_{D}r}}{r}\,.$
(41)

The analogous formula for the static screening of a Gaussian charge is derived
in B.

###  5.1 Damped-dynamic screening

In this section, we compute the effective nuclear potential for a light
nucleus moving in the plasma at a constant velocity. This is done by Fourier
transforming Eq. (40). The velocity of the nucleus is assumed to be the most
probable velocity given by a Boltzmann distribution

$\beta_{\text{N}}=\sqrt{\frac{2T}{m_{N}}}\,.$ (42)

Since the poles of the Eq. (38) can be solved analytically, ideally, one would
perform contour integration to get the position space field. Due to the
intricacy of these poles $\omega_{n}(\boldsymbol{k})$, we find it insightful
to look at the field in a series expansion around velocities of the light
nuclei smaller than the thermal velocity of electrons and positrons, and large
damping.

$\frac{(\boldsymbol{k}\cdot\boldsymbol{\beta}_{\text{N}})^{2}}{\omega_{p}^{2}}\ll\frac{\boldsymbol{k}^{2}}{m_{D}^{2}}\ll\frac{\kappa^{2}}{\omega_{p}^{2}}\,.$
(43)

This expansion is useful during BBN since the temperature is much lower than
the mass of light nuclei and the damping rate $\kappa$ is approximately twice
the Debye mass $m_{D}$, as seen in Fig. 2. Applying this expansion to Eq. (40)
we obtain at first order the potential

$\phi(t,\boldsymbol{x})=\phi_{\text{stat}}(t,\boldsymbol{x})\\\
-Ze\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}e^{i\boldsymbol{k}\cdot(\boldsymbol{x}-\boldsymbol{\beta_{\text{N}}}t)}\frac{i\boldsymbol{k}\cdot\boldsymbol{\beta_{\text{N}}}m_{D}^{4}(\frac{\boldsymbol{k}^{2}}{m_{D}^{2}}-\frac{\kappa^{2}}{\omega_{p}^{2}})}{\boldsymbol{k}^{2}(\boldsymbol{k}^{2}+m_{D}^{2})^{2}\kappa}e^{-\boldsymbol{k}^{2}\frac{R^{2}}{4}}\,.$
(44)

First, we focus on the second term evaluating for a point charge $R\rightarrow
0$

$\Delta\phi(t,\boldsymbol{x})=-Ze\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}e^{i\boldsymbol{k}\cdot(\boldsymbol{x}-\boldsymbol{\beta_{\text{N}}}t)}\frac{i\boldsymbol{k}\cdot\boldsymbol{\beta_{\text{N}}}m_{D}^{4}(\frac{\boldsymbol{k}^{2}}{m_{D}^{2}}-\frac{\kappa^{2}}{\omega_{p}^{2}})}{\boldsymbol{k}^{2}(\boldsymbol{k}^{2}+m_{D}^{2})^{2}\kappa}\,.$
(45)

Numerical calculations are performed using the full Gaussian charge
distribution. We focus on the second term in Eq. (44) since
$\phi_{\text{stat}}$ is the standard static screening result. In order to
perform this integration we first note that this expression can be re-written
as the laboratory time derivative

$\Delta\phi(t,\boldsymbol{x})=Ze\frac{d}{dt}\left[\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}e^{i|\boldsymbol{k}||\boldsymbol{x}-\boldsymbol{\beta_{\text{N}}}t|\cos(\theta)}\frac{m_{D}^{4}(\frac{\boldsymbol{k}^{2}}{m_{D}^{2}}-\frac{\kappa^{2}}{\omega_{p}^{2}})}{\boldsymbol{k}^{2}(\boldsymbol{k}^{2}+m_{D}^{2})^{2}\kappa}\right]\,,$
(46)

where the dot products were replaced by their angular dependence. The angular
integration can be evaluated in spherical coordinates, see B for details,

$\Delta\phi(t,\boldsymbol{x})=2Ze\frac{d}{dt}\left[\int\frac{d\boldsymbol{k}}{(2\pi)^{2}}\frac{\sin(|\boldsymbol{k}||\boldsymbol{x}-\boldsymbol{\beta}_{N}t|)}{|\boldsymbol{k}||\boldsymbol{x}-\boldsymbol{\beta}_{N}t|}\,\frac{m_{D}^{4}(\frac{\boldsymbol{k}^{2}}{m_{D}^{2}}-\frac{\kappa^{2}}{\omega_{p}^{2}})}{\boldsymbol{k}^{2}(\boldsymbol{k}^{2}+m_{D}^{2})^{2}\kappa}\right]\,,$
(47)

finally, this can be integrated over the radial wavevector to find

$\Delta\phi(t,\boldsymbol{x})=\frac{Ze}{4\pi}\frac{m_{D}^{2}}{\kappa}\frac{d}{dt}\Bigg{[}\left(\frac{1+\nu_{\tau}^{2}}{2m_{D}}+\frac{\nu_{\tau}^{2}}{m_{D}^{2}|\boldsymbol{x}-\boldsymbol{\beta}_{N}t|}\right)e^{-m_{D}|\boldsymbol{x}-\boldsymbol{\beta}_{N}t|}\\\
-\frac{\nu_{\tau}^{2}}{m_{D}^{2}|\boldsymbol{x}-\boldsymbol{\beta}_{N}t|}\Bigg{]}\,,$
(48)

here we introduce the ratio of the damping rate to the rate of oscillations in
the plasma $\nu_{\tau}=\kappa/\omega_{p}$. We can apply the time derivative;
note that

$\frac{d}{dt}|\boldsymbol{x}-\boldsymbol{\beta}_{N}t|=-\frac{\boldsymbol{\beta}_{N}\cdot(\boldsymbol{x}-\boldsymbol{\beta}_{N}t)}{|\boldsymbol{x}-\boldsymbol{\beta}_{N}t|}=-\beta_{N}\cos(\psi)\,,$
(49)

where $\psi$ is the angle between $\boldsymbol{x}-\boldsymbol{\beta}_{N}t$ and
$\boldsymbol{\beta}_{N}$. After taking the derivative and using the above
expression Eq. (49) one finds

$\Delta\phi(t,\boldsymbol{x})=\frac{Ze\beta_{N}\cos(\psi)m_{D}^{2}}{4\pi\varepsilon_{0}\kappa}\Bigg{[}\left(\frac{\nu_{\tau}^{2}}{m_{D}^{2}r(t)^{2}}+\frac{\nu_{\tau}^{2}}{m_{D}r(t)}+\frac{1+\nu_{\tau}^{2}}{2}\right)e^{-m_{D}r(t)}\\\
-\frac{\nu_{\tau}^{2}}{m_{D}^{2}r(t)^{2}}\Bigg{]}\,,$ (50)

with $r(t)=|\boldsymbol{x}-\boldsymbol{\beta}_{N}t|$. This expression is valid
for large damping and slow motion of the nucleus or if the velocity of the
nuclei is small. A similar result valid at large distances, which only
includes the last term, was derived previously in Ref. [24] in the context of
dusty (complex) plasmas. The static screening potential $\phi_{\text{stat}}$
must be added to this contribution to get the full potential

$\phi(t,\boldsymbol{x})=\phi_{\text{stat}}(t,\boldsymbol{x})+\Delta\phi(t,\boldsymbol{x})\,.$
(51)

For large distances and large $\nu_{\tau}$, the last term in the second line
is dominant, indicating that the overall potential would be over-damped. In
this regime, the potential is heavily screened in the forward direction and
unscreened in the backward direction relative to the motion of the nucleus. As
$\nu_{\tau}$ becomes small, the $1/2$ in the first portion of the third term,
proportional to $m_{D}^{2}/\kappa$, dominates. This flips the sign of the
damped-dynamic screening contribution causing a wake potential to form behind
the nuclei. This shift indicates the change from damped to undamped screening
where Eq. (50) is no longer valid.

We note in passing that using a Schwinger parameterization similar to the one
used to derive Eq. (68)

$\frac{1}{(k^{2}+m_{D}^{2})^{2}}=\int_{0}^{\infty}ds\,s\,e^{-s(k^{2}+m_{D}^{2})}\,,$
(52)

it is possible to derive an expression for damped-dynamic screening for a
Gaussian charge distribution, which is a trivial extension of Eq. (50).

Figure 4: A comparison of the following screening models plotted along the
direction of motion of a nucleus
$\boldsymbol{r}\cdot\hat{\boldsymbol{\beta}_{\text{N}}}$: static screening
(black), dynamic screening (red dotted)[2], damped-dynamic screening (blue
dashed), and the approximate analytic solution of Eq. (51) (orange dashed). A
black arrow indicates the direction of motion of the nucleus
$\hat{\boldsymbol{\beta}_{\text{N}}}$.

Figure 4 demonstrates that the damped-dynamic response in the analytic
approximation Eq. (50) (shown as orange dashed line) is sufficient to
approximate the full numerical solution (blue dashed line) found by numerical
integration of Eq. (38). The temperature $T=100$ keV, above our upper limit of
BBN temperatures, is chosen to relate to the dynamic screening result found in
Ref. [2]. Our analytic solution differs from the numerical result in Fig. 4 of
Ref. [2] by a factor of $\sqrt{2}$ and is horizontally flipped. This
reflection is due to a difference in convention in the permittivity, as seen
in Eq. (40). We can see that at large distances, dynamic screening is slightly
stronger than damped screening, as expected. At short distances relevant for
thermonuclear reaction rates, damped and undamped screening are very similar.

Dynamic screening in both the damped and undamped cases predicts less
screening behind and more in front of the moving nucleus, as compared to
static screening. This effect was previously observed for subsonic screening
in electron-ion-dust plasmas [24, 25, 26]. As a result, a negative
polarization charge builds up in front of the nucleus. The small negative
potential in front alters the potential energy between light nuclei, possibly
changing the equilibrium distribution of light elements. This effect is much
larger in the undamped case and is known in some cases to lead to the
formation of dust crystals [37]. We discuss possible implications for the
early universe structure formation in the conclusions.

Figure 5: The analytic result from Eq. (50) is compared to the full numerical
(dashed lines) result plotted along the direction of motion of a nucleus in
the temperature range, $20-90\,$keV. At low temperatures, the full numerical
solution returns to static screening (gray solid line). At high temperatures
near the beginning of BBN, the Eq. (50) approaches the full numerical
solution.

In Fig. (5), we present a comparison of the approximate analytic result of Eq.
(50), shown as solid lines ranging from blue to red, to the full numerical
solution of Eq. (38), shown as dashed lines, for damped-dynamic screening.
Interestingly, as the temperature decreases, the screening effect increases in
contrast to the usual Debye-Hückel screening behavior. This is because in the
large damping limit, we have terms in Eq. (50) proportional to
$\kappa/m_{D}^{2}$. Since the debye mass decreases with temperature but
$\kappa$ remains relatively constant during BBN Fig. (2), one finds an overall
increase in long-distance screening. The expansion used to derive the analytic
solution Eq. (43) assumes large temperatures and breaks down at temperatures
below $60\,$keV. At low temperatures, the numerical result matches the static
screening result; at $T=20\,$kev, the numerical solution indicated by a dashed
blue line traces static screening at the same temperature, shown in solid
gray. The analytic solution trends to the numerical result at higher
temperatures, within $\sim 10\%$ for $T=70\,$keV. At temperatures above
$200\,$keV, the non-relativistic approximation should be re-evaluated.

###  5.2 Limitations of linear response: toy model

We now consider the equilibrium phase space distribution in the presence of a
‘strong’ potential. In covariant form and not neglecting the magnitude of the
4-potential $A^{\mu}$ as compared to the energy-momentum $p^{\mu}$ of the
particle in statistical ensemble, one obtains [38, 39]

$f^{(eq)}(x,p)=e^{-u_{\mu}(p^{\mu}+qA_{(\text{eq})}^{\mu}(x))/T}\,,$ (53)

where $u^{\mu}=(1,0,0,0)$ is the global velocity of the plasma in its rest
frame. The density according to this distribution is given by the usual
expression but altered by the exponential factor in potential

$n_{eq}=\int\frac{dp^{3}}{(2\pi)^{3}p^{0}}p^{0}f^{(eq)}(x,p)=n_{eq}e^{-u^{\mu}qA_{(\text{eq})}^{\mu}(x)/T}\,.$
(54)

Inserting the equilibrium density into the Poisson equation

$-\nabla^{2}\phi_{(\text{eq})}(x)=\rho_{\mathrm{tot}}(x)=\rho_{\mathrm{ext}}(x)+\rho_{\mathrm{ind}}(x)\,,$
(55)

gives the nonlinear Poisson-Boltzmann equation in equilibrium

$-\nabla^{2}\phi_{(\text{eq})}(x)+4en_{(\text{eq})}\sinh\left[e\phi_{(\text{eq})}(x)/T\right]=\rho_{\mathrm{ext}}(x)\,.$
(56)

We can re-scale the potential with temperature to rewrite as

$-\nabla^{2}\Phi_{(\text{eq})}(x)+m_{D}^{2}\sinh\left[\Phi_{(\text{eq})}(x)\right]=e\rho_{\mathrm{ext}}(x)/T\,,$
(57)

where $\Phi_{(\text{eq})}(x)=e\phi_{(\text{eq})}(x)/T$. This equation has a
well-known solution for an infinite sheet, where the external charge
$\rho_{\mathrm{ext}}$ becomes a fixed boundary condition. We will set the
potential at the sheet $\phi(0)=\phi_{0}$ and solve the Poisson-Boltzmann
equation in the region to the right of the sheet ($x>0$)

$-\frac{d^{2}}{dx^{2}}\Phi_{(\text{eq})}(x)+m_{D}^{2}\sinh\left[\Phi_{(\text{eq})}(x)\right]=0\,.$
(58)

The solution to this problem can be found for example in Ref. [40]

$\Phi_{(\text{eq})}(x)=\frac{e\phi_{(\text{eq})}(x)}{T}=4\tanh^{-1}\left[\tanh\left(\frac{e\phi_{0}}{4T}\right)e^{-m_{D}x}\right]\,,$
(59)

which we can compare to the usual Debye screening for a charged plane

$\phi_{(\text{eq})}(x)=\phi_{0}e^{-m_{D}x}\,.$ (60)

Figure 6: Comparison of nonlinear screening to Debye linearized screening
(solid top line) in the range $50\leq T\leq 400$ keV; the BBN temperature
range is shaded. In this toy model example $Z=2$ was used in Eq. (61).

We contrast linear screening Eq. (60) and nonlinear screening Eq. (59) in Fig.
(6). The boundary value of the potential $\phi_{0}$ is chosen to be the
maximum value of the potential of a Gaussian charge distribution

$\phi_{0}=\frac{Ze}{8\pi^{3/2}\varepsilon_{0}R}\,.$ (61)

Figure 6 shows that nonlinear screening and linear screening predict different
decay rates for the potential $\phi$ at short distances but have the same
long-range decay rate. Since the screening near the origin is larger in the
nonlinear case, the long-distance potential $\phi$ does not approach the weak
field limit but is offset. As Fig. (6) demonstrates this is especially true in
the BBN temperature range. All solutions have an exponential slope $-m_{D}$
indicating linear or weak field behavior at large distances.

The important prediction of the nonlinear screening in Fig. (6) is a rapid
drop of the potential at short distances, thus a strong screening effect near
the classical turning point of thermal reacting ions, where quantum tunneling
for thermonuclear reactions will be important. This effect increases
nonlinearly with $Z$, see Eq. (59). Note that the use of relativistic Fermi
statistics may be required since the tunneling regime involves potentials that
are comparable to and larger than the electron mass. We expect that nonlinear
screening effects could lead to drastic changes in the Coulomb penetration
factor for $Z>1$ and, in turn, light element abundances.

To describe nonlinear screening beyond the presented toy model, in the context
of light element screening during BBN, we need to solve Eq. (57) with a
realistic spherical charge distribution. That way the difference in asymptotic
behavior will be captured for the solution at large distances. It turns out
that even the simple first look at the 3D-static asymptotic Debye-type
solution requires formidable numerical effort, let alone merging the short-
range strong field limit with damped-dynamic linear response theory as
developed here. Therefore the strong field screening will be addressed
separately, opening the path to the computation of screened-BBN reaction
rates.

##  6 Conclusion and discussion

This study sheds light on the properties and dynamics of the electron-positron
plasma in the early Universe. In particular, we examined the damping rate and
the electron-positron to baryon density ratio and their potential implications
for Big Bang Nucleosynthesis (BBN) through screening within linear response
theory.

Our results confirm that BBN occurred in a dense $e^{-}e^{+}\gamma$ QED plasma
interspersed with the heavy charged dust of protons and light nuclei. The
polarization cloud surrounding an ionized nucleus reduces the Coulomb barrier
needed to be overcome for a thermonuclear reaction to proceed. This increase
in the penetration probability enhances the thermonuclear reaction rates.

We calculated the screening effect on the internuclear Coulomb potential in
the electron-positron plasma, accounting for the thermal motion of the nuclei
and collisional damping. Our use of the current conserving collisional linear
response framework for relativistic plasmas detailed in [22] allowed us to
present an approximate analytic formula in Eq. (50) to describe screening
effects which were previously found only numerically. Our analytic formula can
be readily used to estimate the effect of screening on thermonuclear
reactions. Our results are illustrated in Fig. 5.

Common theoretical BBN models assume that the Universe was essentially void of
anything except for reacting light nucleons and electrons needed to keep the
local baryon density charge-neutral, a situation similar to the laboratory
environment in which empirical nuclear reaction rates are obtained. Newer
models, along with ours, address the BBN in dense electron-positron plasma
using linear response approach [1, 2]. Since our results for damped-dynamic
screening are very similar to dynamic screening found in [2] we expect that
our work will not change existing predictions [2], who obtained numerically
the modification of Coulomb potential in linear response approximation but
without damping. Here, we were able to obtain these results analytically,
which helped us realize the limitations of the linear response approximation
at short distances relevant for thermonuclear reactions.

The prior and present study of BBN within the linear response approach has not
fully recognized these limitations before: We believe that to model the BBN
reaction network processes in dense primordial plasma, we need to obtain the
strong screening effect when the potential strength is close to the thermal
energy $|\phi|/T\approx 1$. This requires developing novel theoretical tools
beyond the scope of this work. We showed, see Fig. (6), that for the simple
one-dimensional toy model, strong field screening could alter particle nuclear
reaction rates, with importance increasing with $Z$. Since this will impact
the BBN element production we will return to strong screening at short
internuclear distances in follow-up work. Consideration of thermonuclear
reaction rate enhancement estimation will follow the study of the short
distance screening potential in the strong field limit.

We currently calculate the damping strength to be much larger than the
screening mass. This implies another necessary improvement in evaluating the
strength of damping, which itself depends on plasma properties through the
screening of scattering potential acting between charged plasma components. In
forthcoming work, we will propose an approach where damping and plasma
properties are determined in a self-consistent manner. We expect a self-
consistent calculation to yield a smaller value as the damping effect reduces
the magnitude of in plasma scattering cross-sections. The behavior of the
plasma can alter significantly if the damping parameter becomes small enough
that plasma oscillations become undamped.

Our study mirrors the theoretical findings of planetary and space dusty plasma
theory [23, 24, 25, 26]. The large distance behavior of the damped-dynamic BBN
screened potential in Eq. (50) matches that of slowly moving dust particles in
plasma [24]. Dusty plasma theory studies a number of effects that are not
currently included in BBN plasma studies, including dust charging, dust
acoustic waves, dust instabilities, and structure formation (dust crystals)
[25]. We expect that these results can be ported to the nuclear light element
dust dynamics in the primordial $e^{-}e^{+}\gamma$ QED plasma. This
interdisciplinary connection, which has not been recognized previously, could
have substantial implications for our understanding of the evolution of matter
in the early Universe, implying in particular primordial matter clumping.

It is compelling to consider that this polarization binding effect could lead
to the early seeding of matter structure in the Universe. This happens due to
the existence of an attractive portion of the internuclear (charged dust)
potential which favors a specific average separation between charged dust
particles. Therefore as the density of the dust in the Universe decreases due
to Universe expansion, density inhomogeneity formation can be expected: The
energetically favorable configuration consists of domains with ‘optimal’
density of nuclear dust and other domains with reduced density or even
entirely devoid of dust. Such inhomogeneities could influence the production
of heavier isotopes and elements in the late stage of the BBN.

Dust crystals are known to form in polarized plasma [37, 41]. In our case
nuclear ‘dust’ clustering driven by in plasma polarization forces may provide
a natural explanation for the recently discovered primordial galaxy
accompanied by a large central black hole [18] existent when the Universe was
barely 500 million years old. This is in our opinion suggestive of primordial
matter concentration rather than evidence for other explanations, such as dark
matter stars [20] or a much longer Universe lifespan [42] which could impact
other predictions of the $\Lambda$CDM model [19]. We believe that developing
novel BBN plasma tools necessary for a screened nuclear reaction network, and
a study of associated matter clumping in plasma driven by polarization forces
could contribute to an improved understanding of JWST-HST spectacular
observations.

Many of our results can be ported to other types of plasma than considered
here. An immediate interdisciplinary application for this work could be the
environment of stellar and ‘laboratory’ plasma fusion with significant effects
expected for the proton-Boron environment of intense current interest [43,
44]. Another vibrant field of current research in acute need of a better
understanding is the production of light nuclei in relativistic heavy ion
collisions. These are emerging from high-temperature hadronic fireballs
comprising predominantly of light pions and a few heavy ‘dust’ nucleons at an
ambient temperature vastly exceeding the nuclear binding. Their relatively
abundant formation lacks a convincing explanation [45]. In the same
experimental relativistic heavy ion collision environment, the new strongly
interacting quark-gluon plasma is formed initially. The gluon and light quark
$u,d,s$ pair-plasma is seeded with heavy $c,\bar{c},b,\bar{b}$ ‘charged dust’
[46]. This situation is also reminiscent of the environment explored here.
However, additional complexity is introduced by the non-Abelian nature of
strong interaction.

In summary, our study has provided advances in understanding electron-positron
plasma in the early Universe and its potential influence on the BBN. We
discussed relation to other domains of physics where similar methods were used
and presented several interdisciplinary directions for future research.
Additionally, we have clarified steps going beyond the linear response method,
which are necessary to understand the implications of our findings for BBN
thermonuclear reaction rates and the evolution of the early Universe.

## Acknowledgements

This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors. However, we thank Tamás Biró
for his hospitality at the PP2023 Margaret Island Symposium where this work
was presented. An allocation of computer time was used from the High
Performance Computing (HPC) resources supported by the University of Arizona
TRIF, UITS, and Research, Innovation, and Impact (RII) and maintained by the
UArizona Research Technologies department.

##  Appendix A External charge distribution

Here we define the free charge density used to describe light nuclei in BBN.
We wish to model a nucleus moving at constant velocity
$\boldsymbol{\beta}_{N}$. For simplicity, we model the charge distribution as
a Gaussian in all directions

$\rho_{\text{ext}}(t,\boldsymbol{x})=\frac{Ze\gamma}{\pi^{3/2}R^{3}}e^{-\frac{1}{R^{2}}(\boldsymbol{x}-\boldsymbol{\beta}_{N}t)^{2}}\,,$
(62)

Since we work in the non-relativistic limit we assume that the Lorentz factor
$\gamma\approx 1$. The normalization is chosen in such a way that

$\int\rho_{\text{ext}}(t,\boldsymbol{x})d^{3}\boldsymbol{x}=Ze\,$ (63)

is the total charge of the nucleus. The Gaussian radius parameter $R$ is
related to the mean squared radius of the nucleus, $\langle r^{2}\rangle$ at
rest, by

$\langle r^{2}\rangle=\frac{1}{Ze}\int
r^{2}\rho_{\text{ext}}(\boldsymbol{x})d^{3}\boldsymbol{x}=\frac{3}{2}R^{2}\,.$
(64)

Experimental measurements of the mean square radius can be found in Ref. [47].
The Fourier transformed charge distribution is

$\widetilde{\rho}_{\text{ext}}(\omega,\boldsymbol{k})=2\pi
Ze\,e^{-\boldsymbol{k}^{2}\frac{R^{2}}{4}}\delta(\omega-\boldsymbol{k}\cdot\boldsymbol{\beta}_{N})\,,$
(65)

where $\delta$ is the Dirac delta function in one dimension.

##  Appendix B Static screening for a Gaussian charge distribution

Starting from Eq. (40) we can take the limit $\beta_{N}\rightarrow 0$ to find
the static screening potential for a Gaussian charge distribution

$\phi_{\text{stat}}(t,\boldsymbol{x})=Ze\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}e^{i\boldsymbol{k}\cdot\boldsymbol{x}}\frac{e^{-\boldsymbol{k}^{2}\frac{R^{2}}{4}}}{\boldsymbol{k}^{2}+m_{D}^{2}}\,.$
(66)

First, we perform the angular integration to find

$\phi_{\text{stat}}(t,\boldsymbol{x})=Ze\int_{0}^{\infty}\frac{k^{2}d{k}}{(2\pi)^{2}}\frac{2\sin{({k}|\boldsymbol{x}|)}}{{k}|\boldsymbol{x}|}\frac{e^{-{k}^{2}\frac{R^{2}}{4}}}{{k}^{2}+m_{D}^{2}}\,.$
(67)

This integral can be evaluated using a Schwinger parameterization

$\phi_{\text{stat}}(t,\boldsymbol{x})=Ze\int_{0}^{\infty}ds\int_{0}^{\infty}\frac{k^{2}d{k}}{(2\pi)^{2}}\frac{2\sin{({k}|\boldsymbol{x}|)}}{{k}|\boldsymbol{x}|}e^{-{k}^{2}\frac{R^{2}}{4}-s({{k}^{2}+m_{D}^{2}})}\,,$
(68)

where we first perform the integration over $k$. The $k$ integration is
calculated by writing the $\sin$ in its exponential form and splitting it into
two terms with opposite signs in the exponents. These two terms can be
recombined into one Fourier transform from $-\infty$ to $+\infty$. Finally,
the integration over $s$ results in the position space static screening
potential

$\phi_{\text{stat}}(t,\boldsymbol{x})=\frac{Ze}{4\pi\varepsilon_{0}}\frac{e^{m_{D}^{2}R^{2}/4}}{r}\Bigg{(}\frac{\text{Erf}\left(\frac{r}{R}-\frac{m_{D}R}{2}\right)}{2}e^{-m_{D}r}\\\
+\frac{\text{Erf}\left(\frac{r}{R}+\frac{m_{D}R}{2}\right)}{2}e^{m_{D}r}-\sinh(m_{D}r)\Bigg{)}\,,$
(69)

where we have reintroduced the vacuum permittivity $\varepsilon_{0}$. Taking
the limit of this expression for $R\rightarrow 0$ yields the usual static
screening result [5]

$\lim_{R\rightarrow
0}\phi_{\text{stat}}(t,\boldsymbol{x})=\frac{Ze}{4\pi\varepsilon_{0}}\frac{e^{-m_{D}r}}{r}\,.$
(70)

## References

  * [1] B. Wang, C. A. Bertulani and A. B. Balantekin, “Electron screening and its effects on Big-Bang nucleosynthesis,”Phys. Rev. C 83, 018801 (2011) doi:10.1103/PhysRevC.83.018801 [arXiv:1010.1565 [astro-ph.CO]]
  * [2] E. Hwang, D. Jang, K. Park, M. Kusakabe, T. Kajino, A. B. Balantekin, T. Maruyama, C. M. Ryu and M. K. Cheoun, “Dynamical screening effects on Big-Bang nucleosynthesis,”JCAP 11 (2021), 017 doi:10.1088/1475-7516/2021/11/017 [arXiv:2102.09801 [nucl-th]]
  * [3] J. Rafelski, J. Birrell, A. Steinmetz and C. T. Yang, “A short survey of matter-antimatter evolution in the primordial Universe,”Universe 9, 309 (2023) doi:10.3390/Universe9070309 [arXiv:2305.09055 [hep-th]]
  * [4] E. E. Salpeter, “Electron screening and thermonuclear reactions,”Austral. J. Phys. 7 (1954), 373-388 doi:10.1071/PH540373
  * [5] P. Debye and E. Hückel, “Zur Theorie der Elektrolyte,”Physikalische Zeitschrift 24 185–206 (1923).
  * [6] E. E. Salpeter and H. M. van Horn, “Nuclear Reaction Rates at High Densities,”Astrophys. J., 155, 183 (1969) doi:10.1086/149858
  * [7] M. A. Famiano, A. B. Balantekin and T. Kajino, “Low-lying Resonances and Relativistic Screening in Big-Bang Nucleosynthesis,”Phys. Rev. C 93 (2016) no.4, 045804 doi:10.1103/PhysRevC.93.045804 [arXiv:1603.03137 [astro-ph.CO]]
  * [8] C. Carraro, A. Schafer and S. E. Koonin, “Dynamic Screening of Thermonuclear Reactions,”Astrophys. J., 331, 565 (1988) doi:10.1086/166582
  * [9] A. V. Gruzinov, “Dynamic screening and thermonuclear reaction rates,”Astrophys. J. 496 (1998), 503 doi:10.1086/305349 [arXiv:astro-ph/9702064 [astro-ph]]
  * [10] M. Opher and R. Opher, “Dynamic screening in thermonuclear reactions,” Astrophys. J. 535 (2000), 473-474 doi:10.1086/308808 [arXiv:astro-ph/9908218 [astro-ph]].
  * [11] X. Yao, T. Mehen and B. Müller, “Dynamical screening of $\alpha$-$\alpha$ resonant scattering and thermal nuclear scattering rate in a plasma,” Phys. Rev. D 95 (2017) no.11, 116002 doi:10.1103/PhysRevD.95.116002 [arXiv:1609.00383 [nucl-th]].
  * [12] N. Sasankan, A. Kedia, M. Kusakabe and G. J. Mathews, ‘Analysis of the Multi-component Relativistic Boltzmann Equation for Electron Scattering in Big Bang Nucleosynthesis,” Phys. Rev. D 101 (2020) no.12, 123532 doi:10.1103/PhysRevD.101.123532 [arXiv:1911.07334 [nucl-th]].
  * [13] A. Kedia, N. Sasankan, G. J. Mathews and M. Kusakabe, “Simulations of Multi-Component Relativistic Thermalization,” Phys. Rev. E 103 (2021), 032101 doi:10.1103/PhysRevE.103.032101 [arXiv:2004.13186 [cond-mat.stat-mech]].
  * [14] C. Pitrou, A. Coc, J. P. Uzan and E. Vangioni, “Precision big bang nucleosynthesis with improved Helium-4 predictions,”Phys. Rept. 754, 1-66 (2018) doi:10.1016/j.physrep.2018.04.005 [arXiv:1801.08023 [astro-ph.CO]]
  * [15] C. Pitrou, A. Coc, J. P. Uzan and E. Vangioni, “Resolving conclusions about the early Universe requires accurate nuclear measurements,“Nature Rev. Phys. 3, 231 (2021) doi:10.1038/s42254-021-00294-6 [arXiv:2104.11148 [astro-ph.CO]]
  * [16] C. A. Bertulani, F. W. Hall and B. I. Santoyo, “Big Bang nucleosynthesis as a probe of new physics,”EPJ Web Conf. 275, 01003 (2023) doi:10.1051/epjconf/202327501003 [arXiv:2210.04071 [nucl-th]]
  * [17] U. G. Meißner, B. C. Metsch and H. Meyer, “The electromagnetic fine-structure constant in primordial nucleosynthesis revisited,”[arXiv:2305.15849 [hep-th]]
  * [18] P. A. Haro et al, “Spectroscopic verification of very luminous galaxy candidates in the early universe,”[arXiv:2303.08918 [astro-ph.GA]]
  * [19] N. Sabti, J. B. Muñoz and M. Kamionkowski, “Insights from HST into Ultra-Massive Galaxies and Early-Universe Cosmology,”[arXiv:2305.07049 [astro-ph.CO]].
  * [20] Cosmin Ilie, Jillian Paulin, and Katherine Freese, “Supermassive Dark Star candidates seen by JWST,” PNAS 120, e2305762120 (2023) doi:10.1073/pnas.2305762120
  * [21] P. L. Bhatnagar, E. P. Gross and M. Krook, “A Model for Collision Processes in Gases. 1. Small Amplitude Processes in Charged and Neutral One-Component Systems,”Phys. Rev. 94, 511 (1954) doi:10.1103/PhysRev.94.511
  * [22] M. Formanek, C. Grayson, J. Rafelski and B. Müller, “Current-conserving relativistic linear response for collisional plasmas,”Annals Phys. 434 (2021), 168605 doi:10.1016/j.aop.2021.168605 [arXiv:2105.07897 [physics.plasm-ph]]
  * [23] D. Montgomery, G. Joyce, and R. Sugihara, “Inverse third power law for the shielding of test particles,”J. Plasma Phys. 4, 1-13 (1970)
  * [24] L. Stenflo, M. Y. Yu, and P. K. Shukla “Shielding of a slow test charge in a collisional plasma”Phys. Fluids 16, 450 (1973) doi:10.1063/1.1694361
  * [25] P. K. Shukla and A. A. Mamun “Introduction to Dusty Plasma Physics”Plasma Phys. Control. Fusion 44 395 (2002) doi:10.1088/0741-3335/44/3/701
  * [26] Martin Lampe, Glenn Joyce, Gurudas Ganguli and Valeriy Gavrishchaka “Interactions between dust grains in a dusty plasma”Physics of Plasmas 7, 3851 (2000) doi:10.1063/1.1288910
  * [27] C. Grayson, M. Formanek, J. Rafelski and B. Mueller, “Dynamic magnetic response of the quark-gluon plasma to electromagnetic fields,”Phys. Rev. D 106 (2022) no.1, 014011 doi:10.1103/PhysRevD.106.014011 [arXiv:2204.14186 [hep-ph]]
  * [28] J. Birrell, C. T. Yang, P. Chen and J. Rafelski, “Relic neutrinos: Physically consistent treatment of effective number of neutrinos and neutrino mass,”Phys. Rev. D 89, 023008 (2014) doi:10.1103/PhysRevD.89.023008
  * [29] R. L. Workman et al. [Particle Data Group], “Review of Particle Physics,”PTEP 2022, 083C01 (2022) doi:10.1093/ptep/ptac097
  * [30] J. Letessier and J. Rafelski, Hadrons and quark-gluon plasma, Cambridge University Press (2002, 2022) doi:10.1017/9781009290753 (Open Access)
  * [31] I. Kuznetsova, D. Habs and J. Rafelski, “Thermal reaction processes in a relativistic QED plasma drop,”Phys. Rev. D 81, 053007 (2010) doi:10.1103/PhysRevD.81.053007
  * [32] I. Kuznetsova and J. Rafelski, “Electron-Positron Plasma Drop Formed by Ultra-Intense Laser Pulses,”Phys. Rev. D 85, 085014 (2012) doi:10.1103/PhysRevD.85.085014
  * [33] M. B. Kislinger and P. D. Morley, “Collective Phenomena in Gauge Theories. 1. The Plasmon Effect for Yang-Mills Fields,”Phys. Rev. D 13, 2765 (1976) doi:10.1103/PhysRevD.13.2765
  * [34] R. Starke and G. A. H. Schober, “Relativistic covariance of Ohm’s law,”Int. J. Mod. Phys. D 25, 1640010 (2016) doi:10.1142/S0218271816400101 [arXiv:1409.3723 [math-ph]]
  * [35] D. Melrose, Quantum Plasmadynamics: Unmagnetized Plasmas, Lect. Notes Phys. 735 (Springer, New York, 2008) doi:10.1007/978-0-387-73903-8
  * [36] A. Steinmetz, C. T. Yang, and J. Rafelski, “Matter-antimatter origin of cosmic magnetism,”[arXiv:2308.14818 [hep-ph]]
  * [37] P. K. Shukla and N. N. Rao “Coulomb crystallization in colloidal plasmas with streaming ions and dust grains”Physics of Plasmas 3, 1770 (1996) doi: 10.1063/1.871695
  * [38] R. Hakim ”Relativistic Perfect Gas in an External Force Field,”Phys. Rev. 162, 128 (1967) doi:10.1103/PhysRev.162.128
  * [39] S. R. De Groot, W. A. Van Leeuwen and C. G. Van Weert, Relativistic Kinetic Theory. Principles and Applications, North-Holland Publishing Company (1980)
  * [40] G. Gouy “Sur la constitution de la charge électrique à la surface d’un électrolyte,”J. Phys. Theor. Appl., 9 , 457 (1910) doi:10.1051/jphystap:019100090045700
  * [41] H. Thomas, G. E. Morfill, V. Demmel, J. Goree, B. Feuerbacher and D. Möhlmann, “Plasma Crystal: Coulomb Crystallization in a Dusty Plasma,”Phys. Rev. Lett. 73, 652 (1994) doi:10.1103/PhysRevLett.73.652
  * [42] R. Gupta, “JWST early Universe observations and $\Lambda$CDM cosmology” MNRAS stad2032, (2023) doi:10.1093/mnras/stad2032
  * [43] C. Labaune, C. Baccou, S. Depierreux, C. Goyon, G. Loisel, V. Yahia and J. Rafelski, “Fusion reactions initiated by laser-accelerated particle beams in a laser-produced plasma,”Nature Commun. 4, 2506 (2013) doi:10.1038/ncomms3506 [arXiv:1310.2002 [physics.plasm-ph]]
  * [44] Daniele Margarone, Julien Bonvalet, Lorenzo Giuffrida, Alessio Morace, Vasiliki Kantarelou, Marco Tosca, Didier Raffestin, Philippe Nicolai, Antonino Picciotto, Yuki Abe, Yasunobu Arikawa, Shinsuke Fujioka, Yuji Fukuda, Yasuhiro Kuramitsu, Hideaki Habara and Dimitri Batani, “In-Target Proton–Boron Nuclear Fusion Using a PW-Class Laser,”MDPI Applied Sciences , 12(3), [e1444]. doi:10.3390/app12031444
  * [45] N. Buyukcizmeci, T. Reichert, A. S. Botvina and M. Bleicher, “Nucleosynthesis of light nuclei and hypernuclei in central Au+Au collisions at $\sqrt{s_{NN}}$=3 GeV,”[arXiv:2306.17145 [nucl-th]]
  * [46] J. Zhao, K. Zhou, S. Chen and P. Zhuang, “Heavy flavors under extreme conditions in high energy nuclear collisions,”Prog. Part. Nucl. Phys. 114, 103801 (2020) doi:10.1016/j.ppnp.2020.103801 [arXiv:2005.08277 [nucl-th]]
  * [47] H. De Vries, C. W. De Jager and C. De Vries, “Nuclear charge and magnetization density distribution parameters from elastic electron scattering,”Atom. Data Nucl. Data Tabl. 36, 495 (1987) doi:10.1016/0092-640X(87)90013-1