INVESTIGATION OF THE TRIPLE-ALPHA REACTION IN A FULL THREE-BODY
APPROACH
By
Ngoc Bich Nguyen

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics − Doctor of Physiology
2013

ABSTRACT
INVESTIGATION OF THE TRIPLE-ALPHA REACTION IN A FULL
THREE-BODY APPROACH
By
Ngoc Bich Nguyen
We have developed a new three-body method to compute the triple-alpha reaction rate
at low temperatures where measurements are impossible and many numerical attempts have
failed before. In this work, the triple-alpha is modeled as a three-body Borromean system
in hyperspherical harmonics coordinates. In the low temperature region, the triple-alpha
proceeds through a quadrupole transition from the 0+ continuum to the 2+ bound state
1
in 12 C. The 2+ bound state is obtained by solving a set of coupled channels equations in
1
hyper-radius coordinates for negative energy with a boundary condition that requires the
wavefunction to go to zero at large distances. The same approach can not be applied to the
0+ continuum state because it requires an exact boundary condition for the three charged
particles. We therefore combine the R-matrix expansion, the R-matrix propagation method,
and the screening technique in the hyperspherical harmonics basis to obtain a numerically
stable three-body continuum wavefunction for the α + α + α system. We employ the AliBodmer potential for the alpha-alpha interaction which reproduces the low energy phase
shifts as well as the 0+ resonance of 8 Be. We add a three-body force to fit experimental
data. Both the 2+ bound state and the 0+ resonant state in 12 C are well reproduced in our
1
2
framework. We find a dominant triangle three-alpha configuration for the Hoyle resonance
by studying the density distribution function. The resonant and non-resonant continuum
states of 12 C(0+ ) are obtained simultaneously, allowing us to include these two processes on
the same footing. Long range Coulomb interactions show important effects especially in the

low temperature regime. We also present a detailed convergence study of the triple-alpha
reaction rate with respect to the screening radius and the size of the model space. The
new rate agrees with the NACRE rate for temperatures greater than 0.07 GK, but a large
enhancement at lower temperatures is found (≈ 1012 at 0.02 GK). We observe a transition
between dominance of the non-resonant versus the resonant triple-alpha capture process
around 0.06 GK. Our results are then compared to previous calculations where additional
approximations were made. We also explore the impact in astrophysics namely, in stellar
evolution, helium accreting white dwarfs, and helium accreting neutron stars.
The presence of very narrow two-body and three-body resonances in addition to the
strong, long-range Coulomb interaction makes the triple-alpha problem very challenging.
Our framework which is the combination of various methods is a new approach to overcome
the well known difficulty of the three charged-particle system. This method is particularly
suited to the very low energy regime where measurements are impossible. In addition, the
method we developed opens new opportunities in addressing three-body, low-energy reactions
in other fields (e.g., atomic and molecular physics) where all three particles have charge.

ACKNOWLEDGMENTS

I am heartily thankful to my supervisor, Prof. Filomena Nunes, whose support and guidance have had great impact on both my academic and personal life. I am grateful for her
enthusiasm and patience in the beginning of my journey with nuclear reaction theory, her
positive attitude and encouragement throughout the years, and her kindness and the care
she gave me during difficult times of my life with health issues. I feel very fortunate to have
her as my advisor.
I would like to express my gratitude to Prof. Ian Thompson for his expertise on the
project and his tremendous help with developing the computer codes. I would like to thank
Prof. Edward Brown for many helpful discussions on the astrophysics aspect of the project.
Special thanks to my guidance committee members Prof. Pawel Danielewicz, Prof. Hendrik
Schatz, Prof. Bhanu Mahanti, Prof. Edward Brown for their comments and helpful suggestions. Their guidance has been invaluable and I owe them my deepest appreciation. I would
also like to thank Prof. Vladimir Zelevinsky, who have led me to Nuclear Physics, for his
inspiration and excellent lectures. A sincere thank to Prof. Ron Johnson from whom I have
learned a lot about reaction theory during his visits.
The Department of Physics and Astronomy at Michigan State University, the National
Superconducting Cyclotron Laboratory, and especially the theory group are sincerely acknowledged for providing the financial, academic and technical support, and a great working
environment. I would also like to acknowledge the National Science Foundation and the
Department of Energy for this support.
I thank my friends and colleagues from Vietnam, United States and other parts of the
world for many belated memories we have shared and their company, emotional support
iv

during my long journey. Special thank to my first and current office mate Brent Barker for
many sincere conversations, to my dearest friends Seungyun Yook and Winifred Nwaefido
for their kindness and great friendship. I would like to thank my group members Neelam
Upadhyay, Muslema Pervin and Luke Titus for many good memories and fun time we have
together. I also want to thank the Association of Vietnamese Students & Scholars for their
strong support and help at the very beginning of my time at Michigan State University.
Last but not least, I am greatly indebted to my family, my parents Chin Nguyen and
Duc Nguyen, my sister Nga Nguyen for their unconditional love and support. Without them
I would not have been able to make it to this point of my life. My heartfelt thanks are
addressed to my husband David Lincoln and his family, especially Lou and Diane Lincoln
whose love and encouragement have guided me stepping through all the obstacles of my
journey. I hope that this work makes you all proud.

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 History of the triple-alpha reaction . . . . . . . . . . . . . . . . .
1.2 The triple-alpha reaction rate and its current status . . . . . . . .
1.3 Three-body reaction methods in solving the triple-alpha problem
1.4 Motivation for this work . . . . . . . . . . . . . . . . . . . . . . .
1.5 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.

1
1
4
11
13
13

Chapter 2 The Hyperspherical Harmonics R-matrix (HHR) method . . .
2.1 Hyperspherical Harmonics method . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Hyperspherical coordinates . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Three-body bound states . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2.1 Asymptotic behavior of the three-body Coulomb interaction
2.1.3 Three-body continuum states . . . . . . . . . . . . . . . . . . . . . .
2.1.4 Properties of the three-alpha wavefunction . . . . . . . . . . . . . . .
2.2 R-matrix method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 R-matrix expansion . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 R-matrix propagation . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Screening technique . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 HHR3a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15
15
16
18
21
24
26
29
30
33
39
42

Chapter 3 Three-body reaction rate . . . . . . . . . . . . . . .
3.1 Mawell-Boltzmann energy distribution for a three-body system
3.2 The reaction rate formalism for three-body capture . . . . . .
3.3 Three-body quadrupole transition function . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

45
46
48
52

Chapter 4 Results . . . . . . . . . . . . . . . . . . .
4.1 Interactions . . . . . . . . . . . . . . . . . . . .
4.2 The 2+ bound state . . . . . . . . . . . . . . .
1
4.3 The 0+ Hoyle state . . . . . . . . . . . . . . .
2
4.4 Testing the new code HHR3a . . . . . . . . . .
4.5 Convergence and uncertainty . . . . . . . . . . .
4.5.1 The impact of using screening potentials
4.5.2 Convergence study . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

58
58
62
67
72
74
75
78

vi

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

4.6

4.5.3 Interaction uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.4 Rate uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 5 Discussion . . . . . . . .
5.1 Comparison with other methods
5.2 Long-range Coulomb effects . .
5.3 Reaction dynamics . . . . . . .
5.4 Astrophysical S-factor . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

83
84
87

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

91
91
94
96
99

Chapter 6 Astrophysics Applications . .
6.1 Evolution of single stars . . . . . . . .
6.2 Binary stars . . . . . . . . . . . . . . .
6.2.1 Helium accreting white dwarfs .
6.2.2 Helium accreting neutron stars

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

104
106
113
114
118

Chapter 7 Summary and outlook . . . . . . . . . . . . . . . . . . . . . . . . . 121
7.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
7.2 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Appendix A: Detailed Derivations of the R-matrix propagation method . . . . . .
A.1
Deriving Eq. (2.58) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2
Deriving Eq. (2.61) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3
Deriving Eq. (2.63), Eq. (2.64), Eq. (2.65), and Eq. (2.66) . . . . . .
Appendix B: Detailed derivation of the quadrupole transition function . . . . . . .
B.1
Deriving Eq. (3.32) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.2
Deriving Eq. (3.42) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Appendix C: Inputs for numerical calculations . . . . . . . . . . . . . . . . . . . .
C.1
Bound state calculation . . . . . . . . . . . . . . . . . . . . . . . . .
C.2
Continuum state calculation . . . . . . . . . . . . . . . . . . . . . . .
Appendix D: Highlights on the (d,p) reaction projects . . . . . . . . . . . . . . . .
D.1
Common physical quantities in (d,p) studies . . . . . . . . . . . . . .
D.2
The adiabatic distorted wave approximation . . . . . . . . . . . . . .
D.2.1
Zero-range adiabatic wave approximation (ZR-ADWA) . . .
D.2.2
Finite-range adiabatic wave approximation (FR-ADWA) . .
D.3
Finite-range effects in (d,p) reactions . . . . . . . . . . . . . . . . . .
D.3.1
Finite-range effects in the potentials . . . . . . . . . . . . .
D.3.2
Finite-range effects in the transfer cross sections . . . . . . .
D.4
Testing the dispersive optical model for (d,p) reactions . . . . . . . .
D.4.1
Details of the calculation . . . . . . . . . . . . . . . . . . . .
D.4.2
Transfer cross section . . . . . . . . . . . . . . . . . . . . .
D.4.3
Spectroscopic factors . . . . . . . . . . . . . . . . . . . . . .
D.5
Application of the finite-range adiabatic distorted wave approximation
(FR-ADWA) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
vii

127
128
128
131
132
135
135
138
141
141
143
145
147
150
152
153
156
157
159
167
168
170
175
178

D.6

Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180

Bibliography . . . . . . . . . . . . . . .

viii

. . . . . . . . . . . . . . . . . . . 183

LIST OF TABLES

Table 4.1

The three-body interactions to reproduce the experimental energies
[69] of the 2+ bound state and 0+ resonant state in 12 C. . . . . . .
1
2

61

Comparison of the three-body binding energy and the rms radius of
the 2+ bound state of 12 C obtained from our three-body calculation
1
with other studies: a three-body hyperspherical adiabatic method
[26], a microscopic Fermionic Molecular Dynamics method [31] and
experimental data [69]. . . . . . . . . . . . . . . . . . . . . . . . . .

64

Table 4.3

Sources of uncertainty in calculating the triple-alpha reaction rate .

86

Table 4.4

Our triple-alpha reaction rate (in cm6 s−1 mol−2 ) after being normalized to NACRE for T ≥ 0.5 GK (the factor of 2 uncertainty is listed
for temperature below 0.08 GK above which the triple-alpha rate
agrees with NACRE). . . . . . . . . . . . . . . . . . . . . . . . . . .

88

Table 4.2

Table 6.1

Temperature sensitivity of triple-alpha rate . . . . . . . . . . . . . . 105

Table D.1

The finite-range effects on the deuteron distorting potential are
sented as the changes in the diffuseness aR and aI of the real
imaginary parts, the depths of the real and imaginary parts V
Ws , as well as the corresponding radii rR and rI . . . . . . . .

Table D.2

Finite-range effects in transfer cross sections. Comparisions are performed with respect to the zero-range Johnson and Soper calculation
(cases with no existing data are presented by *). . . . . . . . . . . 160

Table D.3

Comparison of the DOM overlap function, corrected for non-locality
ϕDOM with the Woods-Saxon single-particle wavefunction ϕWS . In
this table we include: n (principal number), j (the angular momenta
of the valence orbital), Sn (the separation energy), Rrms (the rootmean-square radius of the valence orbital), and bnlj (the modulus of
the single-particle ANC). . . . . . . . . . . . . . . . . . . . . . . . . 169

Table D.4

Comparision of the spectroscopic factors obtained from the FR-ADWA
analysis using different models for optical potentials and overlap functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
ix

preand
and
. . . 158

Table D.5

ANCs extracted for both the ground state and the first excited state
of the 14 C(d,p)15 C reaction. . . . . . . . . . . . . . . . . . . . . . . 178

x

LIST OF FIGURES

Figure 1.1

The triple-alpha reaction mechanism: The energy-level diagram of
12 C and 8 Be relative to the three-alpha threshold are presented. Two
alpha particles combine to form 8 Be which can then capture another
alpha particle to produce 12 C at the 0+ excited state (or the Hoyle
2
resonance). Here, 12 C is likely to decay back into three alpha particles again but it can also decay to the ground state of 12 C through
two successive gamma ray emissions or the emission of one e+ e− pair.
“For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation.”

3

Figure 1.2

Triple-alpha Jacobi coordinates in CDCC study. . . . . . . . . . . .

6

Figure 2.1

Three ways to form a Jacobi coordinate set for a system of three
particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

Figure 2.2

Coordinates defined for Jacobi set i . . . . . . . . . . . . . . . . . .

16

Figure 2.3

Hyper-radial behavior of the potential couplings Vγγ (ρ) in Eq. (2.17).
The black solid curve represents the diagonal coupling and the others
demonstrate the off-diagonal couplings. . . . . . . . . . . . . . . . .

24

