DYNAMIC SCREENING CORRECTION FOR SOLAR p–p REACTION RATES

In this section, we derive an expression for calculating nuclear reaction rates following the treatment of Clayton (1968). The reaction rate per unit volume between particles of types α and β is a product of the number densities of the particles, n_α and n_β, and the average value of the product of the relative velocity v times the cross section σ,

$$\begin{eqnarray} r_{\alpha \beta } & = & \frac{n_{\alpha }n_{\beta }}{1+\delta _{\alpha \beta }}\langle \sigma v \rangle _{\alpha \beta } \nonumber \\ & = & \frac{n_{\alpha }n_{\beta }}{1+\delta _{\alpha \beta }} \int _0^\infty \psi (E)v(E)\sigma (E)dE, \end{eqnarray} \tag{ 1 }$$

where ψ(E) is the relative velocity distribution and δ_αβ accounts for reactions of like particles. Because we will be dealing with the correction due to dynamic screening (a ratio between the unscreened and screened reaction rates), we can ignore the density factor and focus on the reaction rate per pair of particles,

$$\begin{eqnarray} \lambda &=& \langle \sigma v \rangle _{\alpha \beta } \nonumber \\ &=& f(\mu,T) \int _0^\infty E\;{\rm {exp}}\left(-\frac{E}{k_BT}\right) \sigma (E)dE, \end{eqnarray} \tag{ 2 }$$

where

$$\begin{equation} f(\mu,T) = \sqrt{8 \over {\pi \mu }}\:\left(\frac{1}{k_B T}\right)^{3/2}, \end{equation} \tag{ 3 }$$

μ is the reduced mass of the pair, and the Maxwell–Boltzmann distribution is used for ψ(E). The cross section σ(E) can be defined as a product of three separate energy-dependent factors

$$\begin{equation} \sigma (E) = \frac{S(E)}{E} \; {\rm {exp}}\left(\frac{-b}{\sqrt{E}}\right), \end{equation} \tag{ 4 }$$

where b = 31.28Z_αZ_βA^1/2 keV^1/2, with Z_α and Z_β being the charges of the interacting ions and A is the reduced atomic weight. The exponential factor in this expression comes from the barrier penetration probability, the inverse energy dependence comes from the quantum–mechanical interaction between the two particles, and S(E) contains the intrinsically nuclear parts of the probability for a nuclear reaction to occur. With this substitution for σ(E), Equation (2) can be re-written as

$$\begin{equation} \lambda = f(\mu,T) \int _0^\infty S(E)\:{\rm {exp}}\left(-{{E} \over {k_B T}} - {{b} \over {\sqrt{E}}} \right) dE. \end{equation} \tag{ 5 }$$

In the non-resonant reaction case, S(E) is slowly varying with E, so we can treat it as a constant S₀ evaluated at the energy where exp(−E/k_BT − b/E^1/2) is maximum. Then the reaction rate per pair of particles (without screening) can be computed as

$$\begin{equation} \lambda = f(\mu,T) S_0 \int _0^\infty {\rm {exp}}\left(-{{E} \over {k_B T}} - {{b} \over {\sqrt{E}}} \right) dE. \end{equation} \tag{ 6 }$$

Salpeter (1954) developed a treatment to include the effect of static electron screening on nuclear reaction rates. Here we summarize his method which we will use in Section 2 as the inspiration for our calculation of the dynamic screening correction.

We begin by writing the total interaction energy as a combination of the bare Coulomb potential and a contribution from the plasma:

$$\begin{equation} U_{\rm {total}}(r) = \frac{Z_1Z_2e^2}{r} + U(r). \end{equation} \tag{ 7 }$$

Then consider a case in which the classical impact parameter r_c is very small compared with the charge cloud radius R_D and the nuclear radius r_n is much smaller than r_c. Then the barrier penetration factor for r_n < r < r_c depends only on the expression

$$\begin{equation} E - U(r) - \frac{Z_1Z_2e^2}{r}. \end{equation} \tag{ 8 }$$

For distances larger than r_c, the barrier penetration factor hardly depends on the potential. U(r) must be small for distances greater than R_D and approach a constant value U₀ of the order of magnitude of Z₁Z₂e²/R_D for small r. Then,

$$\begin{equation} \frac{r_{\rm {c}}}{R_D} \approx \frac{U_0}{E_{\rm {max}}}\ll 1, \end{equation} \tag{ 9 }$$

