PP Chains#
There are multiple proton-proton chains.
PP-I involves: \(p\), \(d\), \({}^3\mathrm{He}\), \({}^4\mathrm{He}\)
PP-II adds: \({}^7\mathrm{Li}\), \({}^7\mathrm{Be}\)
PP-III adds: \({}^8\mathrm{Be}\), \({}^8\mathrm{B}\)
Here we can explore them all.
Creating a network#
We’ll use pynucastro to create a network with all of these nuclei.
import pynucastro as pyna
net = pyna.network_helper(["p", "d", "he3", "he4",
"li7", "be7", "be8", "b8"])
warning: He4 was not able to be linked in TabularLibrary
warning: d was not able to be linked in TabularLibrary
warning: B8 was not able to be linked in TabularLibrary
warning: Li7 was not able to be linked in TabularLibrary
warning: p was not able to be linked in TabularLibrary
warning: He3 was not able to be linked in TabularLibrary
warning: Be8 was not able to be linked in TabularLibrary
warning: Be7 was not able to be linked in TabularLibrary
Note
ReacLib does not provide the rate for \({}^8\mathrm{Be} \rightarrow \alpha + \alpha\) because it is so short lived. It essentially will break apart into two \(\alpha\) immediately. It instead has a rate that combines:
\[{}^8\mathrm{B} \rightarrow {}^8\mathrm{Be} + e^+ + \nu_e\]
\[{}^8\mathrm{Be} \rightarrow \alpha + \alpha\]
Approximating composition in Sun’s core#
We need a reasonably accurate estimate for the composition in the core of the Sun. We’ll find this by integrating an initial H/He mix for 4.5 billion years–this should find an equilibrium
tmax = 4.5e9 * 365.25 * 24 * 3600
tmax
1.420092e+17
net.write_network("pp.py")
import pp
/opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/pynucastro/rates/derived_rate.py:120: UserWarning: He3 partition function is not supported by tables: set log_pf = 0.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set log_pf = 0.0 by default'))
/opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/pynucastro/rates/derived_rate.py:120: UserWarning: d partition function is not supported by tables: set log_pf = 0.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set log_pf = 0.0 by default'))
/opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/pynucastro/rates/derived_rate.py:120: UserWarning: Be7 partition function is not supported by tables: set log_pf = 0.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set log_pf = 0.0 by default'))
/opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/pynucastro/rates/derived_rate.py:120: UserWarning: B8 partition function is not supported by tables: set log_pf = 0.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set log_pf = 0.0 by default'))
/opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/pynucastro/rates/derived_rate.py:120: UserWarning: Li7 partition function is not supported by tables: set log_pf = 0.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set log_pf = 0.0 by default'))
We’ll integrate this using solve_ivp
import numpy as np
from scipy.integrate import solve_ivp
X0 = np.zeros(pp.nnuc)
X0[pp.jp] = 0.72
X0[pp.jhe4] = 0.28
Y0 = X0 / pp.A
rho = 150
T = 1.5e7
sol = solve_ivp(pp.rhs, [0, tmax], Y0, method="BDF", jac=pp.jacobian,
dense_output=True, args=(rho, T), rtol=1.e-7, atol=1.e-12)
---------------------------------------------------------------------------
TypingError Traceback (most recent call last)
Cell In[8], line 1
----> 1 sol = solve_ivp(pp.rhs, [0, tmax], Y0, method="BDF", jac=pp.jacobian,
2 dense_output=True, args=(rho, T), rtol=1.e-7, atol=1.e-12)
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/scipy/integrate/_ivp/ivp.py:623, in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)
620 if method in METHODS:
621 method = METHODS[method]
--> 623 solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
625 if t_eval is None:
626 ts = [t0]
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/scipy/integrate/_ivp/bdf.py:206, in BDF.__init__(self, fun, t0, y0, t_bound, max_step, rtol, atol, jac, jac_sparsity, vectorized, first_step, **extraneous)
204 self.max_step = validate_max_step(max_step)
205 self.rtol, self.atol = validate_tol(rtol, atol, self.n)
--> 206 f = self.fun(self.t, self.y)
207 if first_step is None:
208 self.h_abs = select_initial_step(self.fun, self.t, self.y,
209 t_bound, max_step, f,
210 self.direction, 1,
211 self.rtol, self.atol)
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/scipy/integrate/_ivp/base.py:158, in OdeSolver.__init__.<locals>.fun(t, y)
156 def fun(t, y):
157 self.nfev += 1
--> 158 return self.fun_single(t, y)
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/scipy/integrate/_ivp/base.py:24, in check_arguments.<locals>.fun_wrapped(t, y)
23 def fun_wrapped(t, y):
---> 24 return np.asarray(fun(t, y), dtype=dtype)
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/scipy/integrate/_ivp/ivp.py:595, in solve_ivp.<locals>.fun(t, x, fun)
594 def fun(t, x, fun=fun):
--> 595 return fun(t, x, *args)
File ~/work/stars/stars/content/h-burning/pp.py:733, in rhs(t, Y, rho, T, screen_func)
732 def rhs(t, Y, rho, T, screen_func=None):
--> 733 return rhs_eq(t, Y, rho, T, screen_func)
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/numba/core/dispatcher.py:424, in _DispatcherBase._compile_for_args(self, *args, **kws)
420 msg = (f"{str(e).rstrip()} \n\nThis error may have been caused "
421 f"by the following argument(s):\n{args_str}\n")
422 e.patch_message(msg)
--> 424 error_rewrite(e, 'typing')
425 except errors.UnsupportedError as e:
426 # Something unsupported is present in the user code, add help info
427 error_rewrite(e, 'unsupported_error')
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/numba/core/dispatcher.py:365, in _DispatcherBase._compile_for_args.<locals>.error_rewrite(e, issue_type)
363 raise e
364 else:
--> 365 raise e.with_traceback(None)
TypingError: Failed in nopython mode pipeline (step: nopython frontend)
NameError: name 'log_scor_He2_He4' is not defined
During: Pass nopython_type_inference
Now let’s get the composition at the final time
X_final = sol.y[:, -1] * pp.A
X_final
array([2.81318714e-01, 2.04411172e-18, 8.77279202e-06, 7.18672513e-01,
2.03687202e-15, 9.99888210e-12, 4.53092291e-54, 1.95666221e-21])
and create a pynucastro Composition from this
comp = pyna.Composition(net.unique_nuclei)
comp.set_array(X_final)
print(comp)
X(p) : 0.28131871449282164
X(d) : 2.0441117249462732e-18
X(He3) : 8.772792020130378e-06
X(He4) : 0.7186725127051565
X(Li7) : 2.0368720231888285e-15
X(Be7) : 9.998882102842114e-12
X(Be8) : 4.530922913017406e-54
X(B8) : 1.956662207786176e-21
Importance of PP-I vs. PP-II#
We can evaluate the rates with this composition
r_sun = net.evaluate_rates(rho, T, comp)
/raid/zingale/development/pynucastro/pynucastro/rates/derived_rate.py:126: UserWarning: He3 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/raid/zingale/development/pynucastro/pynucastro/rates/derived_rate.py:126: UserWarning: d partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/raid/zingale/development/pynucastro/pynucastro/rates/derived_rate.py:126: UserWarning: Be7 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/raid/zingale/development/pynucastro/pynucastro/rates/derived_rate.py:126: UserWarning: B8 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/raid/zingale/development/pynucastro/pynucastro/rates/derived_rate.py:126: UserWarning: Li7 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
Let’s look at just the rates that consume \({}^3\mathrm{He}\)
for r in r_sun:
if pyna.Nucleus("he3") in r.reactants:
print(f"{str(r):30}: {r_sun[r]:12.6g}")
He3 + p ⟶ He4 + e⁺ + 𝜈 : 2.73268e-26
He3 + He4 ⟶ Be7 + 𝛾 : 1.97861e-19
He3 + H2 ⟶ p + He4 : 9.85879e-23
He3 + He3 ⟶ p + p + He4 : 1.42325e-19
Be7 + He3 ⟶ p + p + He4 + He4 : 2.64203e-39
He3 ⟶ p + H2 : 0
Notice that the rate of \({}^3\mathrm{He}({}^3\mathrm{He},pp){}^4\mathrm{He}\) is only slightly slower than \({}^4\mathrm{He}({}^3\mathrm{He},\gamma){}^7\mathrm{Be}\) at this point.
r33 = net.get_rate_by_name("he3(he3,pp)he4")
r34 = net.get_rate_by_name("he4(he3,g)be7")
ratio = r_sun[r33] / r_sun[r34]
ratio
np.float64(0.7193207261903833)
fig = net.plot(rho, T, comp,
curved_edges=True,
always_show_alpha=True,
always_show_p=True,
ydot_cutoff_value=1.e-25)
Now let’s consider a bit hotter. We’ll keep the same composition
T = 3.5e7
fig = net.plot(rho, T, comp,
curved_edges=True,
always_show_alpha=True,
always_show_p=True,
ydot_cutoff_value=1.e-25)
r_hot = net.evaluate_rates(rho, T, comp)
ratio_hot = r_hot[r33] / r_hot[r34]
ratio_hot
np.float64(0.41694689409397745)
Now we see that \({}^4\mathrm{He}({}^3\mathrm{He},\gamma){}^7\mathrm{Be}\) is much faster.
Overall, as the temperature increases, we see that the \({}^4\mathrm{He}({}^3\mathrm{He},\gamma){}^7\mathrm{Be}\) rate is more active.
Note
These rates are strongly dependent on the amount of \({}^3\mathrm{He}\) present, so we should really compute the equilibrium abundance before making the comparisons.