The hyper-radial scattering wavefunction corresponding to the incomi i
ing channel of K i lx ly = 000 and the outgoing channel of Klx ly = 000
in the hyperspherical expansion (Eq. (2.28)) is calculated at E = 0.2
MeV for two different R-matrix box sizes: 50 fm (red) and 100 fm
(blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

Diagram presenting the sectorization in hyper-radius ρ of the Rmatrix propagation method. The subscripts L and R indicate the
evaluations at the left and right sector boundary. The R matrix R
of a sector is calculated at the right boundary. . . . . . . . . . . . .

34

Schematic diagram of the propagation method in obtaining the threebody scattering wavefunction . . . . . . . . . . . . . . . . . . . . . .

37

Figure 2.4

Figure 2.5

Figure 2.6

xi

Figure 2.7

Figure 2.8

The quadrupole transition strength function dB(E2)/dE is calculated using the global (solid) and local (dashed) representations of
the wavefunction propagation method. . . . . . . . . . . . . . . . . .

39

The hyper-radial scattering wavefunction corresponding to the incomi i
ing channel of K i lx ly = 000 and the outgoing channel of Klx ly = 000
in the hyperspherical expansion (Eq. (2.28)) is calculated at E = 0.5
MeV for both cases: with screening (red solid) and without screening
(blue dashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

Figure 2.9

Relative deviation of the quadrupole transition strength dB(E2)/dE
(denoted X in the plot) when diagonal Coulomb screening is introduced. 41

Figure 2.10

Our three-body method is divided into four steps in which we employ the R-matrix expansion, R-matrix propagation and screening
technique in Hyperspherical Harmonics coordinates. . . . . . . . . .

42

Dependence of the three-body binding energy of the 2+ bound state
1
12 C on the depth of a three-body force V (ρ) = V e−ρ2 /ρ2 . The
0
of
0
3b
maximum hyper-angular momentum Kmax is taken as 20 for which
the results are converged. . . . . . . . . . . . . . . . . . . . . . . . .

60

Figure 4.1

Figure 4.2

Dependence of the Hoyle resonant energy of 12 C on the depth of a
2

2

three-body force V3b (ρ) = V0 e−ρ /ρ0 . The maximum hyper-angular
momentum Kmax is taken as 26 for which the results are converged.

61

Dependence of the three-body binding energy of the 2+ bound state
1
of 12 C on the maximum hyper-angular momentum Kmax . . . . . . .

63

Dependence of the rrms radius of the 2+ bound state of 12 C on the
1
maximum hyper-angular momentum Kmax . . . . . . . . . . . . . . .

64

The rrms radius as a function of the maximum radius of the box in
which the 12 C(2+ ) bound state wavefunction is calculated. . . . . .
1

65

Figure 4.6

The density distribution for the 2+ bound state of 12 C. . . . . . . .
1

66

Figure 4.7

The dominant three-alpha cluster in the 12 C(2+ ) bound state: three
1
alpha particles touch each other and form a equilateral triangle. . .

66

Figure 4.3

Figure 4.4

Figure 4.5

xii

Figure 4.8

The quadrupole transition strength dB(E2)/dE as a function the
three-alpha relative energy E. . . . . . . . . . . . . . . . . . . . . .

68

The quadrupole transition strength function dB(E2)/dE around the
Hoyle resonant energy is calculated with different values of hypermomentum Kmax . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

Dependence of the three-body energy of the 0+ Hoyle state of 12 C on
2
the maximum hyper-angular momentum Kmax . . . . . . . . . . . . .

70

Figure 4.11

The density distribution for the 0+ resonance state of 12 C. . . . . .
2

71

Figure 4.12

Three triple-alpha configurations of the 0+ resonance in 12 C: the
2
prolate triangle (a), the oblate triangle (b) and the equilateral triangle
(c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

A numerical solution (dashed) of Eq. (2.30) at E = 0.5 MeV for the
0+ scattering wavefunction of the triple-alpha system when only the
diagonal Coulomb couplings are taken into account is compared with
the analytic solution (solid). . . . . . . . . . . . . . . . . . . . . . .

73

The triple-alpha reaction rate is computed for the diagonal Coulomb
case using the numerical solution of Fig. 4.13 (dashed) and compared
to the analytic solution (solid). . . . . . . . . . . . . . . . . . . . . .

74

Figure 4.9

Figure 4.10

Figure 4.13

Figure 4.14

Figure 4.15

The reaction rate as a function of temperature for the fictitious threealpha system in which only Coulomb interactions are included: Kmax =
20 (solid) and Kmax = 26 (dot-dashed). . . . . . . . . . . . . . . . . 75

Figure 4.16

The quadrupole transition strength dB(E2)/dE at E = 0.01 MeV
as a function of model space Kmax when both nuclear and Coulomb
interaction (a) and only Coulomb interaction (b) are present in the
three-alpha system. . . . . . . . . . . . . . . . . . . . . . . . . . . .

76

Relative differences between screening and no screening calculations
for the fictitious three-alpha scattering system in which only Coulomb
interactions are introduced. The screening calculation is performed
for Kmax = 26 and ρscreen = 800 fm. Calculations without screening
are performed for Kmax = 20 (solid) and Kmax = 26 (dashed). . . .

77

Figure 4.17

xiii

Figure 4.18

The quadrupole transition strength are computed at low energies using three different values of the diffuseness in screening potentials:
ascreen = 10 fm (solid), ascreen = 20 fm (dashed), ascreen = 30 fm
(dot-dashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

Figure 4.19

Relative differences between the triple-alpha rate calculated at several
screening radii and its converged value. The reaction rates calculated
using ρscreen = 600 fm (solid), ρscreen = 800 fm (dashed), ρscreen =
1000 fm (dot-dashed) are compared with the converged rate in Table 4.4. 79

Figure 4.20

Relative differences between the triple-alpha rate calculated at several
maximum hyper-momentum Kmax and its converged value. The reaction rates calculated using Kmax = 20 (solid), Kmax = 22 (dashed),
Kmax = 24 (dot-dashed), Kmax = 26 (dot-dot-dashed) are compared
with the converged rate in Table 4.4. . . . . . . . . . . . . . . . . .

81

Convergence of the quadrupole transition strength dB(E2)/dE with
hyper-momentum Kmax for E = 0.01 MeV. . . . . . . . . . . . . . .

81

Convergence of the quadrupole transition strength dB(E2)/dE with
njac, number of Gauss-Jacobi quadrature grid point, for a given low
energy of 0.01 MeV . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

Convergence of the quadrupole transition strength dB(E2)/dE with
matching radius ρa for a given low energy of 0.01 MeV . . . . . . .

83

The sensitivity of the triple-alpha rate to the interactions: comparison
between α−α Ali-Bodmer+3-body force (solid) and α−α interaction
as in [15]+3-body force (dashed) . . . . . . . . . . . . . . . . . . . .

85

The photo-disintegration cross section is fitted to the Breit-Wigner
shape to obtained the gamma decay width. . . . . . . . . . . . . . .

87

The integrand f (E) and its energy cumulative sum Sm(E) at a very
low temperature of 0.01 GK . . . . . . . . . . . . . . . . . . . . . .

89

The integrand f (E) and its energy cumulative sum Sm(E) at T = 0.1
GK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89

Range of energy relevant for the triple-alpha reaction rate at a given
temperature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

Figure 4.21

Figure 4.22

Figure 4.23

Figure 4.24

Figure 4.25

Figure 4.26

Figure 4.27

Figure 4.28

xiv

Figure 5.1

Different evaluations of the triple-alpha reaction rate: comparing
the Hyperspherical Harmonic R-matrix method (solid) with NACRE
(dotted), CDCC (dashed) and the three-body Breit Wigner (dotdashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

Figure 5.2

The triple-alpha reaction rate: comparing the Hyperspherical Harmonic R-matrix method (solid) with NACRE (dotted), Fynbo (dashed),
and CF88 (dot-dashed). . . . . . . . . . . . . . . . . . . . . . . . . . 94

Figure 5.3

The long-range Coulomb effects are shown in four different calculations: only diagonal Coulomb couplings (dotted), only Coulomb
couplings (dot-dashed), both nuclear and Coulomb interactions with
off-diagonal Coulomb couplings up to 30 fm (dashed) and a fully converged calculation with off-diagonal Coulomb couplings up to 800 fm
(solid) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

Spatial distribution of function fγ i (r, R) for E = 0.05 MeV and γ i =
{0, 0, 0}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

Spatial distribution of function fγ i (r, R) for E = 0.38 MeV and γ i =
{0, 0, 0}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

Figure 5.4

Figure 5.5

Figure 5.6

Spatial distribution of function fγ i (r, R) for E = 0.5 MeV γ i = {0, 0, 0}. 98

Figure 5.7

2
Fitting Eγ σγ to an exponential behavior A e E for energies below
the resonance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

Figure 5.8

Comparison between the energy dependence of the cross section σγ
and the three-body astrophysical S-factor S3b . . . . . . . . . . . . . 101

Figure 5.9

Comparison between the energy dependence of the cross section σγ
and the astrophysical S-factor S0 for Kmax = 0. . . . . . . . . . . . 102

Figure 6.1

Temperature dependence for the HHR rate (solid) and the NACRE
rate (dashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Figure 6.2

Standard evolution curve of 1M stars (surface luminosity as a function of surface temperature): 1 - main sequence, 2 - red giant, 3 core helium flashes, 4 - thermal pulses, 5 - planetary nebula, and 6 white dwarf. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

−2πB
√

xv

Figure 6.3

Full evolutionary track (luminosity vs. surface temperature) from
main sequence to white dwarf of a one solar-mass (1M ) star for HHR
(solid), NACRE (dashed), Fynbo (dotted), and CF88 (dot-dashed).

110

Figure 6.4

Evolutionary track (luminosity vs. surface temperature) from main
sequence to the end of asymptotic giant branch of a one solar-mass
(1M ) star for HHR (solid), NACRE (dashed), Fynbo (dotted), and
CF88 (dot-dashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

Figure 6.5

Surface luminosity as a function of age: HHR rate (solid) and NACRE
rate (dashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

Figure 6.6

Full evolutionary track (luminosity vs. surface temperature) from
main sequence to white dwarf of a 0.8 solar-mass (0.8M ) star for
the HHR rate (solid) and the NACRE rate (dashed). . . . . . . . . . 113

Figure 6.7

Binary system and its Roche lobes . . . . . . . . . . . . . . . . . . . 113

Figure 6.8

Helium ignition curves for the NACRE rate (blue solid) and the HHR
rate (red dashed). Calculations are performed for helium accreting
white dwarf with initial mass of 1M and several values of accretion
rate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

Figure 6.9

Helium ignition curves for the NACRE rate (blue solid) and the HHR
rate (red dashed). Calculations are performed for helium accreting
white dwarf with the accretion rate of 5 × 10−10 M yr−1 and two
different initial masses: 0.734M and 1.315M . . . . . . . . . . . . 117

Figure 6.10

The light curve (luminosity vs. time) for helium accreting a neutron
star at accretion rate of 1.75×10−10 M yr−1 : Comparing the CF88
rate (dashed) and the HHR rate (solid). . . . . . . . . . . . . . . . . 119

Figure D.1

(d,p) transfer reaction: (a) incoming channel, (b) outgoing channel . 151

Figure D.2

Finite-range effects in sub-Coulomb (d,p) reactions: (a) 48 Ca(d,p)49 Ca(g.s.)
Ed = 2 MeV (data from [124]) and (b) 208 Pb(d,p)209 Pb(g.s.) Ed = 8
MeV (data from [125]). . . . . . . . . . . . . . . . . . . . . . . . . 163

Figure D.3

Finite-range effects in (d,p) reactions at energies slightly above the
Coulomb barrier: (a) 86 Kr(d,p)87 Kr(g.s.) Ed = 11 MeV (data from
[126]) and (b) 208 Pb(d,p)209 Pb(g.s.) Ed = 20 MeV (data from [127]). 164

xvi

Figure D.4

Finite-range effects in (d,p) reactions at high energies: (a) 12 C(d,p)13 C(g.s.)
Ed = 56 MeV (data from [128]) and (b) 48 Ca(d,p)49 Ca(g.s.) Ed = 56
MeV (data from [129]). . . . . . . . . . . . . . . . . . . . . . . . . 165

Figure D.5

Systematic finite-range effect as a function of beam energy. The effect
from the adiabatic potentials is presented as open symbols while that
from the evaluation of the matrix elements is plotted as filled symbols.166

Figure D.6

The square of the single-neutron overlap functions in the WoodsSaxon (solid) and the DOM (dashed) model are compared for 49 Ca.

170

Figure D.7

DOM and CH89 optical potentials are compared for n-132 Sn and p132 Sn at 4.7 MeV: (a) the real part and (b) the imaginary part. . . . 171

Figure D.8

48 Ca(d,p)49 Ca

Figure D.9

48 Ca(d,p)49 Ca

at Ed = 2 MeV (angular distributions are normalized
to the data at backward angles). . . . . . . . . . . . . . . . . . . . . 172
at Ed = 19.3 MeV (angular distributions are normalized to the data at the peak). . . . . . . . . . . . . . . . . . . . . . . 173

Figure D.10 48 Ca(d,p)49 Ca at Ed = 56 MeV (angular distributions are normalized
to the data at forward angles). . . . . . . . . . . . . . . . . . . . . . 173
Figure D.11 132 Sn(d,p)133 Sn at Ed = 9.46 MeV (angular distributions are normalized at the peak of the experimental cross section). . . . . . . . 174
Figure D.12 208 Pb(d,p)209 Pb at Ed = 20 MeV (angular distributions are normalized at the peak of the experimental data). . . . . . . . . . . . . . . 174
Figure D.13 208 Pb(d,p)209 Pb at Ed = 8 MeV (angular distributions are normalized at the peak of the experimental data). . . . . . . . . . . . . . . 177
Figure D.14 Angular distribution for the ground state in the 14 C(d,p)15 C reaction. 179

xvii

Chapter 1
Introduction

1.1

History of the triple-alpha reaction

In the early 1950s, physicists were puzzled by the mystery of nucleosynthesis beyond mass
8 since stable isotopes of mass number 5 and 8 can not be found in nature. To explain the
puzzle, Salpeter proposed the theory that stable 12 C is formed by alpha capture on 8 Be in
equilibrium, and heavier nuclei are then synthesized from 12 C [1]. In his theory, two alpha
particles first fuse to create an intermediate particle 8 Be which is unstable by only 92 keV
against alpha decay. The ground state of 8 Be has a very short lifetime of 2.6 × 10−16 sec
which corresponds to a narrow resonance with a width of 2.5 eV [2]. Although short, this
lifetime is sufficient for a small amount of 8 Be in equilibrium to build up. For every 109
alpha particles, one finds approximately 1 8 Be nucleus [3]. This modest concentration of 8 Be
can then capture another alpha particle to form 12 C. However, the corresponding reaction
rate was too small to account for the observed abundances of 12 C in stars. Because both 4 He
and 8 Be nuclei have total spins J π = 0+ , the alpha capture on 8 Be would highly increase its
rate if the reaction could proceed through a 12 C s-wave resonance. Based on that argument,
1

Fred Hoyle predicted the existence of a 0+ resonance in the 12 C energy spectrum, near the
three-alpha threshold, to significantly enhance the reaction rate [4]. This resonant state and
its properties were confirmed by experimentalists [5, 6] shortly thereafter. The resonance was
then named after Hoyle and located as a second excited state of 12 C at 0.38 MeV above the
three-alpha threshold. Fig. 1.1 illustrates the energy-level diagram of 12 C and 8 Be relative to
the three-alpha threshold. The Hoyle resonant state may sequentially decay back into three
alpha particles through the intermediate step of the 8 Be ground state. It may also deexcite
to the ground state of 12 C through two successive gamma ray emissions. While direct decay
to the ground state of 12 C through a photon emission is forbidden by angular momentum
conservation, an emission of a e+ e− pair is possible. The probability for the photon decay
from the Hoyle resonant state to the ground state is much smaller than the probability of
three-alpha breakup but significantly larger than the e+ e− pair decay [6]. A large amount
of 12 C could therefore be produced in stars through this mechanism.
Thus, the two-step (sequential) mechanism of the triple-alpha capture,

α + α → 8 Be(0+ )
1
8 Be(0+ ) + α
1

→ 12 C(0+ ) ,
2

was accepted to explain the mystery of nucleosynthesis and the abundances of 12 C in helium
burning stars. The capture probability is highly enhanced by the narrow resonances of 8 Be
and 12 C. The accuracy of this helium burning rate has a big impact on the late stages of
stellar evolution. It was shown that there is a strong dependence between the triple-alpha
reaction rate and the production of nuclei 26 Al, 44 Ti, 60 Fe in supernova explosions [7]. Over
the range of twice its experimental uncertainty (∼ 12%), one finds that the production of

2

Figure 1.1: The triple-alpha reaction mechanism: The energy-level diagram of 12 C and 8 Be
relative to the three-alpha threshold are presented. Two alpha particles combine to form
8 Be which can then capture another alpha particle to produce 12 C at the 0+ excited state
2
(or the Hoyle resonance). Here, 12 C is likely to decay back into three alpha particles again
but it can also decay to the ground state of 12 C through two successive gamma ray emissions
or the emission of one e+ e− pair. “For interpretation of the references to color in this and
all other figures, the reader is referred to the electronic version of this dissertation.”
60 Fe, 26 Al

and 44 Ti varies by a factor of 7.8, 2.8 and 1.6 respectively [7]. The uncertainty in

the triple-alpha reaction rate also leads to an ambiguity in determining the Fe core in core
collapse supernovae. For an uncertainty of 10% in C/O ratio, there would be an intrinsic
uncertainty of 0.2M in the Fe core mass for a 25M star [7]. The uncertainty in determining
the triple-alpha rate affects the predictive power of stellar evolution models. Because of its
important impact in astrophysics as well as astronomy, many efforts have been done to better
understand this reaction.

3

1.2

The triple-alpha reaction rate and its current status

In order for the triple-alpha reaction to happen as a two step process, the stellar environment
temperature must be high enough so the alpha particles can reach the resonant energies. At
low temperatures, we expect the three-alpha direct capture to become dominant. Even
though the current experimental uncertainty of the triple-alpha resonant rate is 12% [7], the
theory of this process is well established. The first evaluation of the triple-alpha resonant
reaction rate was credited to Caughlan and Fowler [8] in which a Breit-Wigner shape [9] for
the two-body resonant cross sections of α-α and α-8 Be was used:

Γi (E)Γf (E)
π (2J + 1)(1 + δ12 )
.
σBW (E) = 2
k (2j1 + 1)(2j2 + 1) (E − Er )2 + Γ(E)2 /4

(1.1)

In Eq. (1.1), the momentum k and the energy E of the two-body system are related by
k=

√
2µE/ , where µ is the reduced mass of the two particles and

is the Planck constant.

ji are the spins of projectile and target, and J is the spin of the resonance. The Kronecker
delta δ12 accounts for identical particles in the entrance channel. Γi (Γf ) is the resonance
partial width of the entrance (exit) channel, Γ is the total resonance width, and Er is the
resonant energy. The triple-alpha reaction rate was calculated in two steps, through the
formation of the intermediate 8 Be(0+ ) resonance. One first evaluates the reaction rate for
1
alpha capture by 8 Be by integrating the α-8 Be cross section over its corresponding MaxwellBoltzmann distribution:

8
8π
NA σv α Be = NA 2
µ 8

α Be

µα8 Be

3/2

2πkB T

∞
0

4

E

−
σα8 Be (E ; E) e kB T E dE ,

(1.2)

where E and E are the energies of α-α and α-8 Be systems, respectively; kB is the Boltzmann
constant and T is the temperature at which the triple-alpha reaction rate is evaluated. The
triple-alpha reaction rate is then calculated by the Maxwell-Boltzmann averaging of the α8 Be

capture rate weighted with the α-α fusion cross section which is proportional to the

probability of the 8 Be formation:

8π
2
NA σv ααα = 3NA 2
µαα

3/2 ∞
− E
8
µαα
σαα (E)
e kB T NA σv α Be E dE . (1.3)
8
2πkB T
0 Γα ( Be, E)

Eq. (1.3) clearly indicates the sequential mechanism of the triple-alpha reaction used in this
theory. It is important to note that Caughlan and Fowler considered only the Hoyle resonance when calculating the triple-alpha reaction rate. Ten years later, this reaction rate was
improved by Angulo [10] when taking into account the contribution of the 2+ resonance.
2
Astrophysicists refer to this new rate as the NACRE (Nuclear Astrophysics Compilation of
REaction Rate) reaction rate. The 12 C(2+ ) resonance contributes significantly to the triple2
alpha rate at T > 1 GK but has little impact at low temperatures [11]. Uncertainties in
the properties of this 2+ state have been recently addressed [12]. The work in [13] found a
2
0+ broad resonance around 3 MeV above the three-alpha threshold that does not influence
3
the rate at T ∼ 0.1 GK and is much less dominant than the 2+ resonance at T ∼ 1 GK.
2
These high lying resonances leads to variation of the triple-alpha rate in [8], [10] and [13].
Throughout this work we will use NACRE [10] as our reference, but comparisons with other
rates [8, 13] are presented in [14]. It is important to note that none of these resonances
contribute to the triple-alpha rate at low temperatures where the measurements are impossible. To determine the rate at T < 0.1 GK, a well founded theory for the three-alpha direct
(non-resonant) capture is crucial.

5

Figure 1.2: Triple-alpha Jacobi coordinates in CDCC study.
The non-resonant triple-alpha rate tabulated in NACRE is obtained by a simple extrapolation to low energies using the tail of the Breit-Wigner cross sections in the sequential
model as in Eq. (1.3) and Eq. (1.2). The validity of this approach is doubtful since at low
energy the three alpha particles do not have access to the intermediate resonances. The
question of how much the direct capture process contributes to the triple-alpha rate at low
temperature was not addressed in [10]. Ogata et al. [15] was the first group who attempted
to tackle this problem without using an ambiguous extrapolation. In order to describe the
three-alpha system, they used the Jacobi coordinates r and R as shown in Fig. 1.2, which
represent the relative radius between α1 and α2 and the distance from their center of mass
to α3 , respectively. In their studies, the continuum discretized coupled channel (CDCC)
method [16, 17] was employed to solve a three-body scattering equation:

[Tr + TR + v(r) + v(r1 ) + v(r2 ) − E]Ψ(r, R) = 0 ,

6

(1.4)

where E is the total relative energy of the three alpha system, Tr and TR are the kinetic
operators and v is the alpha-alpha potential. The continuum states of the α1 -α2 subsystem
were first discretized into momentum bins [ki , ki+1 ] with corresponding bin wavefunctions
ui (r):
ˆ
ui (r) =
ˆ

ki+1
ki

fi (k) u(k, r) dk/

ki+1
ki

1/2
2

|fi (k)| dk

,

(1.5)

where the two-body scattering wavefunction u(k, r) is averaged over the momentum bin
[ki , ki+1 ] using the weight function fi (k). The bin wavefunctions ui (r) form a discrete
ˆ
orthonormal basis

u∗ (r) uj (r) dr = δij . The total three-body wavefunction Ψ(r, R) was
ˆi
ˆ

then expanded in terms of the products between the bin wavefunction of the 8 Be continuum
states and its relative motion with respect to the third alpha particle:

+
Ψ0
ˆ ,E (r, R)
ki
0

=

(i )
imax
ˆ
ui (r) χi 0 (R)
ˆ
2 1
1
.
ˆ ˆ
π 32π 2 ki Ki
r
R
0 0 i=1

(1.6)

(i )
In Eq. (1.6) the bin wavefunction ui (r) describes the α1 -α2 continuum while χi 0 (R) deˆ
ˆ

scribes the motion of momentum bin ui (r) relative to α3 . The subscript i0 indicates the
ˆ
ˆ
incident channel, k is the average bin momentum of the α1 -α2 subsystem corresponding to
ˆ
the average bin energy ˆ12 , and K is the relative momentum between the continuum bin
and the third alpha. Using the above wavefunction decomposition technique, the three-body
Schrodinger equation becomes a set of coupled channels equations:

(i )
[TR + Vii (R) − (E − ε12,i )] χi 0 (R) = −
ˆ
ˆ

(i )
Vii (R) χ 0 (R) .
ˆ
i

i =i

7

(1.7)

The coupling potential Vii (R) is defined as:

Vii (R) =

u (r)
ˆ
ui (r)
ˆ
|v(r1 ) + v(r2 )| i
r.
r
r

(1.8)

In [15], the authors considered the quadrapole transition from the 0+ continuum states to
the 2+ bound state of 12 C to be the main contribution to the triple-alpha reaction rate.
1
They adopted a microscopic 2+ bound state wavefunction [18, 19] while numerically solving
1
Eq. (1.7) to obtain a description for the 12 C(0+ ) continuum states and from there constructed
the triple-alpha reaction rate. The advantage of this method was that both the resonant and
the non-resonant processes are treated on the same footing. However, the presence of the
long-range Coulomb potentials made it a challenging problem since the truncated CDCC
wavefunction generally does not have the correct asymptotic form for the scattering of three
charged particles [20]. The CDCC results differ from NACRE by 20 orders of magnitude at
T = 0.02 GK. Ogata et al. [15] explained that the Coulomb barrier between the two-alpha
subsystem and the third alpha particle is very much quenched at low energies, making the
non-resonant capture process dominant and enhancing the reaction rate at low temperatures. Such drastic enhancement of the triple-alpha reaction rate can significantly affect
astrophysics. For example, when using the CDCC rate to study the evolution of low and
intermediate mass stars, the red giant phase is strongly suppressed or disappears [21, 22],
which is incompatible with observation. Also, the studies of He ignition in accreting white
dwarfs [23] and accreting neutron stars [24, 25] show that the CDCC rate is barely consistent
with the observations of type Ia supernovae and type I X-ray bursts, respectively.
Another attempt to solve the triple-alpha problem was made by the Madrid-Aarhus collaboration [26, 27]. The idea of their work is to solve the Faddeev equations in coordinate

8

space using the adiabatic hyperspherical expansion method. The hyperspherical coordinates which consist of one radial coordinate ρ and five generalized angles denoted by Ω
were adopted in their studies (details on the hyperspherical coordinates will be presented in
chapter 2) for an adiabatic expansion of the wavefunction:

ΨJM =

1
ρ5/2

fn (ρ) ΦJM (ρ, Ω) .
n

(1.9)

n

Here, the expansion coefficient fn (ρ) is hyper-radial dependent. ΦJM (ρ, Ω) is the solution of
n
an eigenvalue problem with respect to the five dimensional spin-angular part of the Faddeev
operator when considering ρ to be invariant:

(TΩ − λn ) ΦJM +
n,(i)

2m 2
ρ Vi ΦJM = 0 ,
n
2

ΦJM = ΦJM + ΦJM + ΦJM .
n
n,(1)
n,(2)
n,(3)

(1.10)
(1.11)

The index i in ΦJM represents one of the three Jacobi coordinates in the Faddeev repren,(i)
sentation. TΩ is the angular part of the kinetic operator and Vi is the two-body interaction
in Jacobi i in the Faddeev representation. The eigenvalues λn and the eigenfunction ΦJM
n
are hyper-radial ρ dependent. A set of coupled channels equations arising from the adiabatic hyperspherical expansion of the wavefunctions is then solved to obtain different states
in 12 C. The results are used to compute the triple-alpha reaction rate. According to their
studies, the formation of 12 C from three alpha particles can happen in several paths. At
temperatures below 1 GK the quadrapole transition from the 0+ continuum (both resonant
and non-resonant) states to the 2+ bound state dominates. However, contribution from the
1
2+ −→ 0+ transition increases and begins to be dominant at T ≈ 2.5 GK. Other transi1

9

tions from the 2+ continuum states to the 2+ bound state do not contribute significantly
1
to the rate for all temperatures. Although the triple-alpha reaction rate agrees well with
what is obtained in NACRE for T > 0.1 GK, below that it is not computable due to the
numerical limitation of the implementation in [26, 27]. Therefore, in order to estimate the
direct triple-alpha rate at low temperatures assuming it captures to the tail of the Hoyle
state directly, Garrido et al. extrapolated a three-body Breit-Wigner cross section for the
three-alpha capture to low energies [28] in a similar manner as done in the sequential process
[10]. The reaction rate from that work (named as BW(3B)) shows an increase of 7 orders
of magnitude at T = 0.02 GK compared to NACRE due to the three-alpha direct capture.
There are large discrepancies in the results of [10], [15] and [28], demonstrating uncertainties
and ambiguities in our understanding of the low-temperature triple-alpha reaction. We want
to develop a theory that provides a better description of the α + α + α system with the correct Coulomb asymptotic behavior, thereby enabling us to compute the reaction rate at very
low temperatures without extrapolation. It is demonstrated in chapter 6 of this thesis that
the low-temperature triple-alpha rate plays important role in several astrophysical scenarios
such as helium accreting white dwarfs and helium accreting neutron stars. We find that the
triple-alpha reaction rate at low temperatures dictates the helium ignition on the surface of
these objects for very low accretion rates. Since the triple-alpha reaction rate is impossible
to measure at low temperatures and can have dramatic effects in astrophysics, a well founded
theory is critical.

10

1.3

Three-body reaction methods in solving the triplealpha problem

At low temperature, the triple-alpha reaction proceeds through the 0+ continuum states to
the 2+ bound state, which then fully deexcites to form 12 C at ground state. In order to
1
construct a reliable non-resonant reaction rate, we need a good description of both 2+ bound
1
states and 0+ continuum states of 12 C. Recent microscopic theories such as the no-core shell
model [29] and the Green’s function Monte Carlo method [30] struggle to reproduce the
Hoyle state of 12 C by solving the 12-body problem. Even though the fermionic molecular dynamics method produces both the 2+ bound state as well as the 0+ resonance with
1
2
significant triple-alpha configurations [31], a microscopic description for the non-resonant
continuum states is currently not available. Given the dominant triple-alpha structure in
both 12 C(2+ ) and 12 C(0+ ), it is reasonable to construct this as a three-body problem as
1
2
done in [15, 26, 27].
The three-body problem is an important topic not only in nuclear physics but also
in atomic and molecular physics. The first three-body problem was studied for a threenucleon system with no Coulomb interactions in a formalism invented by Faddeev [32].
With Coulomb interactions present, the problem became much more complicated. Theories for the three-body problem have made significant progress over the last few decades.
An accurate solution for a three-body bound state system can now be obtained (e.g. [33]).
However, there are still difficulties remaining for a three-body scattering problem because
no standard boundary condition for this system exist when including the long-range effects
of the Coulomb interaction. Some efforts have been made to solve this problem by using an
approximate boundary condition formulated by Merkuriev [34], but the results were limited
11

to p-d scattering systems [35]. Deltuva et al. [36] have developed a screening and renormalization technique for Coulomb interactions in momentum space which can solve a three-body
scattering problem but for cases where only one pair of particles have charge [36]. Our triplealpha problem is more difficult since all three particles are charged.
In this work, we employ the hyperspherical harmonics (HH) method [33, 37, 38] to tackle
the triple-alpha problem. The HH method is excellent to reproduce the 2+ bound state of
1
12 C.

However, applying it directly to calculate the 0+ continuum states of 12 C is more chal-

lenging due to the strong Coulomb interaction that governs the system at large distances.
Despite that difficulty, using hyperspherical coordinates offers a natural way to express the
asymptotic form for the scattering of this three charged-particle system [20]. We combine
the HH method with the R-matrix method and R-matrix propagation (the new method is
named HHR) to compute the 0+ continuum states of 12 C. We also implement a technique
of screening the off-diagonal Coulomb couplings to ensure numerical stability at very low
energies. The scattering wavefunctions for the triple-alpha system are constructed with the
correct Coulomb asymptotic behavior. Since the resonant and the non-resonant continuum
states come naturally from the solutions of the HHR equations, we are able to treat both
resonant and non-resonant, and direct and sequential capture mechanisms on the same footing. The HHR method enables us to obtain numerically stable solutions at low energies and
therefore successfully obtain the triple-alpha reaction rate in the low temperature regime
where many attempts have failed before.

12

1.4

Motivation for this work

The triple-alpha reaction has a long history of investigation and is still a current topic of
interest because it plays a very important role in many astrophysical processes. Although
the theory for the triple-alpha reaction at high-temperatures is well established, our current knowledge for the low temperature regime is very uncertain and ambiguous as seen
through the large deviation among [10], [15] and [28]. The aim of this work is to resolve this
problem by developing a three-body framework that provides an accurate description of the
α + α + α scattering at low temperatures with the correct Coulomb asymptotic behavior.
The non-resonant capture mechanism is therefore taken into account naturally without any
extrapolation and the role of the non-resonant continuum in the triple-alpha rate will be
clarified.

1.5

Outline

In this work we develop a new three-body method for solving the triple-alpha problem. Details on the hyperspherical harmonics R-matrix method are presented in chapter 2. The
hyperspherical coordinates will be defined. The R-matrix expansion and the R-matrix propagation method in the hyperspherical harmonic (HH) representation are introduced to tackle
the difficulty of three charged particles in the asymptotic region. We also discuss the screening technique for the off-diagonal Coulomb couplings to ensure the numerical stability of
the problem at low energies. The three-body quadrupole strength function and the reaction
rate formalism are derived in chapter 3. In chapter 4 we present the main results of this
work which include a thorough convergence study of the problem, a closer look at the cluster
structure of the 2+ bound state and the 0+ Hoyle resonant state of 12 C, and the new triple1
2
13

alpha reaction rate in which both the resonant and the non-resonant capture mechanisms
are naturally taken into account. We then discuss the reaction dynamics of the triple-alpha
problem in chapter 5. In this chapter we want to understand the effect of long-range Coulomb
effects on the reaction rate at low temperatures as well as the dynamics of the three alpha
particles at low energy when intermediate resonances are impossible. We also compare our
rate with other studies. In chapter 6 we focus on astrophysics, exploring the impact of the
new rate on stars, helium accreting white dwarfs, and helium accreting neutron stars. Will
it produce drastic changes as the rate in [15]? The summary of our work and conclusions
are drawn in chapter 7.

14

Chapter 2
The Hyperspherical Harmonics
R-matrix (HHR) method

2.1

Hyperspherical Harmonics method

The hyperspherical harmonics (HH) method originated in atomic and molecular physics
[39, 40]. Later, it was extended to few-body systems in nuclear physics by Delves [41, 42]. The
theory of the HH method was well developed for Borromean systems [43, 44, 45]. Borromean
nuclei are defined as systems of three particles which are loosely bound and have no bound
states in any of the two-body subsystems. 11 Li and 6 He which have two neutron weakly
coupled to the core, are typical examples of Borromean nuclei. In this work, we formulate
our problem as a Borromean system of three alpha particles. Many details on the HH method
can be found in the literature (e.g., [33], [37], [38]). In the current chapter we will present
an overview of the method.

15

Figure 2.1: Three ways to form a Jacobi coordinate set for a system of three particles.

Figure 2.2: Coordinates defined for Jacobi set i

2.1.1

Hyperspherical coordinates

Let us consider a system of three nuclei with masses mi . For a three-body system there are
three ways to form a Jacobi coordinate set (xi , yi ) as presented in Fig. 2.1. For each Jacobi
set i, we define the scaled coordinates (see Fig. 2.2):

xi =

Aj Ak
r ,
Aj + Ak i

yi =

Ai (Aj + Ak )
R .
Ai + Aj + Ak i

(2.1)

Here, ri = ak − aj is the relative radius from particle j to particle k and Ri = acm − ai is
jk
the radius from particle i to the center of mass of the two-body subsystem (jk). The ratios
16

Ai = mi /m are dimensionless and the scaled mass m is usually taken as the nucleon mass.
The hyperspherical coordinates are then defined as functions of xi and yi :

2
ρ2 = x 2 + yi ,
i

x
θi = arctan i .
yi

(2.2)

The hyper-radius ρ is universal in all three Jacobi coordinates. It does not change under
translations, rotations, and permutations and therefore provides information about the size
of the system:
2
rrms =

2
i Ai rrms (i)

A

+

1
< ρ2 > ,
A

(2.3)

where A = Ai + Aj + Ak ; rrms and rrms (i) are the root-mean-square (rms) radius of the
system and particle i, respectively.
The hyper-angle θi depends on the selected Jacobi set, containing some insights about
the configuration of the system. For example, θ3 ≈ 0 tells us that particle 1 and 2 are much
closer to each other than to particle 3. If θ3 ≈ π/2, we could say that particle 1 and 2
are very far apart and particle 3 lie between them. For the triple-alpha problem in which
three particles are identical, all the Jacobi sets are equivalent, we choose to work with Jacobi
set i = 3 and drop the index i in our equations from now on for convenience. The scaled
coordinates x and y can be reexpressed in terms of the hyperspherical coordinates as

x = ρ sinθ ,
y = ρ cosθ .

17

(2.4)

2.1.2

Three-body bound states

In this section we describe a three-body bound system in the hyperspherical harmonics
representation. Since the triple-alpha involves only spin-zero particles, all the equations in
HH coordinates take this into account because the equations are reduced to much simpler
forms. The general case of non-zero spin particles is presented with details in Chapter 9 of
[46].
We first write down the three-body Schrodinger equation in the scaled Jacobi coordinates:

H3b ΨLM = [T + V ] ΨLM = E ΨLM ,
2

(2.5)

2m

[∆2 + ∆2 ] ,
x
y

(2.6)

Vij + V3b ,

T = −

(2.7)

3

V

=
j>i=1

where, E is the three-body relative energy which is negative for a bound system and positive
for a scattering system; m is the scaled mass;

is the Planck’s constant. The three-body

Hamiltonian H3b is the sum of the kinetic operator T and the potential V . The threebody wavefunction ΨLM describes a state of total spin and angular momentum L and its
projection M .
The hyperspherical expansion separates the radial and angular dependence of the threebody wavefunction. For a bound state, the expansion has the following form:

5

lx ly

ΨLM = ρ− 2

χL x ly (ρ) ϕK (θ) [Ylx ⊗ Yly ]LM ,
Kl
Klx ly

18

(2.8)

where ρ−5/2 cancels out the factor ρ5 in the volume element dV = ρ5 sin2 θcos2 θdρdθdΩx dΩy
to simplify the normalization condition of a three-body bound state wavefunction:

ˆ ˆ
ΨLM (ρ, θ, x, y)

2

χL x ly (ρ)
Kl

dV =

2

dρ = 1

(2.9)

Klx ly

From now on we will drop the L notation in the hyper-radial wavefunction χL l (ρ) for
Klx y
convenience. Eq. (2.8) is written for the case of zero-spin particles as we discussed above (a
lx ly

more general formula can be found in [37]). χKlx ly (ρ), ϕK (θ) and Y are the hyper-radial
functions, the hyper-angular functions, and the spherical harmonics functions respectively;
lx and ly are the orbital angular momenta corresponding to the Jacobi coordinates x and y.
This expansion introduces a new quantum number, the hyper-momentum K, an extended
concept of angular momentum for a three-body system.
We use the definition of hyperspherical coordinates Eq. (2.4) and the hyperspherical
wavefunction expansion Eq. (2.8) to rewrite the kinetic energy in Eq. (2.6) in terms of the
hyperspherical variables ρ and θ:
2

1 ∂
1
∂
ρ5
+ 2 Λ2 ,
5 ∂ρ
2m ρ
∂ρ
ρ
∂
∂
lx (lx + 1) ly (ly + 1)
1
sin2 2θ
−
−
Λ2 =
.
2 2θ ∂θ
∂θ
cos2 θ
sin
sin2 θ

T =−

(2.10)
(2.11)

The hyper-angle operator Λ2 has its eigenvalues as a function of hyper-momentum K and
lx ly

its eigensolutions being the hyper-angular functions ϕK (θ) :
lx ly

lx ly

Λ2 ϕK (θ) = −K(K + 4) ϕK (θ) .

19

(2.12)

lx +1/2,ly +1/2

These eigensolutions can be written in terms of Jacobi polynomials Pn

lx ly

lx ly

ϕK (θ) = NK

lx +1/2,ly +1/2

(sinθ)lx (cosθ)ly Pn

(cos2θ) .

(cos2θ):

(2.13)

lx ly

In order for the hyper-angular functions ϕK (θ) to form an orthonormal set with the weight
factor sin2 θ cos2 θ:
π/2

lx ly

ϕK (θ) ϕ

0

lx ly
(θ) sin2 θ cos2 θ dθ
K

(2.14)

= δKK

lx ly

the normalization coefficients NK ,
lx ly

NK

=

2lx +ly +2n+3

(lx + ly + 2n + 2) n! (lx + ly + n + 1)!
,
[2(n + lx ) + 1]!! [2(n + ly ) + 1]!! π

(2.15)

can be derived with the help from [47]. Using Eq. (2.12) and Eq. (2.8), the three-body
Schrodinger equation for a bound state reduces to a set of coupled channels equations that
only depends on the hyper-radius ρ:
2

2m

d2 ∆K (∆K +1)
−
+E
dρ2
ρ2

χγ (ρ) =

Vγγ (ρ) χγ (ρ) ,

(2.16)

γ

where ∆K = K + 3/2 and γ = {K, lx , ly }. The coupling potentials Vγγ (ρ) are the sum of
interactions defined in Eq. (2.7), integrated over all variables but ρ:
3

Vγγ (ρ) =

ˆ ˆ
Ωγ (θ, x, y)

ˆ ˆ
Vjk + V3b Ωγ (θ, x, y)

,

(2.17)

k>j=1

lx ly

ˆ ˆ
where Ωγ (θ, x, y) = ϕK (θ) [Ylx (ˆ ) ⊗ Yly (ˆ )] are the hyper-harmonic basis functions conx
y
taining all the angular dependence in coordinate space of Eq. (2.8). For a system of three

20

charged particles like the triple-alpha, Vγγ (ρ) are long-range causing significant difficulties
in determining the boundary conditions of the problem. It is thus important to study the
properties of these couplings.

2.1.2.1

Asymptotic behavior of the three-body Coulomb interaction

In the asymptotic region (ρ −→ ∞), the potential couplings Vγγ (ρ) only contains a sum of
three pairwise Coulomb interactions. All the nuclear forces as well as the three-body force
go to zero at large distance:

V

asy
(ρ)
γγ

C(12)
C(23)
C(31)
(ρ) + V
(ρ) + V
(ρ).
γγ
γγ
γγ

=V

(2.18)

As we choose to present all equations in the T-basis (Jacobi coordinates i = 3 (Fig. 2.1c)),
C(12)
(ρ)
γγ

calculating V

is straight forward. At large radius, the Coulomb interaction between

particle 1 and 2 is expressed as V C(12) = Z 2 e2 /r3 where Z is the charge of an alpha
particle; e is the elementary electric charge and r3 is the radius between them. Although
this expression is modified at small radii to take into account the finite size of each alpha
particle (see Sec. 4.1), it is exact for the two-body Coulomb interaction in the asymptotic
region. Using Eq. (2.1) and Eq. (2.4), V

V

C(12)
(ρ)
γγ

=

Z 2 e2
ρ

A
2

=

Z 2 e2
ρ

A
2

γ, 3

ˆ ˆ
Ωγ (θ3 , x3 , y3 )

C(12)
(ρ)
γγ

1
γ ,3
sinθ3

is expressed as:

,

1
ˆ ˆ
Ω∗ (θ3 , x3 , y3 ) sin2 θ3 cos2 θ3 dθ3 dΩx3 dΩy3 ,
sinθ3 γ

21

(2.19)

where A is the atomic number of an alpha particle. We then simplify the integration in
Eq. (2.19) to obtain:

C(12)
(ρ)
V
γγ

Z 12
=

12
Zγγ = Z 2 e2

γγ

ρ
A
2

,

(2.20)
lx ly

ϕK (θ3 )

∗lx ly
1
(θ) sin2 θ3 cos2 θ3 dθ3 ,
ϕ
K
sinθ3
C(23)
(ρ)
γγ

It is more difficult to express the other two components V

(2.21)

C(31)
(ρ)
γγ

and V

in

Jacobi i = 3 because vectors r1 and r2 are not proportional to the scaled coordinate x3 .
The advantage of the HH method is that the hyper-radius ρ is the same for all three Jacobi
sets, leading to a possibility of six-dimensional coordinate rotations. For example, Jacobi
coordinates (x3 , y3 ) can be obtained by rotating Jacobi coordinates (x2 , y2 ) by an angle φ:

x3 = −cosφ x2 + sinφ y2 ,

(2.22)

y3 = −sinφ x2 − cosφ y2 ,

(2.23)

where the rotation angle φ is a function of the mass of the three particles. Based on this
argument, Raynal and Revai developed a way to rotate the hyper-harmonic basis function
ˆ ˆ
Ωγ (θ, x, y) from one coordinate set to another using the multiplying Raynal-Revai coefficients
C(23)
C(31)
(ρ) and V
(ρ)
γγ
γγ

which can be calculated analytically [48]. The potential couplings V

are first computed in Jacobi i = 1 and i = 2 in Fig. 2.1a and Fig. 2.1b respectively and then
rotated back to Jacobi i = 3 using the Raynal-Revai frame-rotation technique. In the end,
the potential couplings Vγγ (ρ) in the asymptotic region can then be expressed as:
Z eff
Vγγ (ρ) =

22

γγ

ρ

,

(2.24)

where Z eff is constant for any given channel γ and γ and has a non-trivial expression:
γγ


eff
Zγγ = Z 2 e2

1
A
γ, 3
γ ,3 +
2
sinθ3

32 23
γα α γ
αα

α, 2

1
α ,2
sinθ2


31 13
γβ β γ

+

β, 1

ββ

where

ij
βγ

1
β ,1  ,
sinθ1

(2.25)

are the Raynal-Revai coefficients [48]; i = 1, 2, 3 represent the three Jacobi

coordinates (xi , yi ) (Fig. 2.1). Each bra-ket matrix element is calculated as in Eq. (2.19).
Eq. (2.24) shows that both the diagonal and the off-diagonal couplings decay slowly as
1/ρ for a three charge system in the HH representation, as seen in [20]. Fig. 2.3 demonstrates
the hyper-radial behavior of the potential couplings Vγγ (ρ). In this graph, couplings are
calculated by fixing γ = (0, 0, 0) and varying γ . Although the diagonal term γ = γ =
(0, 0, 0) is dominant, all the couplings do not vanish at a large hyper-radius ρ = 800 fm.
This demonstrates the difficulty of our problem in the asymptotic region.
For neutral Borromean systems, the long range couplings are not present and the threebody bound state wavefunction decays exponentially at large distances [49, 50]:

ρ→∞
χγ (ρ) −→ e−κρ .

(2.26)

When introducing Coulomb interactions, a simple analytic expression for the asymptotic
bound state wavefunction of three charged particles no longer exists. However, for a wellbound system, imposing the boundary condition that the wavefunction goes to zero at large

23

0.2

Diagonal term Vγγ

Vγγ’

0.1
0

γ’=(0,0,0)
γ’=(2,0,0)
γ’=(4,2,2)
γ’=(4,0,0)
γ’=(6,2,2)
γ’=(6,0,0)

-0.1
-0.2
-0.3
-0.4

100 200 300 400 500 600 700 800

ρ [fm]
Figure 2.3: Hyper-radial behavior of the potential couplings Vγγ (ρ) in Eq. (2.17). The black
solid curve represents the diagonal coupling and the others demonstrate the off-diagonal
couplings.
distances is sufficient for the numerical calculation [37]:

χγ (ρ → ∞) → 0 .

2.1.3

(2.27)

Three-body continuum states

i i
For a three-body scattering wavefunction, an additional set of subscripts K i lx ly is needed to

describe the incoming channels of the system in which three particles come in with direction
k (κ, θκ ) in momentum space:

ΨLM =

i i
K i lx ly

1
(κρ)5/2

i i
K i lx ly
Klx ly

i i
lx ly

lx ly

χKl l (κρ) ϕK (θ) [Ylx ⊗ Yly ]LM ϕ i (θκ ) [Yli ⊗ Yli ]κ , (2.28)
xy
K
x
y LM

24

where κ is the hyper-radial momentum and relates to the three-body energy E through
√
κ = 2mE/ . The (κρ)5/2 factor is introduced to simplify the normalization condition of
the three-body scattering wavefunction:
∞
0

χγγ i (κρ) χγγ i (κ ρ) dρ = δ(κ − κ ) .

(2.29)

i i
where γ = (K, lx , ly ) [γ i = (K i , lx , ly )] are the outgoing [incoming] channels. Using this HH

expansion, we obtain a set of coupled channels equations for a three-body scattering system:
2

d2 ∆K (∆K +1)
−
+E
2m dρ2
ρ2

χγγ i (κρ) =

Vγγ (ρ) χγ γ i (κρ) ,

(2.30)

γ

where Vγγ (ρ) are the potential couplings as defined in Eq. (2.17).
Solving Eq. (2.30) for positive energy E to obtain continuum states for a system of three
charges is numerically challenging since it requires an exact boundary condition. In order
to overcome that difficulty we implement a screening technique in which all the off-diagonal
couplings are screened at large distances. Because the diagonal couplings are dominant in the
asymptotic region (see Fig. 2.3), the screening method can help to achieve numerical stability
and the correct three-body scattering boundary condition while preserving the physics of the
three Coulombian particle problem. More details on the screening method will be presented
in Sec. 2.2.3. In the region where all the off-diagonal couplings vanish, Eq. (2.30) can be
rewritten as:
d2 χγγ i (x)
dx2

+ 1−

2ηγ
∆ (∆ + 1)
− K K
x
x2

χγγ i (x) = 0 ,

(2.31)

where x = κρ and ∆K = K + 3/2. The variable ηγ is an equivalent of the Sommerfeld parameter in two-body systems but is more complex and hyper-momentum K dependent. ηγ
25

eff
eff
is obtained from the parameter Zγγ in Eq. (2.25) through the relationship ηγ = mZγγ / 2 κ.

It is well known that Eq. (2.31) has analytic solutions which are the regular FK+3/2 (ηγ , κρ)
and irregular GK+3/2 (ηγ , κρ) three-body Coulomb functions [46]. Therefore, in the asymptotic region where all the off-diagonal couplings are negligible, the α + α + α scattering
wavefunction behaves as:

ρ→∞

−
+
χHHR (ρ) −→ HK+3/2 (ηγ , κρ) δγγ i − HK+3/2 (ηγ , κρ) Sγγ i .
i
γγ

(2.32)

In this equation, H ± = G ± iF are the Coulomb functions describing the outgoing and
incoming spherical waves and Sγγ i is the scattering matrix [46].

2.1.4

Properties of the three-alpha wavefunction

In this section we want to discuss about the symmetry properties of the three-alpha wavefunction. An alpha particle is considered as a boson of spin zero, thus the wavefunction
must be unchanged for any permutation of these two particles. This selects certain partial
waves contributing to the HH wavefunction expansion (Eq. (2.8) and Eq. (2.28)). Next, we
derive the relation describing this selection rule. Fig. 2.1 represents three different Jacobi
coordinates. These give rise to three representations of the three-body wavefunction:

|Ψ1 = [(lx1 , ly1 )L1 , (s2 , s3 )Sx1 ]J1 , s1 ; J [[t1 ⊗ t2 ] ⊗ t3 ]T ,

(2.33)

|Ψ2 = [(lx2 , ly2 )L2 , (s3 , s1 )Sx2 ]J2 , s2 ; J [[t1 ⊗ t2 ] ⊗ t3 ]T ,

(2.34)

|Ψ3 = [(lx3 , ly3 )L3 , (s1 , s2 )Sx3 ]J3 , s3 ; J [[t1 ⊗ t2 ] ⊗ t3 ]T .

(2.35)

26

We denoted si (ti ) as the spin (isospin) of particle i and J(T ) as the total spin (total isospin)
of the three-body system. lxi , lyi are the the relative orbital angular momenta of the twobody subsystem in Jacobi i and the orbital angular momenta between its center of mass and
particle i respectively. If we interchange the position of particle 1 and particle 2 which are
identical, the total wavefunction in the third Jacobi coordinate becomes:

Ψ3 = [(−lx3 , ly3 )L3 , (s2 , s1 )Sx3 ]J3 , s3 ; J [[t2 ⊗ t1 ] ⊗ t3 ]T .

(2.36)

We can rewrite this new wavefunction in terms of Ψ3 with the aid from [51]:

Ψ3 = (−1)lx3 (−1)Sx3 −s1 −s2 (−1)T21 −t2 −t1
× [(lx3 , ly3 )L3 , (s1 , s2 )Sx3 ]J3 , s3 ; J [[t1 ⊗ t2 ] ⊗ t3 ]T ,
Ψ3 = (−1)lx3 +Sx3 +T21 |Ψ3 .

(2.37)

Since particle 1 and particle 2 are bosons, their spins and isospins are even numbers which
can then be dropped from Eq. (2.37). Because the wavefunction is symmetric under the
exchange of two identical bosons, Eq. (2.37) requires lx3 + Sx3 + T21 to be even. Therefore
we can conclude that for three identical bosons, expressed in Jacobi coordinates 3, only
partial waves satisfying the selection rule lx + Sx + Tx = even will contribute to the total
wavefunction. Here, lx , Sx and Tx are the relative angular momentum, the total spin, and
total isospin of the two-body subsystem in respect to that Jacobi coordinate system. For
the triple-alpha problem, each alpha particle has spin and isospin zero, the selection rule is
thus reduced to lx = even.
Using the same techniques, we could write down the wavefunction Ψ2 in the second

27

Jacobi coordinate after exchange particle 1 and particle 2:

Ψ2 = [(−lx1 , ly1 )L1 , (s3 , s2 )Sx1 ]J1 , s1 ; J [[t2 ⊗ t1 ] ⊗ t3 ]T .

(2.38)

Rearranging the order of momentum couplings we have:

Ψ2 = (−1)lx1 (−1)Sx1 −s3 −s2 (−1)T21 −t2 −t1
× [(lx1 , ly1 )L1 , (s2 , s3 )Sx1 ]J1 , s1 ; J [[t1 ⊗ t2 ] ⊗ t3 ]T ,
Ψ2 = (−1)lx1 +Sx1 +T21 Ψ1 .

(2.39)

Again, spins and isospins of particle 1 and 2 are even, therefore they can be dropped from
Eq. (2.39). From Eq. (2.39), we can generalize a relationship for the wavefunction of a three
identical boson system between any pair of Jacobi sets:

Ψj = Pjk Ψk ,

(2.40)

T +l +S
where the permutation operator Pjk is defined as Pjk = (−1) jk xk xk . For the triple-

alpha problem, this quantity equals to one when applying the even partial wave selection
rule that we discuss above.
Eq. (2.40) indicates that only one of these wavefunctions should be explicitly included in
the three-component Faddeev equations:






Vi
Vi
Ti + Vi − E
  Ψi 

 

 

  Ψj  = 0 .
Vj
Tj + Vj − E
Vj

 

 
Vk
Vk
Tk + Vk − E
Ψk
28

(2.41)

Here, subscript i and the corresponding variables keep their common definitions in Jacobi
set i. With the use of Eq. (2.40), we can effectively reduce three sets of Faddeev equations
to a set of coupled channel equations in one set of Jacobi coordinates:

i|Ti + Vi − E|i i|Ψi + i|Vi |j j|Ψi + i|Vi |k k|Ψi = 0 ,

(2.42)

which can be rewritten as:

[Ti + Vi + i|Vi |j j|i + i|Vi |k k|i − E] Ψi = 0 .

(2.43)

The constraint condition of the partial wave in combination with the symmetry property of
Eq. (2.43) allow us to fully take into account the symmetrization of the system and highly
increase the computational efficiency.

2.2

R-matrix method

In order to obtain the continuum states of a three charged-particle system, we need to solve
Eq. (2.30) for positive values of energy E. Theoretically, we could integrate N linearly independent solutions of Eq. (2.30) from ρ = 0 to a certain large radius ρa where the off-diagonal
couplings Vγγ (ρ) are negligible [43, 52]. Then, a combination of these solutions would be
matched with the asymptotic wavefunction Eq. (2.32) at ρa to construct the physical continuum wave of a three-body system. However for a three charged-particle system, ρa has to be
very large (∼ 1000 fm) due to the long range Coulomb interaction, making it inefficient to
use the direct integration methods. Since we are interested in the low energy region in which
the wavefunction is needed below the Coulomb barrier, the direct integration methods are
29

not numerically stable for solving coupled channels equations with strong repulsive couplings
[46]. In this section, we discuss about the R-matrix method in hyperspherical coordinates
which enables us to overcome these problems.
Wigner and Eisenbud [53] were the first to introduce the R-matrix method into nuclear
physics to study resonances in nuclear scattering. Later on, this method was developed and
extended in detail by Lane and Thomas [54]. The general idea of this method is to use
an orthonormal basis expansion inside an R-matrix box. An R matrix is then constructed
to match to the asymptotic wavefunction outside the box [46]. Even though the R-matrix
method was originally developed for a two-body scattering problem, it has been generalized
to a three-body system in HH coordinates by I.J. Thompson et al. [55].

2.2.1

R-matrix expansion

In the R-matrix expansion, the eigenfunctions of the diagonal parts of the Hamiltonian in
each separate channel γ are used as a basis:
(K +
d2
−
2
2m dρ
2

−

3 )(K
2
ρ2

5
+ 2)

+ Vγγ (ρ) − εn
γ

n
ωγ (ρ) = 0 ,

(2.44)

n
where εn is the eigenenergy corresponding to the eigenfunction ωγ (ρ). These basis functions
γ

are defined to have the same logarithmic derivatives β =

n
dlnωγ (ρ)
dρ

at the R-matrix radius ρm .

n
The constant feature of the logarithmic derivatives ensures that ωγ (ρ) forms an orthonormal

basis set within the range of [0, ρm ].
ρm

ωn (ρ) ωm (ρ) dρ = δnm .
0

30

(2.45)

We then solve the coupled channels equations in the interior region [0, ρm ]:
2

−

(K +
d2
−
2
2m dρ

3 )(K
2
ρ2

+ 5)
2

p
− p φγ (ρ) +

p
γ

Vγγ (ρ)φ (ρ) = 0 ,

(2.46)

γ

p

n
by expanding the eigensolution φγ at eigenenergy ep in terms of the orthonormal basis ωγ (ρ):
N
p

pn

n
cγ ωγ (ρ) .

φγ (ρ) =

(2.47)

n=1
pn

The expansion coefficients cγ can easily be obtained by solving a set of linear equations:

pn

pn
γ

ε n cγ +
γ

n
n
< ωγ |Vγγ |ωγ > c

pn

= ep cγ .

(2.48)

nγ
p

These eigenfunctions φγ (ρ) will also form an orthonormal basis within the R-matrix radius
ρm . The three-body scattering wavefunction χγγ i (ρ) at energy E, a solution of Eq. (2.30)
p

in the whole coordinate space, can be expanded in terms of φγ (ρ):

χγγ i =

p

p

p

A i φγ .
γγ

(2.49)

In our calculation, the logarithmic derivative β is chosen to be zero. This option may
lead to the discontinuity of the derivative of the wavefunction at the R-matrix boundary
when using a truncated basis as discussed in [56]. However, a sufficiently large basis is
employed in our calculation to ensure a smooth behavior of the wavefunction at all energies.
In addition, the zero logarithmic derivative highly reduces the complexity of the R-matrix
propagation method (Sec. 2.2.2), therefore increasing the efficiency and accuracy of the
numerical performance when propagation is needed to a very large radius. All the equations
31

presented in this thesis consider β = 0 while more general formulae can be found in [46].
The R matrix at energy E relates the wavefunction χγγ i (ρ) and its derivative:

ρ Rγγ (E) χ

χγγ i (ρ) =

γ γi

(ρ) .

(2.50)

γ

The constant logarithmic derivative feature enables us to calculate the R matrix Rγγ (E)
p

from the known eigenfunctions φγ at radius ρm :

2

Rγγ (E) =

2mρm

P
p=1

p

p

φγ (ρm ) φ (ρm )
γ
ep − E

.

(2.51)

where P is the number of poles used in our calculation. If we substitute Eq. (2.49) and
Eq. (2.51) into Eq. (2.50), the only missing quantities to calculate the scattering wavefunction
p

are the expansion coefficients A i . As discussed in Sec. 2.1.3, Eq. (2.32) is the asymptotic
γγ
form of a continuum wavefunction for the three alpha system when the off-diagonal couplings
can be neglected. For sufficiently large ρm , the wavefunction will be in the asymptotic regime.
p

The unknown coefficients A i are obtained by performing an asymptotic matching of the
γγ
three-body scattering wavefunction at ρm :
2
1
p
A i=
p
γγ
2m eγ − E

φ

p
γ

−

+

δγ γ i H K +3/2 (κρm ) − Sγ γ i H K +3/2 (κρm ) .

(2.52)

γ

All functions in Eq. (2.52) are evaluated at the R-matrix boundary ρm . The scattering
matrix Sγ γ i and the R matrix are directly related [46]:
−

δγ γ i H −
(κρm ) − ρm Rγ γ i H K +3/2 (κρm )
K +3/2
Sγ γ i =
.
δγ γ i H +
(κρm ) − ρm Rγ γ i H +
(κρm )
K +3/2
K +3/2

32

(2.53)

2.2.2

R-matrix propagation

The triple-alpha system is driven by a strong Coulomb interaction with long-range offdiagonal couplings. Therefore, the R-matrix radius ρm has to be very large to fully capture
the physics of the problem. Solving the coupled channels equations in a large R-matrix box
will lead to numerical instability. Fig. 2.4 presents the scattering wavefunction calculations
at 0.2 MeV using two different R-matrix radius. The red curve corresponds to an R-matrix
box of ρm = 50 fm. The blue curve represents a calculation with ρm = 100 fm. We observe
a numerical instability in the small radius region when using a larger R-matrix box. Large
R-matrix box requires more terms in the R-matrix expansion Eq. (2.49) and Eq. (2.51). The
R-matrix basis wavefunctions around the origin are expected to be infinitesimally small due
to the strong Coulomb barrier. Accumulation error of adding many small numbers could
lead to the instability of the solution at small radius as seen in Fig. 2.4.
In order to overcome this difficulty, we employ the R-matrix propagation technique originally developed by Light and Walker for the atom-molecule scattering problem [57, 58, 59].
First, the set of coupled channels equations (2.30) is solved in a small R-matrix box with
radius ρm . Then the R-matrix at ρm is propagated to ρa

ρm where it becomes safe to

ignore the off-diagonal couplings, and matching to known Coulomb functions of Eq. (2.32)
can be performed. This propagation method is well known for its fast convergence and high
numerical stability in atomic and molecular system. We first rewrite Eq. (2.30) as
d2 χγγ i (ρ)
dρ2

˜
Vγγ χγ γ i (ρ) ,

=
γ

33

(2.54)

10

χ(ρ)

10
10
10

0

-3

ρm=100 fm
ρm=50 fm

-6

-9

-12

10 0

20 40 60 80 100 120 140 160 180 200

ρ [fm]
Figure 2.4: The hyper-radial scattering wavefunction corresponding to the incoming channel
i i
of K i lx ly = 000 and the outgoing channel of Klx ly = 000 in the hyperspherical expansion
(Eq. (2.28)) is calculated at E = 0.2 MeV for two different R-matrix box sizes: 50 fm (red)
and 100 fm (blue).

Figure 2.5: Diagram presenting the sectorization in hyper-radius ρ of the R-matrix propagation method. The subscripts L and R indicate the evaluations at the left and right sector
boundary. The R matrix R of a sector is calculated at the right boundary.

34

where
2m
˜
Vγγ = 2 Vγγ +

∆K (∆K + 1) 2mE
− 2
ρ2

δγγ .

(2.55)

The interval from the R-matrix radius ρm to the matching radius of the three-body asymptotic wavefunction ρa is divided into M sectors, illustrated in Fig. 2.5. We choose the size hp
of sector p to be sufficiently small that the interaction within is considered to be constant
˜
λγ (p) in each channel γ. We diagonalize Vγγ by solving the equation:

˜
˜
(Tp )T V(ρp ) Tp = λ(p)2 ,

(2.56)

where ρp is taken at the center of sector p. The wavefunction and its derivative at the left
and right boundary of sector p can be related through a local diagonal representation of the
propagating functions Gp :



p
p
p
 G1 G2   −χR 

.

=

p
p
p
p
χL
G3 G4
χL


p
χR



(2.57)

The subscripts R and L imply evaluations at the right and left side of the sector boundary,
respectively. In this section, all the equations are written in matrix form for simplicity. The
p

bold letter indicates a matrix quantity, e.g χR is the wavefunction matrix of which each
p
matrix element is (χγγ i )R corresponding to the incoming γ i and outgoing γ channels. The
p

p

˜
propagating functions are expressed as Gi = Tp gi Tp , where gp is a simple function of λγ (p)

35

(see Appendix A for detailed derivation):

p
(g1 )γγ

p

=

p
(g4 )γγ

= δγγ



 1
−

coth|hp λγ |,
|λ |
γ

λ2 > 0
γ


 1


cot|hp λγ |, λ2 ≤ 0
γ
|λγ |


 1
−

csch|hp λγ |, λ2 > 0
γ
|λ |
γ

p

(g2 )γγ = (g3 )γγ = δγγ

 1


csc|hp λγ |,
|λ |
γ

.

λ2 ≤ 0
γ
(2.58)

The R matrix of sector p relates the wavefunction and its derivative at the right sector
boundary through a matrix equation:

p

p

p

χR = ρR Rp χR .

(2.59)

We can obtain a similar expression for sector p−1:

p

p−1

p

χL = ρR Rp−1 χL .

(2.60)

Inserting Eq. (2.59) and Eq. (2.60) into Eq. (2.57), we derive an equation that propagates
the R matrix from sector p − 1 to the next (see Appendix A for detailed derivation):

1
p p
p−1
p
p
Rp = p G2 [G4 − ρR Rp−1 ]−1 G3 − G1 .
ρR

(2.61)

The R-matrix propagation in Eq. (2.61) is the first step in the propagation method to
calculate the three-body scattering wavefunction (see Fig. 2.6). We first propagate the R

36

Figure 2.6: Schematic diagram of the propagation method in obtaining the three-body scattering wavefunction
matrix from a small radius ρm to a large distance ρa where the asymptotic matching can
be performed and the wavefunction (χM ) and its derivative (χM ) at the boundary of the
R
R
final sector M are obtained. These quantities are then used in the global propagation which
will be discussed later to compute the wavefunction (χ1 ) at the R-matrix boundary ρm . We
L
use that to determine the wavefunction expansion coefficient Aγγ i in Eq. (2.49) from which
the three-body scattering wavefunction within the R-matrix radius can be calculated. For
ρ > ρm , its value is obtained by propagating the wavefunction from χ1 to the desired radius
L
instead.
We now discuss the global propagation and the wavefunction propagation which are important in constructing the three-body scattering wavefunction. The global propagation is
designed to propagate the wavefunction from an arbitrary sector p to sector 1 and vice versa
(see Appendix A for details):
37





p
χR

χ1
L







p
G 1


=

p
G2 

p
G3

p
G4



p
−χR

χ1
L



,

(2.62)

p

where the global propagation functions G i are calculated by a set of recursive equations:

p

p

p

p

p−1 −1 p
] G3 ,

p

p

p

p

p−1

[G4 + G 1

p

p−1

− G3

G 1 = G1 − G2 [G4 + G 1

(2.63)

p−1 −1 p−1
] G2 ,

(2.64)

p−1 −1 p
] G3 ,

(2.65)

G 2 = G2 [G4 + G 1
G3 = G3
G4 = G4

p

p−1

p

p−1 −1 p−1
] G2 .

[G4 + G 1

(2.66)

Given the wavefunction χM and its derivative χM at the matching radius ρa , we can
R
R
evaluate its value at the R-matrix boundary ρm by using Eq. (2.62). The wavefunction
is then propagated to larger radius by a sequential application of the local propagation
p

functions Gi in Eq. (2.57):

p

p−1

χL = G1

p−1 −1 p−1
] χL

[G3

p−1

+ G2

p−1

− G1

p−1 −1 p−1
] G4

[G3

p−1

χL

,

(2.67)

p

or a single application of the global propagation functions G i in Eq. (2.62):

p

p

p

χR = G 1 p [G 3 ]−1 χ1 + G 2 − G 1 p [G 3 p ]−1 G 4 p χ1 .
L
L

(2.68)

Fig. 2.7 compares the quadruple transition strength functions dB(E2)/dE calculated using
the global (solid) and local (dashed) representations of the wavefunction propagation method.
Since the two curves are identical, we conclude that these methods are numerically equivalent

38

0

-1

dB(E2)/dE [e fm MeV ]

10

-20

2

4

10
10
10

Global
Local

-40

-60

-80

10 0

0.1

0.2

0.3

0.4

0.5

E [MeV]
Figure 2.7: The quadrupole transition strength function dB(E2)/dE is calculated using the
global (solid) and local (dashed) representations of the wavefunction propagation method.
and provide consistent results.

2.2.3

Screening technique

The R-matrix propagation method is numerically stable when the Coulomb off-diagonal couplings are not introduced. The presence of very narrow two-body and three-body resonances
in addition to the strong, long-range Coulomb interaction in our problem creates numerical
instabilities in the propagated three-body scattering wavefunctions. The blue dashed curved
in Fig. 2.8 represents the first component of the three-alpha scattering wavefunction in the
hyperspherical expansion at 0.5 MeV which is calculated when considering both nuclear and
Coulomb interactions to all orders. The numerical instability occurs around 90 fm which
is also the turning point for this channel. Because the propagation method calculates the
wavefunction at a given sector using the previous sector’s information, the numerical insta-

39

10

χ(ρ)

10
10
10

3

0

-3

screen
no screen

-6

-9

10 0

20 40 60 80 100 120 140 160 180 200

ρ [fm]
Figure 2.8: The hyper-radial scattering wavefunction corresponding to the incoming channel
i i
of K i lx ly = 000 and the outgoing channel of Klx ly = 000 in the hyperspherical expansion
(Eq. (2.28)) is calculated at E = 0.5 MeV for both cases: with screening (red solid) and
without screening (blue dashed).
bility at the turning point region is propagated outward as seen in Fig. 2.8.
We overcome these difficulties by introducing a Woods-Saxon screening factor on the
off-diagonal potentials of Eq. (2.17):

Vγγ =

Vγγ
1 + exp ρ−ρscreen
ascreen

,

(2.69)

where γ = γ . Fig. 2.3 demonstrates that the diagonal couplings Vγγ (ρ) are dominant at
large radius (∼ factor of 10 larger than the off-diagonal couplings). Even though the offdiagonal couplings are small, the fact that they decay slowly as 1/ρ is responsible for the
numerical instability of the wavefunctions as shown above. When the off-diagonal couplings
are screened, we are able to improve the numerical stability and at the same time retain
the long-range asymptotic behavior of the problem. The red solid curve in Fig. 2.8 shows
40

∆X/X

10

24

10

16

10
10

8

0

0.05

0.1

0.15

0.2

E [MeV]
Figure 2.9: Relative deviation of the quadrupole transition strength dB(E2)/dE (denoted
X in the plot) when diagonal Coulomb screening is introduced.
an improvement of the scattering wavefunction when using the screening technique. The
screening radius ρscreen is chosen sufficiently large to gain convergence, yet small enough to
ensure numerical stability. Note that we should not screen the diagonal Coulomb couplings.
The additional diagonal Coulomb screening leads to a large changes in the quadrupole transition strength dB(E2)/dE at low energies (see Fig. 2.9).
Fig. 2.10 summarizes the different steps in our method. We divide the hyper-radial space
into four regions. The Hyperspherical Harmonics R-matrix expansion is first applied in a
small box of radius ρm (∼ 50 fm). We then use the R-matrix propagation method to much
larger distances. All Coulomb couplings are included from ρm to ρscreen (∼ 800 fm). After
that, the off-diagonal couplings are screened up to ρa (∼ 3000 fm) where it is safe to perform
the asymptotic matching with Eq. (2.32).

41

Figure 2.10: Our three-body method is divided into four steps in which we employ the
R-matrix expansion, R-matrix propagation and screening technique in Hyperspherical Harmonics coordinates.

2.2.4

HHR3a

A new code named HHR3a is generated to solve Eq. (2.16) for bound states and Eq. (2.30)
for continuum states of the triple-alpha system. It is developed from the program FaCE
[33] and STURMXX [60] which are originally designed by I. J. Thompson et al. for a
core + n + n system. In FaCE, these coupled equations are solved by expanding the hyperradial wavefunction in terms of the Laguerre basis introduced in [61] :

χγ (ρ) =

an,γ Rn (ρ, ρ0 ) ,

(2.70)

n

Rn (ρ, ρ0 ) =

ρ5/2 ρ−3
0

1/2
n!
L5 (z)exp(−z/2) ,
n
(n + 5)!

(2.71)

where z = ρ/ρ0 . ρ0 is a scaling radius and Ln (z) is the associated Laguerre polynomial
of the order n = 0, 1, 2, ... This basis forms an orthonormal set with respect to the weight
function ρ5 :
∞
0

Rn (ρ, ρ0 )Rn (ρ, ρ0 )ρ5 dρ = δn,n .

42

(2.72)

Using this basis, we can rewrite the coupled equation Eq. (2.16) in matrix form:

Hv = Ev.

(2.73)

Then, a diagonalization method is used to obtain eigenenergies E and eigenvectors v from
this set of linear equations. The boundary condition implemented in FaCE for bound state
calculations is that χγ (ρ) vanishes at both ρ = 0 and ρ −→ ∞. This is ensured by the
properties of the Laguerre basis in Eq. (2.71).
The hyperspherical harmonics R-matrix (HHR) method is implemented in STURMXX,
used to solve the coupled equations Eq. (2.30) for positive energy E. First, the coupling
matrix from FaCE is fed to STURMXX which then uses this information to construct an
R-matrix expansion basis. The coupled equations Eq. (2.30) are solved in a small R-matrix
box by expanding its solutions in terms of that orthonormal basis set built on the diagonal
couplings and an R matrix is calculated at the boundary of the box. We then propagate
this R matrix to larger radius by the use of the R-matrix propagation method (see Sec. 2.2
for details). At a large enough distance, the asymptotic matching is performed to obtain
the three-body continuum wavefunction. STURMXX operates very well for a core + n +
n system, providing fast converged and highly stable solutions of the scattering problem.
Our current problem, involving three identical alpha particles, has different symmetries and
is more difficult because of the long-range Coulomb interaction. Therefore, considerable
modifications are introduced to the new code HHR3a in order to tackle this problem. First
of all, HHR3a is corrected for the symmetry properties given by the three identical bosons.
As we discuss in Sec. 2.1.4, only partial waves with lx = even will contribute to the total
wavefunction. Here, lx are the relative angular momentum the two identical alpha particles

43

being interchanged. This constraint is applied for the wavefunction in each of the three Jacobi
coordinates. Because the wavefunction is symmetric for any pair of interchanged particles,
we are able to reduce three Faddeev components to a set of coupled channel equations
in one Jacobi coordinate. We thus highly increase the computational efficiency and fully
take into account the symmetrization of the system. In addition, the correct Coulomb
asymptotic wavefunctions Eq. (2.32) are used instead of plane waves in STURMXX, since
our problem involves charged particles. The propagation method is put to work in our
new version of the code. Because for the triple-alpha problem, the R-matrix propagation
method is needed for very large radii, a new technique of screening the off-diagonal couplings
(see Sec. 2.2.3) is added to the HHR3a code in order to improve the numerical stability of
the calculation. When the off-diagonal couplings die off we can perform a matching of the
scattering wavefunction with the correct asymptotic behavior Eq. (2.32). We also implement
the local representation of the wavefunction propagation in HHR3a to test the validity of
this method.

44

Chapter 3
Three-body reaction rate
Nuclear burning is the main source of energy for a star to maintain its active existence. It is
not only the Q-value of a given nuclear reaction but also its reaction rate that is important
to determine the amount of energy produced in a star. The nuclear reaction rate is defined
as the number of reactions per unit time and unit volume. This quantity is well understood
for a two-body reaction. It certainly depends on the nuclear cross section σ, which measures
the probability for this reaction to occur with respect to any pair of interacting particles.
The reaction rate per particle pair is defined as σv. However in a stellar plasma, the relative
velocity v is not the same for all pairs but follows a thermal distribution, namely the MaxwellBoltzmann distribution P (v):

P (v) dv = P (E) dE
2

mab 3/2 − mab v
e 2kT 4πv 2 dv ,
=
2πkT
√ −E
2
1
=√
E e kT dE ,
π (kT )3/2

45

(3.1)

where mab is the reduced mass of the two interacting particles a and b; k is the Boltzmann
constant and T is the stellar temperature. The reaction rate per pair of interacting particles
is obtained by averaging the cross section over the thermal distribution of relative velocities
vP (v):
∞

σv =

∞

v P (v) σ(v)dv =
0

=

v P (E) σ(E)dE ,
0

1/2
1
8
πmab
(kT )3/2

∞

E

−
E σ(E) e kT dE .

(3.2)

0

At high temperatures, Eq. (3.2) can be used to evaluate the triple-alpha rate since this
reaction is considered as two two-body problems [10]. In the low temperature regime where
the three-alpha direct capture mechanism is dominant, Eq. (3.2) is no longer applicable. In
this chapter we will present a detailed derivation of the three-body reaction rate formula.

3.1

Mawell-Boltzmann energy distribution for a threebody system

Since the astrophysical reaction rate is calculated by integrating the energy dependent rate
over the Mawell-Boltzmann energy distribution, it is important to derive the formula of this
distribution for a three-body system. We define v1 as the relative velocity of particle 1
and 2, and v2 as the relative velocity between particle 3 and the center of mass of system
12. The probability of finding particle 1, 2, 3 with these relative velocities in the range of
(v1 , v1 + dv1 ) and (v2 , v2 + dv2 ) can be expressed as:
2
mµ1 3 − mµ1 v1
2
P (v1 ) dv1 P (v2 ) dv2 =
e 2kT
2πkT

46

2
3
mµ2 2 − mµ2 v2
e 2kT dv1 dv2 ,
2πkT

(3.3)

where k = 8.617 × 105 eV K−1 is Boltzmann constant and T is the temperature. µ1 and µ2
m

are defined in terms of the ratio between the mass of particle i and the scaled mass, Ai = mi :
A1 A2
,
A1 + A2
A (A + A2 )
µ2 = 3 1
.
A1 + A2 + A3

µ1 =

(3.4)
(3.5)

The conjugate Jacobi momenta corresponding to x and y coordinates are defined as:

p
mµ v
√
kx = √ 1 = √ 1 1 = m µ1 v 1 ,
µ1
µ1
p
mµ v
√
ky = √ 2 = √ 2 2 = m µ2 v 2 .
µ2
µ2

(3.6)
(3.7)

Using Eq. (3.6) and Eq. (3.7) we can rewrite Eq .(3.3) as:
2

2

1 − kx − ky
1
e 2mkT e 2mkT dkx dky .
P (v1 ) dv1 P (v2 ) dv2 =
(2πkT )3 m3

(3.8)

2
2
The fact that kx + ky = κ2 allows us to simplify the above equation to:

P (v1 ) dv1 P (v2 ) dv2 =

1
1 − κ2
e 2mkT κ5 sin2 θκ cos2 θκ dκ dθκ dΩkx dΩky .
3 m3
(2πkT )

(3.9)

Performing the integration of Eq. (3.9) over all the angles, we obtain the Maxwell Boltzmann
distribution for the three-body system:

P (E) =

1
−E
e kT E 2 .
2(kT )3

47

(3.10)

Comparing to Eq. (3.1) which is obtained for a two-body system, the three-body MaxwellBoltzmann distribution has very different energy and temperature dependence. In the next
section we will apply this new distribution to derive a general formula for the three-body
reaction rate.

3.2

The reaction rate formalism for three-body capture

We first consider a reaction in which three particles 1, 2, 3 fuse to create two new particles
4, 5:
1+2+3

4 + 5.

(3.11)

The ratio of the reaction rates between the direct and inverse process is defined as in [62]:

123
g g
= 4 5
45
g1 g2 g3

3

2
A4 A5
A1 A2 A3

3

2π 2 2 1 + ∆123 − Q
e kT ,
mkT
1 + δ45

(3.12)

where ∆123 = δ12 + δ23 + δ31 + 2δ123 to account for the effects introduced by the number of
identical particles in the reactions; the Q-value of the reaction is defined as Q = E123 − E45 ;
gi is the degenerate factor for particle i which can be calculated from its total spin Ji :
gi = 2Ji +1; , k, Ai and m have their usual meanings as defined before; T is the temperature
of the environment.
The three-body 123 and two-body 45 reaction rates are defined using Eq. (3.10) and

48

Eq. (3.1) respectively for averaging over the Maxwell-Boltzmann distribution:

123 =

1
2(kT )3

45 =

2
8
πm45

R123 (E123 ) e
1

E
− 123
kT

2
E123 dE123 ,

(3.13)

E45

1

−
σ(E45 ) e kT E45 dE45 .

3
(kT ) 2

(3.14)

In these equations, E123 and E45 are the three-body and the two-body relative energies,
respectively. σ(E45 ) is the reaction cross section of system (4, 5) which is equivalent to
m m

R123 (E123 ) in the three-body case and m45 = m 4 5 is the reduced mass of a two-body
4 +m5
system (4, 5).
Using the Q-value definition, we rewrite Eq. (3.13) in terms of an integration over the
two-body energy E45 :

123 =

E45

Q
1
−
e kT
2(kT )3

−
2
R123 (E123 ) e kT E123 dE45 .

(3.15)

Inserting Eq. (3.15) and Eq. (3.14) into Eq. (3.12) and after some simple algebraic reduction,
we obtain:
E45

−
2
R123 (E123 ) e kT E123 dE45
E45

−
σ(E45 ) e kT E45 dE45

3

3

g g
= 4 5
g1 g2 g3

2 1 + ∆123
A4 A5
A1 A2 A3
1 + δ45

2π 2 2
mkT

1

2
32
.
πm45

(3.16)
In order for Eq. (3.16) to be valid with all reactions as described by Eq. (3.11), the integrands
themselves have to satisfy:

g g
2
R123 (E123 ) E123 = 4 5
g1 g2 g3

3

3

2 1 + ∆123
A4 A5
A1 A2 A3
1 + δ45

2π 2 2
mkT

1

2
32
σ(E45 )E45 .
πm45

(3.17)

49

This formula is correct when all 1, 2, 3, 4, and 5 are particles. In our triple-alpha reaction, one
of the products is a photon, so the formula for this case is slightly different. We now need to
transform the cross section σ45 for the case of two real particles into the photo-disintegration
cross section where one of these particles is a photon. Let us consider a two-body process:

4+5

γ + 0.

(3.18)

The ratio between the direct and inverse cross section is defined in [9]:
2
Eγ
gγ g0
σ(E45 )
(1 + δ45 ) ,
=
σγ (Eγ )
g4 g5 2 m45 c2 E45

(3.19)

where all the variables in this equation have the same meaning as described before.
Inserting Eq. (3.19) into Eq. (3.17) we have:

gγ g0
2
R123 (E123 ) E123 =
g1 g2 g3

3

3

2
A4 A5
(1 + ∆123 )
A1 A2 A3

2π 2 2
mkT

1

2
Eγ
2
32
σγ (Eγ ) .
πm45
2m45 c2
(3.20)

Expressing m45 in terms of A4 and A5 , we could further reduce Eq. (3.20) to a simpler form:
3

R123 (E123 ) = 2
c

2
gγ g0 Eγ
(1 + ∆123 )
σγ (Eγ ) .
3
2
g1 g2 g3 E123
3 (µ µ ) 2
m
1 2

8π

(3.21)

Eq. (3.21) is the energy dependent rate for the fusion of three particles 1, 2, 3 into particle
0 with the emission of a photon γ:

1+2+3

50

γ + 0.

(3.22)

We denote E = E123 for convenience and then integrate R123 (E) over the three-body
Maxwell-Boltzmann distribution Eq. (3.10) to obtain the astrophysical three-body reaction
rate:
2
R123 = NA

∞

1
2(kT )3

0

E

−
e kT E 2 R123 (E) dE .

(3.23)

Inserting Eq. (3.21) into Eq. (3.23) and using the expression for the reaction Q-value (Q =
E − Eγ ), we have:

2
R123 = (1 + ∆123 ) NA

Q 3
1
−
e kT 2
2(kT )3
c

8π
m3 (µ1 µ2 )

gγ g0
3 g g g
2 1 2 3

∞
|Q|

Eγ

−
2
e kT Eγ σγ (Eγ ) dEγ .

(3.24)
Since the photon has only 2 polarization directions, we have gγ = 2. We replace (1 + ∆123 )
by p! where p is the number of identical particles. Eq. (3.24) is now reduced to:

R123 =

2
p! NA

3

c2

8π
m3 (µ1 µ2 )

Q
1
g0
−
e kT
3 g g g
(kT )3
2 1 2 3

∞

e
|Q|

−

Eγ
kT

2
Eγ σγ (Eγ ) dEγ .

(3.25)

Eq. (3.25) relates the radiative capture reaction rate to the photo-disintegration cross section
σγ which depends on the electromagnetic transition strength dB(Eλ)/dE through which the
reaction occurs [11, 27, 46, 63]:

σγ =

(2π)3 (λ + 1)
λ[(2λ + 1)!!]2

Eγ 2λ−1 dB(Eλ)
.
c
dE

(3.26)

The missing ingredient to calculate the three-body reaction rate (Eq. (3.25)) is the electromagnetic transition strength dB(Eλ)/dE. In the next section, we will discuss how to
calculate this quantity for a three-body system in the HH representation.

51

3.3

Three-body quadrupole transition function

At low energies, the triple-alpha reaction proceeds through a quadrupole transition from the
0+ continuum to the 2+ bound state in 12 C. The quadrupole transition operator E2 for any
1
three-body system has the general form in coordinate space:

a
a
a
E2m = e[Z1 a2 Y2m (ˆ1 ) + Z2 a2 Y2m (ˆ2 ) + Z3 a2 Y2m (ˆ3 )] ,
3
2
1

(3.27)

Where, Zi is the charge of particle i and Y is the spherical harmonic function and m =
−2, −1, 0, 1, 2 are the projections of angular momentum l = 2; a1 , a2 , and a3 are the distances
from the center of mass of the three-body system to each particle as seen in Fig. 2.2. Vectors
ai can be expressed in terms of the scaled Jacobi coordinates (x, y):
√

a1
a2
a3

√
µ1
µ2
=−
x−
y,
A1
A1 + A2
√
√
µ1
µ2
x−
y,
=
A2
A1 + A2
√
µ2
=
y.
A3

(3.28)
(3.29)
(3.30)

For a given vector relationship a = a1 + a2 , we have the following properties for spherical
functions as proven in [51]:
1

al C

a
lm (ˆ )

=
λ,µ

2 l−λ λ
2l!
a1 a2 Cl−λ,m−µ (ˆ1 ) Cλ,µ (ˆ2 ) l − λ m − µ λ µ|l m ,
a
a
2λ!2(l − λ)!

(3.31)
−1/2
Ylm .
where Clm = 2l+1
4π

Applying Eq. (3.31) to Eq. (3.28), Eq. (3.29) and Eq. (3.30) we can express the transition
operator E2 in terms of hyperspherical coordinates (ρ, θ). For the triple-alpha system in
52

which three particles are identical: A1 = A2 = A3 = A and Z1 = Z2 = Z3 = Z, the
quadruple transition operator reduces to much simpler form (see Appendix B for details):

E2m =

eZ
(ρ sinθ)2 Y2m (ˆ ) + (ρ cosθ)2 Y2m (ˆ ) .
x
y
A

(3.32)

A quadrupole transition from the initial continuum state of total spin L and momentum κ
to the final bound state of spin L is characterized by the strength function:

dB(E2, L → L )
=
dκ dΩκ
5

L M |E2m | LM ; κ

2

,

(3.33)

mM M

where dΩκ = sin2 θκ cos2 θκ dθκ dΩkx dΩky . Because our system is composed of spin zero
5
particles, the total angular momentum L also indicates the total spin. The total B(E2) is
computed by integrating over all the continuum states |LM, κ :

B(E2, L → L ) =
√

Applying the relationship κ =

2mE

B(E2, L → L ) =

dB(E2, L → L ) 5
κ dκ dΩκ .
5
dκ dΩκ
5

(3.34)

to Eq. (3.34) we obtain:

dB(E2, L → L ) E 2
dκ dΩκ
2
5

2m 3
2

dΩκ dE .
5

(3.35)

Using Eq. (3.33) and comparing Eq. (3.35) with the following equation:

B(E2, L → L ) =

dB(E2, L → L )
dE ,
dE

53

(3.36)

we obtain the energy dependent strength function:
dB(E2, L → L )
E2
=
dE
2

2m 3

2

L M |E2m | LM ; κ

2

dΩκ .
5

(3.37)

mM M

We use the hyperspherical expansion of the wavefunctions in Eq. (2.8) and Eq. (2.28) to
expand the squared modulus in Eq. (3.37):
i i
lx ly
ϕ i
K

L M |(E2)m | LM ; κ =
i i
K i lx ly Klx ly K lx ly

∗

Yl i ⊗ Yl i
x
y LM κ

i i
× K lx ly L M |(E2)m | Klx ly K i lx ly LM s ,

(3.38)

where the subscript s denotes the spatial part (ρ, θ, Ω5 ) in the wavefunction expansions.
The integration in Eq. (3.37) is then simplified using the following properties of spherical
lx ly

harmonics functions Ylm and hyper-angular functions ϕK :

dΩkx Y ∗ i (Ωkx ) Yli mi (Ωkx ) = δli li δmi mi ,
i
lx mx
x 1x x 1x
1x 1x

(3.39)

dΩky Y ∗ i (Ωky ) Yli mi (Ωky ) = δli li δmi mi ,
i
lx mx
y 1y
y 1y
1y 1y

(3.40)

i i
lx ly

i i
lx ly

sin2 θκ cos2 θκ dθκ (ϕ i (θκ ))∗ ϕ i (θκ ) = δK i ,K i .
K
K
1

(3.41)

1

The integration in Eq. (3.37) is reduced to (see Appendix B for details):

L M |(E2)m | LM ; κ

2

dΩκ
5
2

i i
K lx ly L M |E2m | Klx ly K i lx ly LM s .

=
i i
K i lx ly Klx ly K lx ly

54

(3.42)

Inserting Eq. (3.42) into Eq. (3.37) and employing the Wigner-Eckart theorem from [51] we
obtain:
dB(E2, L → L )
E2
=
dE
2

ˆ2
2m 3 L
2
ˆ
L2
2

i i
K lx ly L ||E2m || Klx ly K i lx ly L s

×

,

(3.43)

i i
K i lx ly Klx ly K lx ly

ˆ
where L2 = 2L + 1. Let’s first calculate the matrix element in Eq. (3.43) using Eq. (3.32):

eZ
i i
K lx ly L (ρ sinθ)2 Y2 (ˆ ) Klx ly K i lx ly L s
x
A
eZ
i i
+
K lx ly L (ρ cosθ)2 Y2 (ˆ ) Klx ly K i lx ly L s .
y
A

i i
K lx ly L ||E2m || Klx ly K i lx ly L s =

(3.44)

We denote M.E.1 and M.E.2 as the first and second term on the r.h.s of Eq. (3.44). M.E.1
can be factorized into two terms of which one contains the hyper-variable dependence and
the other is angular momentum dependent:

M.E.1 =

eZ
i i
K lx ly (ρ sinθ)2 Klx ly K i lx ly
A

lx ly L ||Y2 (ˆ)|| lx ly L .
x

(3.45)

The first bra-ket term on the r.h.s of Eq .(3.45) is calculated by taking an integral over the
hyper-radial and hyper-angular parts of the wavefunctions:

K lx ly

(ρ sinθ)2

i i
Klx ly K i lx ly

55

=

2mE −5/4
2

Iρ I1,θ ,

(3.46)

where we define:

Iρ =
I1,θ =

i i
K i lx ly

dρ χ∗ l l (ρ) ρ2 χKl l (ρ) ,
K xy
xy

(3.47)

lx ly
lx ly
(θ))∗ sin2 θ ϕK (θ) .
K

sin2 θ cos2 θ dθ (ϕ

(3.48)

The second bra-ket term on the r.h.s of Eq. (3.45) can be explicitly calculated using angular
momentum algebra [51] as:

lx ly L ||Y2 (ˆ)|| lx ly L = δ(ly , ly )
x




L 2   lx 2 lx 

.


0 0 0
x lx ly
(3.49)



 L
5 ˆˆ ˆ
2+ly +L+2lx
Llx lx (−1)

4π
 l

Performing a similar calculation for the second term on the r.h.s of of Eq. (3.44) (M.E.2),
we obtain the final result for a three-body quadrupole transition strength function:
2

dB(E2, L → L )
=
dE

m 1 ˆ2
√ L
2 2 E

eZ
A

2

Iρ I1,θ A1 + I2,θ A2

,

i i
K i lx ly Klx ly K lx ly

(3.50)
where

Iρ =

i i
K i lx ly

dρ χ∗ l l (ρ) ρ2 χKl l (ρ),
K xy
xy

(3.51)

I1,θ =

sin2 θ cos2 θ dθ (ϕ

lx ly
(θ))∗
K

sin2 θ ϕK (θ),

I2,θ =

sin2 θ cos2 θ dθ (ϕ

lx ly
(θ))∗
K

cos2 θ ϕK (θ),

56

lx ly

lx ly

(3.52)
(3.53)

with coefficients:

A1 = δ(ly , ly )



 L
5 ˆˆ
2+ly +L+2lx
lx lx (−1)

4π
 l

x

A2 = δ(lx , lx )

L
lx



 L
5 ˆˆ
2+lx +L +ly +ly
ly l (−1)

4π y
 l

y



2   lx



ly
0


L 2 



ly lx


2 lx 
(3.54)
,
0 0

ly 2 ly 
 . (3.55)
0 0 0

The quadrupole transition strength in Eq. (3.50) is used to obtain the triple-alpha reaction
rate. Calculation of Eq.( 3.50) requires correct solutions of the 2+ bound state and the 0+
1
continuum states of 12 C. The convergence study will focus of these quantities of our interest.

57

Chapter 4
Results

4.1

Interactions

Although there has been progress in many body techniques as discussed in Sec. 1.3, the nonresonant continuum description of 12 C at a microscopic 12 body level is not yet attainable.
Since the alpha particle is the strongest bound nucleus among few-nucleon systems, it is
reasonable to formulate the 12-body problem of 12 C as an effective three-body system. This
argument is supported by the dominant triple-alpha structure found in the microscopic
calculation of both the 2+ bound state and the Hoyle resonant state of 12 C [31]. In order to
1
obtain a good description for this three-body system, a well-constrained effective two-body
interaction for any pair of alpha particles is necessary. Ali and Bodmer [64] had developed
a phenomenological alpha-alpha potential for s, d and g waves by fitting the experimental
phase shifts at low energies. We employ the effective alpha-alpha interaction from [65] which
modified the Ali-Bodmer potential to exactly reproduced the 8 Be s-wave resonance:

2
2
2
2
ˆ
ˆ
Vαα = (125Pl=0 + 20Pl=2 ) e−r /1.53 − 30.18 e−r /2.85 .

58

(4.1)

This l-dependent potential has a very strong repulsive core in the s-wave to simulate the
Pauli exclusion principle in the alpha-alpha system. It reproduces successfully the low energy phase shifts as well as the 0+ resonant state of 8 Be at 0.092 MeV.
Another interaction we investigated is a local Gaussian alpha-alpha interaction independent of energy and angular momentum [66]. Even though this interaction fits the scattering
data very well, it creates deep s-wave states which are forbidden and need to be projected
out of the model space. Although the handling of forbidden states in computing the 12 C(2+ )
1
bound state was possible, this additional difficulty was hard to overcome within the threebody continuum.
In order to describe the triple-alpha system, a three-body force V3b is introduced to the
sum of three pairwise alpha-alpha interactions Vαα . It is important to understand the physical interpretation of an effective three-body force used in our calculation. The three-body
force was first considered in nuclear physics to explain the discrepancy between the measured
binding energy of the triton and the theoretical calculations using different realistic nucleonnucleon interactions [67, 68]. Although the triton is one of the simplest systems in nuclear
physics, treating nucleons as fundamental particles that interact via a two-body potential
gave rise to such discrepancy. An effective three-body force was then needed to account for
this simplification. Our triple-alpha system is largely reduced since each alpha particle is
composed of four nucleons. We introduce a three-body force in our model to compensate for
the reduction of the model space, allowing us to reproduce the experimental observable of a
2 2
three-alpha system. In this work, a three-body interaction V3b (ρ) = V0 e−ρ /ρ0 is added to

the Hamiltonian to reproduce the experimental bound state energy of the 2+ state as well
1
as the Hoyle resonant state 0+ of 12 C (similarly to [65]). Since the effective three-body force
2
originates from the fact that the alpha particle is not a fundamental particle, we expect the
59

-2.835

E(21 ) [MeV]

-2.85

+

-2.865

-2.875 MeV

-2.88

-2.895
-2.91
-16.1 -16.05

-16

-15.95 -15.9 -15.85 -15.8

V0 [MeV]
Figure 4.1: Dependence of the three-body binding energy of the 2+ bound state of 12 C on the
1
2

2

depth of a three-body force V3b (ρ) = V0 e−ρ /ρ0 . The maximum hyper-angular momentum
Kmax is taken as 20 for which the results are converged.
range of this interaction to be approximately the overall size of the three-alpha system. We
therefore take ρ0 = 6 fm as the position where the three alpha particles touch each other
and adjust the depth V0 of the three-body interaction to measurements.
The blue solid curve in Fig. 4.1 represents the dependence of the 12 C(2+ ) binding
1
energy on the three-body potential depth V0 for a fixed alpha-alpha interaction. The red
dash line is the experimental value −2.875 MeV of this state [69]. As can be seen, we need
V0 = −15.94 MeV to reproduce the measured excitation energy of the 2+ bound state. The
1
same story applies to the Hoyle state. Fig. 4.2 shows a dependence of the resonant energy
on the depth of a three-body force V0 (blue curve). The red curve is again the experimental
energy (0.3795 MeV) [69]. A depth of V0 = −19.46 MeV reproduces the measured Hoyle
state. This interaction is then used to calculate all the 0+ continuum states. In Table 4.1,
we summarize the parameters for the three-body interactions used in our calculations.
60

0.395

+

E(02 ) [MeV]

0.39
0.385

0.3795 MeV

0.38
0.375
0.37
-19.5

-19.4

-19.45

-19.35

-19.3

V0 [MeV]
Figure 4.2: Dependence of the Hoyle resonant energy of 12 C on the depth of a three-body
2 2
force V3b (ρ) = V0 e−ρ /ρ0 . The maximum hyper-angular momentum Kmax is taken as 26 for
which the results are converged.

Table 4.1: The three-body interactions to reproduce the experimental energies [69] of the
2+ bound state and 0+ resonant state in 12 C.
1
2

JΠ

V0 (MeV) ρ0 (fm) E(MeV)

2+
1

-15.94

6

-2.875

0+
2

-19.46

6

0.380

In addition to the nuclear interaction we include the Coulomb potential:

Coul
Vαα (r)

=

Z 2 e2

×






3
2

−

r2
2
2r
Coul

 1



1
rCoul

r ≤ rCoul
(4.2)
r ≥ rCoul .

r

61

Here, Z is the charge of an alpha particle, r is the distance between two alpha particles, and
rCoul is the Coulomb radius which is taken as twice the alpha particle radius rCoul = 2.94
fm.

4.2

The 2+ bound state
1

A three-body bound state is obtained by solving the set of coupled channels equations
(Eq. (2.16)) with a boundary condition that the wavefunction goes to zero at large ρ [33].
These coupled equations are solved by expanding the hyper-radial wavefunction in terms of
the Laguerre basis [61] which forms a complete and orthogonal set (Sec. 2.2.4). The matrix diagonalization is then applied to obtain eigenenergies and eigenvectors. We choose the
scaling radius ρ0 in Eq. (2.71) for the Laguerre polynomial to be 0.3 fm. The standard mass
m is taken as the nucleon mass (939 MeV). The hyper-angular integrations are performed
using Gauss-Jacobi quadrature on a grid with 70 points. The hyper-radial integrations are
performed using Gauss-Laguerre quadrature on a grid with 180 points. 40 basis polynomials
are used in the hyper-radial wave expansion Eq. (2.70). A actual input file with all the
necessary variables for a bound state calculation can be seen in Appendix C.
Fig. 4.3 shows the convergence of the three-body binding energy of the 2+ bound state
1
of 12 C with the size of the model space represented here by the maximum hyper-angular
momentum Kmax . For a given Kmax we can calculate how many channels are included in
the wavefunction expansion (Eq. (2.8)), e.g., Kmax = 20 produces 36 channels in the expansion. From Fig. 4.3, the bound state energy starts converging at Kmax = 12. Therefore,
energy convergence is guaranteed by choosing Kmax = 20. The bound state energy of the
2+ state of 12 C converges to the experimental value E = −2.875 MeV [69] with respect to

62

-2.6

+

E(21 ) [MeV]

-2.65
-2.7
-2.75
-2.8
-2.85
4

6

8

10

12

Kmax

14

16

18

20

Figure 4.3: Dependence of the three-body binding energy of the 2+ bound state of 12 C on
1
the maximum hyper-angular momentum Kmax .
the three-alpha threshold.
The convergence of the rms radius of the 12 C(2+ ) bound state with respect to Kmax is
1
shown in Fig. 4.4. Again, the result converges very well for Kmax = 20. Our final value for
the rms radius is 2.459 fm which is close to 2.50 fm obtained in a microscopic calculation
done by Chernykh et al.[31]. This agreement validates the three-body approximation for the
2+ bound state.
In Table 4.2, we compare the three-body binding energy and the rms radius of the 2+
1
bound state of 12 C obtained from our three-body calculation with experimental data and
other studies. Our results reproduce the experimental data and agree very well with an
independent three-body study using the hyperspherical adiabatic method [26]. The microscopic calculations in [31] overestimate the binding energy of the 12 C(2+ ) state but provide
1
a similar rms radius.

63

2.46

+

rrms(21 ) [fm]

2.465

2.455

2.45
6

8

10

12

14

16

18

20

Kmax
Figure 4.4: Dependence of the rrms radius of the 2+ bound state of 12 C on the maximum
1
hyper-angular momentum Kmax .
Table 4.2: Comparison of the three-body binding energy and the rms radius of the 2+
1
bound state of 12 C obtained from our three-body calculation with other studies: a threebody hyperspherical adiabatic method [26], a microscopic Fermionic Molecular Dynamics
method [31] and experimental data [69].

This work

[26]

[31]

Exp [69]

E (MeV)

-2.875

-2.875

-3.74

-2.875

rrms (fm)

2.459

2.45

2.50

In Fig. 4.5, we plot the rrms radius as a function of the maximum radius ρmax of the
calculated 2+ bound state wavefunction. As can be seen, the rrms no longer depends on the
1
maximum radius of the bound state wavefunction for ρmax > 20 fm. All the calculations
through out this work employ a bound state wavefunction up to 200 fm, thus the correct

64

+

rrms(21 ) [fm]

2.46

2.459

2.458

2.457

2.456

20

30

ρmax [fm]

40

50

Figure 4.5: The rrms radius as a function of the maximum radius of the box in which the
12 C(2+ ) bound state wavefunction is calculated.
1
contribution from the long-range part of the wavefunction to the rate at low temperatures
is ensured.
Given that the wavefunction of the 12 C(2+ ) bound state is available in our work, it
1
is interesting to look at the cluster structure of this state. We next construct the density
distribution function from the 2+ wavefunction to study its structure:
1

P (r, R) =

|Ψ(r, R)|2 dˆ dR .
r ˆ

(4.3)

Here r is the radius between the two alpha particles and R is the distance from its center of
mass to the third alpha. In Fig. 4.6 we present the density distribution P (r, R) for the 2+
1
bound state of 12 C in contour color style. We observe a maximum in the density distribution
at r ∼ R ∼ 3 fm. This indicates that the dominant configuration for the 12 C(2+ ) bound
1
state is a nearly equilateral triangle in which each pair of particles is ∼ 3 fm apart. Since an
65

8
0 .0
0 .0
0 .0
0 .0
0 .0
0 .0
0 .0
0 .0
0 .0

R [fm ]

6

4

7 5
6 6
5 7
4 8
3 9
3 0
2 1
1 2
0 3

2

0
0

2

4

6

8

r [fm ]
Figure 4.6: The density distribution for the 2+ bound state of 12 C.
1

Figure 4.7: The dominant three-alpha cluster in the 12 C(2+ ) bound state: three alpha
1
particles touch each other and form a equilateral triangle.

66

alpha particle has a radius of 1.47 fm, three alpha particles almost touch each other in this
dominant configuration. This configuration is described in Fig. 4.7.

4.3

The 0+ Hoyle state
2

The three alpha particle system is driven by a strong Coulomb interaction of which the
off-diagonal couplings are long-range in the HH representation. Therefore, solving the coupled channels equations (Eq. (2.30)) for positive energy E is not a trivial task. The HHR
method which combines the R-matrix expansion and R-matrix propagation in the HH basis
is employed to overcome this difficulty. Details of this method are presented in Sec. 2.2.
The numerical calculations are performed using our HHR3a code (see Sec. 2.2.4). We use
an R-matrix box size of 50 fm with 50 poles included in the R-matrix expansion and a radial step size of 0.05 fm. The hyper-angular integrations are performed using Gauss-Jacobi
quadrature on a grid with 100 points. As mentioned before, the fixed logarithmic derivative
in the R-matrix calculation is selected to be zero. Our converged calculation is obtained with
the maximum hyper-momentum Kmax = 26, the screening radius ρscreen = 800 fm and the
screening difusseness ascreen = 10 fm. The convergence of these quantities are thoroughly
tested. The R matrix is propagated out to 3000 fm, sufficiently large to perform the Coulomb
asymptotic matching. A detailed input file for the calculation of 0+ continuum states can
be viewed in Appendix C.
The 0+ continuum states of 12 C are calculated for E = 0.01 − 0.5 MeV. The quadrupole
transition strength Eq. (3.37) is then constructed using the 2+ bound state wavefunction and
1
the 0+ resonant and non-resonant continuum states. Therefore, in our calculations we treat
the resonant and non-resonant process on the same footing. Fig. 4.8 shows the quadrupole

67

10

4

10
10

2

-1

dB(E2)/dE [e fm MeV ]

10

10
10
10
10

0

-10
-20
-30
-40
-50
-60

10 0

0.1

0.2

0.3

0.4

0.5

E [MeV]
Figure 4.8: The quadrupole transition strength dB(E2)/dE as a function the three-alpha
relative energy E.
transition strength as a function of the relative energy for the three interacting alpha particles. The curve peaks at the measured Hoyle resonance energy (E = 0.38 MeV). As we
expect, the strength function decreases as it approaches the lower energy regime. Below 0.05
MeV there is a sharp reduction of the transition strength below which the formation of 12 C
becomes unlikely.
It is important to discuss our method to extract the Hoyle resonant energy as well as
its convergence. We repeated the calculation of the quadrupole transition strength function around the 0+ resonant energy in Fig. 4.8 for different values of Kmax . The results
2
are presented in Fig. 4.9. As can be verified, the position of the peak in the quadrupole
transition strength curve is sensitive to the size of our model space. As Kmax increases the
discrepancy between the two consecutive curves is reduced which indicates a converged behavior. For each Kmax , the resonant energy is obtained from the maximum of the transition

68

0

Kmax=12
Kmax=14
Kmax=16
Kmax=18
Kmax=20
Kmax=22
Kmax=24
Kmax=26
Kmax=28

-1

dB(E2)/dE [e fm MeV ]

10

-2

2

4

10
10

-4

10

-6

10

-8

0.36 0.38

0.4

0.42 0.44 0.46 0.48

0.5

E [MeV]
Figure 4.9: The quadrupole transition strength function dB(E2)/dE around the Hoyle resonant energy is calculated with different values of hyper-momentum Kmax .
strength function and plotted in Fig. 4.10. A clearer convergence pattern of the 12 C(0+ )
2
Hoyle resonant energy with respect to the size of the model space is depicted in this graph.
As expected, we find that the convergence for the Hoyle resonant state is slower than for the
bound state. The result begins to converge at Kmax = 26 and this is the point we choose to
fit the three-body force to reproduce the experimental energy of 0.38 MeV [69] with respect
to the three-alpha threshold. One might wonder why can we use the phase-shift analysis
to extract the resonant energy instead? Because the problem involves so many channels, it
is not possible for us to continuously link the eigenphases at various energies. In addition
one cannot associate the eigenphases with certain hyper-spherical harmonics channels. A
specific algorithm would need to be developed to tackle those challenges.
In order to estimate the uncertainty in extracting the Hoyle resonant energy, the data in
Fig. 4.10 is fitted to an exponential function y = A + Be−Cx . When Kmax goes to infinity

69

0.43

+

E(02 ) [MeV]

0.42
0.41
0.4
0.39
0.38
0.37
12

14

16

18

20

22

24

26

28

30

Kmax
Figure 4.10: Dependence of the three-body energy of the 0+ Hoyle state of 12 C on the
2
maximum hyper-angular momentum Kmax .
this function approaches its converged value which is then compared to the value of the
energy at Kmax = 26 to obtain an uncertainty of 4%.
By using the wavefunction at the Hoyle resonance energy we are able to construct
the density distribution function as in Eq. (4.3) for the 0+ state in 12 C. This quantity il2
lustrates the spatial structure of the Hoyle state, a topic of tremendous interest for many
years [31, 70, 71, 72, 73]. Fig. 4.11 depicts the density distribution in a color contour plot.
The dominant configuration for the Hoyle state, corresponding to the biggest maxima in the
density plot, is the prolate triangle in which two alpha particles are near each other (∼ 3.7
fm) and further away from the third alpha (∼ 5.2 fm). This result is in agreement with the
finding in [65, 73] and the alpha cluster configuration for this case is presented in Fig. 4.12a.
Thus, a 8 Be-α like structure dominates the 12 C Hoyle resonant state as expected. We also
observe two smaller maxima for the density distribution function. One of them supports

70

1 0

0 .0
0 .0
0 .0
0 .0
0 .0
0 .0
0 .0
0 .0
0 .0

8

R [fm ]

6

2 7
2 4
2 1
1 8
1 5
1 2
0 9
0 6
0 3

4
2
0
0

2

4

6

8

1 0

r [fm ]
Figure 4.11: The density distribution for the 0+ resonance state of 12 C.
2

Figure 4.12: Three triple-alpha configurations of the 0+ resonance in 12 C: the prolate triangle
2
(a), the oblate triangle (b) and the equilateral triangle (c).

71

an oblate triangle configuration of three alpha particles, where two alpha particles are ∼ 7
fm apart and its center of mass is ∼ 1.5 fm from the third alpha. This configuration forms
an almost chain-like structure as depicted in Fig. 4.12b. An ab initio calculation using lattice effective field theory method [72] also supports the existence of this three-alpha cluster
structure for the Hoyle state. The other small maximum indicates an almost equilateral
triangle with a distance of ∼ 3 fm between any pair of alpha particles. This is the same
structure obtained for the 2+ bound state and is presented in Fig. 4.12c. The weights of
1
these two latter configurations are about two times smaller than that of the prolate triangle
configuration. In conclusion, we find three different configurations for the triple-alpha in
which the prolate triangle is dominant, but the oblate triangle (almost chain-like) and the
equilateral triangle also contribute significantly to the structure of the Hoyle state of 12 C.
Note that our results are limited to the relative magnitudes of r and R only. The average
angle between these two vectors requires further investigation. Fig. 4.12 does not consider
the angular orientation of the two radius vectors.

4.4

Testing the new code HHR3a

HHR3a is developed from the programs FaCE [33] and STURMXX [60], which are originally
designed for a core + n + n problem, to introduce the correct symmetry and the right boundary conditions for the three identical charged-particle system. It is important to ensure that
the new code operates properly. When only the diagonal Coulomb couplings are considered,
Eq. (2.30) becomes Eq. (2.31) which has analytic solutions as seen in Sec. 2.1.3. We therefore can use this case to test our numerical method. In Fig. 4.13, we compare the numerical
solution of Eq. (2.31) calculated at E = 0.5 MeV for the 0+ scattering wavefunction of the

72

2
1.5
1

χ(ρ)

0.5
0
-0.5
-1

analytic
HHR

-1.5
-2
0

50

100

150

200

250

300

ρ [fm]
Figure 4.13: A numerical solution (dashed) of Eq. (2.30) at E = 0.5 MeV for the 0+ scattering
wavefunction of the triple-alpha system when only the diagonal Coulomb couplings are taken
into account is compared with the analytic solution (solid).
three-alpha system with its analytic solution expressed in terms of the confluent hypergeometric function of the first kind [74]. We observe a fine agreement between the two curves
indicating the correct behavior of the codes in the absence of off-diagonal couplings. In this
graph we only present one component of the wavefunction in the hyperspherical expansion of
which the incoming channel and the outgoing channel are characterized by the set of quantum numbers Klx ly = 000, however all other components show equal level of agreement.
We then perform a calculation for the triple-alpha rate when only the diagonal Coulomb
couplings are taken into account by employing the 0+ continuum wavefunctions above, obtained both numerically and analytically. The results are compared and plotted in Fig. 4.14.
We find that the reaction rate calculated from our numerical solutions produced by HHR3a
agrees very well with that obtain using the analytic wavefunctions. Therefore, we conclude
that the new code is able to reproduce the analytic results for a simplified three-alpha prob73

-15

6 -1

-2

〈Rααα〉 [cm s mol ]

10

-30

10

10

analytic
HHR

-45

-60

10

0.01

0.1

1

T [GK]
Figure 4.14: The triple-alpha reaction rate is computed for the diagonal Coulomb case using
the numerical solution of Fig. 4.13 (dashed) and compared to the analytic solution (solid).
lem in which only the diagonal Coulomb interaction is present.
We also test the new code for a realistic case when both nuclear and Coulomb interactions
are included. The HHR3a code produces the Hoyle resonant state at the right energy and
the results agree well with an independent study in [26] given the same interactions. We
also reproduce the results from [75] for the bound states as well as the resonant states of
12 C(0+ ).

4.5

These reinforce the validity of our new code.

Convergence and uncertainty

It is important for us to carefully study the convergence properties of our method, given
the well known difficulty of the three charged-particle problem. Since the 2+ bound state
wavefunction is easily obtained and fully converged for the effective three-alpha system when
both the nuclear and Coulomb interactions are present, in all of our tests presented in this
74

6 -1

-2

〈Rααα〉 [cm s mol ]

10
10
10
10

-20

-30

-40

Kmax=20
Kmax=26

-50

-60

10 0.01

0.10

1.00

T [GK]
Figure 4.15: The reaction rate as a function of temperature for the fictitious three-alpha
system in which only Coulomb interactions are included: Kmax = 20 (solid) and Kmax = 26
(dot-dashed).
section we will fix the bound state wavefunction and only study the convergence of the
solutions of the coupled channels equations (Eq. (2.30)) for scattering states.

4.5.1

The impact of using screening potentials

With the help of the R-matrix propagation method we are able to perform calculations out
to a very large radius and obtain stable results with Coulomb interactions only. Fig. 4.15
presents the calculated reaction rate for a fictitious three-alpha system in which only the
Coulomb interactions play a role. We perform two calculations of the rate with different
sizes of the model space Kmax = 20 and Kmax = 26 in order to study the convergence of
this system when just long-range couplings are introduced. It is important to note that the
bound state wavefunction which is employed to compute the rate is obtained using a realistic
triple-alpha model where the nuclear and Coulomb interactions are fully coupled. For the
75

4

-1

dB(E2)/dE [e fm MeV ]

10

-48 Coulomb + Nuclear

-51

2

10
10

Coulomb Only
(b)

(a)

10

no screen
screen

-54

10

10

10

-57

10

-48

-51

-54

-57

10 12 14 16 18 20 22 24 26 10 12 14 16 18 20 22 24 26

Kmax

Kmax

Figure 4.16: The quadrupole transition strength dB(E2)/dE at E = 0.01 MeV as a function
of model space Kmax when both nuclear and Coulomb interaction (a) and only Coulomb
interaction (b) are present in the three-alpha system.
0+ continuum, the R matrix is propagated out to a very large radius (3000 fm) in order
for the rate to be fully converged. This demonstrates the challenge of our problem. The
two calculations in Fig. 4.15 are different by 0.3% indicating a very good convergence of the
triple-alpha rate with respect to Kmax when only the Coulomb interaction is considered. It
is important to note that the results converge faster when no nuclear interaction is taken
into account. In fact, we only need Kmax = 20 for this case to converge while Kmax = 26 is
required for the realistic triple-alpha system.
When nuclear interactions are introduced, the very narrow resonances occur in both the
two-body and three-body systems, and lead to a degradation of our propagation technique,
increasing numerical instabilities, especially in the turning point region. This is confirmed
in Fig. 4.16 in which the quadrupole transition strength dB(E2)/dE at E = 0.01 MeV from
the 0+ continuum states to the 2+ bound state of 12 C is plotted as a function of Kmax .
1

76

|Rnoscreen-Rscreen|/Rnoscreen

0.5

0.45

Kmax=20
Kmax=26

0.4

0.35

0.01

0.1

1

T [GK]
Figure 4.17: Relative differences between screening and no screening calculations for the
fictitious three-alpha scattering system in which only Coulomb interactions are introduced.
The screening calculation is performed for Kmax = 26 and ρscreen = 800 fm. Calculations
without screening are performed for Kmax = 20 (solid) and Kmax = 26 (dashed).
dB(E2)/dE on the left panel (Fig. 4.16a, blue-circle curve) is computed for a real triplealpha system where both nuclear and Coulomb interactions are present and no screening
potential is used. Unlike the Coulomb only case (Fig. 4.16b), the result does not converge
as we increase the size of our model space.
As mentioned before, we tackle this problem by screening the off-diagonal couplings by
a Woods-Saxon multiplying factor [1 + exp((ρ − ρscreen )/ascreen )]−1 . Note that this only
affects the Coulomb part of the couplings. The red-diamond curve in Fig. 4.16a shows a
much more stable behavior of the quadrupole transition strength dB(E2)/dE when screening is introduced. However, it is necessary for us to ensure that important physics is not
left out. We therefore compare in Fig. 4.17 screening and no screening calculations for a
scattering system of three alpha particles when only Coulomb interactions are included. For
this case, there are no resonances in either the two-body or three-body system, thus the

77

propagation technique is stable and we are able to determine the impact of introducing
screening potentials. We perform a calculation of the triple-alpha rate using the screening
technique with ρscreen = 800 fm. This calculation fully converges at Kmax = 26. This result
is then compared to calculations with no screening, presented as the solid and dashed curves
for Kmax = 20 and Kmax = 26, respectively. Uncertainty of the reaction rate due to the
screening technique ranges from 35% to 39% for temperatures within 0.01 − 1 GK. These
are very small numbers given the many orders of magnitude involved in the problem.

4.5.2

Convergence study

It was shown in Fig. 4.16a that when the screening technique is not used, the quadrupole
transition strength dB(E2)/dE at 0.01 MeV diverges as the size of the model space increases.
However if one screens the off-diagonal couplings at large radii, the result becomes more stable and a converged triple-alpha rate can be obtained (see red-diamond line in Fig. 4.16a).
The convergence study of our triple-alpha rate is performed carefully.
Although using screening improves the numerical stability and the convergence of the result, unfortunately it also introduces an approximation into our calculation. It is important
to verify that the triple-alpha rate converges with screening radius and we capture the important physics of our problem. In Fig. 4.18 we present three calculations of the quadrupole
transition strength dB(E2)/dE at low energies by varying the diffuseness ascreen of the
Woods-saxon screening potential from 10 fm to 30 fm while fixing the screening radius at
800 fm. The results are not sensitive to this parameter, we thus fix the diffuseness of 10 fm
throughout this work.
In order to study the convergence of the triple-alpha rate with respect to the screening
radius ρscreen , we perform calculations of the triple-alpha rate using three different screening
78

-20

10

-22

ascreen=10
ascreen=20
ascreen=30

2

4

-1

dB(E2)/dE [e fm MeV ]

10

10

-24

10

-26

0.06

0.07

0.08

0.09

0.1

E [MeV]
Figure 4.18: The quadrupole transition strength are computed at low energies using three
different values of the diffuseness in screening potentials: ascreen = 10 fm (solid), ascreen = 20
fm (dashed), ascreen = 30 fm (dot-dashed).

|Rααα-R(ρscreen)|/Rααα

1
0.8
0.6

ρscreen=600
ρscreen=800
ρscreen=1000

0.4
0.2
0
0.01

0.1

1

T [GK]
Figure 4.19: Relative differences between the triple-alpha rate calculated at several screening
radii and its converged value. The reaction rates calculated using ρscreen = 600 fm (solid),
ρscreen = 800 fm (dashed), ρscreen = 1000 fm (dot-dashed) are compared with the converged
rate in Table 4.4.

79

radii: 600, 800 and 1000 fm. These results are plotted in comparison with the converged rate
in Table 4.4. As seen from Fig. 4.19, the triple-alpha reaction rate in general converges for
ρscreen ≥ 800 fm with a small uncertainty of 5% for temperature within 0.02 − 1 GK. Only
at T = 0.01 GK, the uncertainty gets much larger (∼ 50%). At those low temperatures,
a larger screening radius is needed because the capture happens at much larger distances.
We will see in chapter 6 that the triple-alpha reaction rate at T = 0.01 GK does not have
significant impact on many astrophysical phenomena.
Fig. 4.20 shows the convergence of the triple-alpha reaction rate with hyper-momentum
Kmax . We perform four calculations of the rate as a function of temperature by varying
Kmax from 20 to 26. The results are plotted in comparison with the converged rate in Table 4.4. Kmax = 20 is not large enough to reproduce the experimental 0+ Hoyle resonant
2
state, therefore the reaction rate for this case is larger than other curves for T > 0.1 GK.
The uncertainty of the triple-alpha rate is larger at high temperatures due to the sensitivity
of the Hoyle resonant energy with Kmax . The uncertainty is less than 20% for Kmax = 26.
In Fig. 4.21 we look at the quadrupole strength function, which is a principal ingredient
for calculating the reaction rate, at a specific energy of E = 0.01 MeV. If dB(E2)/dE converges, then the cross section and the reaction rate will also converge. Since the dB(E2)/dE
has an exponential behavior, namely y = A − Be−Cx as a function of hyper-momentum
Kmax , we are able to extrapolate the converged value of dB(E2)/dE when Kmax → ∞. For
this very low energy, the extrapolated value of the rate differs from that at Kmax = 26 by
only 4%. Our calculation is therefore very well converged with the size of the model space.
Although the screening radius ρscreen and the size of the model space Kmax are the major
sources of uncertainties in our calculation due to the difficulty of the long-range behavior in a
three charged-particle system, we have fully tested the convergence of our results with other
80

|Rααα-R(Kmax)|/Rααα

0.01
1.4

0.1
Kmax=20

Kmax=20
Kmax=22
Kmax=24
K+max=26

1.2
1
0.8

1
150
100
50
0

0.6
0.4
0.2
0
0.01

0.1

1

T [GK]
Figure 4.20: Relative differences between the triple-alpha rate calculated at several maximum
hyper-momentum Kmax and its converged value. The reaction rates calculated using Kmax =
20 (solid), Kmax = 22 (dashed), Kmax = 24 (dot-dashed), Kmax = 26 (dot-dot-dashed) are
compared with the converged rate in Table 4.4.

dB(E2)/dE [arb. unit]

-C*x

y=A-Be
HHR

0.3

0.25

0.2
10

12

14

16

18

20

22

24

26

28

Kmax
Figure 4.21: Convergence of the quadrupole transition strength dB(E2)/dE with hypermomentum Kmax for E = 0.01 MeV.

81

-71

-1

dB(E2)/dE [e fm MeV ]

3.2×10

-71

2

4

3.1×10
3.0×10
2.9×10
2.8×10

-71

-71

-71

60

80

100 120 140 160 180 200

njac
Figure 4.22: Convergence of the quadrupole transition strength dB(E2)/dE with njac, number of Gauss-Jacobi quadrature grid point, for a given low energy of 0.01 MeV
variables in the numerical method for completeness. In Fig. 4.22, we study the convergence
of the quadrupole transition strength dB(E2)/dE with njac, the size of the Gauss-Jacobi
quadrature grid for the hyper-angular integrations, for a given low energy of 0.01 MeV. The
result converges very well for njac ≥ 100. As mention earlier we use njac=100 in our calculation of the triple-alpha rate.
It is important that the matching radius ρa is large enough so that we can match the
interior wavefunction with the diagonal Coulomb wavefunction in the asymptotic region.
We perform a calculation of the quadrupole strength dB(E2)/dE at very low energy (0.01
MeV) for which case the matching radius is expected to be very large. Different values of
the matching radius are used to compute dB(E2)/dE. The results are plotted in Fig. 4.23
and it is clear that our calculation converges with the matching radius for ρa ≥ 2000 fm.
We use a matching radius of 3000 fm in our calculation of the triple-alpha rate.

82

-71

-1

dB(E2)/dE [e fm MeV ]

3.06×10

-71

2

4

3.05×10
3.04×10
3.03×10

-71

-71

-71

3.02×10 1500

2000

2500

3000

ρa [fm]

3500

4000

Figure 4.23: Convergence of the quadrupole transition strength dB(E2)/dE with matching
radius ρa for a given low energy of 0.01 MeV

4.5.3

Interaction uncertainty

Since the alpha-alpha interaction is phenomenological based, there are ambiguities that need
to be considered. In order to study the sensitivity of the triple-alpha rate to the nuclear
interaction and the three-body force, we consider employing the same Hamiltonian as [15]
to calculate the 0+ continuum. In [15] the alpha-alpha interaction is given as:

2
2
2
2
Vαα (r) = 100.0 e−r /1.00 − 30.35 e−r /2.13 .

(4.4)

Unlike in [15], we find that a three-body interaction is needed to reproduce the relevant Hoyle
state. As seen from Fig. 4.24, the new interaction produces a triple-alpha reaction rate at
low temperatures 4 orders of magnitude higher than that obtained with the Ali-Bodmer
interaction [64, 65]. While the Ali-Bodmer potential reproduces the α-α phase shifts, the

83

interaction in [15] does not. It therefore provides an upper limit for the error associated with
Vαα ambiguities.
Another issue concerning the uncertainty of the alpha-alpha interactions is as discussed
in [75]. It was stated in [75] that additional narrow resonances are produced when using
the shallow potentials as compared to the deep alpha-alpha interaction. This may lead to
unexpected contributions to the triple-alpha reaction rate. As seen from Fig. 4.8, no other
peaks in addition to the Hoyle state is found in the quadrupole transition strength function
plot in the low energy regime. Therefore the enhancement in our low-temperature triplealpha rate purely comes from the three-alpha non-resonant contribution. The difference
in the spectroscopic properties between the deep and shallow alpha-alpha interactions may
have larger impact in the evaluation of the triple-alpha rate at high temperatures where high
lying resonances play an important role. However, we expect this effect to be too small to be
detected in our temperature region of interest, given the many orders of magnitude involved
in the problem.

4.5.4

Rate uncertainty

Our final rate is normalized to the standard rate from NACRE for temperatures T ≥ 0.5
GK, where the Hoyle resonance completely dominates. At this temperature, the triple-alpha
reaction proceeds through the very narrow Hoyle resonance and its rate is proportional to
the gamma decay width Γγ [9], known experimentally Γγ = 3.7 × 10−9 MeV [69]. In our
calculation we fit the photo-disintegration cross section σγ to a Breit-Wigner shape:

σγ =

Γααα Γγ
1 π 2 c2
,
2 (E − E )2 + Γ2 /4
2 Eγ
R

84

(4.5)

-10

6 -1

-2

〈R ααα〉 [cm s mol ]

10
10
10
10
10

-20
-30
-40
-50

Ali-Bodmer V3b=-19.46 r3b=6
Ogata V3b=-13.426 r3b=7

-60

10 0.01

0.10

1.00

T [GK]
Figure 4.24: The sensitivity of the triple-alpha rate to the interactions: comparison between
α−α Ali-Bodmer+3-body force (solid) and α−α interaction as in [15]+3-body force (dashed)
where Γααα and Γγ are the particle and gamma decay widths respectively; Γ = Γααα + Γγ
is the total width.

As shown in Fig. 4.25, we fit our calculated cross sections to a Breit-Wigner shape
and obtain a value of ∼ 10−9 MeV for the partial decay width, which is the same order
of magnitude as experiment. In practice, a factor of 2 is needed for normalizing to the
NACRE rate, counting for both the uncertainty of the decay width and the convergence of
the problem. In Table 4.3, we summarize different sources of uncertainty in calculating
the triple-alpha reaction rate. At low temperatures the uncertainty mostly comes from the
screening technique. Errors coming from Kmax increase for higher temperatures when the
contribution from the Hoyle resonance becomes important because such a narrow resonance is
very sensitive to the maximum hyper-momentum. Table 4.3 shows that the total uncertainty

85

Table 4.3: Sources of uncertainty in calculating the triple-alpha reaction rate

T (GK)

Uncertainty
Screening ρscreen Kmax

Total

0.01

38%

50%

16%

65%

0.02

36%

4%

2%

36%

0.03

32%

1%

2%

32%

0.04

31%

2%

2%

31%

0.05

37%

2%

2%

37%

0.06

39%

2%

3%

39%

0.07

39%

0.1%

25%

46%

0.08

39%

0.2%

24%

46%

0.09

38%

0.2%

23%

44%

0.1

38%

0.2%

22%

44%

> 0.1

36%

0.2%

16%

39%

of the triple-alpha rate is much smaller than the normalization factor, we therefore consider
this factor of 2 to be our conservative error of this problem.

86

-2

HHR
BW fit

10

-4

2

Eγ σγ / πh c

2 2

10

10

-6

0.379295

0.3793

0.379305

0.37931

0.379315

E [MeV]
Figure 4.25: The photo-disintegration cross section is fitted to the Breit-Wigner shape to
obtained the gamma decay width.

4.6

Rate

With the 2+ bound state and the 0+ continuum states of 12 C obtained in Sec. 4.2 and
1
Sec. 4.3, we are able to compute the triple-alpha reaction rate at different temperatures.
Since both the resonant and the non-resonant continuum 0+ states come naturally from
solving the coupled channels equation (Eq. (2.30)), the HHR method enables us to treat the
resonant and non-resonant processes on the same footing. We present our new triple-alpha
rate in Table 4.4 for temperatures ranging from 0.01 GK to 1.0 GK beyond which we need
to consider other high lying resonances.
The triple-alpha reaction proceeds primarily through either the resonant or the nonresonant path depending on the temperature of the stellar environment. In order to have a
full understanding of the triple-alpha reaction mechanism, we estimate the energy range relevant for the reaction rate at a given temperature. For convenience we denote the integrand
87

Table 4.4: Our triple-alpha reaction rate (in cm6 s−1 mol−2 ) after being normalized to
NACRE for T ≥ 0.5 GK (the factor of 2 uncertainty is listed for temperature below 0.08
GK above which the triple-alpha rate agrees with NACRE).

T (GK)
Rααα
0.010
8.47x10−53

Uncertainty
x2

T (GK)
0.15

Rααα
1.53x10−18

0.015

2.11x10−47

x2

0.20

9.90x10−16

0.020

2.86x10−44

x2

0.25

4.13x10−14

0.025

4.44x10−42

x2

0.30

4.50x10−13

0.030

2.08x10−40

x2

0.35

2.31x10−12

0.04

5.72x10−38

x2

0.4

7.44x10−12

0.05

3.11x10−36

x2

0.5

3.44x10−11

0.06

6.79x10−35

x2

0.6

8.63x10−11

0.07

4.18x10−32

x2

0.7

1.55x10−10

0.08

7.12x10−29

0.8

2.28x10−10

0.09

2.26x10−26

0.9

2.95x10−10

0.10

2.19x10−24

1.0

3.51x10−10

of Eq. (3.25) as function f (E):
E

−
f (E) = (E − Q)2 σγ (E) e kB T .

(4.6)

At a given temperature, there is a certain range of energy that dominates the integrand
f (E). In Fig. 4.26 we plot in dashed-blue the integrand f as a function of the threealpha relative energy E at low temperature (0.01 GK). The energy cumulative sum of this
88

[arb. unit]

10
10
10
10
10

-100

f(E)
Sm(E)

-150

-200

-250

-300

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

E [MeV]
Figure 4.26: The integrand f (E) and its energy cumulative sum Sm(E) at a very low
temperature of 0.01 GK

[arb. unit]

10
10

-60

10
10

-45

f(E)
Sm(E)

-75

-90

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

E [MeV]
Figure 4.27: The integrand f (E) and its energy cumulative sum Sm(E) at T = 0.1 GK
E
function (Sm(E) = 0.01 f (E )dE ) is presented in red. This gives us an estimation for the

energy range that is relevant to calculate the triple-alpha rate at a given temperature. From
89

Fig. 4.26, one concludes that the contribution from the Hoyle resonance to the triple-alpha
rate is negligible at T = 0.01 GK.
A similar calculation is presented in Fig. 4.27 but for a higher temperature T = 0.1

Figure 4.28: Range of energy relevant for the triple-alpha reaction rate at a given temperature.

GK. The energy cumulative sum of function f (E) in this case indicates that the dominant
contribution to the triple-alpha rate comes indeed from the Hoyle state.
We summarize our results in Fig. 4.28. The shaded area represents the energy range
relevant for the evaluation of the triple-alpha reaction rate at a given temperature. As we
can see the resonant energy at 0.38 MeV completely dominates the integration for T ≥ 0.07
GK. This is valid for temperatures below 1 GK because when T > 1 GK additional high lying
resonances contribute significant. For T = 0.06 − 0.07 GK, there is a competition between
the resonant and the non-resonant processes, marking the transition region between those.
The non-resonant capture mechanism dominates for T < 0.06 GK.

90

Chapter 5
Discussion
The HHR method which combines the R-matrix expansion and the R-matrix propagation
in the hyperspherical harmonics representation is developed to solve an exact three-body
problem. We have succeeded in providing the triple-alpha reaction rate for T < 1 GK,
including temperatures below 0.1 GK. In this section we compare our result with other
studies ([10], [15], [28]) and discuss different aspects of the reaction mechanism.

5.1

Comparison with other methods

Fig. 5.1 presents our results (solid line) in comparison with other studies: NACRE, CDCC
and BW(3B). Although the new rate (HHR) is slightly reduced below 0.07 GK, it is significantly enhanced for T < 0.06 GK as compared to NACRE. We also obtained a very different
temperature dependence in this region. In NACRE, the triple-alpha reaction is assumed as
a two-step process for all temperatures. An enhancement of the rate at low temperatures
in our calculation is attributed to the non-resonant three-body direct capture which is not
considered in NACRE. The results assuming an extrapolation of a three-body Breit-Wigner

91

6 -1

-2

〈R ααα〉 [cm s mol ]

10
10
10
10
10

-10
-20
-30

HHR
NACRE
CDCC
BW(3B)

-40
-50
-60

10 0.01

0.10

1.00

T [GK]
Figure 5.1: Different evaluations of the triple-alpha reaction rate: comparing the Hyperspherical Harmonic R-matrix method (solid) with NACRE (dotted), CDCC (dashed) and
the three-body Breit Wigner (dot-dashed).
cross section to low energies BW(3B) (dot-dashed line) [28] have a similar behavior but
the reaction rate increases to a lesser extent. Although the non-resonant treatment in [28]
is crude by assuming the low temperature triple-alpha reaction is driven by the tail of a
three-alpha Breit-Wigner resonance, without the intermediate 8 Be resonant state, that work
shares the same finding as ours, namely that an enhancement of the rate is expected in the
low temperature regime due to the direct capture mechanism. The CDCC calculation in [15]
(long-dashed) produces a large enhancement of the triple-alpha rate for temperatures up to
0.2 GK when comparing with NACRE. That is the first attempt to include the non-resonant
contribution in the triple-alpha reaction rate at low temperatures. However, the effect in
[15] is much stronger than what is seen in our studies. In [15], the three-body wavefunction
is expanded in terms of the continuum states in the two-body subsystem (8 Be), approximately reducing a three-body problem to a 2 + 1 body problem. In practise, truncation in
92

maximum excitation energy and angular momentum of the two-body subsystem is necessary
for numerical performance. Our hypothesis is that this leads to the incorrect asymptotic
form for the three-body scattering wavefunction when all particles have charge [20]. The
CDCC calculation in [15] are combined with a 2+ bound state from a microscopic cluster
1
model (not the same framework). However, we do not expect this would introduce the many
orders of magnitude difference that we observe. Our HHR method approaches the triplealpha problem in a different way, expanding the total three-body wavefunction in terms of
hyperspherical harmonics, equally considering all the possible configurations of the three alpha particles, and highly improving the asymptotic treatment of the three charged-particle
system as compared to the CDCC method. We know that the results in [15] overestimate
the enhancement of the rate leading to a failure of the description of the red giant phase in
the evolution of low mass stars [21].
There exists a kink in the HHR curve around T ≈ 0.06 GK as seen in Fig. 5.1. This marks
the transition between the resonant and the non-resonant processes for which the temperature dependences are different. Above T ≈ 0.06 GK the resonant process dominates, while
below there is mostly direct non-resonant capture. This agrees with the finding in Fig. 4.28.
Although we consider the NACRE rate as our reference, there exist other rates that
are commonly used in astrophysics such as CF88 from [8] and Fynbo from [13]. The theory behind these rates are the same, assuming the triple-alpha reaction proceeds through
a sequential two-step process for all energies. However, the treatment of the high energy
resonances of 12 C causes deviation between them. In Fig. 5.2 we compare the HHR rate
with others namely NACRE, Fynbo and CF88. It is shown that the difference in the rates
typically used in astrophysics is insignificant to the large enhancement that we obtain in the
HHR rate for T < 0.06 GK. The result is published in [14].
93

-2

10

〈R ααα〉 [cm s mol

10

6 -1

10

10
10
10

0

-10
-20
-30

HHR
NACRE
Fynbo
CF88

-40
-50
-60

10 0.01

0.10

1.00

T [GK]
Figure 5.2: The triple-alpha reaction rate: comparing the Hyperspherical Harmonic R-matrix
method (solid) with NACRE (dotted), Fynbo (dashed), and CF88 (dot-dashed).

5.2

Long-range Coulomb effects

Since we obtain a large enhancement for the triple-alpha reaction rate at low temperatures,
it is important to isolate the source of this effect. We therefore perform calculations of the
12 C(0+ )

C
continuum states for cases in which only the diagonal Coulomb couplings Vγγ or the

full Coulomb couplings V C are included. The results are then compared to calculations inγγ

cluding both nuclear and Coulomb interactions (V C +V N )γγ with the off-diagonal Coulomb
couplings up to 30 fm and 800 fm. The reaction rate for each case is then constructed by
using those scattering wavefunctions and fixing the 12 C(2+ ) bound state. Fig. 5.3 shows the
1
results of these calculations. When only the diagonal Coulomb couplings are present (dotted
line), we are able to obtain an analytic solution of Eq. (2.32) (Sec. 2.1.3). The inclusion
of the off-diagonal Coulomb couplings (dot-dashed line) significantly increases the reaction
rate at low temperatures. When both nuclear and Coulomb couplings are fully included
94

-2
6 -1

〈R ααα〉 [cm s mol ]

10
10

-30

10
10

-15

-45

Coul Diag
Coul Only
Nu+Coul (30 fm)
Nu+Coul (800 fm)

-60

0.01

0.10

1.00

T [GK]
Figure 5.3: The long-range Coulomb effects are shown in four different calculations: only
diagonal Coulomb couplings (dotted), only Coulomb couplings (dot-dashed), both nuclear
and Coulomb interactions with off-diagonal Coulomb couplings up to 30 fm (dashed) and a
fully converged calculation with off-diagonal Coulomb couplings up to 800 fm (solid)
in our calculation (dashed and solid lines), we observe an increase in the reaction rate at
high temperatures due to the resonant contribution. Comparison between the solid and the
dashed curves confirms that long-range off-diagonal Coulomb couplings drive the increase
in the triple-alpha reaction rate at low temperatures. Fig. 5.3 indicates that the effect of
off-diagonal long-range couplings are relatively small at high temperatures but very important in the low temperature regime. About 10 orders of magnitude enhancement in the rate
is found at T = 0.01 GK due to these effects, demonstrating the importance of including
Coulomb correctly.

95

5.3

Reaction dynamics

In order to understand the mechanism for 12 C production in more detail, we rewrite Eq. (3.37)
as
dB(E2)
∼
dE

2

fγ i (r, R) dr dR

,

(5.1)

γi

where r is the radius between two alpha particles and R is the distance from their center
i i
of mass to the third alpha particle. γ i = {K i , lx , ly } indicates an incoming channel in the

hyperspherical wave expansion for a scattering state. The spatial distribution of function
fγ i (r, R) at different three-body kinetic energies E contains information about the dynamics of the triple-alpha reaction. Fig. 5.4, 5.5, and 5.6 illustrate the spatial distributions of
fγ i (r, R) at a very low energy E = 0.05 MeV, at the resonant energy E = 0.38 MeV, and
at an energy well above the resonance E = 0.5 MeV, respectively. We just present here
i
the distribution functions corresponding to the first incoming channel (K i = 0, lx = 0, and
i
ly = 0) which is the dominant contribution to the quadrupole strength function dB(E2)/dE.

Other channels exhibit the same trends.
At low energy E = 0.05 MeV, the spatial distribution of function fγ i (r, R) shown in
Fig. 5.4 has a different symmetry in comparison with the higher energies (Fig. 5.5, 5.6). We
observe comparable contributions coming from two different triple-alpha configurations: the
prolate triangle and the oblate triangle as shown in Fig. 4.12a and Fig. 4.12b, respectively.
There is a large cancellation between these two contributions resulting in a small value of
dB(E2)/dE at low energies (for example, the two contributions are ∼ 10−12 but their sum
is ∼ 10−16 ). Nevertheless these cancellations are well within the numerical accuracy of our
computations.
Fig. 5.5 and Fig. 5.6 illustrate the spatial distribution of function fγ i (r, R) at the res96

4 5

- 1 4

-1 0

3 5

-1 0

3 0

R [fm ]

- 1 3

-1 0

4 0

1 0

2 5

1 0

- 1 5
- 1 5
- 1 4
- 1 3

1 0

2 0
1 5
1 0
5
0
0

5

1 0

1 5

2 0

2 5

3 0

3 5

4 0

4 5

r [fm ]
Figure 5.4: Spatial distribution of function fγ i (r, R) for E = 0.05 MeV and γ i = {0, 0, 0}.

- 1 0 0
- 5 0
- 1 0
1 0
5 0
1 0 0

1 2
1 0

R [fm ]

8
6
4

C
2
0
0

B
2

4

6

8

1 0

1 2

r [fm ]
Figure 5.5: Spatial distribution of function fγ i (r, R) for E = 0.38 MeV and γ i = {0, 0, 0}.

97

- 4

12

-10
- 5

-5x10

10

- 5

-10
- 5

R [fm]

8

10
-5

5x10

6

- 4

10

4

C

2
0

B
0

2

4

6

8

10

12

r [fm]

Figure 5.6: Spatial distribution of function fγ i (r, R) for E = 0.5 MeV γ i = {0, 0, 0}.
onant energy E = 0.38 MeV and higher E = 0.5 MeV. These two cases share the same
symmetry. The dominant contribution to the quadrupole strength function dB(E2)/dE
comes from the region within the smallest contour in Fig. 5.5 and Fig. 5.6 that contains
both maxima B and C. The maximum C in the spatial distribution of function fγ i (r, R) is
caused by the triple-alpha equilateral triangle configuration in both the 2+ bound state and
1
the Hoyle resonant state (see the density distribution for each state in Fig. 4.6 and Fig. 4.11).
The contribution to the maximum B mostly comes from the three-alpha oblate configuration which only appears in the Hoyle resonant state. Even though the three-alpha prolate
configuration dominates the Hoyle state’s structure, it does not contribute significantly to
the strength function because the 2+ bound state wavefunction is zero in that region.
1

98

5.4

Astrophysical S-factor

The astrophysical S-factor is often introduced to facilitate the calculation of the reaction
rate at low temperatures in a two-body reaction in which the reaction rate is defined as:

σv =

1/2
8
1
πmab
(kT )3/2

∞

E

−
E σ(E) e kT dE .

(5.2)

0

The cross section σ(E) is expressed in terms of the astrophysical S-factor S(E):

σ(E) =

1 −2πη
e
S(E) ,
E

(5.3)

where η is the two-body Sommerfeld parameter measuring the strength of the Coulomb
Z Z e2
barrier η = a b

mab 1/2
.
2E

Eq. (5.3) removes explicitly the penetration factor, resulting

on an astrophysical S-factor S(E) with weak energy dependence. This improves the quality
of the extrapolation of the cross section to low energies where measurements are not possible
or a simplification in calculating the reaction rate [9, 46]. If this property holds for a threebody S-factor, we will be able to estimate the low-temperature triple-alpha rate in a simpler
way. As derived in Sec. 3.2, the three-body reaction rate has the formula:
3

2
R123 = p! NA 2
c

8π
m3 (µ

1 µ2 )

g0
1
3 g g g (kT )3
2 1 2 3

∞
0

E

−
2
e kT Eγ σγ (Eγ ) dE .

(5.4)

We apply the same analogy in the two-body problem to obtain the relationship between the
three-body astrophysical S-factor S3b (E) and the photo-disintegration cross section σγ :

σγ (Eγ ) =

1 −2πζ
e
S3b (E) ,
2
Eγ

99

(5.5)

10

-12

HHR
fit
-15

2

Eγ σγ

10
10

-18

10

-21

0.2

0.21

0.22

0.23

0.24

0.25

E [MeV]
2
Figure 5.7: Fitting Eγ σγ to an exponential behavior A e
nance.

−2πB
√
E

for energies below the reso-

where ζ is an equivalent to the Sommerfeld parameter in the two-body case.
For a three-body system, the picture of Coulomb interaction is more complex. As seen
from Eq. (2.25) in Sec. 2.1.2, the strength of the Coulomb barrier in a three-body system is
not a constant as in the two-body case but rather depends on the incoming and outgoing
channels in the HH representation. We thus do not have an exact expression for ζ except
√
the fact that ζ should be proportional to 1/ E because the three-body Coulomb couplings
2
also vary as 1/ρ in our method. We estimate ζ by fitting Eγ σγ to an exponential behavior

Ae

−2πB
√
E

for energies below the Hoyle resonance. A and B are fitting parameters. As seen

from Fig. 5.7, the fitting procedure is performed for the energy range (0.2-0.25) MeV, below
which we are not able to obtain any reasonable fit. The result from this fit is then used to
construct the three-body S-factor as in Eq. (5.5). In Fig. 5.8, we plot the energy dependence
of both the cross section σγ and the three-body astrophysical S-factor S3b . A significant
100

E [MeV]

S3b(E)

0.1

0.2

0.3

0.4

0.5

0.2

0.3

0.4

0.5

20

10

10

15

10

2

σγ [fm ]

10

10
10
10

0

-10
-20
-30

10 0.1

E [MeV]
Figure 5.8: Comparison between the energy dependence of the cross section σγ and the
three-body astrophysical S-factor S3b .
reduction in the energy dependence is observed for the S-factor in the low energy range
(0.15-0.33) MeV. For energies below 0.15 MeV, the astrophysical S-factor strongly varies
with energy. This is not fixable by fitting ζ to the lower energy range since no reasonable
fit can be obtained. In conclusion, we are not able to construct a three-body astrophysical
S-factor for the triple-alpha system at very low energies that has weak energy dependence.
For the intermediate energies (0.15-0.33) MeV, a desired three-body S-factor in the same
two-body analogy is possible.
As stated above, the parameter ζ varies with the HH channels causing difficulty in extracting the three-body Sfactor. In order to verify that statement, we perform similar calculations
but only include the first incoming and outgoing channels (0,0,0). This allows us to uniquely
define ζ:
ζ=

eff
Z00

101

m
,
2E

(5.6)

E [MeV]
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

S0 (E)

10
10

2

6

10
10

σγ [fm ]

9

10
10

3

0

-25

-50

10

-75

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

E [MeV]
Figure 5.9: Comparison between the energy dependence of the cross section σγ and the
astrophysical S-factor S0 for Kmax = 0.
eff
where Z00 is calculated from Eq. (2.25). We therefore expect to obtain an S-factor with

much weaker energy dependence for this case. The results are presented in Fig. 5.9. As
you can see, the extracted S-factor is almost constant with relative energy, reinforcing our
conclusion on the difficulty in obtaining a three-body astrophysical S-factor. We also explore
different fitting functions which are commonly used in an effort to extract a desired threebody S-factor. For example, we slightly modify the energy dependence of ζ:

σγ (Eγ ) =

1 − 2πB
e E C S3b (E) ,
2
Eγ

(5.7)

or introduce a polynomial into the S-factor formula:

σγ (Eγ ) =

√
1 − 2πB
E (A0 + A1 ∗ E + A2 ∗ E 2 ) .
e
2
Eγ

102

(5.8)

However, we continue to find difficulty in obtaining a three-body S-factor with weak energy
dependence using these simple expressions.

103

Chapter 6
Astrophysics Applications
As mentioned in chapter 1, the triple-alpha reaction is one of the key reactions in astrophysics. It is responsible for the creation of 12 C and provides a mechanism to overcome the
difficulty in nucleosynthesis of heavy elements due to the instability of nuclei with mass 5
and 8. Helium burning is a critical source of energy for horizontal branch stars, a significant
phase in the life of a low-mass star. Its rate determines the tip of the red giant branch in
a normal star’s evolution. Competition between the triple-alpha rate and the 12 C(α,γ) rate
determines the C/O ratio in the core of a star with mass ranging from 0.4M

to 2M

[9].

Helium burning can also be found in helium accreting white dwarfs and helium accreting neutron stars in compact binary systems. If the accretion rate is sufficiently low, the evolution
of these binary systems could lead to very rare energetic events such as type .Ia supernovae
and superbursts of which many related studies are based on computational simulations. It
is critical that the triple-alpha rate and other inputs in these models are well constrained.
We obtain a large enhancement of the triple-alpha reaction rate at T < 0.06 GK as
compared to NACRE. The new rate also has a very different temperature dependence. In Table. 6.1, we present the temperature dependence of different rates by evaluating d ln Rααα /d ln T .
104

Table 6.1: Temperature sensitivity of triple-alpha rate

T (GK) d ln Rααα /d ln T

0.01
0.02
0.04
0.08
0.16

HHR
34.1
23.3
18.5
51.7
24.4

NACRE
56.5
45.5
47.7
48.3
24.4

dLog〈Rααα〉/dLogT

70
60

HHR
NACRE

50
40
30
20
10
0
7.00

7.50

8.00

8.50

9.00

Log T
Figure 6.1: Temperature dependence for the HHR rate (solid) and the NACRE rate (dashed).
Fig. 6.1 plots the temperature dependence of the HHR rate (solid line) and the NACRE rate
(dashed line). As we can see, the HHR rate has a much weaker temperature dependence than
the NACRE rate for temperatures below 0.06 GK. Beyond this point, the two rates have the
same temperature sensitivity. Since our rate significantly differs from the standard NACRE
rate at low temperatures, we want to investigate its influence on important astrophysical
events.

105

6.1

Evolution of single stars

In this section, we describe a general picture of the evolution of a normal star with solar
metallicity Z

and initial mass ranging from 0.4M

to 2M

(M

is the solar mass). For

heavier stars, we do not expect our rate would have significant impact because helium burns
at higher temperatures [9]. In Fig. 6.2 we present the standard evolution curve of a 1M
star. Stars are first formed from an interstellar gas cloud which is mainly composed of hydrogen and helium. In the first stage of their lives as premain-sequence stars, they undergo
gravitational contraction. This leads to the increase of density, causing the temperature in
the star to rise. The gas cloud is hot enough that some primordial nuclei start to burn. The
nuclear energy generated at this stage is small but sufficient to halt the gravitational collapse
of the central part of the cloud, forming the core of the star. When the core temperature
reaches ∼ 0.01 GK, hydrogen burning starts. It becomes the only energy supply of the core
when the star approaches its main-sequence stage (label “1” in Fig. 6.2). Depending on the
mass of the star, hydrogen fusion can proceed via the pp chain or the CNO cycle. Stars with
M

1.5M

burn hydrogen at higher temperature through the CNO cycle, whereas lighter

stars do it through the pp chain. The energy produced in the core is carried outward and
increases the surface luminosity.
A significant amount of helium is accumulated in the core due to hydrogen fusion. The
star eventually runs out of hydrogen in the core. It starts burning hydrogen in a shell surrounding the inert helium core. The star begins to leave the main sequence. Because there
is no energy generated in the core, the helium core contracts again due to gravity. This increases the temperature in the core as well as in the hydrogen burning shell. More energy is
produced in the shell causing a rapid increase in the surface luminosity. The star’s envelope

106

Figure 6.2: Standard evolution curve of 1M stars (surface luminosity as a function of
surface temperature): 1 - main sequence, 2 - red giant, 3 - core helium flashes, 4 - thermal
pulses, 5 - planetary nebula, and 6 - white dwarf.
becomes fully convective and expands dramatically. The drastic expansion together with the
star’s surface reddening indicates that the star enters a new phase in its evolution process,
the so-call red giant stage (label “2” in Fig. 6.2). The maximum luminosity is achieved when
the star gets to the tip of the red giant branch.
The gravitational contraction of the core leads to an increase in the central temperature
and density. Around 0.1 GK, the star starts burning helium in the core. Since the central
density is too high, matter becomes electron degenerate and therefore does not respond to
the increasing temperature in the core due to helium fusion. The degenerate pressure prevents the star’s surface expansion while the temperature keeps rising. More and more energy
is generated inside the star. Eventually the rising temperature produces a thermonuclear
runaway, or a core helium flash (label “3” in Fig. 6.2). This lifts the degeneracy and the star
is able to expand and cool off to a stable state. The star settles at the horizontal branch

107

and burns the remaining core helium in a stable manner.
When helium in the core runs out, gravitational contraction starts again and heats up
the stellar interior. The increasing temperature ignites the helium shell next to the carbonoxygen core. The hydrogen shell next to the helium burning region also begins to burn.
The star enters the asymptotic giant branch. The hydrogen shell burning is so active that
more and more helium produced by hydrogen fusion in this shell is added to the helium
zone, causing its density and temperature to increase. However, the inner helium shell is
such a thin layer that energy generated can not be carried outward fast enough by radiative
diffusion. As a result, a thermonuclear runaway occurs. This unstable helium shell burning
is also called helium shell flash or thermal pulse (label “4” in Fig. 6.2).
Eventually, hydrogen and helium in the shell are exhausted, and the star leaves the
asymptotic giant branch. Due to a strong stellar wind, most of the star’s envelope is ejected.
The star enters the planetary nebula stage (label “5” in Fig. 6.2). The surface temperature
of the star keeps increasing because hotter layers are now exposed.
When there is no hydrogen envelope left and the hydrogen shell burning completely dies
out, the star begins to lose its luminosity (label “6” in Fig. 6.2). The star will end its life
as a white dwarf which is mainly composed of carbon and oxygen. The electron degeneracy
pressure prevents the white dwarf from gravitational collapse. It eventually cools by radiating away the thermal energy.
It is clear that the triple-alpha reaction plays an important role in the evolution of a star.
Changes in the rate can lead to a very different evolution of the star. When Ogata et al.
[15] first claimed a drastic enhancement in the reaction rate due to the non-resonant contribution, an investigation of the impact of that finding on the evolution of a single star was
performed immediately thereafter. In [21] Dotter et al. compared the evolutionary tracks
108

obtained using the standard NACRE rate and the rate from [15] (denoted as OKK in [21]
and CDCC in our work). Using the OKK (CDCC) rate severely shortens the red giant phase.
This is because the drastic enhancement of the triple-alpha rate in OKK (CDCC) results in
core helium burning at much lower luminosity and temperature. Since the predictions from
the stellar evolution model using the NACRE rate agree very well with observations of red
giant branch (RGB) stars [76], Dotter et al. [21] conclude that the OKK (CDCC) result is
inconsistent with observations.
In this work we obtain a new triple-alpha rate that is significantly larger than the standard NACRE rate for T < 0.06 GK. The rate enhancement is less than that of Ref. [15].
It is important to explore the impact of our new rate on the stellar evolution. We repeat
the calculation in [21] for a 1M

star with standard composition of 0.02 metallicity (Z )

and 0.7 hydrogen using our HHR rate. A code named mesa (Modules for Experiments in
Stellar Astrophysics) is used for the calculations [77, 78]. The result is then compared to
those using NACRE rate.
Fig. 6.3 presents the full evolutionary track of a 1M

star calculated using our HHR

rate (solid curve) and the NACRE rate (dashed curve). We also include the results using the
Fynbo (dotted curve) and the CF88 (dot-dashed) rates for completeness. The star evolves
from the main sequence through the red giant phase, the asymptotic giant branch, the planetary nebula stage and ends its life as a white dwarf. Fig. 6.4 blows up a small portion of the
track in which the star goes from the main sequence to the end of thermal pulsation in the
asymptotic giant branch. Because the HHR rate is normalized to the NACRE rate at high
temperature, the difference in the evolutionary tracks reflects the impact of our rate in the
low temperature region. It is clear that both curves are identical from the main sequence
to the end of the red giant branch. This is not the case for the reaction rate from [15], for
109

Figure 6.3: Full evolutionary track (luminosity vs. surface temperature) from main sequence
to white dwarf of a one solar-mass (1M ) star for HHR (solid), NACRE (dashed), Fynbo
(dotted), and CF88 (dot-dashed).

Figure 6.4: Evolutionary track (luminosity vs. surface temperature) from main sequence to
the end of asymptotic giant branch of a one solar-mass (1M ) star for HHR (solid), NACRE
(dashed), Fynbo (dotted), and CF88 (dot-dashed).

110

which the red giant phase either severely shortens or completely vanishes. We observe a
small difference beginning at the thermal pulsation phase in the asymptotic giant branch, as
the enhanced HHR rate at T < 0.06 GK tends to produce a slightly larger convective zone
during helium shell flashes. The difference at the onset of hydrogen and helium thermal shell
flashes causes the two tracks to depart from each other during the planetary nebula phase.
However those differences are probably smaller than the inherent uncertainties in using a
one-dimensional code to simulate this highly unstable stage. The two tracks rejoin when the
star cools off and dies as a white dwarf. The final white dwarfs obtained in both calculations
(HHR and NACRE) has the same carbon-oxygen (C/O) ratio of 1 : 2.08.
The Fynbo and CF88 tracks are identical but noticeably differ from the NACRE track
in the asymptotic giant branch. Because these rates are relatively consistent at low temperatures, the above differences must come from the treatment of the high lying resonances in
12 C.

Using the Fynbo or the CF88 rate, the final white dwarf of a 1M

star obtains a C/O

ratio of 1 : 2.37 which is ∼ 12% smaller than that from NACRE. This issue, which does
not arise from the low-temperature behavior, must be addressed because many astrophysical
scenarios burn helium at higher temperatures.
Fig. 6.5 plots the surface luminosity as a function of the stellar age for our HHR rate
(solid curve) and the NACRE rate (dashed curve). The peak of the surface luminosity when
using the HHR rate occurs at about a million years earlier than for the results using the
NACRE rate. This is because the HHR rate is much larger than the NACRE rate for temperatures below 0.06 GK. At this low temperature, the star with HHR rate will burn helium
faster than that with the NACRE rate.
For completeness, we also perform comparisons between HHR and NACRE for other
stars with mass M = 0.08M , 1.25M

and 1.5M . The evolutionary track for a 1.25M
111

Figure 6.5: Surface luminosity as a function of age: HHR rate (solid) and NACRE rate
(dashed).
star exhibits similar behavior as seen in the 1M
tionary track of a 0.8M

case. In Fig. 6.6, we show the full evolu-

star calculated using our HHR rate (solid curve) and the NACRE

rate (dashed curve). For this case, helium burning happens at temperatures above 0.07 GK,
where our rate and NACRE are the same, and we therefore obtain identical evolutionary
tracks. The same result is obtained for a 1.5M star. In conclusion, a large enhancement in
the HHR rate at T < 0.06 GK does not produce severe changes in the evolutionary tracks of
normal stars with solar metallicity, contrary to what happens when the CDCC rate is used.
Our rate preserves the evolution of stars from main sequence to the red giant branch. Only
small differences are obtained during the thermal pulsation phase for stars of mass 1M and
1.25M .

112

Figure 6.6: Full evolutionary track (luminosity vs. surface temperature) from main sequence
to white dwarf of a 0.8 solar-mass (0.8M ) star for the HHR rate (solid) and the NACRE
rate (dashed).

Figure 6.7: Binary system and its Roche lobes

6.2

Binary stars

A binary star system is defined as a pair of gravitationally bound stars. In a close binary
system, the evolution of one star can strongly affect the other. A binary system is characterized by the Roche lobe, which forms a teardrop-shaped region in space (see Fig. 6.7).

113

It is defined as the equipotential surface of the binary system due to the gravitational and
centrifugal forces when considering a rotating coordinate frame of the center of mass. Materials confined in this region are bound to the corresponding stars by gravity. The point of
contact between the two Roche lobes is called the inner Lagrangian point. If one star fills
its Roche lobe, material from that star can flow through the inner Lagrangian point to the
other companion. This gives rise to some interesting astrophysical phenomena such as type
.Ia supernovae in helium accreting white dwarfs and superbursts in helium accreting neutron
stars. In this work we explore the influence of the new triple-alpha rate in these types of
astrophysical events.

6.2.1

Helium accreting white dwarfs

We consider a binary system consisting of a helium donor and a carbon-oxygen white dwarf
accretor which is often called AM Canum Venaticorum (AM CVn) binary [79]. Helium from
the donor is transferred onto the white dwarf through the inner Lagrangian point. It is
accumulated on the white dwarf’s surface. Due to the intensive gravitational contraction,
the helium layer is strongly compressed. As a result, the temperature at the bottom of the
helium layer keeps rising as more and more helium is added to the white dwarf’s surface.
When the temperature is high enough, helium burning occurs. The bottom layer density
is so high that matter becomes electron degenerate and therefore a thermonuclear runaway
ensues. During this process a very bright outburst is produced. Depending on the mass of
the accreted layer when the explosion occurs, these phenomena are classified as helium nova
or type .Ia supernovae.
The ignition of helium at the base of the accreted layer is closely related to the nuclear

114

heating time scale:
τn =

Cp T

,

(6.1)

ααα

where T is the star’s temperature and Cp is the specific heat calculated at constant pressure.
The triple-alpha generation rate ααα is given as:

Q 2
ρ
ααα = NA
6

Y 3
fscr Rααα .
4

(6.2)

In this equation, NA is Avogadro number; Q is the Q-value of the triple-alpha reaction; ρ
is the density of the star; Y is the helium mass fraction; fscr is the screening factor and
Rααα is the triple-alpha reaction rate.
The work in [80] found that the condition for helium ignition to occur is when τn = 106
yrs. Solving Eq. (6.1), one obtains a set of ignition points (ρ, T ) which show the density
and temperature of the star when the helium flash triggers. It is clear from Eq. (6.1) and
Eq. (6.2) that the helium ignition depends on the triple-alpha reaction rate. It is interesting
to see if our new rate produces significant changes of helium ignition in these dynamical
events.
For uncertainty estimation, two values of the nuclear heating time scale τn = 105 yrs and
τn = 107 yrs are used to produce the helium ignition curve. Calculations are performed with
both NACRE and HHR rates. It is important to mention that the screening factor fscr in
Eq. (6.2) is taken from mesa [81]. This screening factor is constructed as if the triple-alpha
reaction is a two-step process, which is not a valid assumption at low temperatures. Given
that the direct triple-alpha screening description is not yet available, using a consistent
screening factor from mesa in both calculations will provide some insight of the impact
that the HHR rate may have on helium accreting white dwarfs. In Fig. 6.8, we plot the
115

Figure 6.8: Helium ignition curves for the NACRE rate (blue solid) and the HHR rate (red
dashed). Calculations are performed for helium accreting white dwarf with initial mass of
1M and several values of accretion rate.
helium ignition curves calculated using the NACRE rate (blue solid curve) and the HHR
rate (red dashed curve). The helium shell flashes are expected to occur in the shaded region
for each case [23]. At high temperatures, the HHR and the NACRE rates are identical
resulting in a very good agreement between the two corresponding ignition curves. These
curves diverge as the ignition temperature decreases. As seen from Fig. 6.8, the differences
are significant between the NACRE and HHR calculations. At the same temperature, the
helium ignition in the HHR calculation occurs at much lower density as compared to the
NACRE calculation. When does this difference become important? In order to answer that
question, we show in Fig. 6.8 the evolutionary tracks at the base of the helium accreted
layer, for three different accretion rates dM/dt (M yr−1 ) = 10−8 , 10−9 and 5 × 10−10 . The
ignition of helium in all cases falls within the predicted ignition curves. At high accretion

116

rate dM/dt = 10−8 M yr−1 , the NACRE and HHR evolutionary tracks are almost identical.
Helium ignition in this case occurs at higher temperature where the NACRE and our rate
are the same. As expected, as the accretion rate decreases, the ignition happens at lower
temperature. The evolutionary tracks are significantly different between NACRE and HHR.
For the lowest accretion rate in which we see the most difference between our rate and
NACRE, the ignition occurs when the mass of the helium accreted layer is about 0.315M
for HHR rate and 0.360M

for NACRE rate. These events fall into the type .Ia supernova

category.
In Fig. 6.9 we study the dependence of helium ignition on the initial mass of the white

Figure 6.9: Helium ignition curves for the NACRE rate (blue solid) and the HHR rate (red
dashed). Calculations are performed for helium accreting white dwarf with the accretion
rate of 5 × 10−10 M yr−1 and two different initial masses: 0.734M and 1.315M .

dwarf accretor. Two different initial masses of 0.734M

and 1.315M

are employed. The

calculations are performed using a low accretion rate of 5 × 10−10 M yr−1 where the HHR
117

rate is likely to differ from the NACRE rate. The evolutionary tracks at the bottom of the
helium accreted layer and the ignition curves are plotted in the same graph (Fig. 6.9) similarly
to Fig. 6.8. As expected, we obtain significant differences in the evolutionary tracks between
the NACRE and the HHR calculations for both masses. The helium ignitions happen at
much lower density when using the HHR rate. It is also seen that the ignition temperature
is higher for a white dwarf with larger initial mass. This is because a stronger gravitational
force is exerted on the accreted materials causing the bottom layer to heat up faster when
the white dwarf is larger.
Type .Ia supernovae were predicted by theoretical calculations [82]. They should last for
about 10 days and have one-tenth the luminosity of a normal type Ia supernova, hence the
name. We find that our triple-alpha rate produces very different results from NACRE for
this type of explosive scenarios. However, type .Ia supernovae are very rare events that are
predicted to occur once every 5000-15000 yr in a galaxy of 1011 M

[82]. There are so far

only two observed candidates for type .Ia supernovae: SN 2002bj [83] and 2010X [84]. The
properties of those have not been fully understood. In addition, the current one-dimensional
code is not designed to simulate such energetic explosion. Therefore comparison between our
calculated type .Ia supernova and observations are not yet available. We hope our results
would motivate astronomers to collect more data of these phenomena.

6.2.2

Helium accreting neutron stars

In a binary system consisting of a neutron star and a helium companion, helium is transfered
from the donor and accretes onto the neutron star’s surface. This accreted layer is strongly
compressed by the gravitational force causing the temperature at its base to rise. At some
point helium ignites, and the thin shell instability [85] enhanced by electron degeneracy
118

Luminosity (1033 erg s−1 )

7
6
5
4
3
2
1
0

CF88
HHR

0

10

20

30

40

Time (years)

50

60

70

Figure 6.10: The light curve (luminosity vs. time) for helium accreting a neutron star at
accretion rate of 1.75 × 10−10 M yr−1 : Comparing the CF88 rate (dashed) and the HHR
rate (solid).
leads to a thermonuclear runaway. The neutron star releases an extensive amount of energy.
During this process an outburst of X-rays is produced. This phenomena is called a type I
X-ray burst. Depending on the duration, bursts are classified as normal, intermediate or
superbursts. Note that many superburst theories favor the ignition of a thick carbon layer
instead [86].
We explore the impact of our new rate on helium accreting neutron stars. We expect
to see the effect of our results most strongly at low accretion rates because helium burning
then occurs at lower temperatures, where the HHR rate is different from NACRE and other
rate predictions. Simulations are performed by L. Keek [87] for a helium accreting neutron
star, at an accretion rate of 1.75 × 10−10 M

yr−1 and a low crustal heating of 0.056

MeV/nucleon. We compare results obtained with our HHR rate and the CF88 rate [8]. The
results are presented in Fig. 6.10. While the CF88 calculation exhibits a sudden increase
in the luminosity once the X-ray burst ignites, we observe a slow rise of the light curve
in the HHR calculation for half of the burst recurrence time. In addition, the recurrence

119

time obtained from the HHR calculation is substantially shorter than that from the CF88
calculation. We attain an ignition mass of 0.704 × 10−8 M
10−8 M

for the HHR rate and 1.156 ×

for the NACRE rate. Based on the thermal timescale of the ignited helium layer,

the duration of this burst would fall into the superburst category. In the HHR calculation,
helium is ignited first due to a large enhancement in the rate at low temperatures, leading
to a smaller amount of helium piled up on the neutron star’s surface. As a result, the burst
is less energetic than that in the CF88 calculation, which explains a slower increase in the
luminosity of the light curve (Fig. 6.10). Our preliminary results demonstrate the large
impact of the new triple-alpha rate on this burst.
Because of its long recurrence time, the number of observed instances for this burst is
very low. In our model, we use approximately the same low accretion rate as the system
with the observed superburst 0614+091 [88]. Because the event is too energetic, we are
not able to simulate the entire burst. Therefore, the peak luminosity of this burst is not
obtained. In addition, the start of the 0614+091 [88] superburst was not observed. For
these two reasons, no comparison between the simulated and the observed superburst can
be made at this stage. We certainly hope more data of this event would be collected in the
future. Also, the upcoming advanced telescopes may be able to detect the slow rise of the
light curve produced by HHR.

120

Chapter 7
Summary and outlook

7.1

Summary

We have successfully calculated the triple-alpha reaction rate in the low temperature regime
T < 0.1 GK, where many numerical attempts have failed before. In this work, the triplealpha reaction is modeled as a three-body problem. We use the hyperspherical harmonics
(HH) method to tackle this problem. In the low temperature region, the triple-alpha reaction
proceeds through a quadrupole transition from the 0+ continuum to the 2+ bound state in
1
12 C.

The 2+ bound state is obtained by solving a set of coupled channels equations in
1

hyper-radius coordinates for negative energy, under the condition that the wavefunction
goes to zero at large distances. The same approach can not be applied to the 0+ continuum
states because they require an exact boundary condition for the three charged particles. We
overcome this difficulty by combining the R-matrix expansion and the R-matrix propagation
method in the hyperspherical harmonics basis. We first solve a set of coupled channels
equations for a positive energy in a small R-matrix box with fixed logarithmic derivatives.
The R-matrix at the boundary of the box is then propagated to a sufficiently large radius,
121

where the Coulomb asymptotic matching can be performed. We also implement a technique
of screening the off-diagonal Coulomb couplings to ensure numerical stability at very low
energies. The screening technique allows us to obtain analytic solutions for the Coulomb
asymptotic wavefunction, which are then used in the matching procedure. We also study
the symmetry properties of our problem. Since an alpha particle can be considered as a
boson of spin zero, the wavefunction must be unchanged for any permutation of the alphaalpha pair. This only allows even partial waves to contribute to the total wavefunction.
The symmetrization also enables us to reduce three Faddeev components to a set of coupled
channel equations in the hyper-radius formed from one pair of Jacobi coordinates. This
strongly increases the computational efficiency.
For the alpha-alpha interaction we use the Ali-Bodmer potential, which reproduces the
low energy phase shifts as well as the 0+ resonance of 8 Be. A three-body force is necessary to
fit experimental data. We obtain good convergence for the binding energy and the root-meansquare radius of the 2+ bound state in 12 C. Our binding energy agrees with measurements
1
and the root-mean-square radius is compatible with a microscopic study, supporting the
correct structure for the bound state wavefunction. In the HHR framework, we are not
only able to calculate the 0+ non-resonant continuum states in 12 C, but also reproduce
its well known Hoyle state. Therefore, this method allows us to include the resonant and
non-resonant process on the same footing. The density distribution function for the Hoyle
state is constructed to study its structure. We find the prolate triangle configuration of the
triple-alpha to be dominant.
When both the 2+ bound state and the 0+ continuum states of 12 C are obtained, we
1
are able to construct the quadrupole transition function and the cross section, from which
we determine the triple-alpha reaction rate. A careful assessment of different sources of
122

uncertainty for this rate, arising from the numerical calculation is also provided. We first
check our method for the pure Coulomb case where no resonances are present in either the
two-body or three-body systems. The propagation method for this case is stable without
screening. An relatively small uncertainty of 35% due to the screening technique is introduced
in the calculation of the reaction rate at low temperatures. A thorough convergence study
is then performed for the realistic triple-alpha problem. Our reaction rate converges well
with screening radius ρscreen and the size of our model space Kmax . We estimate an overall
uncertainty of a factor of 2 for our triple-alpha rate, due to the screening approximation and
the resonant width estimation.
Our triple-alpha reaction rate obtained by the HHR method agrees with NACRE at high
temperatures but increases significantly at temperatures below 0.06 GK, accompanied by a
very different temperature dependence. The enhancement of the HHR rate at low energy is
not as strong as that seen in the CDCC result of [15]. We find a dominant contribution of
the non-resonant continuum to the triple-alpha rate at low temperatures. Our study shows
that this non-resonant contribution is driven by the long-range Coulomb effects. In this work
we include the Coulomb couplings to all orders, up to sufficiently large radius (800 fm). The
transition between the resonant and non-resonant process occurs around 0.06 GK.
Our triple-alpha rate differs strongly from previous predictions. While the CDCC rate
does not reproduce the red giant phase in the evolution of low-mass stars, our HHR rate does.
The evolutionary tracks are very similar between the HHR and the NACRE calculations.
Using our new rate, stars still undergo the red giant phase, contrary to the results using the
CDCC rate. We observe small differences in the thermal pulsation and the planetary nebula
stages. However, we expect these differences to be smaller than the inherent uncertainties
in the astrophysical model. In conclusion, our rate does not drastically change the evolution
123

of low-mass stars, opposed to the CDCC rate. The implication of our new rate is also
investigated for other astrophysical scenarios that burn helium at lower temperatures. We
look at helium accreting white dwarfs and helium accreting neutron stars with small accretion
rates, which can give rise to very energetic events like type .Ia supernovae and superbursts.
Helium burning in these scenarios occurs at lower temperatures where the HHR rate and
the NACRE rate are very different. We find that helium ignition condition are very different
between the HHR and the NACRE rates. At a given temperature, the helium shell flash
ignites at much lower density when the HHR rate is used. As a result, the ignition mass in
the HHR calculation is smaller then that in the NACRE calculation. There is also significant
difference in the light curve of a helium accreting neutron star between the HHR rate and the
NACRE rate. The recurrence time of the burst is much shorter using the HHR rate. Because
the current one-dimensional codes are not designed for such energetic events like type .Ia
supernovae and superbursts, where our rate demonstrates large impact, full simulations of
these explosive scenarios are not possible. We should also note that, until today, the number
of observed instances for such events is very low, therefore comparison with our calculation
is not yet available at this stage. Our results serve as an additional stimulus to collect more
data on these types of phenomena.

7.2

Outlook

The presence of very narrow two-body and three-body resonances, in addition to the strong,
long range Coulomb interaction makes the triple-alpha problem very challenging. Our HHR
framework which is the combination of various methods is a new approach to overcome the
well known difficulty of the three charged-particle system. This method allows us to approach

124

the very low energy regime where measurements are impossible without using extrapolation
and results obtained have an impact in astrophysics. Our success in solving the triple-alpha
problem opens the opportunity for further applications of the method to other three-body
systems.
Within the triple-alpha system, there is still room for further developments. Throughout
this work, we employ a fixed zero logarithmic derivative for the R-matrix method due to its
mathematical simplicity. One might extend the R-matrix propagation method for arbitrary
values of the logarithmic derivative β. This could improve the numerical stability of the
method and free us from using the screening technique. Because the aim of this work is
to calculate the triple-alpha rate at low temperatures, our convergence study is currently
performed with respect to this quantity. One can further explore the convergence of the threebody scattering wavefunction from phase-shift and scattering matrix analysies. Because the
problem involves so many channels, it is not trivial to continuously link the eigenphases
at various energies and associate the eigenphases with certain hyper-spherical harmonics
channels. A specific algorithm would need to be developed. The scattering matrix and
phase-shift analysies can also improve the accuracy of extracting the properties of the Hoyle
resonant state of 12 C. Cluster structures of the 2+ bound state and the 0+ resonant state
1
2
can be further investigated to obtain information about the angular orientation between the
three alpha particles. A recent ab initio lattice calculation [89] finds a “bent-arm” triangle
configuration for the Hoyle state when taking into account the effect of rotational excitation.
It will be interesting to compare that, with results obtained from our three-body framework.
Although specifically designed for the triple-alpha problem, the HHR method can be
applied to other three-body systems. It has a large application for many important reactions
in astrophysics. For example, one can use this method to directly calculate the production
125

rate of elements that are relevant in particular astrophysical environments such as 6 He and
9 Be.

The study of dineutron decay in nuclear reactions (such as 16 Be in [90] ) can be studied

using the same framework. Because these three-body systems do not involve the Coulomb
interaction, the calculations are expected to converge much faster. Also important, it opens
new opportunities in addressing three-body, low energy reactions for particles with charge
where the physics is more challenging. For example, the diproton decay experiment of 6 Be
[91] can be analyzed using this method. The HHR method can be extended to any threebody problem in atomic and molecular physics. The work in [92] employs a three-body
approach to calculate the ground state energy of the two electron atoms and ions such as
negatively charged hydrogen, helium, and positively charged lithium. The results can be
further improved by using our method.

126

APPENDICES

127

Appendix A
Detailed Derivations of the R-matrix
propagation method

A.1

Deriving Eq. (2.58)

In this section, we will present the detailed derivation of the propagating functions gi in
Eq. (2.58). This is the key element in the R-matrix propagation method. First, we insert
Eq. (2.56) into Eq. (2.54) for sector p to obtain:
p

d2 χ i (ρ)
γγ
dρ2

p
αγ

p

˜
Tγα λ2 (p) T
α

=

χ

p
(ρ) .
γ γi

(A.1)

αγ

Because T is a unitary matrix which means

F

p
(ρ)
αγ i

γ

˜p T p
T
αγ

γβ

= δαβ , we can define:

˜p χp (ρ) ,
T
αγ γ γ i

=
γ

p
χ i
γγ

p

=

Tγα F
α

128

p
(ρ) .
αγ i

(A.2)

Inserting Eq. (A.2) into Eq. (A.1), we obtain a simpler equation:
p
(ρ)
αγ i
dρ2

d2 F

p
(ρ) .
αγ i

= λ2 (p) F
α

(A.3)

The general solutions of Eq. (A.1) are well known:

p
F i (ρ)
αγ

=


 p

A sinh(|λ (p)|ρ) + B p cosh(|λ (p)|ρ),
 α
α
α
α

 p
A sin(|λ (p)|ρ) + B p cos(|λ (p)|ρ),
 α
α
α
α

λ2 (p) > 0
α
λ2 (p)
α

.

(A.4)

≤0

From now on we will drop the index p for convenience and keep in mind that these are the
solutions for sector p. The constants Aα and Bα are determined by the boundary conditions
of sector p at ρL and ρR (see Fig. 2.5):










 FR  g1 g2   −FR 

,

=
g3 g4
FL
FL

(A.5)

where gi are the diagonal matrices. Eq. (A.5) for F is equivalent to Eq. (2.57) for χ due to
the relationship in Eq. (A.2). We rewrite Eq. (A.5) in the following form:

(Fαγ i )R = −(g1 )αα (Fαγ i )R + (g2 )αα (Fαγ i )L ,
(Fαγ i )L = −(g3 )αα (Fαγ i )R + (g4 )αα (Fαγ i )L .

(A.6)

For simplicity, we assume ρL = 0 and ρR = ρR − ρL = hp without losing the generality of
the problem.

129

For λ2 > 0 :
α
We first evaluate Eq. (A.4) at the two boundaries of sector p and then plug the results into
Eq. (A.6):

Aα sinh|λα hp | + Bα cosh|λα hp | =
− (g1 )αα Aα |λα | cosh|λα hp | + Bα |λα | sinh|λα hp | + (g2 )αα Aα ,

(A.7)

Bα = −(g3 )αα Aα |λα | cosh|λα hp | + Bα |λα | sinh|λα hp | + (g4 )αα Aα .

(A.8)

In order for Eq. (A.7) and Eq. (A.8) to be satisfied with all values of Aα and Bα , we have
to constrain (gi )αα by the following equations:

(g1 )αα |λα | cosh|λα hp | − (g2 )αα + sinh|λα hp | = 0 ,

(A.9)

(g1 )αα |λα | sinh|λα hp | + cosh|λα hp | = 0 ,

(A.10)

(g3 )αα |λα | cosh|λα hp | − (g4 )αα = 0 ,

(A.11)

(g3 )αα |λα | sinh|λα hp | + 1 = 0 .

(A.12)

Solving these equations we obtain the sector propagating functions gi :
coth|λα hp |
,
|λα |
1
p
p
(g2 )αα = (g3 )αα = −δαα
.
|λα | sinh|λα hp |
p

p

(g1 )αα = (g4 )αα = −δαα

130

(A.13)
(A.14)

For λ2 ≤ 0 :
α
In a very similar manner we obtain the sector propagating functions gi for this case:
cot|λα hp |
,
|λα |
1
p
p
(g2 )αα = (g3 )αα = δαα
.
|λα | sin|λα hp |
p

p

(g1 )αα = (g4 )αα = δαα

(A.15)
(A.16)

The results in Eq. (A.13), Eq. (A.14), Eq. (A.15) and Eq. (A.16) are briefly presented in
Eq. (2.58) in Sec. 2.2.2

A.2

Deriving Eq. (2.61)

In this section, we derive the propagating equation of the R-matrix: Eq. (2.61). In order to
that we first insert Eq. (2.59) and Eq. (2.60) into Eq. (2.57):

p

p

p

p

p

p

ρR Rp χR = −G1 χR + G2 χL ,
p−1

ρR

p

p

p

p

(A.17)
p

Rp−1 χL = −G3 χR + G4 χL ,

(A.18)

Arranging Eq. (A.17) and Eq. (A.18) we obtain:

p

p

p

p

p

(−G1 + ρR Rp ) χR = G2 χL ,
p

p

p

p−1

G3 χR = (G4 − ρR

p

p

Rp−1 ) χL ,

p−1

We first multiply both sides of Eq. (A.20) with (G4 − ρR
p

(A.19)
(A.20)

Rp−1 )−1 to the left and then

with G2 also to the left. Afterward, we compare the new equation with Eq. (A.19) and

131

obtain:
p

p

p

p

p−1

(−G1 + ρR Rp ) = G2 (G4 − ρR

p

Rp−1 )−1 G3 .

(A.21)

Eq. (A.21) reduces to Eq. (2.61) in Sec. 2.2.2, after trivial manipulation.

A.3

Deriving Eq. (2.63), Eq. (2.64), Eq. (2.65), and
Eq. (2.66)

This section presents the detailed derivation of several equations namely Eq. (2.63), Eq. (2.64),
Eq. (2.65), and Eq. (2.66) which result in the propagation of the three-body scattering wavefunctions. We first rewrite Eq. (2.62):

p

p

p

(A.22)

p

p

p

p

(A.23)

χR = −G 1 χR + G 2 χ1 ,
L
χ1 = −G 3 χR + G 4 χ1 ,
L
L

We obtain a similar set of equations for sector p − 1:

p−1

χR

p−1

= −G 1

p−1

χ1 = −G 3
L

p−1

χR
p−1

χR

p−1

+ G2

p−1

+ G4

χ1 ,
L

(A.24)

χ1 ,
L

(A.25)

p−1

Because the wavefunctions are continuous at the sector boundaries, i.e χR
p−1

χR

p

= χL and

p

= χL , we can express the Eq. (2.57) in the following form:

p

p

p

p

p−1

χR = −G1 χR + G2 χR
p−1

χR

p

p

p

p−1

= −G3 χR + G4 χR
132

,

(A.26)
.

(A.27)

p

p−1 −1
]

We multiply to the left of both sides of Eq. (A.22) and Eq. (A.24) with [G 2 ]−1 and [G 2
respectively and then subtract the results:

p

p

p−1 −1
]

[G 2 ]−1 χR − [G 2

p−1

p−1 −1
]

p−1

p−1

G1

χR

p−1 −1 p
p−1
p−1 p−1
] G4 +[G 2 ]−1 G 1 )χR

p

p

= [G 2

p

− [G 2 ]−1 G 1 χR .

−([G 2

χR

(A.28)

Inserting Eq. (A.27) into Eq. (A.28) we obtain:

p

p

[G 2 ]−1 χR = ([G 2

p−1 −1 p
p
p p
] G3 +[G 2 ]−1 G 1 )χR

. (A.29)

p

Multiplying G 2 to the left of both sides of Eq. (A.29) and comparing with Eq. A.26 we have:

p

p

p−1 −1
]

G3 + [G 2 ]−1 G 1 ) ,

p

p

p−1 −1
]

G4 + [G 2

p

G2 = G 2 ([G 2

p

p

p

p−1 −1
]

p−1 −1
]

G1 = G 2 ([G 2

p

G2

p−1

G1

(A.30)
).

(A.31)

Eq. (A.31) is equivalent to Eq. (2.64):

p

p

G 2 = G2 [G4 + G 1

p−1

.

(A.32)

p

(A.33)

Inserting Eq. (A.32) into Eq. (A.30) we arrive at Eq. (2.63):

p

p

p

p

p−1 −1
]

G 1 = G1 − G2 [G4 + G 1

G3 .

Comparing Eq. (A.23) and Eq. (A.25) we have:

p

p−1 −1
]

χ1 = [G 4 − G 4
L

p

p

p−1

(G 3 χR − G 3

133

p−1

χR

).

(A.34)

Inserting Eq. (A.34) into Eq. (A.24) we obtain:

p−1

χR

p−1

= G2

p

p−1 −1 p p
] G 3 χR

[G 4 − G 4

p−1

− (G 1

p−1

+ G2

p

p−1 −1 p−1
p−1
] G 3 ) χR

[G 4 − G 4

. (A.35)

Comparing Eq. (A.35) and Eq. (A.27) we have:

p

p

p−1

[G 4 − G 4

p

p−1

− G2

G3 = −G 2

G4 = −G 1

p−1 −1
]

p−1

p

p

G3 ,

(A.36)

p−1 −1
]

[G 4 − G 4

p−1

G3

.

(A.37)

Eq. (A.37) reduces to Eq. (2.66), after trivial algebra :

p

p−1

G4 = G4

p−1

− G3

p

p−1 −1
]

[G4 + G 1

p−1

G2

.

(A.38)

Inserting Eq. (A.38) into Eq. (A.36) we obtain Eq. (2.65):

p

p−1

G3 = G3

p

p−1 −1
]

[G4 + G 1

134

p

G3 .

(A.39)

Appendix B
Detailed derivation of the quadrupole
transition function

B.1

Deriving Eq. (3.32)

In this section, we derive the quadrupole transition operator E2m (Eq. (3.32)) in the hyperspherical harmonics representation. This is an important step in constructing the three-body
transition strength function. First, we apply Eq. (3.31) to Eq. (3.28) and obtain:

5 2
a C (ˆ ) ,
a
4π 1 2m 1
1
√
√
2−λ
λ
2
µ1
µ2
4!
−
x
−
y
2λ!2(2 − λ)!
A1
A1 + A2

Z1 a2 Y2m (ˆ1 ) = Z1
a
1
= Z1

5
4π

λ,µ

× C2−λ,m−µ (ˆ ) Cλ,µ (ˆ ) 2 − λ m − µ λ µ|2 m .
x
y

135

(B.1)

We do the same with Eq. (3.29):

5 2
a
a C (ˆ ) ,
4π 2 2m 2
1
√
√
2−λ
λ
2
µ1
µ2
4!
x
−
y
2λ!2(2 − λ)!
A2
A1 + A2

a
Z2 a2 Y2m (ˆ2 ) = Z2
2
5
4π

= Z2

λ,µ

× C2−λ,m−µ (ˆ ) Cλ,µ (ˆ ) 2 − λ m − µ λ µ|2 m .
x
y

(B.2)

Adding Eq. (B.1) and Eq. (B.2), we arrive at:

Z1 a2 Y2m (ˆ1 ) + Z2 a2 Y2m (ˆ2 )
a
a
1
2
=

5
4π

√

1

λ,µ

2 √
4!
( µ1 x)2−λ
2λ!2(2 − λ)!

−

× C2−λ,m−µ (ˆ ) Cλ,µ (ˆ ) 2 − λ m − µ λ µ|2 m
x
y

λ
µ2
y
A1 + A2

1 2−λ
Z1 −
+ Z2
A1

1 2−λ
.
A2
(B.3)

Due to the properties of angular momentum, λ can only be 0,1,2. Firstly, we calculate the
r.h.s of Eq. (B.3) for λ = 0:

r.h.s =

5
4π

√
0
µ2
−
y
A1 + A2

1

4! 2 √
( µ1 x)2
0!4!

× C2m (ˆ ) C00 (ˆ ) 2 m 0 0|2 m
x
y
= µ1 x2 Y2m (ˆ )
x

Z1
Z
+ 2
2
A1 A2
2

.

136

Z1 −

1 2
+ Z2
A1

1 2
A2
(B.4)

For λ = 1, the r.h.s of Eq. (B.3) becomes:

5
4π

r.h.s =
×

1

4! 2 √
( µ1 x)1
2!2!

√
1
µ2
−
y
A1 + A2

Z1 −

1 1
+ Z2
A1

1 1
A2

1 m − µ 1 µ|2 m C1 m−µ (ˆ ) C1µ (ˆ )
x
y
µ

√
µ1 µ2
6π
5 A1 + A2

=2

Z1
Z
− 2
A1 A2

1 m − µ 1 µ|2 m Y1 m−µ (ˆ ) Y1µ (ˆ ) . (B.5)
x
y

xy
µ

We do the same thing for the case of λ = 2

r.h.s =

5
4π

1
√
2
µ2
4! 2 √
0 −
y
( µ1 x)
4!0!
A1 + A2

× C00 (ˆ ) C2m (ˆ ) 0 0 2 m|2 m
x
y
=

µ2
y2
2
(A1 + A2 )

1 0
Z1 −
+ Z2
A1

1 0
A2

Y2m (ˆ ) [Z1 + Z2 ] .
y

(B.6)

Putting Eq. (B.4), Eq. (B.5), and Eq. (B.6) together we have:

Z1 a2 Y2m (ˆ1 ) + Z2 a2 Y2m (ˆ2 ) + Z3 a2 Y2m (ˆ3 )
a
a
a
1
2
3
Z
µ2
Z1
y 2 [Z1 + Z2 ] Y2m (ˆ )
y
+ 2 Y2m (ˆ ) +
x
2
2
(A1 + A2 )2
A1 A2
√
µ1 µ2
6π
Z1
Z
− 2 xy
1 m − µ 1 µ|2 m Y1 m−µ (ˆ ) Y1µ (ˆ )
x
y
5 A1 + A2 A1 A2
µ

= µ 1 x2
+2
+ Z3

µ2 2
y Y2m (ˆ ) .
y
A2
3

(B.7)

137

Inserting Eq. (B.7) into Eq. (3.27) we obtain an expression for the quadrupole transition
operator in a three-body hyperspherical approach:

E2m = eµ1
+ 2e

Z
Z3
Z1 + Z2
Z1
+ 2 x2 Y2m (ˆ ) + eµ2
x
+
y 2 Y2m (ˆ )
y
2
2
2
A1 A2
A3 (A1 + A2 )2
√
µ1 µ2
6π
Z1
Z
− 2 xy
1 m − µ 1 µ|2 m Y1 m−µ (ˆ )Y1µ (ˆ ) . (B.8)
x
y
5 A1 + A2 A1 A2
µ

For the triple-α problem we have A1 = A2 = A3 = A (4) and Z1 = Z2 = Z3 = Z (2), the
quadrupole transition operator formulae can be simplified:

eZ
A
eZ
=
A

E2m =

B.2

x2 Y2m (ˆ ) + y 2 Y2m (ˆ ) ,
x
y
(ρ sinθ)2 Y2m (ˆ ) + (ρ cosθ)2 Y2m (ˆ ) .
x
y

(B.9)

Deriving Eq. (3.42)

This section present the missing steps in deriving Eq. (3.42). We first take the squared
modulus of both sides of Eq. (3.38) we have:

L M |E2m | LM ; κ
∗
i i
lx ly
ϕ i (θκ )
K

i i
l1x l1y

ϕ i
K

1

2

=
Klx ly K lx ly K1 l1x l1y K l l K i li li K i li li
x y 1 1x 1y
1 1x 1y
∗
κ
κ

(θκ )

Yl i ⊗ Yl i
y LM
x

Yl i ⊗ Yl i
1y
1x

LM

i i
× K lx ly L M |(E2)m | Klx ly K i lx ly LM s
i i i
× K1 l1x l1y L M |(E2)m | |K1 l1x l1y K1 l1x l1y LM ∗ .
s

138

(B.10)

In order to simplify Eq. (B.10), we first consider its angular part in momentum space. The
coupling between two angular momentums can be expressed in term of Clebsch-Gordan
coefficients and spherical harmonics functions:
∗

κ

Yl i ⊗ Yl i
x
y LM

κ

Yl i ⊗ Yl i
1x
1y LM
i
i
i
i
lx mi ly mi |LM ∗ l1x mi l1y mi |LM
x
y
1y
1x

=
mi mi mi mi
x y 1x 1y

× Y ∗ i (Ωkx ) Y ∗ i (Ωky ) Yli mi (Ωkx ) Yli mi (Ωky ) .
i
i
lx mx
ly my
1x 1x
1y 1y

(B.11)

Using Eq. (3.39), Eq. (3.40) and Eq. (3.41) to integrate Eq. (B.10) over dΩκ we obtain:
r

L M |E2m | LM ; κ

2

dΩκ =
r
Klx ly K lx ly K1 l1x l1y K l l K i li li
xy
1 1x 1y







i
i
lx mi ly mi |LM
x
y

mi mi
x y

2

i i i
 K lx ly L M |E2m | Klx ly K lx ly LM s

i i
× K1 l1x l1y L M |E2m | K1 l1x l1y K i lx ly LM ∗ .
s

(B.12)

There exists a sum rule for the Clebsch-Gordan coefficients

i
i
lx mi ly mi |LM
x
y
mi mi
x y

139

2

= 1,

(B.13)

which reduces Eq. (B.12) to:

L M |E2m | LM ; κ

2

dΩκ =
r
Klx ly K lx ly K1 l1x l1y K l l K i li li
xy
1 1x 1y

i i
× K lx ly L M |E2m | Klx ly K i lx ly LM s
i i
× K1 l1x l1y L M |E2m | K1 l1x l1y K i lx ly LM ∗ .
s

(B.14)

Applying further algebraic simplification we arrive at Eq. (3.42):

L M |E2m | LM ; κ

2

dΩκ
5
2

i i
K lx ly L M |E2m | Klx ly K i lx ly LM s .

=
i i
K i lx ly Klx ly K lx ly

140

(B.15)

Appendix C
Inputs for numerical calculations

C.1

Bound state calculation

&fname
n f i l e =’2 bs ’
d e s c= ’ 2 bs : t r i p l e −alpha ’ /
&s c a l e
amn=939. hc =197.3/
&n u c l e i
name= ’ alpha ’ , ’ alpha ’ , ’ alpha ’
mass= 4 4 4
z= 2 2 2
r a d i u s =1.47 1 . 4 7 1 . 4 7 /
&i d e n t i c a l
i d=F , F , T, i s o =0.0 , 0 . 0 , 0 . 0 /
&t o t a l
ngt = 1 , g t o t ( 1 ) = 2 . 0 , g p a r i t y (1)=+1 /
&p a r t i c l e s
ns (1)=1 , s p i n (1 ,1)= 0 . 0 , p a r i t y ( 1 , 1 ) = 1 , e n e r g y ( 1 , 1 ) = 0 . 0 ,
ns (2)=1 , s p i n (1 ,2)= 0 . 0 , p a r i t y ( 1 , 2 ) = 1 , e n e r g y ( 1 , 2 ) = 0 . 0 ,
ns (3)=1 , s p i n (1 ,3)= 0 . 0 , p a r i t y ( 1 , 3 ) = 1 , e n e r g y ( 1 , 3 ) = 0 . 0 /
&em
c o r e k ( 1 ) = 0 . 0 d e f ( 2 , 1 ) = 0 . 0 Qmom( 1 ) = 0 . 0 Mmom( 1 ) = 0 . 0
c o r e k ( 2 ) = 0 . 0 d e f ( 2 , 2 ) = 0 . 0 Qmom( 2 ) = 0 . 0 Mmom( 2 ) = 0 . 0
c o r e k ( 3 ) = 0 . 0 d e f ( 2 , 3 ) = 0 . 0 Qmom( 3 ) = 0 . 0 Mmom( 3 ) = 0 . 0 /
&waves
141

auto (1)=T, kmaxa (1)=10 , lxmax (1)=20 , lymax (1)=20 ,
auto (2)=T, kmaxa (2)=10 , lxmax (2)=20 , lymax (2)=20 ,
auto (3)=T, kmaxa (3)=10 , lxmax (3)=20 , lymax (3)=20/
&poten
d e t a i l ( 1 ) =’2 alpha ’
rcoul (1)=2.94
typc ( 1 ) =’gau ’
ps ( 1 , 1 ) = 1 2 5 . 0 ps ( 2 , 1 ) = 1 . 5 3 ps (3 ,1)= −30.18 ps ( 4 , 1 ) = 2 . 8 5 ,
pd (1 ,1)=20 pd ( 2 , 1 ) = 1 . 5 3 pd (3 ,1)= −30.18 pd ( 4 , 1 ) = 2 . 8 5 ,
pa (1 ,1)= −30.18 pa ( 2 , 1 ) = 2 . 8 5 ,
t y p s o (1)= ’ nul ’
t y p s s (1)= ’ nul ’
t y p t ( 1 ) =’ nul ’
d e t a i l ( 2 ) =’2 alpha ’
rcoul (2)=2.94
typc ( 2 ) =’gau ’
ps ( 1 , 2 ) = 1 2 5 . 0 ps ( 2 , 2 ) = 1 . 5 3 ps (3 ,2)= −30.18 ps ( 4 , 2 ) = 2 . 8 5 ,
pd (1 ,2)=20 pd ( 2 , 2 ) = 1 . 5 3 pd (3 ,2)= −30.18 pd ( 4 , 2 ) = 2 . 8 5 ,
pa (1 ,2)= −30.18 pa ( 2 , 2 ) = 2 . 8 5 ,
t y p s o (2)= ’ nul ’
t y p s s (2)= ’ nul ’
t y p t ( 2 ) =’ nul ’
d e t a i l ( 3 ) = ’ 2 alpha ’
rcoul (3)=2.94
typc ( 3 ) =’gau ’
ps ( 1 , 3 ) = 1 2 5 . 0 ps ( 2 , 3 ) = 1 . 5 3 ps (3 ,3)= −30.18 ps ( 4 , 3 ) = 2 . 8 5 ,
pd (1 ,3)=20 pd ( 2 , 3 ) = 1 . 5 3 pd (3 ,3)= −30.18 pd ( 4 , 3 ) = 2 . 8 5 ,
pa (1 ,3)= −30.18 pa ( 2 , 3 ) = 2 . 8 5 ,
t y p s o (3)= ’ nul ’
t y p s s (3)= ’ nul ’
t y p t ( 3 ) =’ nul ’ /
&pot3b typ3b =’gau ’ , s3b (1)= −15.94 , r3b ( 1 ) = 6 . 0 /
&gamso gamso1 =0 ,0 ,0 , gamso2 =0 ,0 ,0 /
&g r i d s r r =0.3 n l a g =180 n j a c =70 / / Methods :
&t r a c e
/
&continuum docont=F
/
&b 2 s t a t e s n 2 s t a t e s =0 /
&s o l v e eimin = −5.0 , eimax =0, eqn =’T’ nbmax=40 meigs=1 /
EOF
#t i m e r
f a c e 1 3 c 3 < data . $$ > 2 bs . out

142

C.2

Continuum state calculation

&fname
n f i l e =’0 exc ’
d e s c= ’ 0 exc : t r i p l e −alpha ’ /
&s c a l e
amn=939. hc =197.3/
&n u c l e i
name= ’ alpha ’ , ’ alpha ’ , ’ alpha ’
mass= 4 4 4
z= 2 2 2
r a d i u s =1.47 1 . 4 7 1 . 4 7 /
&i d e n t i c a l
i d=F , F , T, i s o =0.0 , 0 . 0 , 0 . 0 /
&t o t a l
ngt = 1 , g t o t ( 1 ) = 0 . 0 , g p a r i t y (1)=+1 /
&p a r t i c l e s
ns (1)=1 , s p i n (1 ,1)= 0 . 0 , p a r i t y ( 1 , 1 ) = 1 , e n e r g y ( 1 , 1 ) = 0 . 0 ,
ns (2)=1 , s p i n (1 ,2)= 0 . 0 , p a r i t y ( 1 , 2 ) = 1 , e n e r g y ( 1 , 2 ) = 0 . 0 ,
ns (3)=1 , s p i n (1 ,3)= 0 . 0 , p a r i t y ( 1 , 3 ) = 1 , e n e r g y ( 1 , 3 ) = 0 . 0 /
&em
c o r e k ( 1 ) = 0 . 0 d e f ( 2 , 1 ) = 0 . 0 Qmom( 1 ) = 0 . 0 Mmom( 1 ) = 0 . 0
c o r e k ( 2 ) = 0 . 0 d e f ( 2 , 2 ) = 0 . 0 Qmom( 2 ) = 0 . 0 Mmom( 2 ) = 0 . 0
c o r e k ( 3 ) = 0 . 0 d e f ( 2 , 3 ) = 0 . 0 Qmom( 3 ) = 0 . 0 Mmom( 3 ) = 0 . 0 /
&waves
auto (1)=T, kmaxa (1)=26 , lxmax (1)=20 ,
auto (2)=T, kmaxa (2)=26 , lxmax (2)=20 ,
auto (3)=T, kmaxa (3)=26 , lxmax (3)=20/
&poten
d e t a i l ( 1 ) =’2 alpha ’
rcoul (1)=2.94
typc ( 1 ) =’gau ’
ps ( 1 , 1 ) = 1 2 5 . 0 ps ( 2 , 1 ) = 1 . 5 3 ps (3 ,1)= −30.18 ps ( 4 , 1 ) = 2 . 8 5 ,
pd (1 ,1)=20 pd ( 2 , 1 ) = 1 . 5 3 pd (3 ,1)= −30.18 pd ( 4 , 1 ) = 2 . 8 5 ,
pa (1 ,1)= −30.18 pa ( 2 , 1 ) = 2 . 8 5 ,
t y p s o (1)= ’ nul ’
t y p s s (1)= ’ nul ’
t y p t ( 1 ) =’ nul ’
d e t a i l ( 2 ) =’2 alpha ’
rcoul (2)=2.94
typc ( 2 ) =’gau ’
ps ( 1 , 2 ) = 1 2 5 . 0 ps ( 2 , 2 ) = 1 . 5 3 ps (3 ,2)= −30.18 ps ( 4 , 2 ) = 2 . 8 5 ,
pd (1 ,2)=20 pd ( 2 , 2 ) = 1 . 5 3 pd (3 ,2)= −30.18 pd ( 4 , 2 ) = 2 . 8 5 ,
pa (1 ,2)= −30.18 pa ( 2 , 2 ) = 2 . 8 5 ,
143

t y p s o (2)= ’ nul ’
t y p s s (2)= ’ nul ’
t y p t ( 2 ) =’ nul ’
d e t a i l ( 3 ) = ’ 2 alpha ’
rcoul (3)=2.94
typc ( 3 ) =’gau ’
ps ( 1 , 3 ) = 1 2 5 . 0 ps ( 2 , 3 ) = 1 . 5 3 ps (3 ,3)= −30.18 ps ( 4 , 3 ) = 2 . 8 5 ,
pd (1 ,3)=20 pd ( 2 , 3 ) = 1 . 5 3 pd (3 ,3)= −30.18 pd ( 4 , 3 ) = 2 . 8 5 ,
pa (1 ,3)= −30.18 pa ( 2 , 3 ) = 2 . 8 5 ,
t y p s o (3)= ’ nul ’
t y p s s (3)= ’ nul ’
t y p t ( 3 ) =’ nul ’ /
&pot3b typ3b =’gau ’ , s3b (1)= −19.46 , r3b ( 1 ) = 6 . 0 /
&gamso gamso1 =0 ,0 ,0 , gamso2 =0 ,0 ,0 /
&g r i d s r r =0.3 n l a g =3000 n j a c =100 r s c r e e n =800 a s c r e e n =10 / Methods :
&t r a c e
/
&continuum docont=F
/
&b 2 s t a t e s n 2 s t a t e s =0 /
&s o l v e eimin = −6.0 , eimax =0, eqn =’T’ nbmax=0 meigs=0 /
EOF
#t i m e r
f a c e 1 3 c 3 c < data . $$>0 e x c f a c e . out
$mefile
FPOT = ’ 0 exc ’
$end
$sturmin
Eo ( 1 ) = 0 . 0 0 egsmin =−6.0 egsmax =3.0 sturm (1)= ’F ’ ngsmax=0
s t e p i = 0 . 0 5 , rmax=50 , s t e p s =0.25 , eeigmax =5000.0 l p r w f=3 wfrmax=200
nsturm =50 , b u t t l e =5 r a f i n =3000 meigs=0 maxnop=100
ebeg =0.3525 e s t e p =0.0005 emax=0.5 wfs=T numax=5 s t r f u n=F
n o g s c a t=F e l i n e a r=T betalw=1e −10
f l w p =’/tmp/0 exc . lwra ’ , f l w p r =’/tmp/0 exc . lwrmats ’
$end
$bstate
k p o l e = 1 , emin = 1 . 0 , emax=0.0 ne = 30 numaxgs=10
$end
’EOF’
#t i m e r
sturmxx96 < data . $$ > 0 exc . out

144

Appendix D
Highlights on the (d,p) reaction
projects
Transfer reactions have been considered as an effective tool to extract the single-particle
structure and the properties of nuclei across the nuclear chart [93, 94, 95, 96]. Systematic
studies on (d,p) reactions have provided a better understanding about the properties of nuclear interactions. For example, the work [96] on Sn isotopes using single-nucleon transfer
reactions discovers the reduction of the spin-orbit strength when approaching the neutron
drip line.
Transfer reactions also play an important role in astrophysics since they provide an indirect method to measure the (n, γ) rates which are crucial for the r-process nucleosynthesis.
Neutron capture reactions (n, γ) on unstable nuclei are not possible to measure in the laboratory because both the beam and the target are unstable. Recently, the neutron transfer
(d,p) inverse kinematics reactions on unstable targets has become possible due to the development of rare isotope facilities and detector technical advances [97, 98, 99]. Studies in

145

[100, 101] show that the neutron capture rates can be derived indirectly from (d,p) reactions.
Significant effort has been invested in reducing the theoretical uncertainty of this indirect
method [102, 103].
In the late 70s the single-particle structure of almost all stable nuclei has been successfully mapped out using (d,p) transfer reactions. Nuclear science has now shifted toward
unstable rare isotopes, and this imposes severe challenges. The inverse kinematics technique is employed to study the structure and properties of exotic nuclei. Pioneering studies
[97, 98, 99, 104, 105, 106, 107] showed promising results in applying this technique to unstable nuclei. With the development of new-generation rare-isotope facilities where beam rates
will be enhanced, our possibility to have a detailed understanding of the exotic nuclei and
the underlying nucleon-nucleon force is highly increased. In parallel with the experimental
development, a well founded theory for (d,p) transfer reactions is crucial.
The traditional method to study transfer reactions is the distorted wave Born approximation (DWBA) [108, 109, 110] in which the T-matrix is truncated after the first term in
the Born series. This method does not take into account deuteron break-up since it uses the
deuteron-target elastic scattering wavefunction as an input [111]. At energies above 2.23 MeV
(deuteron binding energy), deuteron tends to breakup as an intermediate step in transfer
processes. For these cases, the break-up effects can be significant and the DWBA method is
likely to fail to describe the data. In order to tackle this shortcoming of the DWBA method,
Johnson and Soper [112] introduced the adiabatic distorted wave approximation (ADWA)
framework which provides a practical prescription for including the deuteron break-up effects. In [112] the zero-range approximation is made for the n-p interaction, and this reduces
the transfer matrix element to a form similar to that in the DWBA method. However, in
the ADWA method the incoming wave is the adiabatic three-body wavefunction instead of
146

the elastic deuteron wavefunction. Using this method to analyze single-nucleon transfer reactions one often finds improvement in the description of experimental data [113]. Another
formal development on the ADWA method was later introduced by Johnson and Tandy [114]
to take into account the finite-range effects of the n-p interaction. However, the numerical
implementation and application has not been systematically explored.
In the past, the zero-range (ZR-ADWA) approximation was used to perform (d,p) calculations due to computational limitations. The finite-range effects were then corrected by the
local energy approximation (LEA) procedure which was developed by Buttle and Goldfarb
[115]. The resulting codes persist to this day. It is not clear that this approximated correction still holds at high deuteron beam energies. In this appendix we will address this issue
and further explore the implication of the FR-ADWA method.
Appendix D is organized into six sections. In section D.1, we discuss the physical quantities that are commonly investigated in (d,p) studies. Section D.2 presents the theory of
the ADWA method. We then summarize our results on the study of the finite-range effects
in (d,p) reactions in section D.3. The dispersive optical model is tested for (d,p) reactions
using the FR-ADWA framework in section D.4. An application of the finite-range adiabatic
distorted wave approximation (FR-ADWA) method for extracting the normalization coefficients from the 14 C(d,p)15 C reaction is presented thereafter in section D.5. In the last
section (D.6) of this appendix, we draw the general conclusions from our (d,p) studies.

D.1

Common physical quantities in (d,p) studies

Deuteron transfer reactions are considered as one of the most effective methods to study the
single-particle structure of nuclei. A quantity of relevance for this process, in which nucleus

147

B is formed by the transfer of a neutron from deuteron projectile to the target nucleus A
(see Fig. D.1) is the spectroscopic factor, defined as the square of the norm of the overlap
lj

function φI :I (r):
A B
IA I
Slsj B =

lj

|φI :I (r)|2 dr ,
A B

(D.1)

where,
lj

φI :I (r) = ΦA (ξA )|ΦB (ξA , r) .
IA
IB
A B

(D.2)

Here ΦA and ΦB are the many-body wave functions of nucleus A and B, respectively. IA
IB
IA
(IB ) is the spin state of nucleus A (B) and ξA denotes the internal A coordinates. lsj refers
to the quantum state of the transfered neutron. r is defined as the vector radius from the
IA I
target nucleus A to the neutron. The spectroscopic factor Slsj B therefore measures the

probability to form a nucleus B in quantum state IB having a single-particle structure with
core nucleus A in quantum state IA and valence neutron n in quantum state lsj. Since at
large radius (r −→ ∞) the overlap function goes to zero, contribution to the integration
in Eq. (D.1) comes mostly from the nuclear interior (r < rN ). The spectroscopic factor is
dominated by the radial behavior of the overlap function within the nucleus, thus probing
reactions at intermediate and high energies.
Another structure quantity relevant in peripheral (d,p) reactions is the asymptotic norIA I
malization coefficient (ANC) Clsj B which determines the strength of the tail of the overlap

function:
r→∞ IA I
lj
φI :I (r) −→ Clsj B W
(−2kB r) ,
−η,l+ 1
A B

(D.3)

2

where, W

−η,l+ 1
2

is the Whittaker function; the Coulomb Sommerfeld parameter η is taken

as zero in our study because the transfered particle is neutral; kB relates to the binding

148

energy B of nucleus B through kB =

√

2mB B / ; mB is the mass of the nucleus and

is

Planck’s constant. The ANC reflects the peripheral properties of the bound state, therefore
it is the relevant ingredient in evaluating reactions occurring at low energies.
lj

In order to obtain the overlap function φI :I (r), one needs to solve a many-body probA B
lem which is highly complex and numerically difficult. Traditionally in (d,p) studies one
replaces the many-body overlap function by a single-particle overlap function:

IA I (sp) 1/2
lj(sp)
lj
] ϕ(lsj)I (r) ,
φI :I (r) ≈ φI :I (r) = [Slsj B
B
A B
A B

(D.4)

where ϕ(lsj)I (r) is the normalized single-particle bound state wave function of nucleus B.
B

IA I (sp)
is
Since ϕ(lsj)I (r) is normalized to one, the single-particle spectroscopic factor Slsj B
B
IA I
equal to the many-body spectroscopic factor Slsj B in Eq. (D.1) if the bound state wave

function and the overlap function have similar radial behavior in both the interior and
exterior of the nucleus. The asymptotic behavior of the bound state wave function is given
by:
r→∞ I

B
(−2kB r) ,
ϕ(lsj)I (r) −→ blsj W
−η,l+ 1
B

(D.5)

2

I

B
where blsj is the single-particle ANC which defines the amplitude of the tail of the neutron

bound state wave function. From Eq. (D.3), Eq. (D.4), and Eq. (D.5) we can derive a
relationship between the single-particle spectroscopic factor, the ANC, and the single-particle
ANC:
IA I (sp)
Slsj B

=

IA I
(Clsj B )2
IB
(blsj )2

.

(D.6)

In practice, an experimental spectroscopic factor is determined by the normalization of the
calculated cross section to the experimental cross section at the first peak of the angular

149

distribution:
S exp =

dσ
dΩ exp
dσ
dΩ th

.

(D.7)

While S exp is not an observable, dσ/dΩ is.
Because ANC is a peripheral quantity, probing the region outside the nuclear radius
rN , it is best extracted using reactions at very low energies. This quantity should not
depend on the choice of the neutron-target binding potential. The spectroscopic factor on
contrary contains the nuclear interior information, therefore is commonly determined from
high-energy reactions. As indicated in Eq. (D.6), the spectroscopic factor is sensitive to the
assumptions made on the geometry of the neutron-target binding interaction, which leads
to large ambiguities. The work in [116] provides a method to extract both the ANC and the
spectroscopic factor in a more consistent and less ambiguous way. In our (d,p) study, we
will often refer to these two quantities.

D.2

The adiabatic distorted wave approximation

Let us consider a (d,p) transfer reaction: d + A −→ B + p. The incoming and outgoing
channel for this reaction are presented in Fig. D.1. In the adiabatic model of [112, 114], a
three-body problem is considered for the n + p + A system in which the deuteron scattering
wavefunction Ψ(+) (r, R) in the incident channel is obtained by solving the inhomogeneous
differential equation:

[E + i − Tr − TR − UnA − UpA − Vnp ]Ψ(+) (r, R) = ı φd (r) exp(ıKd .R) ,

150

(D.8)

Figure D.1: (d,p) transfer reaction: (a) incoming channel, (b) outgoing channel
where
2

Tr = −

2µnp
2

TR = −

2µdA

2
r
2
R

(D.9)
(D.10)

In Eq. (D.8,D.9,D.10), we define Tr and TR as the kinetic energy operators describing the
n-p relative motion and the motion of its center of mass relative to the target. µnp and µdA
are the reduced mass of the n-p and d-A systems, respectively. r = rp − rn [R = (rn + rp )/2]
is the relative coordinate (center-of-mass coordinate) of the n-p system with rp and rn being
the neutron and proton coordinates relative to the center of mass of target A. We take
UnA (rn ), UpA (rp ) and Vnp (r) to be the neutron-target, proton-target and neutron-proton
interactions, respectively. The term proportional to ı on the r.h.s of Eq. (D.8) ensures that
there is an incoming wave only in the deuteron channel.
In this approach, the exact T matrix is:

151

(−)

T =< φnA χpB |Vnp + UpA − UpB |Ψ(+) > ,

(D.11)

(−)

where φnA is the neutron-target bound state wavefunction, and χpB is the proton scattering
distorted wave in the outgoing channel. Note that the T matrix in Eq. (D.11) is exact; the
three-body deuteron wavefunction Ψ(+) (r, R) of Eq. (D.8) is used instead of the deuteron
elastic scattering wave as in DWBA. In Eq. (D.11) the remnant term UpA − UpB is often
small and thus can be neglected except for light nuclei. If that is valid, then it is noted from
Eq. (D.11) that the exact solution of the three-body wavefunction Ψ(+) (r, R) in Eq. (D.8)
is only necessary within the range of the n-p interaction. This important realization in both
[112] and [114] allows the construction of an adiabatic method for (d,p) reactions.

D.2.1

Zero-range adiabatic wave approximation (ZR-ADWA)

This method is developed by Johnson and Soper [112] in which the three-body wavefunction
Ψ(+) (r, R) is expanded in term of the eigenstates of the n-p Hamiltonian:

Ψ(+) (r, R) = φd (r)χd (R) +

(+)

dkφk (r)χk (R) ,

(D.12)

where the wave function φ(r) describes the two nucleon motion, while χ(R) is the wave
function of the n + p center of mass and target A system. The scattering states φk (r)
of the n-p system and the deuteron bound state φd (r) form a complete set since they are
eigensolutions of the n-p Schrodinger equation:

(Tr + Vnp )φd (r) = −εd φd (r) ,
(Tr + Vnp )φk (r) = +εk φk (r) .

152

(D.13)

In Eq. (D.13), εd and εk are the bound and scattering eigenenergies for the n-p system. In
ZR-ADWA, the deuteron excitation energies εk + εd is assumed to be much smaller than the
deuteron kinetic energy E. Inserting Eq. (D.12) into Eq. (D.8) in the zero-range limit of the
n-p interaction and using the adiabatic assumption, we arrive at an optical model equation:

(E + d − TR − UJS (R))χJS (R) = 0 ,
d

(D.14)

where the effective potential UJS does not describe deuteron elastic scattering, but rather
incorporates deuteron breakup effects within the range of Vnp :

UJS (R) = UnA (R) + UpA (R) .

(D.15)

Within this model, the transfer amplitude reduces to

TJS = D0

dR φ∗ (R)χ∗ (R)χJS (R)
nA
pB
d

(D.16)

where D0 is the zero-range constant of the deuteron:

D0 =

D.2.2

dr Vnp (r) φd (r) .

(D.17)

Finite-range adiabatic wave approximation (FR-ADWA)

This method is developed by Johnson and Tandy [114] in order to take into account the
finite-range effects of the n-p interaction. Instead of using a set of Hnp eigenstates for
the wave function expansion as done in the ZR-ADWA method, they employ the Weinberg
representation. Since the exact three-body wavefunction is only needed within the range of
153

Vnp , the Weinberg basis is an excellent choice because they form a complete square integrable
set over this range. The Weinberg basis is defined as:

(Tr + αi Vnp )φi (r) = −εd φi (r) ,

(D.18)

where d is the deuteron binding energy and αi are the eigenvalues which are real and
increase with i = 1, 2, ... Here φi (r) denote the Weinberg states and satisfy the orthonormality
relation:
φi |Vnp |φj = −δij .

(D.19)

The three-body wavefunction is then expanded as:
∞

Ψ(+) (r, R) =

φi (r)χi (R).

(D.20)

i=1

Inserting this expansion into the three-body equation Eq. (D.8), one arrives at a set of
coupled channel equations:

¯
[E + i − TR − Uii (R)]χi (R) = ı δi1 Nd exp(ıKd .R) +

¯
Uij (R)χj (R).

(D.21)

j=i

¯
In Eq. (D.21), one defines the coupling potentials Uij (R) as:

¯
Uij (R) = Uij (R) + βij (αj − 1) ,

154

(D.22)

where Uij (R) and βij are derived from the nucleon optical potentials UnA and UpA and the
n-p interaction Vnp , in a complex folding procedure:

Uij (R) = − φi |Vnp (UnA + UpA )|φj ,
2
βij = φi |Vnp |φj .

(D.23)

The normalization coefficient Nd appearing on the r.h.s of Eq. (D.21) is computed as:

Nd = − φ1 |Vnp |φd .

(D.24)

It is not trivial to solve Eq. (D.21) exactly as seen in [117]. However if only the first term
in the Weinberg expansion Eq. (D.20) is needed to describe the three-body wave function in
the range of Vnp , Eq. (D.21) can be reduced to a single-channel optical-model like equation:

(E + d + iε − TR − U11 (R))χJT (R) = iεNd exp(ıKd .R) .
1

(D.25)

In this case, the effective deuteron potential U11 (R) is calculated in a more elaborate manner
in comparison to that in the ZR-ADWA method:

U11 (R) = − φ1 (r)|Vnp (UnA + UpA )|φ1 (r) .

(D.26)

The first Weinberg state φ1 is different from the deuteron ground state wavefunction φd by
a normalization coefficient Nd : | φd >= Nd | φ1 >. Therefore one can rewrite the potential

155

U11 (R) in terms of φd :

U11 (R) ≡ UJT (R) =

φd (r)|Vnp (UnA + UpA )|φd (r)
.
φd (r)|Vnp |φd (r)

(D.27)

When taking into account the finite-range effects of the n-p interaction, the transfer amplitude Eq. (D.11) now has a complex 6-dimensional integral formula:

TJT =
=

D.3

(−)

φnA χpB |Vnp + UpA − UpB |φ1 χJT (R)
1
(−)

φnA χpB |Vnp + UpA − UpB |φd χJT (R)/Nd .
1

(D.28)

Finite-range effects in (d,p) reactions

In order to explore the finite-range effects we perform a systematic study of 26 (d,p) reactions on stable targets [118], involving nuclei with masses ranging A = 12−208 and deuteron
energies Ed = 2 − 70 MeV within the framework of ADWA. For all calculations throughout
this appendix, we use the Reid soft-core potential [119] as the n-p interaction. The nucleon optical potentials are taken from the Chapel Hill global parameterization [120]. The
neutron-target final bound state is produced by a Woods-Saxon potential with the standard
radius r = 1.25 fm and diffuseness a = 0.65 fm. The depth of this potential is adjusted to
reproduce the known neutron separation energy. We employ a subroutine contained in the
code TWOFNR [121] to generate the adiabatic potential for the deuteron in the incident
channel. The cross section for each case is obtained using the code FRESCO [122].
The finite-range effects introduced by the n-p interaction are presented in both the adiabatic deuteron potential Eq. (D.27) and the evaluation of the transfer matrix element

156

Eq. (D.28). By comparing UJT (Eq. (D.27)) and UJS (Eq. (D.15)), TJT (Eq. (D.28)) and
TJS (Eq. (D.16)) we can disentangle the separate effects of the finite-range. We also study
the accuracy of the LEA method.

D.3.1

Finite-range effects in the potentials

Before looking at the changes in cross section due to the finite-range effect, we analyze
the effect on the adiabatic potentials by comparing UJS and UJT . In order to simplify
the comparison, the real part and imaginary part of the potential are fitted to a volume
and surface Woods-Saxon shape, respectively. The finite-range effects in the potentials
are obtained by comparing the Woods-Saxon parameters of UJS and UJT . The results
are presented in Table D.1. We also compare our results with an approximate method
developed by Wales and Johnson [123] to estimate the finite-range effects in the deuteron
potential (denoted UW J ). In Table D.1, WJ and JT represent the percentage difference in
the potential of the Wales and Johnson method and finite-range adiabatic approach, relative
to the zero-range adiabatic prescription. The main difference between the finite-range and
zero-range potentials in the adiabatic method is a constant increase in diffuseness and a
slight systematic decrease in radius. In the Wales and Johnson study, although the radius is
fixed, the main features of the finite-range effects in potentials are captured. The diffuseness
also increases but to a lesser extent. We find that the differences in the depths of the real
and imaginary parts of the potentials in the adiabatic approach are more subtle and vary
case to case while Wales and Johson obtains a small systematic decrease in the depth of the
potential.

157

Table D.1: The finite-range effects on the deuteron distorting potential are presented as the
changes in the diffuseness aR and aI of the real and imaginary parts, the depths of the real
and imaginary parts V and Ws , as well as the corresponding radii rR and rI .
target

parameter

WJ

JT

all

∆aR

+4%

+7%

∆aI

+3%

+8-9%

∆V

-5.6%

-1.98%

∆rR

0%

-1.25%

∆Ws

-4.6%

-4.52%

∆rI

0%

+0.72%

∆V

-2.1%

-0.04%

∆rR

0%

-0.93%

∆Ws

-3.7%

+1.6%

∆rI

0%

-0.97%

∆V

-0.7%

+0.06%

∆rR

0%

-0.35%

∆Ws

-3.3%

+1.2%

∆rI

0%

-0.35%

12 C

48 Ca

208 Pb

158

D.3.2

Finite-range effects in the transfer cross sections

In this section, we study the differences in the cross section for several model calculations
to disentangle the finite-range effects in the deuteron potential and in the evaluation of
the transfer matrix element. All the cross section calculations are performed by the code
FRESCO [122] with the interaction inputs generated by the code TWOFNR [121]. We
consider four different types of calculation for each reaction: UJS deuteron potential and
zero-range transfer matrix Eq. (D.16) (ZR-JS, the reference), UJS deuteron potential and
finite-range transfer matrix Eq. (D.28) (FR-JS), UJT deuteron potential and finite-range
transfer matrix Eq. (D.16) (FR-JT) and a zero-range calculation for both deuteron potential
and transfer matrix with a local energy approximation (LEA) finite-range correction [115].
Since LEA method has been widely used in the past, it is important to study the validity of
this approximation at different deuteron beam energies.
We summarize our results for all 26 reactions in Table D.2. The last four columns
present different aspects of the finite-range effects. ∆(LEA) shows the effect of finite-range
introduced by LEA to the evaluation of the transfer matrix. ∆(FR-JS) represents the effect
of finite-range only in the calculation of the transfer matrix. ∆(FR-JT) shows the total
finite-range effects in both the deuteron adiabatic potential and the evaluation of the transfer
matrix. ∆(JS-JT) is the finite-range effects coming only from the adiabatic potential. All the
percentage differences are evaluated at the first peak of the angular distribution for energies
above the Coulomb threshold and at the backward angles otherwise. The peak angles are
also presented in Table D.2. In all cases, we compute the percentage differences of cross
sections relative to the zero-range Johnson and Soper method in which UJS and Eq. (D.16)
are used.

159

Table D.2: Finite-range effects in transfer cross sections. Comparisions are performed with
respect to the zero-range Johnson and Soper calculation (cases with no existing data are
presented by *).
Target

Ed

θ

∆(LEA)

∆(FR-JS)

∆(FR-JT)

∆(JT-JS)

12 C

4

25

+5.6%

+5.5%

+4.5%

-1.0%

12 C

12

13

+2.6%

+2.9%

-1.5%

-4.3%

12 C

19.6

10

+11%

+13%

+7.7%

-4.2%

12 C

56

6

-37%

-27%

-36%

-12%

48 Ca

2

180

+6.5%

+6.3%

+2.6%

-3.5%

48 Ca

13

12

+4.9%

+3.8%

-2.8%

-6.2%

48 Ca

19

8

+5.0%

+4.0%

-0.30%

-4.1%

48 Ca*

30

4

+7.3%

+4.8%

-2.3%

48 Ca*

40

0

-5.4%

-5.9%

-10%

48 Ca*

50

0

-1.9%

-19%

-18%

48 Ca

56

0

-5.2%

-6.5%

-24%

-18.6%

69 Ga

12

14

+4.3%

+4.7%

-1.1%

-5.49%

86 Kr

11

25

+4.8%

+5.5%

-0.40%

-5.63%

The results in Table D.2 can be divided into three groups. The first group is the subCoulomb reactions which occur at very low energies under the Coulomb barrier. In this
case, we find the finite-range effects are only a few percent and mostly come from the evaluation of the transfer matrix. The contribution from the adiabatic potential is small. In
160

Table D.2 (cont’d)
Target

Ed

θ

∆(LEA)

∆(FR-JS)

∆(FR-JT)

∆(JT-JS)

90 Zr

2.7

138

+6.2%

+7.3%

+5.5%

-1.7%

90 Zr

11

26

+5.4%

+5.0%

-0.90%

-5.6%

124 Sn

5.6

175

+6.1%

+11%

+7.5%

-2.8%

124 Sn

33.3

0

+2.9%

+4.6%

0%

-4.4%

124 Sn*

40

12

-1.1%

-2.4%

-1.4%

124 Sn*

50

11

-3.9%

-4.3%

-0.44%

124 Sn*

60

9

-11%

-30%

-21%

124 Sn*

70

0

+5.1%

-29%

-44%

-21%

208 Pb

8

180

+6.1%

+7.2%

+6.1%

-0.96%

208 Pb

12

98

+5.7%

+8.8%

+2.2%

-6.1%

208 Pb*

20

30

+4.5%

-2.3%

-6.6%

208 Pb*

40

9

+1.4%

-6.9%

-8.1%

208 Pb*

60

0

+0.14%

-8.8%

-9.0%

208 Pb*

80

0

-62%

-86%

-63%

161

Fig. D.2 we plot two examples:(a) 48 Ca(d,p)49 Ca at Ed = 2 MeV and b) 208 Pb(d,p)209 Pb
at Ed = 8 MeV. The angular distributions of these two reactions peak at backward angles as
expected. Since extracting the spectroscopic information is not the purpose of this work, the
data presented here is just to indicate that the ingredients of our model are realistic. The
finite-range effects are 3% in 48 Ca and 6% in 208 Pb. The LEA corrections to the finite-range
effects for sub-Coulomb reactions are below 5%. The LEA method is able to capture most
of the finite-range effects in this low energy regime.
We then look at the reactions that happen at intermediate deuteron beam energies (1020 MeV). In Fig. D.3 we present the angular distributions for 86 Kr(d,p)87 Kr at Ed = 11 MeV
(a) and 208 Pb(d,p)209 Pb at Ed = 20 MeV (b). The finite-range effect from the deuteron potential reduces the cross section while that from the evaluation of the transfer matrix element
increases the cross section. Although the two separate finite-range effects are significantly
large, the overall effect is small for these cases due to cancellation. For example, we obtain
the overall finite-range effect of −1% for the 69 Ga at Ed = 12 MeV, 0.4% for the 86 Kr at
Ed = 11 MeV and 6% in 208 Pb at Ed = 20 MeV. In this energy regime, the LEA method
starts to deviate from the full finite-range calculation.
Finally we study the reactions at high deuteron beam energies (50-80 MeV). There
are only two studied cases in which the experimental data at these energies are available:
12 C(d,p)13 C

at Ed = 56 MeV and 48 Ca(d,p)49 Ca at Ed = 56 MeV. The angular distribu-

tions for these two reactions are plotted in Fig. D.4. In order to ensure that our results are
not biased by low-mass systems, we perform the same calculations for 124 Sn and 208 Pb at
high deuteron energies. The dominant contribution to the finite-range effects in this case
comes from the adiabatic potential which largely reduces the cross section. We find that the
finite-range effects from the evaluation of the transfer matrix is in the same direction but
162

dσ/dΩ (mb/srad)

0.5
0.4

ZR
LEA
FR-JS
FR-JT

0.3
0.2
(a)

0.1
0.0
0

20 40 60 80 100 120 140 160 180

θ (degrees)

dσ/dΩ (mb/srad)

0.8

0.6

ZR
LEA
FR-JS
FR-JT

0.4

0.2
(b)
0.0
0

20 40 60 80 100 120 140 160 180

θ (degrees)
Figure D.2: Finite-range effects in sub-Coulomb (d,p) reactions: (a) 48 Ca(d,p)49 Ca(g.s.)
Ed = 2 MeV (data from [124]) and (b) 208 Pb(d,p)209 Pb(g.s.) Ed = 8 MeV (data from [125]).

163

dσ/dΩ (mb/srad)

15
ZR
LEA
FR-JS
FR-JT

10

5

0
0

(a)

20

40

60

80

100

θ (degrees)

dσ/dΩ (mb/srad)

6.0
ZR
LEA
FR-JS
FR-JT

4.0

2.0

0.0
0

(b)

50

100

150

θ (degrees)
Figure D.3: Finite-range effects in (d,p) reactions at energies slightly above the Coulomb barrier: (a) 86 Kr(d,p)87 Kr(g.s.) Ed = 11 MeV (data from [126]) and (b) 208 Pb(d,p)209 Pb(g.s.)
Ed = 20 MeV (data from [127]).

164

dσ/dΩ (mb/srad)

15

10

ZR
LEA
FR-JS
FR-JT

5
(a)
0
0

20

40

θ (degrees)

80

60

dσ/dΩ (mb/srad)

8.0
ZR
LEA
FR-JS
FR-JT

6.0
4.0
2.0
0.0
0

(b)
10

20

30

40

50

θ (degrees)
Figure D.4: Finite-range effects in (d,p) reactions at high energies: (a) 12 C(d,p)13 C(g.s.)
Ed = 56 MeV (data from [128]) and (b) 48 Ca(d,p)49 Ca(g.s.) Ed = 56 MeV (data from
[129]).

165

10

∆(%)

-10
-30

12

C
Ca
124
Sn
208
Pb
48

-50
-70
0

20

40

60

80

Ed(MeV)
Figure D.5: Systematic finite-range effect as a function of beam energy. The effect from
the adiabatic potentials is presented as open symbols while that from the evaluation of the
matrix elements is plotted as filled symbols.
smaller. The overall finite-range effects are important for (d,p) reactions at high deuteron
beam energies and should not be neglected. The effects are about −36% for 12 C at 56 MeV
beam, −24% for 48 Ca at 56 MeV beam and −43.5% for 124 Sn at 70 MeV beam. We find
that the LEA method is no longer accurate to estimate the finite-range correction in this
high energy regime. The results from the LEA method largely deviate from and sometimes
are in the opposite direction of the full finite-range calculation.
In Fig. D.5 we plot the two separate finite-range effects as a function of the beam energy for four different (d,p) reactions. The solid symbols represent the finite-range effects
attributed to the evaluation of the transfer matrix element while the open symbols indicate
the finite-range effects coming from the deuteron potential. This figure summarizes all the
main features of our study. We can see that the finite-range effects are less than 10% for
reactions at low energies. These effects become more important as the deuteron beam energy
166

increases. We estimate the transition at which finite-range effects cannot be neglected to be
around 20 MeV/u for lighter systems (A < 50) and 30 MeV/u for the heavy systems. We
also obtain similar results when performing a systematic study of (d,p) reactions in which
the outgoing target populates excited states [130, 131].

D.4

Testing the dispersive optical model for (d,p) reactions

In the conventional approach to study (d,p) reactions, a global deuteron optical potential
fitted to elastic scattering is employed for the deuteron-target interaction in the entrance
channel. Similarly, a global proton optical potential fitted to p+A elastic scattering is used
for the proton-target interaction in the exit channel. In our framework (the FR-ADWA
method), the deuteron adiabatic potential is constructed from the neutron and proton optical potentials (Eq. (D.27)). Although these nucleon optical potentials are less ambiguous
than the deuteron optical potentials in the DWBA framework and the deuteron binding
is well constrained, there remains concern about the effective interaction used for the final
neutron bound state. Currently VnA is modeled by an attractive single-particle potential
with Woods-Saxon shape in the FR-ADWA method. We employ standard parameters for
the radius r = 1.25 fm and the diffuseness a = 0.65 fm. The depth of this effective potential
is then adjusted to reproduce the experimental binding energy of the final nucleus. However,
the procedure of choosing the neutron binding potential is not connected to the nucleon
optical potentials mentioned above. The disconnect between these two types of interaction
is unsatisfying.
Recently, a group of scientists at Washington University have successfully implemented
167

the dispersive optical model (DOM) [132, 133, 134] which naturally connects the continuum
and bound state information of the nucleon-target interaction, thus overcoming the inconsistency of the standard applications of the FR-ADWA method. As opposed to the global
optical model where only elastic scattering data are fitted, the DOM method fits simultaneously the nucleon elastic, total and reaction cross section data, as well as bound-state
properties extracted from (e, e p) experiments. The nonlocality correction is also implemented in the fit. The study in [134] observes a nucleon asymmetry dependence of the DOM
potentials, opening a possibility for an extrapolation to exotic nuclei. Because of its very
promising features, we want to explore the validity of the DOM method for (d,p) transfer
reactions using the FR-ADWA framework [135]. In this work we study several cases of doubly magic nuclei for which the experimental data exist: 40 Ca(d,p)41 Ca at Ed = 20 and 56
MeV, 48 Ca(d,p)49 Ca at Ed = 2, 13, 19.3 and 56 MeV, 132 Sn(d,p)133 Sn at Ed = 9.46 MeV,
and 208 Pb(d,p)209 Pb at Ed = 8 and 20 MeV. We used the Reid potential for Vnp [119] and
the DOM potentials for all nucleon-target interactions. The results are also compared with
those produced using the CH89 nucleon optical potential [120] and the standard neutron
single-particle binding interaction.

D.4.1

Details of the calculation

In this work we employ the DOM fits for the Ca, Sn and Pb isotopes [134]. They include
both nucleon scattering and bound state data for Z = 20, 28 and N = 28 nuclei in the fit
for the Ca isotopes. The DOM fit for Sn isotopes includes the proton and 112−124 Sn elastic
scattering data and the neutron elastic-scattering data on 116 Sn, 118 Sn, 120 Sn, and 124 Sn.
Data for nuclei with Z = 82 are used for the DOM fit of Pb isotopes. For this study, the
DOM group provides us all optical potentials necessary to describe A(d, p)B reaction in the
168

Table D.3: Comparison of the DOM overlap function, corrected for non-locality ϕDOM with
the Woods-Saxon single-particle wavefunction ϕWS . In this table we include: n (principal
number), j (the angular momenta of the valence orbital), Sn (the separation energy), Rrms
(the root-mean-square radius of the valence orbital), and bnlj (the modulus of the singleparticle ANC).

Nucleus
41 Ca
49 Ca
133 Sn
209 Pb

Overlap n j
ϕWS
0f7/2
ϕDOM
ϕWS
1p3/2
ϕDOM
ϕWS
1f7/2
ϕDOM
ϕWS
1g9/2
ϕDOM

Sn [MeV] Rrms [f m] |bnlj | [fm−1/2 ]
3.985
2.285
8.362
3.965
2.261
4.606
5.818
5.146
4.820
6.098
6.080
0.844
2.469
6.513
0.831
6.498
1.650
3.936
6.746
1.827

FR-ADWA framework, namely UnA , UpA , UpB , and the neutron binding interaction VnA .
In Table D.3 we present the properties of the final neutron-target bound state. The DOM
potentials generated for Ca isotopes have a radius parameter of 1.18 fm while that of the
standard Woods-Saxon (WS) is 1.25 fm. For the 41 Ca case, the nonlocality correction in the
DOM cancels this difference causing the DOM overlap function ϕDOM to be almost identical
to the Woods-Saxon overlap function ϕWS , thus the root-mean-square radius Rrms from the
two models are exactly the same. The same results are not seen in other cases where a node
is present in the overlap function (n = 0). We find that the nonlocality correction is more
pronounced when the number of nodes in the overlap function differs from zero, making the
DOM rms radius larger than its WS counterpart. An illustration of this effect is plotted for
the 49 Ca in Fig. D.6. Here, the nonlocality increases the single-nucleon overlap function at
large radii.
In Fig. D.7 we present the difference between the DOM and CH89 nucleon optical
potential evaluated at 4.7 MeV for 132 Sn. For both the real part and the imaginary part of

169

WS
DOM

-3

|ϕ(r)| [fm ]

0.3

2

0.2

0.1

0
0

10

5

15

r [fm]
Figure D.6: The square of the single-neutron overlap functions in the Woods-Saxon (solid)
and the DOM (dashed) model are compared for 49 Ca.
the potentials, we find a reduction (enhancement) in the strength of the neutron (proton)
potential in the DOM as compared to the CH89. While the radius parameters are the same
for the real part of the neutron and proton potentials in CH89, this is not the case in the
DOM model. In this model, the real proton potential has a larger radius than that for the
neutron. This is due to the difference between the neutron and proton surface absorption.
As indicated in [134], for Sn DOM potentials the proton imaginary surface term has a linear
asymmetry dependence while this dependency is insignificant for the neutron.

D.4.2

Transfer cross section

For each reaction we perform three calculations using the FR-ADWA framework but different
inputs: (i) the nucleon optical potentials from CH89 and the standard Woods-Saxon (WS)
overlap function ϕWS for the neutron bound state (CH89+WS); (ii) DOM optical potentials

170

0

V [MeV]

(a)

-20
Vn DOM
Vn CH89
Vp DOM
Vp CH89

-40

-60
0

10

5

15

r [fm]

0
(b)

W [MeV]

-5

-10
Wn DOM
Wn CH89
Wp DOM
Wp CH89

-15

-20
0

10

5

15

r [fm]
Figure D.7: DOM and CH89 optical potentials are compared for n-132 Sn and p-132 Sn at
4.7 MeV: (a) the real part and (b) the imaginary part.

171

dσ/dΩ [mb/srad]

0.8

0.6

EXP
CH89+WS
DOM+WS
DOM

0.4

0.2

0.0
0

50

100

θ [deg]

150

Figure D.8: 48 Ca(d,p)49 Ca at Ed = 2 MeV (angular distributions are normalized to the
data at backward angles).
and the WS neutron overlap function ϕW S (DOM+WS); and (iii) both the optical potentials
and the neutron overlap function corrected for nonlocality ϕDOM from DOM (DOM).
In Fig. D.8, D.9, and D.10 we plot the angular distributions for the 48 Ca(d,p)49 Ca reaction at three different deuteron beam energies Ed = 2, 19.3 and 56 MeV respectively. The
calculated cross sections are normalized to data at the first peak of the angular distribution
for energies above the Coulomb threshold and at the backward angles otherwise. It is important to note that the DOM is able to reproduce well the shape of the angular distribution
for these cases. Although there are differences in the shape of the neutron overlap functions,
we do not obtain significant changes in the shape of the angular distribution between the
three calculations: CH89+WS, DOM+WS and DOM. We observe the same behavior for the
40 Ca(d,p)41 Ca

reaction in which DOM performs as well as CH89+WS.

Fig. D.11 and Fig. D.12 display the results for the 132 Sn(d,p)133 Sn and 208 Pb(d,p)209 Pb

172

dσ/dΩ [mb/srad]

40

30
EXP
CH89+WS
DOM+WS
DOM

20

10

0
0

20

40

60

θ [deg]

80

Figure D.9: 48 Ca(d,p)49 Ca at Ed = 19.3 MeV (angular distributions are normalized to the
data at the peak).

dσ/dΩ [mb/srad]

8

6
EXP
CH89+WS
DOM+WS
DOM

4

2

0
0

20

θ [deg]

40

60

Figure D.10: 48 Ca(d,p)49 Ca at Ed = 56 MeV (angular distributions are normalized to the
data at forward angles).

173

dσ/dΩ [mb/srad]

15

10

EXP
CH89+WS
DOM+WS
DOM

5

0
0

20

40

60

θ [deg]

80

100

Figure D.11: 132 Sn(d,p)133 Sn at Ed = 9.46 MeV (angular distributions are normalized at
the peak of the experimental cross section).

dσ/dΩ [mb/srad]

6
EXP
CH89+WS
DOM+WS
DOM

4

2

0
0

50

100

θ [deg]

150

Figure D.12: 208 Pb(d,p)209 Pb at Ed = 20 MeV (angular distributions are normalized at the
peak of the experimental data).

174

reactions. The calculated cross sections are also normalized to data in each graph. We observe a considerable change in shape of the angular distribution between the DOM and
CH89 for these two cases. The diffraction pattern is shifted toward smaller angles when
the interactions from the DOM are used. Since DOM+WS and DOM provide the same
angular distribution after normalization to data, the modification in the shape of the angular distribution is not caused by the neutron overlap function but rather by the optical
potentials. The real part of DOM potentials have a larger radius than that of the CH89,
causing the angular distribution to peak at more forward angles. The larger radius of the
DOM neutron overlap function only changes the magnitude of the calculated cross section.
The DOM method improves the angular distribution for the 132 Sn(d,p)133 Sn reaction but
deviates further from data for 208 Pb(d,p)209 Pb case.

D.4.3

Spectroscopic factors

In this section we study the impact of using the DOM potentials in extracting spectroscopic
factors, as in Eq. (D.7). In Table D.4 we present the spectroscopic factors extracted from
CH89+WS, DOM+WS and DOM calculations.
The conventional approach CH89+WS which employs the CH89 optical potentials and
the standard WS for the neutron overlap function, provides a wide range of spectroscopic
factors depending on beam energy. One can see that the spectroscopic factor varies from
0.77 to 1.1 within the energy range of 2 − 56 MeV for the 48 Ca(d,p)49 Ca reaction. A weaker
energy dependence of the spectroscopic factor is obtained for the DOM calculations with
the exception of the 208 Pb(d,p)209 Pb reaction. The comparison between CH89+WS and
DOM+WS shows a significant reduction of the spectroscopic factor. This is due to the
optical potentials alone because the neutron overlap functions are the same (WS) in the
175

two calculations. Using the overlap function from the DOM instead of the standard WS
can further reduce the spectroscopic factors for all cases except the 40 Ca(d,p)41 Ca reaction.
This is expected: while the DOM and WS overlap functions are the same for 41 Ca, this is
not the case for 49 Ca, 133 Sn and 209 Pb. There, the nonlocality in the DOM causes a shift
of the overlap functions toward the surface region which leads to larger cross sections and
smaller spectroscopic factors. The results using DOM inputs are in better agreement with
those extracted from knockout measurements (e, e p) [141] except for the Pb cases.
Through out this work, the 208 Pb(d,p)209 Pb reaction is proved to be the bad represenTable D.4: Comparision of the spectroscopic factors obtained from the FR-ADWA analysis
using different models for optical potentials and overlap functions.
Nucleus
41 Ca

49 Ca
133 Sn
209 Pb

Ed data CH89+WS
20 [136]
0.96
56 [137]
0.88
2
[124]
0.94
0.82
13 [138]
0.77
19.3 [139]
56 [129]
1.1
9.46 [140]
1.1
8
[125]
1.7
20 [127]
0.89

DOM+WS
0.85
0.73
0.72
0.67
0.68
0.70
1.0
1.5
0.61

DOM
0.86
0.74
0.66
0.61
0.62
0.62
0.72
1.2
0.51

tative case for the DOM study. We obtain a large discrepancy in extracting spectroscopic
factors at Ed = 8 MeV and Ed = 20 MeV. The same behavior is seen for the extracted
ANCs in this case. In addition, the spectroscopic factor for the sub-Coulomb reaction is
much larger than unity although the shape of the angular distribution is reproduced (see
Fig. D.13). These two problems appear not only in the DOM calculation but also in the
traditional CH89+WS approach. In order to investigate the possible cause for these issues, we look into the reaction mechanism. Coupled-channel Born approximation (CCBA)
176

dσ/dΩ [mb/srad]

0.8
0.6

EXP
CH89+WS
DOM+WS
DOM

0.4
0.2
0.0
0

50

100

θ [degree]

150

Figure D.13: 208 Pb(d,p)209 Pb at Ed = 8 MeV (angular distributions are normalized at the
peak of the experimental data).
calculations including the low-lying 3− and 2+ states in 208 Pb are performed using the
deuteron global optical potentials [142] to explore if this problem comes from target excitation. Although this CCBA approach is not meant to extract the spectroscopic factor, a
better understanding of the reaction mechanism can be obtained. The results show a strong
effect of target excitation for 208 Pb(d,p)209 Pb at 20 MeV. The effect is rather weak at the low
energy of 8 MeV. We can conclude that at least for high deuteron energies, the mechanism of
the 208 Pb(d,p)209 Pb reaction can not be described by the present ADWA framework. The
unrealistically large spectroscopic factor at Ed = 8 MeV could indicate the failure of the
present ADWA method at this low energy. It may also question the optical potentials or the
overlap functions in general. Unfortunately we can not verify that by comparing with the
exact Faddeev calculation because the reaction involves a very strong Coulomb field. The
exact Faddeev calculation is presently limited for nuclei with Z ≤ 20 [143].

177

Table D.5: ANCs extracted for both the ground state and the first excited state of the
14 C(d,p)15 C reaction.

State

Method

C2 (fm−1 )
lj
Experiment

Uncertainty
Potentials Peripherality

Total

s1/2

FR-ADWA

1.64 ± 0.26

10%

7%

11%

16%

d5/2

FR-ADWA

(3.55 ± 0.43) 10−3

10%

6%

4%

12%

D.5

Application of the finite-range adiabatic distorted
wave approximation (FR-ADWA)

As an illustration of the application of the FR-ADWA method, we extract the normalization
coefficient (ANC) from the 14 C(d,p)15 C reaction data [144]. This is our work in collaboration
with a group of physicists at Texas A&M University and experimentalists at Czech Academy
of Sciences who provide a measurement of 14 C(d,p)15 C reaction at 17.06 MeV deuteron beam
for both the ground state and the first excited state of 15 C. The single nucleon transfer reaction 14 C(d,p)15 C allows an indirect measurement of the 14 C(n,γ)15 C reaction which in some
scenarios, is very important in the nucleosynthesis of heavy nuclei. A direct measurement of
the 14 C(n,γ)15 C reaction cross section is difficult because both neutron and 14 C are unstable and radioactive. At low energies relevant for astrophysics where 14 C(d,p)15 C reaction is
peripheral, the ANC can be extracted and is used to estimate 14 C(n,γ)15 C cross section.
In this work, we employ the FR-ADWA method to extract the ANC from these new measurements. We first calculate the experimental spectroscopic factor as defined in Eq. (D.7).

178

60
EXP
FR-ADWA

dσ/dΩ [mb/srad]

50
40
30
20
10
0

10

20

30

40

θ [deg]

50

60

70

80

Figure D.14: Angular distribution for the ground state in the 14 C(d,p)15 C reaction.
Eq. (D.6) is then used to extract the nuclear ANC. In order to obtain an accurate ANC, it
is important to make sure that the reaction is peripheral. We check the peripherality of the
reaction by cutting 3 fm of the interior, leading to 11% change in the cross section for the
ground state and 4% for the first excited state of 15 C. In order to estimate the uncertainty
due to the optical potentials we perform FR-ADWA calculations using both CH89 [145] and
KD [146] optical potentials. This leads to 7% (6%) standard deviation of the ANC for the
ground state (first excited state). We also include 10% experimental systematic error in
evaluating the total uncertainty of the calculation. In Table. D.5, we present the square of
the ANC and its associate sources of uncertainty for both the ground state and the first
excited state. The total uncertainty of the square of the ANC is 16% (12%) for the ground
state (excited state) obtained by adding in quadrature all the listed uncertainties.
2
The square of the ANC C01/2 = 1.64 ± 0.26 fm−1 obtained from our analysis shows a

very good agreement with the results extracted from the Coulomb dissociation of 15 C on
179

208 Pb

[100] and also agrees with the measurement of a direct radiative capture in [147].

In Fig. D.14 we plot the 14 C(d,p)15 C angular distributions for the transition to the
ground state of 15 C. The FR-ADWA calculation when rescaled by the appropriate normalization factors reproduces well the first peak of the experimental angular distribution which
is important in extracting the ANC. There is a discrepancy between the experimental data
and the FR-ADWA calculation at large angles, but we expect this contribution is insignificant compared to the large cross section at the forward peak. Note that the cross section at
the first peak is 20 times greater then at the first minimum.

D.6

Conclusion

We explore the effect of finite-range introduced by the n-p interaction by performing a
systematic study on (d,p) reactions. The adiabatic approach for zero-range [112] and finiterange [114] calculations is used to evaluate this effect. We take into account the contributions
from both adiabatic distorting potential and the transfer matrix evaluation. Our findings
show significant reduction in the cross section due to the finite-range effect for reactions at
high energy (Ed > 30 MeV/u). We also find that the LEA method which is traditionally
applied to compensate for the finite-range effect is no longer valid for heavy systems and
high energy reactions. Our results clarify the role of finite-range effects in (d,p) reaction.
It is one of the many efforts in reducing the uncertainty in (d,p) studies. For example,
the work in [148] evaluates the uncertainties in the optical potentials, [149] considers the
effects of coupling to excited states of the target, ambiguities due to the single particle wave
functions are studied in [116, 150], and recently new ways of calculating overlap functions
are introduced by Timofeyuk and Barbieri [151, 152]. We looked at a different aspect of

180

the reaction theory, namely: the finite-range effect. We now have a better understanding of
different sources of uncertainty in our theory for transfer reactions. Applying the finite-range
adiabatic method to extract the ANC from 14 C(d,p)15 C data, we obtain consistent results
with other independent studies [100, 147].
We then use the finite-range adiabatic framework to test the performance of the dispersive
optical model potentials and the corresponding overlap functions developed by the physicists
from Washington University, St. Louis. In general, our work shows that DOM can be used
in the study of (d,p) transfer reactions. It produces the angular distributions as well as the
CH89 global optical potentials and the extracted spectroscopic factors have weaker energy
dependence. An interesting feature of this method is that the asymmetry dependence of the
potentials obtained can be used to extrapolate toward more exotic nuclei. The successful
application of this feature to the case of 132 Sn indicates that DOM can provide a promising
path to pin down optical potentials and therefore improves the uncertainty of (d,p) cross
section .

181

BIBLIOGRAPHY

182

BIBLIOGRAPHY
[1] E. E. Salpeter, Astrophys. J. 115, 326 (1952).
[2] D. D. Clayton, Principles of stellar evolution and nucleosynthesis, McGraw-Hill Book
Company, United States of America, 1968.
[3] J. Audouze and S. Vauclair, An introduction to nuclear astrophysics, D. Reidel Publishing Company, Netherlands, 1980.
[4] F. Hoyle, Astrophys. J. Suppl. Ser. 1, 121 (1954).
[5] D. N. F. Dunbar, R. E. Pixley, W. A. Wenzel, and W. Whaling, Phys. Rev. 92, 649
(1953).
[6] C. W. Cook, W. A. Fowler, C. C. Lauritsen, and T. Lauritsen, Phys. Rev. 107, 508
(1957).
[7] S. M. Austin, Nucl. Phys. A 758, 375c (2005).
[8] G. R. Caughlan and W. A. Fowler, At. Nucl. Data Tables 40, 238 (1988).
[9] C. Iliadis, Nuclear Physics of Stars, Germany, 2007.
[10] C. Angulo et al., Nucl. Phys. A 656, 3 (1999).
[11] R. de Diego, E. Garrido, D. V. Fedorov, and A. S. Jensen, Phys. Lett. B 695, 324
(2011).
[12] M. Freer et al., Phys. Rev. C 80, 041303(R) (2009).
[13] H. O. U. Fynbo et al., Nature 433, 136 (2005).
183

[14] N. B. Nguyen, F. M. Nunes, I. J. Thompson, and E. F. Brown, PoS NIC XII, being
published.
[15] K. Ogata, M. Kan, and M. Kamimura, Prog. Theor. Phys. 122, 1055 (2009).
[16] N. Austern, Y. Iseri, M. Kamimura, G. Rawitscher, and M. Yahiro, Phys. Rep. 154,
125 (1987).
[17] M. Yahiro, N. Nakano, Y. Iseri, and M. Kamimura, Prog. Theo. Phys. 67, 1464
(1982).
[18] E. Hiyama, M. Kamimura, T. Motoba, T. Yamada, and Y. Yamamoto, Prog. Theor.
Phys. 97, 881 (1997).
[19] E. Hiyama, Y. Kino, and M. Kamimura, Prog. Part. Nucl. Phys. 51, 223 (2003).
[20] V. Vasilevsky, A. V. Nesterov, F. Arickx, and J. Broeckhove, Phys. Rev. C 63, 034606
(2001).
[21] A. Dotter and B. Paxton, Astron. Astrophys. 507, 1617 (2009).
[22] T. Suda, R. Hirschi, and M. Fujimoto, Astrophys. J. 741, 61 (2011).
[23] M. Saruwatari and M. Hashinoto, Prog. Theor. Phys. 124, 925 (2010).
[24] Y. Matsuo et al., Prog. Theor. Phys. 126, 1177 (2011).
[25] F. Peng, and C. D. Ott, Astrophys. J. 725, 309 (2010).
[26] R. Alvarez-Rodriguez, E. Garrido, A. S. Jensen, D. V. Fedorov, and H. O. U. Fynbo,
Eur. Phys. J. A 31, 303 (2007).
[27] R. de Diego, E. Garrido, D. V. Fedorov, and A. S. Jensen, Eur. Phys. Lett. 90, 52001
(2010).
[28] E. Garrido, R. de Diego, D. V. Fedorov, and A. S. Jensen, Eur. Phys. J. A 47, 102
(2011).
[29] P. Maris, J. P. Vary, and A. M. Shirokov, Phys. Rev. C 79, 014308 (2009).

184

[30] S.C. Pieper, Nucl. Phys. A 751, 516 (2005).
[31] M. Chernykh, H. Feldmeier, T. Neff, P. von Neumann-Cosel, and A. Richter, Phys.
Rev. Lett. 98, 032501 (2007).
[32] L. D. Faddeev, Zh. Eksp. Teor. Fiz. 39, 1459 (1960); Sov. Phys. JETP 12, 1014
(1961).
[33] I. J. Thompson, F. M. Nunes, and B. V. Danilin, Comp. Phys. Comm. 161, 87 (2004).
[34] S. P. Merkuriev, Ann. Phys. 130, 395 (1980).
[35] A. Kievsky, M. Viviani, S. Rosati, Phys. Rev. C 56, 2987 (1997).
[36] A. Deltuva, A. C. Fonseca, P. U. Sauer, Phys. Rev. C 71, 054005 (2005).
[37] F. M. Nunes, J. A. Christley, I. J. Thompson, R. C. Johnson, and V. D. Efros, Nucl.
Phys. A 609, 43 (1996).
[38] F. M. Nunes, Core excitation in few-body systems: Application to light exotic nuclei,
PhD thesis, University of Surrey, England, 1995.
[39] T. H. Gronwall, Phys. Rev. 51, 655 (1937).
[40] C. D. Lin, Phys. Rep. 257, 1 (1995).
[41] L. M. Delves, Nucl. Phys. 9, 391 (1959).
[42] L. M. Delves, Nucl. Phys. 20, 268 (1962).
[43] M. Zhukov, B. Danilin, D. Fedorov, J. Bang, I. Thompson, and J. Vaagen, Phys. Rep.
231, 151 (1993).
[44] T. Tarutina, I. Thompson, and J. Tostevin, Nucl. Phys. A 733, 53 (2004).
[45] P. Descouvemont, C. Daniel, and D. Baye, Phys. Rev. C 67, 044309 (2003).
[46] I. J. Thompson and F. M. Nunes, Nuclear Reactions for Astrophysics, Cambridge
University Press, United Kingdom, 2009.
185

[47] M. Abramowitz and I. A. Stegun, Handbook of mathematical functions, Dover Publications, 1972.
[48] J. Raynal and J. Revai, Nuovo Cimento A 68, 612 (1970).
[49] S. P. Merkuriev, C. Gignoux and A. Laverne, Ann. Phys. (NY) 99, 1 (1976).
[50] S. P. Merkuriev, Sov. J. Nucl. Phys. 19, 447 (1974).
[51] D. M. Brink and G. R. Satchler, Angular Momentum, Oxford Science Publications,
1994.
[52] B. V. Danilin and M. V. Zhukov, Phys. At. Nucl. 56, 460 (1993).
[53] E. P. Wigner, L. Eisenbud, Phys. Rev. 72, 29 (1947).
[54] A. M. Lane, R. G. Thomas, Rev. Mod. Phys. 30, 257 (1958).
[55] I. J. Thompson, B. V. Danilin, V. D. Efros, J. S. Vaagen, J. M. Bang and M. V.
Zhukov, Phys. Rev. C 61, 24318 (2000).
[56] P. Descouvemont and D. Baye, Rep. Prog. Phys. 73, 036301 (2010).
[57] J. C. Light, R. B. Walker, J. Chem. Phys. 65, 4272 (1976).
[58] E. B. Stechel, R. B. Walker, J. C. Light, J. Chem. Phys. 69, 3518 (1978).
[59] B. I. Schneider, R. B. Walker, J. Chem. Phys. 70, 2466 (1979).
[60] I. J. Thompson, http://www.ianthompson.org/computation.htm, version 92.
[61] G. Erens, J. L. Visschers, and R. van Wageningen, Ann. Phys. 67, 461 (1971).
[62] W. A. Fowler, G. E. Caughlin, and B. A. Zimmerman, Annu. Rev. Astron. Astrophys.
5, 525 (1967).
[63] C. Forss´n, N. B. Shul’gina, and M. V. Zhukov, Phys. Rev. C 67, 045801 (2003).
e
[64] S. Ali and A. R. Bodmer, Nucl. Phys. 80, 99 (1966).
186

[65] D. V. Fedorov, A. S. Jensen, Phys. Lett. B 389, 631 (1996).
[66] B. Buck, H. Friedrich and C. Wheatley, Nucl. Phys. A 275, 246 (1977).
[67] J. L. Friar, B. F. Gibson, and G. L. Payne, Phys. Rev. C 37, 2869 (1988).
[68] T. Y. Saito and I. R. Afnan, Phys. Rev. C 50, 2756 (1994).
[69] F. Ajzenberg-Selove and J. H. Kelley, Nucl. Phys. A 506, 1 (1990).
[70] G. S. Anagnostatos, Phys. Rev. C 51, 152 (1995).
[71] T. Neff and H. Feldmeier, Nucl. Phys. A 738, 357 (2004).
[72] E. Epelbaum, H. Krebs, D. Lee, and Ulf-G. Meissner, Phys. Rev. Lett. 106, 192501
(2011).
[73] V. Vasilevsky, F. Arickx, W. Vanroose, and J. Broeckhove, Phys. Rev. C 85, 034318
(2012).
[74] http://mathworld.wolfram.com/CoulombWaveFunction.html.
[75] P. Descouvemont, J. Phys. G 37, 064010 (2010).
[76] M. Salaris, S. Cassisi, and A. Weiss, PASP, 114, 375 (2002).
[77] B. Paxton, L. Bildsten, A. Dotter, F. Herwig, P. Lesaffre, and F. Timmes, Astrophys.
J. Suppl. Ser. 192, 3 (2011).
[78] http://mesa.sourceforge.net/ (version 4442).
[79] K. J. Shen and L. Bildsten, Astrophys. J. 699, 1365 (2009).
[80] K. Nomoto, Astrophys. J. 253, 798 (1982).
[81] H. Graboske et al., Astrophys. J. 181, 457 (1973).
[82] L. Bildsten, K. J. Shen, Astrophys. J. 662, L95 (2007).

187

[83] D. Poznanski et al., Science, 327, 58 (2010)
[84] M. M. Kasliwal et al., Astrophys. J. Lett, 723, L98 (2010)
[85] C. J. Hansen and H. M. van Horn, Astrophys. J. 195, 735 (1975).
[86] L. Keek, J. J. M. In’t Zand, PoS INTEGRAL08, 032 (2008).
[87] L. Keek, private communication, 2012.
[88] E. Kuulkers et al., A&A, 514, A65 (2010).
[89] E. Epelbaum, H. Krebs, T. Lahde, D. Lee, and Ulf-G. Meissner, arxiv:1208.1328v2,
2012.
[90] A. Spyrou et al., Phys. Rev. Lett 108, 102501 (2012).
[91] I. A. Egorova et al., Phys. Rev. Lett 109, 202502 (2012).
[92] Md. A. Khan, Few-body Syst. 52, 53 (2012).
[93] H. Niewodniczanski, J. Nurzynski, A. Strzalkowski and G. R. Satchler, Phys. Rev.
146, 799 (1966).
[94] W. W. Daehnick and Y.S. Park, Phys. Rev. 180, 1062 (1969).
[95] D. E. Bainum, R. W. Finlay, J. Rapaport, J. D. Carlson and W. G. Love, Phys. Rev.
C 16, 1377 (1977).
[96] J. P. Schiffer et al., Phys. Rev. Lett 92, 162501 (2004).
[97] K.L. Jones et al., Phys. Rev. C 70, 067602 (2004).
[98] J. S. Thomas et al., Phys. Rev. C 71, 021302(R) (2005).
[99] J. S. Thomas et al., Phys. Rev. C 76, 044302 (2007).
[100] N. C. Summers and F. M. Nunes, Phys. Rev. C 78, 011601(R) (2008).

188

[101] R. L. Kozub et al., Phys. Rev. Lett. 109, 172501 (2012).
[102] D. Y. Pang, F. M. Nunes, and A. M. Mukhamedzhanov, Phys. Rev. C 75, 024601
(2007).
[103] A. M. Mukhamedzhanov, F. M. Nunes, and P. Mohr, Phys. Rev. C 77, 051601(R)
(2008).
[104] A. H. Wuosmaa et al., Phys. Lett. 94, 082502 (2005).
[105] A. H. Wuosmaa et al., Phys. Rev. C 72, 061301 (2005).
[106] H. B. Jeppesen et al., Phys. Lett. B 642, 449 (2006).
[107] W. N. Catford et al., J. Phys. G 31, S1655 (2005).
[108] N. Austern, Direct Nuclear Reaction Theories, John Wiley & Sons, New York, 1970.
[109] G. R. Satchler, Direct Nuclear Reactions, Oxford University Press, Oxford, 1983.
[110] N. K. Glendenning, Direct Nuclear Reactions, World Scientific Publishing, 2004.
[111] P. G. Roos et al., Nucl. Phys. A 225, 187 (1975).
[112] R. C. Johnson and P. J. R. Soper, Phys. Rev. C 1, 976 (1970).
[113] R. C. Johnson, AIP Conf. Conf. Procs. 79, 128 (2005).
[114] R. C. Johnson and P. C. Tandy, Nucl. Phys. A 235, 56 (1974).
[115] P. J. A. Buttle and L. J. B. Goldfarb, Proc. Phys. Soc. 83, 701 (1964).
[116] A. M. Mukhamedzhanov and F. M. Nunes, Phys. Rev. C 72, 017602 (2005).
[117] A. Laid, J. A. Tostevin, and R. C. Johnson, Phys. Rev. C 48, 1307 (1993).
[118] N. B. Nguyen, F. M. Nunes, and R. C. Johnson, Phys. Rev. C 82, 014611 (2010).
[119] V. Reid, Ann. Phys. (N.Y.) 50, 411 (1968).
189

[120] R. L. Varner et al., Phys. Rep. 201, 57 (1991).
[121] M. Igarashi et al., computer program TWOFNR, University of Surrey version, 2008.
[122] I. J. Thompson, Comput. Phys. Rep. 7, 167 (1988).
[123] G. L. Wales and R. C. Johnson, Nucl. Phys. A 274, 168 (1976).
[124] J. Rapaport, A. Sperduto, and M. Salomaa, Nucl. Phys. A 197, 337 (1972).
[125] G. M. Crawley, B. V. N. Rao, and D. L. Powell, Nucl.Phys. A 112, 223(1968).
[126] K. Haravu, C. L. Hollas, P. J. Riley, and W. R. Coker, Phys. Rev.C 1, 938 (1970).
[127] D. G. Kovar, N. Stein, and C. K. Bockelman, Nucl. Phys. A 231, 266 (1974).
[128] K. Hatanaka et al., Nucl. Phys. A 419, 530 (1984).
[129] Y. Uozumi et al., Nucl. Phys. A 576, 123 (1994).
[130] J. Liu, N.B. Nguyen, and F.M. Nunes, internal report, NSCL/MSU August 2012.
[131] F. M. Nunes, J. Liu, N. B. Nguyen, L. Titus, and N. J. Upadhyay, Proceedings of 13th
International Conference on Nuclear Reaction Mechanism, Varenna, June (2012)
[132] R. J. Charity, L. G. Sobotka, and W. H. Dickhoff, Phys. Rev. Lett. 97, 162503 (2006).
[133] R. J. Charity, J. M. Mueller, L. G. Sobotka, and W. H. Dickhoff, Phys. Rev. C 76,
044314 (2007).
[134] J. M. Mueller et al., Phys. Rev. C 83, 064605 (2011).
[135] N. B. Nguyen, S. J. Waldecker, F. M. Nunes, R. J. Charity, and W. H. Dickhoff,
Phys. Rev. C 84, 044611 (2011).
[136] F. J. Eckle el. al., Nucl. Phys. A 506, 159 (1990).
[137] K. Hatanaka et. al., Nucl. Phys. A 419, 530 (1984).

190

[138] W. D. Metz et. al, Phys. Rev. C 12, 827 (1975).
[139] W. D. Metz et. al, Phys. Rev. C 12, 827 (1975).
[140] K. Jones et al., Nature 465, 454 (2010).
[141] L. Lapik´s, Nucl. Phys. A 553, 297c (1993).
a
[142] J.M. Lohr and W. Haeberli, Nucl. Phys. A 232, 381 (1974).
[143] N. J. Upadhyay, A. Deltuva, and F. M. Nunes, Phys. Rev. C 85, 054621 (2012).
[144] A. M. Mukhamedzhanov, V. Burjan, M. Gulino, Z. Hons, V. Kroha, M. McCleskey,
J. Mrazek, N. Nguyen, F. M. Nunes, S. Piskor, S. Romano, M. L. Sergi, C. Spitaleri,
and R. E. Tribble, Phys. Rev. C 84, 024616 (2011).
[145] R. L. Varner et al., Phys. Rep. 201, 57 (1991).
[146] A. J. Koning, J. P. Delaroche, Nucl. Phys. A 713, 231(2003).
[147] R. Reifarth et al., Phys. Rev. C 77, 015804 (2008).
[148] X. D. Liu, M. A. Famiano, W. G. Lynch, M. B. Tsang, and J. A. Tostevin, Phys.
Rev. C 69, 064313 (2004).
[149] F. Delaunay, F. M. Nunes, W. G. Lynch, and M. B. Tsang, Phys. Rev. C 72, 014610
(2005).
[150] D. Y. Pang, F. M. Nunes, and A. M. Mukhamedzhanov, Phys. Rev. C 75, 024601
(2007).
[151] N. K. Timofeyuk, Phys. Rev. Lett. 103, 242501 (2009).
[152] C. Barbieri, Phys. Rev. Lett. 103, 202502 (2009).

191