where E_max is the relative kinetic energy for which the integrand in Equation (2) reaches a sharp maximum. If this inequality is satisfied, U(r) can be replaced by the potential at the origin U₀. By examining expression (8), we can see that the screening potential has effectively increased the kinetic energy by a magnitude of U₀, so the cross section factors for U_total = U_Coulomb + U₀ are equivalent to the unscreened factors with energy E − U₀. Equation (2) can then be replaced by

$$\begin{equation} \lambda = f(\mu,T) \int _0^\infty E\;{\rm {exp}}\left(-\frac{E}{k_BT}\right) \sigma (E-U_0)dE. \end{equation} \tag{ 10 }$$

With the change of variables E' = E − U₀ and the approximation (E' + U₀) ≈ E', the reaction rate per pair of particles becomes

$$\begin{equation} \lambda = f(\mu,T) \int _{-U_0}^\infty E^{\prime }\;{\rm {exp}}\left(-\frac{E^{\prime }+U_0}{k_BT}\right) \sigma (E^{\prime })dE^{\prime }. \end{equation} \tag{ 11 }$$

Because the penetration factor for E' = −U₀ is so small, the lower limit of the integral can be set to zero without significantly changing the value of the integral. We then see that

$$\begin{equation} \lambda _{\rm {screened}}={\rm {exp}}\left(\frac{-U_0}{k_BT}\right)\lambda _{\rm {bare}}, \end{equation} \tag{ 12 }$$

illustrating that the reaction rate for the statically screened potential can be approximated by multiplying the rate for the bare Coulomb potential by exp(−U₀/k_BT).

Salpeter then derives U₀ by solving the Poisson–Boltzmann equation for electrons and ions in a plasma under the condition of weak screening (U_total(r) ≪ k_BT). He arrives at an expression for the screening energy that is equivalent to that of the Debye–Hückel theory of dilute solutions of electrolytes (Debye & Hückel 1923):

$$\begin{equation} U_0 = -\frac{Z_1Z_2e^2}{R_D}. \end{equation} \tag{ 13 }$$

The Debye length, R_D, is the characteristic screening length of a plasma at temperature T with number density n which is defined by

$$\begin{equation} R_D^2 = \frac{\epsilon _0 k_B T}{e^2(n_e+n_iZ_i^2)}. \end{equation} \tag{ 14 }$$

For a neutral proton–electron plasma, the electron number density n_e and ion number density n_i are both n/2, so the Debye length is just

$$\begin{equation} R_D^2 = \frac{ \epsilon _0 k_B T}{ne^2}. \end{equation} \tag{ 15 }$$

Although Salpeter's expression accurately describes the effect of static screening, the issue of dynamic screening in the hot, dense plasma of the solar interior remains an open question. Dynamic screening occurs when the screened interaction energy of a pair of ions depends on the relative velocity of the pair. Most of the ions in the solar plasma are much slower than the electrons and the fastest ions. The thermal ions are therefore not able to rearrange themselves as quickly around individual faster moving ions. Since nuclear reactions require energies several times the average thermal energy, the ions that are able to engage in nuclear reactions in the Sun are the faster moving ions, which are not accompanied by a full static screening cloud.

Salpeter's derivation uses the mean-field approach in which the many-body interactions are reduced to an average interaction that simplifies calculations. This technique is quite useful for calculations describing the average behavior of the plasma. However, dynamic effects for the fast-moving, interacting ions in hot, dense plasma lead to a screened potential that deviates from the average value. Therefore, the mean-field approximation is not appropriate for computing stellar nuclear reaction rates. Instead, we use the molecular-dynamics method of Shaviv & Shaviv (1996) to model the motion of protons and electrons in a plasma under solar conditions in order to investigate dynamic screening in p–p reactions. The advantage of the molecular-dynamics method is that it does not assume a mean field. Nor does it assume a long-time average potential for the scattering of any two charges, which is necessary in the statistical way to solve Poisson's equation to obtain the mean potential in a plasma.

In previous work, Mao et al. (2009) present simulation results for the velocity-dependent screening energy of p–p reactions in a plasma with the temperature and density of the solar core (T = 1.6 × 10⁷ K, ρ = 1.6 × 10⁵ kg m⁻³). They demonstrate that the static screening result does not accurately represent this plasma, and they compute a screening energy that depends on the relative kinetic energy of a pair of interacting ions (see Figure 1). In this paper, we use their simulation results to compute a correction to the solar p–p reaction rate due to the dynamic screening they observe.

