1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
/* Copyright 2021 Neil Zaim
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef PROTON_BORON_FUSION_CROSS_SECTION_H
#define PROTON_BORON_FUSION_CROSS_SECTION_H
#include "Utils/WarpXConst.H"
#include <AMReX_Math.H>
#include <AMReX_REAL.H>
#include <cmath>
/**
* \brief Computes the total proton-boron fusion cross section in the range 0 < E < 3.5 MeV using
* the analytical fits given in W.M. Nevins and R. Swain, Nuclear Fusion, 40, 865 (2000).
* For the record, note that there is a typo in equation (1) of this paper: the total cross section
* should read S(E)/E*exp(-sqrt(E_G/E)) instead of S(E)/E*exp(sqrt(E_G/E)) (minus sign in the
* exponential).
*
* @param[in] E_keV the kinetic energy of the proton-boron pair in its center of mass frame, in
* keV.
* @return The total cross section in barn.
*/
AMREX_GPU_HOST_DEVICE AMREX_INLINE
amrex::ParticleReal ProtonBoronFusionCrossSectionNevins (const amrex::ParticleReal& E_keV)
{
using namespace amrex::literals;
using namespace amrex::Math;
// If kinetic energy is 0, return a 0 cross section and avoid later division by 0.
if (E_keV == 0._prt) {return 0._prt;}
// Fits also use energy in MeV
const amrex::ParticleReal E_MeV = E_keV*1.e-3_prt;
constexpr amrex::ParticleReal joule_to_MeV = 1.e-6_prt/PhysConst::q_e;
// Compute Gamow factor, in MeV
constexpr auto one_pr = 1._prt;
constexpr auto Z_boron = 5._prt;
constexpr amrex::ParticleReal m_boron = 10.7319_prt * PhysConst::m_p;
constexpr amrex::ParticleReal m_reduced = m_boron / (one_pr + m_boron/PhysConst::m_p);
constexpr amrex::ParticleReal gamow_factor = m_reduced / 2._prt *
(PhysConst::q_e*PhysConst::q_e * Z_boron /
(2._prt*PhysConst::ep0*PhysConst::hbar)) *
(PhysConst::q_e*PhysConst::q_e * Z_boron /
(2._prt*PhysConst::ep0*PhysConst::hbar)) *
joule_to_MeV;
// Compute astrophysical factor, in MeV barn, using the fits
constexpr auto E_lim1 = 400._prt; // Limits between the different fit regions
constexpr auto E_lim2 = 642._prt;
amrex::ParticleReal astrophysical_factor;
if (E_keV < E_lim1)
{
constexpr auto C0 = 197._prt;
constexpr auto C1 = 0.24_prt;
constexpr auto C2 = 2.31e-4_prt;
constexpr auto AL = 1.82e4_prt;
constexpr auto EL = 148._prt;
constexpr auto dEL_sq = 2.35_prt*2.35_prt;
astrophysical_factor = C0 + C1*E_keV + C2*powi<2>(E_keV) +
AL/((E_keV - EL)*(E_keV - EL) + dEL_sq);
}
else if (E_keV < E_lim2)
{
constexpr auto D0 = 330._prt;
constexpr auto D1 = 66.1_prt;
constexpr auto D2 = -20.3_prt;
constexpr auto D5 = -1.58_prt;
const amrex::ParticleReal E_norm = (E_keV-400._prt) * 1.e-2_prt;
astrophysical_factor = D0 + D1*E_norm + D2*powi<2>(E_norm) + D5*powi<5>(E_norm);
}
else
{
constexpr auto A0 = 2.57e6_prt;
constexpr auto A1 = 5.67e5_prt;
constexpr auto A2 = 1.34e5_prt;
constexpr auto A3 = 5.68e5_prt;
constexpr auto E0 = 581.3_prt;
constexpr auto E1 = 1083._prt;
constexpr auto E2 = 2405._prt;
constexpr auto E3 = 3344._prt;
constexpr auto dE0_sq = 85.7_prt*85.7_prt;
constexpr auto dE1_sq = 234._prt*234._prt;
constexpr auto dE2_sq = 138._prt*138._prt;
constexpr auto dE3_sq = 309._prt*309._prt;
constexpr auto B = 4.38_prt;
astrophysical_factor = A0 / ((E_keV-E0)*(E_keV-E0) + dE0_sq) +
A1 / ((E_keV-E1)*(E_keV-E1) + dE1_sq) +
A2 / ((E_keV-E2)*(E_keV-E2) + dE2_sq) +
A3 / ((E_keV-E3)*(E_keV-E3) + dE3_sq) + B;
}
// Compute cross section, in barn
return astrophysical_factor/E_MeV*std::exp(-std::sqrt(gamow_factor/E_MeV));
}
/**
* \brief Computes the total proton-boron fusion cross section in the range E > 3.5 MeV using a
* simple power law fit of the data presented in Buck et al., Nuclear Physics A, 398(2), 189-202
* (1983) (data can also be found in the EXFOR database).
*
* @param[in] E_keV the kinetic energy of the proton-boron pair in its center of mass frame, in
* keV.
* @return The total cross section in barn.
*/
AMREX_GPU_HOST_DEVICE AMREX_INLINE
amrex::ParticleReal ProtonBoronFusionCrossSectionBuck (const amrex::ParticleReal& E_keV)
{
using namespace amrex::literals;
constexpr amrex::ParticleReal E_start_fit = 3500._prt; // Fit starts at 3.5 MeV
// cross section at E = E_start_fit, in barn
constexpr amrex::ParticleReal cross_section_start_fit = 0.2168440845211521_prt;
constexpr amrex::ParticleReal slope_fit = -2.661840717596765_prt;
// Compute fitted value
return cross_section_start_fit*std::pow(E_keV/E_start_fit, slope_fit);
}
/**
* \brief Computes the total proton-boron fusion cross section. When E_kin_star < 3.5 MeV, we use
* the analytical fits given in W.M. Nevins and R. Swain, Nuclear Fusion, 40, 865 (2000). When
* E_kin_star > 3.5 MeV, we use a simple power law fit of the data presented in Buck et al.,
* Nuclear Physics A, 398(2), 189-202 (1983). Both fits return the same value for
* E_kin_star = 3.5 MeV.
*
* @param[in] E_kin_star the kinetic energy of the proton-boron pair in its center of mass frame,
* in SI units.
* @return The total cross section in SI units (square meters).
*/
AMREX_GPU_HOST_DEVICE AMREX_INLINE
amrex::ParticleReal ProtonBoronFusionCrossSection (const amrex::ParticleReal& E_kin_star)
{
using namespace amrex::literals;
// Fits use energy in keV
constexpr amrex::ParticleReal joule_to_keV = 1.e-3_prt/PhysConst::q_e;
const amrex::ParticleReal E_keV = E_kin_star*joule_to_keV;
constexpr amrex::ParticleReal E_threshold = 3500._prt;
const amrex::ParticleReal cross_section_b = (E_keV <= E_threshold) ?
ProtonBoronFusionCrossSectionNevins(E_keV) :
ProtonBoronFusionCrossSectionBuck(E_keV);
// Convert cross section to SI units: barn to square meter
constexpr auto barn_to_sqm = 1.e-28_prt;
return cross_section_b*barn_to_sqm;
}
#endif // PROTON_BORON_FUSION_CROSS_SECTION_H
|
