WarpX
ProtonBoronFusionCrossSection.H
Go to the documentation of this file.
1 /* Copyright 2021 Neil Zaim
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_PROTON_BORON_FUSION_CROSS_SECTION_H
9 #define WARPX_PROTON_BORON_FUSION_CROSS_SECTION_H
10 
11 #include "Utils/WarpXConst.H"
12 
13 #include <AMReX_Math.H>
14 #include <AMReX_REAL.H>
15 
16 #include <cmath>
17 
30 amrex::ParticleReal ProtonBoronFusionCrossSectionNevins (const amrex::ParticleReal& E_keV)
31 {
32  using namespace amrex::literals;
33  using namespace amrex::Math;
34 
35  // If kinetic energy is 0, return a 0 cross section and avoid later division by 0.
36  if (E_keV == 0._prt) {return 0._prt;}
37 
38  // Fits also use energy in MeV
39  const amrex::ParticleReal E_MeV = E_keV*1.e-3_prt;
40  constexpr amrex::ParticleReal joule_to_MeV = 1.e-6_prt/PhysConst::q_e;
41 
42  // Compute Gamow factor, in MeV
43  constexpr auto one_pr = 1._prt;
44  constexpr auto Z_boron = 5._prt;
45  constexpr amrex::ParticleReal m_boron = 11.00930536_prt * PhysConst::m_u;
46  constexpr amrex::ParticleReal m_hydrogen = 1.00782503223 * PhysConst::m_u;
47  constexpr amrex::ParticleReal m_reduced = m_boron / (one_pr + m_boron/m_hydrogen);
48  constexpr amrex::ParticleReal gamow_factor = m_reduced / 2._prt *
49  (PhysConst::q_e*PhysConst::q_e * Z_boron /
51  (PhysConst::q_e*PhysConst::q_e * Z_boron /
53  joule_to_MeV;
54 
55  // Compute astrophysical factor, in MeV barn, using the fits
56  constexpr auto E_lim1 = 400._prt; // Limits between the different fit regions
57  constexpr auto E_lim2 = 642._prt;
58  amrex::ParticleReal astrophysical_factor;
59  if (E_keV < E_lim1)
60  {
61  constexpr auto C0 = 197._prt;
62  constexpr auto C1 = 0.24_prt;
63  constexpr auto C2 = 2.31e-4_prt;
64  constexpr auto AL = 1.82e4_prt;
65  constexpr auto EL = 148._prt;
66  constexpr auto dEL_sq = 2.35_prt*2.35_prt;
67  astrophysical_factor = C0 + C1*E_keV + C2*powi<2>(E_keV) +
68  AL/((E_keV - EL)*(E_keV - EL) + dEL_sq);
69  }
70  else if (E_keV < E_lim2)
71  {
72  constexpr auto D0 = 330._prt;
73  constexpr auto D1 = 66.1_prt;
74  constexpr auto D2 = -20.3_prt;
75  constexpr auto D5 = -1.58_prt;
76  const amrex::ParticleReal E_norm = (E_keV-400._prt) * 1.e-2_prt;
77  astrophysical_factor = D0 + D1*E_norm + D2*powi<2>(E_norm) + D5*powi<5>(E_norm);
78  }
79  else
80  {
81  constexpr auto A0 = 2.57e6_prt;
82  constexpr auto A1 = 5.67e5_prt;
83  constexpr auto A2 = 1.34e5_prt;
84  constexpr auto A3 = 5.68e5_prt;
85  constexpr auto E0 = 581.3_prt;
86  constexpr auto E1 = 1083._prt;
87  constexpr auto E2 = 2405._prt;
88  constexpr auto E3 = 3344._prt;
89  constexpr auto dE0_sq = 85.7_prt*85.7_prt;
90  constexpr auto dE1_sq = 234._prt*234._prt;
91  constexpr auto dE2_sq = 138._prt*138._prt;
92  constexpr auto dE3_sq = 309._prt*309._prt;
93  constexpr auto B = 4.38_prt;
94  astrophysical_factor = A0 / ((E_keV-E0)*(E_keV-E0) + dE0_sq) +
95  A1 / ((E_keV-E1)*(E_keV-E1) + dE1_sq) +
96  A2 / ((E_keV-E2)*(E_keV-E2) + dE2_sq) +
97  A3 / ((E_keV-E3)*(E_keV-E3) + dE3_sq) + B;
98  }
99 
100  // Compute cross section, in barn
101  return astrophysical_factor/E_MeV*std::exp(-std::sqrt(gamow_factor/E_MeV));
102 }
103 
114 amrex::ParticleReal ProtonBoronFusionCrossSectionBuck (const amrex::ParticleReal& E_keV)
115 {
116  using namespace amrex::literals;
117 
118  constexpr amrex::ParticleReal E_start_fit = 3500._prt; // Fit starts at 3.5 MeV
119  // cross section at E = E_start_fit, in barn
120  constexpr amrex::ParticleReal cross_section_start_fit = 0.2168440845211521_prt;
121  constexpr amrex::ParticleReal slope_fit = -2.661840717596765_prt;
122 
123  // Compute fitted value
124  return cross_section_start_fit*std::pow(E_keV/E_start_fit, slope_fit);
125 }
126 
139 amrex::ParticleReal ProtonBoronFusionCrossSection (const amrex::ParticleReal& E_kin_star)
140 {
141  using namespace amrex::literals;
142 
143  // Fits use energy in keV
144  constexpr amrex::ParticleReal joule_to_keV = 1.e-3_prt/PhysConst::q_e;
145  const amrex::ParticleReal E_keV = E_kin_star*joule_to_keV;
146  constexpr amrex::ParticleReal E_threshold = 3500._prt;
147 
148  const amrex::ParticleReal cross_section_b = (E_keV <= E_threshold) ?
151 
152  // Convert cross section to SI units: barn to square meter
153  constexpr auto barn_to_sqm = 1.e-28_prt;
154  return cross_section_b*barn_to_sqm;
155 }
156 
157 #endif // WARPX_PROTON_BORON_FUSION_CROSS_SECTION_H
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::ParticleReal ProtonBoronFusionCrossSectionNevins(const amrex::ParticleReal &E_keV)
Computes the total proton-boron fusion cross section in the range 0 < E < 3.5 MeV using the analytica...
Definition: ProtonBoronFusionCrossSection.H:30
AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::ParticleReal ProtonBoronFusionCrossSection(const amrex::ParticleReal &E_kin_star)
Computes the total proton-boron fusion cross section. When E_kin_star < 3.5 MeV, we use the analytica...
Definition: ProtonBoronFusionCrossSection.H:139
AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::ParticleReal ProtonBoronFusionCrossSectionBuck(const amrex::ParticleReal &E_keV)
Computes the total proton-boron fusion cross section in the range E > 3.5 MeV using a simple power la...
Definition: ProtonBoronFusionCrossSection.H:114
static constexpr auto ep0
vacuum permittivity: dielectric permittivity of vacuum [F/m]
Definition: constant.H:46
static constexpr auto hbar
reduced Planck Constant = h / tau [J*s]
Definition: constant.H:59
static constexpr auto m_u
dalton: unified atomic mass unit [kg]
Definition: constant.H:56
static constexpr auto q_e
elementary charge [C]
Definition: constant.H:50