Zoom In Zoom Out Reset image size

How does the dynamic screening energy seen in Figure 1 result in a reaction rate that is so close to the unscreened reaction rate? To answer that question, we compare the components of the integrands from Equations (6) and (20). Both integrands can be written in the general form

$$\begin{equation} F(E) = H(E) J(E), \end{equation} \tag{ 21 }$$

where

$$\begin{equation} H(E) = \frac{E}{G(E)}, \end{equation} \tag{ 22 }$$

$$\begin{equation} J(E) = {\rm {exp}}\left[-\frac{E}{k_BT}-\frac{b}{\sqrt{G(E)}} \right], \end{equation} \tag{ 23 }$$

and

$$\begin{equation} G(E) = E-U. \end{equation} \tag{ 24 }$$

The three cases only differ in the screening energy U which is shown for each case in Table 1.

In Figure 2(a), we see that the dynamic G(E) approaches the unscreened G(E) for high energies. Figure 2(b) shows that this leads to the dynamic H(E) approaching the unscreened H(E) for energies above ∼2 k_BT, while the static H(E) is very different from both the dynamic and bare H(E). Below energies of ∼2 k_BT the dynamic H(E) drops rapidly away from the unscreened value of 1. However, the factor H(E) is multiplied by J(E) in the integrand of the reaction-rate equations. As seen in Figure 2(c), J(E) is very small to zero below energies of ∼2 k_BT, damping out the region of H(E) in which the unscreened and dynamic results diverge. This leads to integrands for the dynamic and unscreened cases that are nearly identical, as seen when H(E) and J(E) are multiplied together to give F(E) in Figure 2(d). The Gamow-peak-like factor J(E) acts as a weighting function to devalue the H(E) contribution of the slow pairs of ions that rarely participate in nuclear reactions. The faster pairs that cause less polarization of the surrounding plasma and therefore see less screening provide the main contribution to the reaction-rate integral.

Zoom In Zoom Out Reset image size

Now that we have defined the integrand F(E) for the dynamic case, we can return to the issue of assumptions and approximations that were justified in the static screening rate derivation and adopted in the dynamic screening rate derivation. Here we assess the validity of these assumptions and approximations in the case of dynamic screening.

The first assumption is that r_c is very small compared with the charge cloud radius R_D and that the nuclear radius r_n is much smaller than r_c, leading to expression (8). Although the dynamic case does not have a traditional static screening cloud, the inequality r_n ≪  r_c ≪  R_D is still satisfied. This can be seen in the definition of R_D as the distance beyond which an appreciable fraction of the nuclear charge is screened by the polarization charge cloud, a distance that is large for the dynamic case.

Next, we examine the inequality in Equation (9), U₀/E_max ≪ 1. In the dynamic case, E_max is the kinetic energy for which the integrand of Equation (20) is maximum. We see from Figure 2(d) and U₀ in Figure 1 that the peak of the integrand F(E) occurs at energies for which U₀/E_max ≪ 1.

Finally, we address the approximation that U(r) can be replaced by the potential at the origin U₀. While this is easy to see for the static case, in the case of dynamic screening we do not have the screening energy in the form U(r, E) to examine this claim. Instead, we define U₀(E) to be the screening energy computed at the turning point of the approaching protons, since this is the relevant U(r, E) for nuclear reaction-rate calculations.

The argument for dynamic screening has been damaged by several derivations of alternate screening formulae (see, for example, Carraro et al. 1988; Opher & Opher 2000; Shaviv & Shaviv 1996; Tsytovich 2000) that were subsequently shown to be incorrect. Here we discuss two examples, Carraro et al. (1988) and Shaviv & Shaviv (1996). Carraro et al. (1988) derived a modified screening potential for fusing ions when the Gamow velocity is greater than the thermal velocity. Brown & Sawyer (1997) showed that including processes of excitation or de-excitation of the plasma in an interaction with one of the fusing ions exactly cancels the dynamic modifications proposed by Carraro et al. (1988). Shaviv & Shaviv (1996) then introduced a factor of 3/2 on the screening energy. They arrived at this result by including the interaction of the screening cloud from each fusing ion in the total interaction potential. Bruggen & Gough (1997) showed that Shaviv & Shaviv (1996) misinterpreted the thermodynamics and used an incorrect potential in the Schrödinger equation for the system. Bahcall et al. (2002) summarize the problems with several different alternative screening formulae. However, finding flaws in these (and other) derivations of analytical expressions for screening deviations does not rule out the effect of dynamic screening. This argument only highlights the difficulty in developing a general analytical formalism to describe dynamic screening.

For the Gibbs distribution, probabilities for momenta and coordinates are independent and cannot influence each other. This leads to the argument that velocities of fusing particles cannot have an effect on screening because the distribution function is not factorable. However, it is clear when examining individual ions that their relative velocity affects how close the ions can be to each other. This argument extends to ions in a plasma, where the configuration of the screening cloud of approaching ions depends on the relative velocity of those ions. Over the whole system, this velocity-dependent effect averages out to the Gibbs distribution, but each screening cloud is not identical to the average configuration of a screening cloud in that system.

Solar nuclear reactions select a biased sample of the ions in the system. These nuclear reactions involve mainly the fastest ions, not a random sample of all ions in the system. Therefore, the velocity distribution must be multiplied by the velocity-dependent screening energy before integration instead of beginning with an average value of the distribution.

Do dynamic screening results imply that higher-ordered terms are required? The screening energy can be expressed as a power series expansion in the plasma coupling parameter Λ:

$$\begin{equation} \frac{U(r)}{k_BT} = -\frac{\Lambda }{x} \left(1 - \exp \left(-x \right) \right) - \Lambda ^2 f(x,\Lambda) \;, \end{equation} \tag{ 25 }$$

where

$$\begin{equation} \Lambda = \frac{Z_1Z_2e^2}{R_Dk_BT} \;, \end{equation} \tag{ 26 }$$

x = r/R_D, and f(x, Λ) is given by DeWitt (1965). The first term reduces to the Debye–Hückel weak screening result shown in Equation (13).

For the temperature and density of our simulations (T = 1.6 × 10⁷ K, ρ = 1.6 × 10⁵ kg m⁻³), Λ = 0.05, so higher-order terms are small. Therefore, if dynamic screening could be described by higher-order terms in the expansion, the effect should be much smaller than the first term. However, the dynamic correction is not just a higher-order term in the expansion. Dynamic effects come from a different approach to determining screening effects. Instead of deriving an expression for screening based on average properties of the system, we examine the formation of the screening clouds themselves and do not average out the velocity-dependent nature of the clouds.

What observational evidence can confirm any screening effect in the nuclear reactions in the Sun? Many early discussions of dynamic screening were motivated by the neutrino problem which has since been resolved with neutrino oscillation theory. In addition, including dynamic screening corrections in models with the Grevesse & Noels (1993) or Grevesse & Sauval (1998) abundances worsened agreement with helioseismic constraints (see, for example, Weiss et al. 2001). However, the solar abundance problem provides renewed motivation for exploring dynamic screening in solar nuclear reactions.

Before 2005, solar models with the latest input physics reproduced the sound-speed profile determined from helioseismic inversions to within 0.4% and also provided good agreement with the seismically inferred convection zone depth and convection zone helium abundance. Then Asplund et al. (2005, 2009), Caffau et al. (2008, 2009), and Ludwig et al. (2009) began using three-dimensional hydrodynamic models of the solar atmosphere with improved input physics and non-local thermodynamic equilibrium effects to determine solar atmospheric abundances. These revised calculations lowered element abundances by up to 1/3. When the lower abundances are incorporated in solar models, the sound-speed profiles, convection zone depths, and convection zone helium abundances give worse agreement with helioseismic constraints than models with the old, higher abundances. Many attempts have been made to improve agreement by adjusting the physics or evolutionary assumptions in solar models (see Basu & Antia 2008; Guzik 2008; Guzik & Mussack 2010). Although these adjustments have shown some improvement, no model using the new lower abundances agrees as well with the helioseismic constraints as the models using the older abundances.

In a forthcoming paper, K. Mussack & J. A. Guzik (2011, in preparation) incorporate the dynamic screening correction shown here for solar p–p reaction rates into solar models with the new lower abundances. They show that including this correction in solar models improves the sound-speed discrepancy in the solar core, as shown in Figure 3. This improvement does not fully reconcile the new abundances with helioseismic constraints, but it is a step in the right direction. Perhaps in combination with other changes, dynamic screening corrections could contribute to a solution to the solar abundance problem.

Zoom In Zoom Out Reset image size