Imaginary-time theory for triple-alpha reaction rate T. Akahori,1 Y. Funaki,2 and K. Yabana3,1 1 Institute of Physics, University of Tsukuba, Tsukuba 305-8571, Japan 2 Nishina Center for Accelerator-Based Science, The Institute of Physical and Chemical Research (RIKEN), Wako 351-0198, Japan 3 Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8571, Japan Usingimaginary-timetheory,itisshownthatthetriple-alphareactionratecanbereliablycalcu- latedwithouttheneedtosolvescatteringproblemsinvolvingthreechargedparticles. Thecalculated reactionrateisfoundtoagreewellwiththeempiricalNACRErate,whichiswidelyadoptedinstellar evolutioncalculations. ThereasonforthisisexplainedusingR-matrixtheory. Extremelyslowcon- vergenceis found to occur when a coupled-channelexpansion is introduced,which helps toexplain thevery different reaction rates obtained using different theoretical approaches. 4 1 0 The triple-alpha reaction is a key process that influ- is both surprising and puzzling. However, two possible 2 ences the production of all heavy elements in the uni- explanations can be put forward. The first is the lack of n verse. Accurate knowledge of the reaction rate is essen- a rigorous scattering theory for three charged particles. a tialforunderstandingstellarevolutionandnucleosynthe- Thesecondisrelatedtothequantum-tunnelingnatureof J sis. Sinceexperimentalmeasurementsarenotfeasiblefor the process: the α-particle travels through the Coulomb 7 this reaction, theoretical evaluation of the reaction rate barrieroveralongdistance,typicallyafewhundredfem- 1 is crucially important. tometers,causingthereactionratetobeextremelysmall, ] In the triple-alpha process, the importance of 12C and and to vary by 60 orders of magnitude within the range h 8Be resonancesis wellrecognized[1, 2]. At hightemper- of astrophysically relevant temperatures. t - ature, the reaction proceeds dominantly through a res- Recently,wehaveproposedanewtheoreticalapproach l c onant 0+ state of 12C at 7.65 MeV, which is known as fordeterminingtheradiativecapturereactionrate,which u the Hoyle state. At lower temperatures, processes that we referto as imaginary-timetheory[18]. Inthis theory, n [ do not involve the Hoyle state become important. An imaginary time is identified with inverse temperature as empiricalreactionrate assumingsuccessivetwo-bodyre- isoftenusedinquantummany-bodytheoriesofnonequi- 1 actions of α-α and α-8Be has been derived [3–5], and is libriumsystems. Arelatedapproachhas been developed v adopted in the NACRE compilation [6] as the standard for the theory ofchemicalreactionrates[19]. Inthis let- 0 9 rate to be used in stellar evolution calculations. ter,wereporttheapplicationofimaginary-timetheoryto 3 However, the validity of the empirical rate formula thedeterminationofthetriple-alphareactionrate. Since 4 shouldbeconfirmedbycalculationsbasedonmicroscopic the theory does not require any scattering problems to . 1 quantumtheories. Severaltheoreticalattempts to calcu- be solved, it is ideally suited to the triple-alpha process, 0 late the rate using quantum theory involving three α- for which there is no formal scattering theory available. 4 particles have recently been undertaken. The first was Indeed,it willbe shownthatthe reactionratecanbe re- 1 conducted by Ogata and coworkers [7], and employed liably calculated without any numerical problems. The : v continuum-discretized coupled-channel (CDCC) theory, calculated rate is found to be virtually identical to the i X whichisawell-establishedtheoryfordirectnuclearreac- empirical NACRE rate, and no enhancement occurs at tions [8]. A surprisingly high value for the reaction rate low temperatures. The reason for this good agreement r a was found at low temperatures, and at T = 0.01 GK it is investigated analytically by combining R-matrix the- waslargerthantheNACREvalueby26ordersofmagni- ory [20] with imaginary-time theory. tude[6]. Soonafterthis reportwaspublished,theconse- The following expressiondescribes the triple-alpha re- quences ofthe new rate for the presentunderstanding of action rate [18], stellar evolutionwereinvestigated[9–11]. It waspointed out that such a high rate would not be compatible with 2πβ~2 3 8π(λ+1) N2hαααi=6·33/2N2 thestandardpictureofstellarevolution;forexample,the A A(cid:18) M (cid:19) ~λ((2λ+1)!!)2 α red giant phase disappears if the rate is adopted [9, 11]. 2λ+1 H −E Following the report by Ogata et al., calculations using × hΦ |M e−βH f PM† |Φ i (1) f λµ (cid:18) ~c (cid:19) λµ f different quantumthree-body approacheshave been car- MXfµ riedout[12–17]. Unfortunately,thereisalargedegreeof scatter in the reported rates at low temperatures, which where β = 1/k T is the inverse temperature, M is the B α vary between the NACRE rate [6] and that determined mass of an α-particle, H is the Hamiltonian for three by the CDCC calculation [7]. In view of the number of α particles, Φ is the wave function for 12C in the final f successful achievements of nuclear three-body reaction state after γ-ray emission, E is the energy of the final f theories,this hugediscrepancyamongthe reportedrates state measured from the three α threshold, and M is f 2 the magnetic quantum number for the final state. M T [GK] λµ is the multipole transition operator for γ-ray emission 1.0 0.1 0.05 0.03 0.02 0.015 0.01 1 with a multipolarity λ. The reactionof three α-particles withtotalangularmomentumJ =0isconsidered,which 10−10 rmax=Rmax=50 fm lsetaadtesotof 1t2hCe aemt 4is.s4i4onMoefVafoγr-rtahyewfiintahlλwa=ve2,fuanncdtiothneΦ2f+. −2mol ] 1100−−2300 rrrmmmaaaxxx===RRRmmmaaaxxx===124000000 fffmmm P is the projectionoperatorwhicheliminates anybound −1s rmax=Rmax=600 fm eigenstates of the three-body Hamiltonian. It can easily 6m 10−40 NACRE be shownthat Eq. (1) exactly coincides with the expres- [ c> 10−50 sion for the triple-alpha reaction rate, Eqs. (2) ∼ (5) of α α [21],byinsertingacompletenessrelationforafunctional α< 10−60 space of three α-particles into Eq. (1). 10−70 Inpractice,toevaluatethereactionrateusingEq.(1), it is necessary to calculate Ψ(β) = e−βHM† Φ . This 10−80 λµ f 100 200 300 400 500 600 700 800 900 10001100 is achieved by evolving the wave function along the β [1/MeV] imaginary-time axis, FIG. 1: (Color online) Calculated triple-alpha reaction rates ∂ fordifferentchoicesofradialcutoffdistances,rmax andRmax. − Ψ(β)=HΨ(β), (2) ∂β The NACRErate is also shown for comparison. starting with the initial wave function, Ψ(β = 0) = M놵Φf. The reaction rate at an inverse temperature V3α(r1,r2,r3)=V3e−µ(r122+r223+r321) with µ=0.15 fm−2. β is then evaluated using the wave function at β/2, The value of V is chosen so that the resonance energy 3 of the Hoyle state, the 0+ state of 12C, is reproduced at 2 5 β H −E β 379.8 keV above the three α threshold. The final wave f hαααi∝ Ψ Ψ . (3) MXfµ(cid:28) (cid:18)2(cid:19)(cid:12)(cid:12)(cid:12)(cid:18) ~c (cid:19) (cid:12)(cid:12)(cid:12) (cid:18)2(cid:19)(cid:29) fuusnincgtiothneΦofrtfhoorg1o2nCalJitπy=co2n+ditaiton4.4m4oMdeelV[2i4s].constructed (cid:12) (cid:12) The numerical calculations are carried out using the In Fig. 1, the calculated triple-alpha reaction rates model space and the three-body Hamiltonian described for different spatial areas specified by Rmax and rmax below. In the calculations, the α-particles are treated are compared. The NACRE rate [6] is also shown for as point particles. The assumption of dominant J = 0 comparison. It can be seen that when Rmax and rmax contribution is expected to be valid below T < 1.0 are larger than 400 fm, a fully converged reaction rate GK [21, 22]. A Jacobi coordinate system is used, de- is obtained in the entire temperature region, and this fined by r =r −r andR=(r +r )/2−r ,where r rate coincides well with the NACRE rate. Calculations 1 2 1 2 3 i (i = 1−3) are the coordinates of the three α-particles. within smaller spatial areas yield a rate which is valid The three-body wave function is expanded in partial only in limited higher-temperature regions. It has been waves, Ψ(β) = (u (R,r,β)/Rr)[Y (Rˆ)Y (rˆ)] . In previously shown that for a two-body radiative capture L L L L J=0 the present worPk, only the L = 0 component is consid- process, the maximum radius required for full conver- ered, since this is expected to be the most important gence roughly corresponds to the exit point of quan- at low temperature. Ogata et al. adopted the same tum tunneling at the Gamow peak energy [18]. Us- model space for their CDCC calculations [7]. The ra- ing Rmax = rmax = 400 fm at T = 0.01 GK, the dial variables R and r are discretized with a grid size Gamow peak energy for the triple-alpha reaction can be ∆R=∆r =0.5 fm, and radial grid points are employed estimated from the Coulomb potential, and is given by up to the maximum values, Rmax and rmax. The differ- E0 ≃3×4e2/400fm=43 keV. ential operators in the Hamiltonian are treated using a AlthoughthecalculatedreactionrateandtheNACRE nine-pointfinite differenceformula. TosolveEq.(2),the rateshowgoodagreementonthelogarithmicscaleinFig. Taylorexpansionmethodisusedforshort-timeevolution 1, the calculated rate is smaller than the NACRE rate with a step size of ∆β =0.004 MeV−1. over the entire temperature region by a factor of up to The Hamiltonian H for the three α-particles is con- sixtimes. Thereasonforthisdifferencewillbediscussed structed as follows. For the potential between two α- later. particles, the Ali-Bodmer potential is used, considering IntheempiricalNACREformula,therearethreetem- only the l = 0 angular momentum channel [23]. The peratureregionsthat aredistinguished by different reac- potential parameter is modified slightly so that it ac- tion mechanisms [25]: T > 0.074 GK dominated by the curately reproduces the resonance corresponding to the Hoyle state process, 0.074 GK > T > 0.028 GK domi- ground state of 8Be at 92.08 keV. A three-body poten- nated by the α-8Be two-body nonresonant process, and tialamongthethreeα-particlesisadded,andisgivenby T < 0.028 GK dominated by a nonresonant process in- 3 T [GK] formula [3, 6] is justified. 1.0 0.1 0.05 0.03 0.02 0.015 0.01 In the empiricalformula, it is assumedthat successive β=100 [1/MeV] <H> 0.6 60 α-α and α-8Be reactions occur, and the reaction cross m] 40 sections are described using Breit-Wigner formulas. It 0.5 R [f 20 β=250 [1/MeV] willbeshownthatitispossibletoderiveaformulaquite 300 similar to the empirical one starting from Eq. (1), by 0 MeV] 0.4 0 20r [fm ]40 60 R [fm] 120000 600 β=700 [1/MeV] aasnsdumthienngatphpatrotxhiemtahtirnege-ibtoudsyinHgaRm-imltaotnriiaxntihseoserpya[r2a0b].le, E [ 0.3 0 m] 400 Theseparabilityassumptionforthethree-bodyHamil- 0 10r0 [fm 2]00 300R [f 200 tonian is written as 0.2 0 0 200 400 600 0.1 r [fm] H =hαα(r)+hαBe(R), (4) 0 where the α-α Hamiltonian, h (r) = T +V (r), has 0 100 200 300 400 500 600 700 800 900 1000 1100 αα r αα β [1/MeV] a resonance at Eαα = 92.08 keV. The normalized wave r function for the resonance is expressed as φαα(r). In FIG.2: (Coloronline)Energyexpectationvalueandvariance r addition, a simple potential model is assumed for the α- asfunctionoftemperature. Theinsetsshowdensitydistribu- tions for threedifferent temperatures. 8Be relative motion, so that hαBe(R) = TR +VαBe(R). ThepotentialV (R)ischosensoastogivearesonance αBe at EαBe =287.7 keV with the normalized wave function r volvingthreeα-particles. AcarefullookatFig.1reveals of φαrBe(R). The Hoyle state is then described by the that the calculated rate curves show changes in slope at wave function product, φαrα(r)φαrBe(R), at the summed exactly the same temperatures. resonance energy of Erαα+ErαBe =379.8 keV. To illustrate this more clearly, Figure 2 shows the The following approximation is then introduced. The energy expectation value and the variance, defined by problemofpotentialscatteringineithertheα-αorα-8Be E¯ =hΨ(β/2)|H|Ψ(β/2)i/hΨ(β/2)|Ψ(β/2)iand(∆E)2 = system is considered. The resonance energy is denoted hΨ(β/2)|(H −E¯)2|Ψ(β/2)i/hΨ(β/2)|Ψ(β/2)i, as a func- as Er and its normalized radial wave function as ur(r). tion of the inverse temperature, respectively. The Using R-matrix theory [20], the radial wave function at insets show the density distribution, ρ(R,r,β/2) = anenergyE aroundaresonancewithanasymptoticform uL=0(R,r,β/2)2/R2r2, for three typical temperatures. ofuE(r)→(2µ/π~2k)1/2sin(kr+δ)canbeapproximated by In the high temperature region T > 0.074 GK, the energy expectation value coincides with the resonance energy of the Hoyle state, E¯ =379.8 keV, indicating the uE(r)=ur(r) L(E,Er,Γr(E)), (5) dominance of the Hoyle state process. In addition, the p density ρ(R,r,β/2) is localized within a small R and r where L(E,E ,Γ ) is a Lorentzian function given by r r region, which is consistent with the resonant picture. In the medium temperature region, 0.074 GK > T > 0.028 1 Γ (E) r L(E,E ,Γ (E))= . (6) GK, most of the density is contained within a small r r r 2π(E−E )2+Γ (E)2/4 r r region, r < 10 fm, whereas it is extended along the R direction. This indicates that two of the α-particles are Here the shift of the resonance energy is ignored. The forming a 8Be resonance, with the third remaining out- energy-dependentwidth Γ (E) is relatedto the widthat r side. In the lowest temperature region, T < 0.028 GK, the resonance energy Γ by Γ (E) = Γ P(E)/P (E ), r r r l l r the density extends in both the R and r directions, in- where P (E) is the penetrability. l dicating the nonresonantcharacterof the reaction. Note Using Eq. (5), any function f(H) of the three-body that the average energy at T =0.01 GK is E¯ = 55 keV, separableHamiltonianH inEq.(4)canbeapproximated which is in reasonable agreement with the Gamow peak as energy of 43 keV estimated from the radial convergence shown in Fig. 1. f(H)=|φαα(r)ihφαα(r)|·|φαBe(R)ihφαBe(R)|· The agreement between the calculated and NACRE r r r r rates, not only in terms of the magnitude, but also × dE dE L(E ,Eαα,Γαα(E )) Z ααZ αBe αα r r αα with regard to the change in reaction mechanism, in- dicates that imaginary-time theory provides quantum- ×L(EαBe,ErαBe,ΓαrBe(EαBe))f(Eαα+EαBe), (7) mechanicalsupport for the conventionaldescription. An analytic investigation based on microscopic three-body where Γαα and ΓαBe are the α-decay widths of the α- r r theorywasnextcarriedouttodetermine ifthe empirical α and α-8Be resonances, respectively. Substituting this 4 into Eq. (1) gives T [GK] 1.0 0.1 0.05 0.03 0.02 0.015 0.01 1 2πβ~2 3 14 ch. hαααi=6·33/2 10−10 20 ch. (cid:18) Mα (cid:19) 100 ch. ×Z dEααZ dEαBeL(Eαα,Erαα,Γαrα(Eαα)) −2mol ] 1100−−2300 11402000000 ccchhh... ×L(EαBe,ErαBe,ΓαrBe(EαBe))e−βEαα−βEαBe 6−1m s 10−40 No C.C. ×Γ (12C;0+) Eαα+EαBe−E(12C;2+) 5, (8) [ c> 10−50 γ 2 (cid:18)Er(12C;0+2)−E(12C;2+)(cid:19) ααα< 10−60 where Γγ(12C;0+2) is the radiative decay width of the 10−70 HoylestateandE(12C;2+)istheexcitationenergyofthe first2+ stateof12C.TherateexpressiongiveninEq.(8) 10−80 100 200 300 400 500 600 700 800 900 10001100 is almost equivalent to the empirical NACRE formula β [1/MeV] [6]. Thus, a formula quite close to the NACRE formula couldbesuccessfullyderived,startingwithamicroscopic FIG. 3: (Color online) Calculated triple-alpha reaction rate three-body Hamiltonian. using coupled-channel expansion for different numbers of However, a question still remains regarding the valid- channels, nmax. ity of assuming that the three-body Hamiltonian is in fact separable. To resolve this, a numericalinvestigation was carried out to determine how much this assumption states associated with low eigenvalues are characterized changes the calculated reaction rate. It was found that by a large α-α separation outside the Coulomb barrier, using the separable Hamiltonian (4) when solving equa- except for the resonant state corresponding to the 8Be tion(2)changedthereactionratebyonlyafactoroftwo ground state, which appears as the 14th eigenstate. In or less. This indicates that if the Hamiltonian is con- the coupled-channel approach, the wave function is ex- structed such that the 8Be ground state and 12C Hoyle panded in the form u(R,r,β) = nvn(R,β)wn(r) and state resonancesare reasonablydescribed, the separabil- the imaginary-timeevolution of vnP(R,β) is calculated in ity approximation does not seriously affect the reaction theformofamatrixdifferentialequation. Itwasfirstnu- rate. merically confirmed that, employing all 1200 eigenstates The largest numerical difference between Eq. (8) and intheexpansion,thecalculatedrateexactlymatchesthe theNACREexpressionisassociatedwiththewidthΓαBe. result shown in Fig. 1, as is expected. However, if the r Inthederivationusedinthepresentpaper,thisquantity number of basis functions in the expansion is truncated, representsthedecaywidthoftheα-8Beresonance,while the results depend on the degree of truncation. The de- in the NACRE derivation it is the particle decay width pendence of the convergence behavior on the number of of the Hoyle state, which may include both α-α and α- basis functions is shown in Fig. 3. In each coupled- 8Bedecay. As previouslystated,thereactionrateshown channel calculation, the strength of the three-body po- in Fig. 1 is smaller than the NACRE rate by a factor of tential among the three α-particles is adjusted so that up to six. It is considered that this is largely due to the the Hoyle state always appears at 379.8 keV. difference in ΓαBe. Employingstates below the 8Be groundstate (n = r max Finally, the reason why different theoretical ap- 14), the reaction rate is much higher than that given by proaches yield such widely different reaction rates at the fully converged calculation, and the difference is 24 low temperatures is considered. In particular, although orders of magnitude at T = 0.01 GK. This can be eas- the CDCC approach[7] employs almost the same model ily understood because low-energy eigenstates of h (r) αα spaceasthatadoptedinthepresentpaper,theresultdif- are characterized by a large separation between two α- fers by 26 order of magnitude at low temperature from particlesoutsidetheCoulombbarrier. TheCoulombbar- both the NACRE rate and the rate obtained here. To rier for the third α-particle is then very small in these clarify the originof this discrepancy, the imaginary-time channels, yielding a high reaction rate. This artificial evolution of Eq. (2) was investigated using the coupled- enhancement of the calculated reaction rate for a small channel method. channelnumberisreducedbyoff-diagonalcouplingterms In the coupled-channelapproach,the eigenvalue prob- asthe number ofchannelsincreases. However,asseenin lemforα-αrelativemotiondescribedbytheHamiltonian Fig. 3, convergence is extremely slow. Even with 400 h (r)isfirstsolved. Discretizingtheradialvariablerin channels,thereactionrateisstill13ordersofmagnitude αα 0.5-fm steps up to 600 fm gives 1200 grid points for this larger than the fully convergedresult. Since the α-α en- coordinate. Diagonalizing the radial Hamiltonian then ergyforthe400thchannelisaround46MeV,whichisfar gives 1200 eigenfunctions, w (r) (n = 1−1200). Eigen- above the Coulomb barrier, the coupling effect relevant n 5 to the slow convergence is not related to physicaldistor- [6] C. Anguloet.al, Nucl. Phys. A656, 3 (1999). tion. It is considered that this indicates a difficulty in [7] K. Ogata, M. Kan, M. Kamimura, Prog. Theor. Phys. numerically expressing exponentially small functions in 122, 1055 (2009). [8] M. Kamimura, M. Yahiro, Y. Iseri, Y. Sakuragi, H. the tunneling process using the basis expansion method. Kameyama, M. Kawai, Prog. Theor. Phys. Suppl. 89, Insummary,imaginary-timetheorywasappliedtode- 1 (1986). termine the triple-alpha reaction rate. Since the the- [9] A. Dotter and B. Paxton, Astronomy and Astrophysics ory does not require solving any scattering problems, it 507, 1617 (2009). is quite suitable for the triple-alpha process. Indeed, a [10] M. Saruwatari and M. Hashimoto, Prog. Theor. Phys. converged reaction rate was obtained without any nu- 124, 925 (2010). merical problems. The calculated rate agreed well with [11] T. Suda, R. Hirschi, M.Y. Fujimoto, Astrophys. J. 741, 61 (2011). the conventional NACRE rate, not only in terms of the [12] E.Garrido,R.deDiego,D.V.Fedorov,andA.S.Jensen, magnitude, but also the critical temperatures where the Eur. Phys. J. A47, 102 (2011). dominantreactionmechanismchanges. Noenhancement [13] S. Ishikawa, Few-Body Syst.(2012). of the rate was found at low temperature. The reason [14] N.B. Nguyen, F.M. Nunes, I.J. Thompson, E.F. Brown, for the good agreement was analytically clarified using Phys. Rev.Lett. 109, 141101 (2012). R-matrix theory. It was found that extremely slow con- [15] S. Ishikawa, Phys. Rev.C87, 055804 (2013). vergence occurs if a coupled-channel expansion of the [16] N.B. Nguyen, F.M. Nunes, and I.J. Thompson, Phys. Rev. C87, 054615 (2013). wave function is used, which helps to explain the very [17] Y. Suzuki,P. Descouvemont,arXiv:1308.4021. different reaction rates obtained using different theoreti- [18] K. Yabana, Y.Funaki,Phys. Rev.C85, 055803 (2012). cal approaches. [19] T. Seideman and W.H. Miller, J. Chem. Phys. 96, 4412 (1992); 97, 2499 (1992). [20] A.M. Lane, R.G. Thomas, Rev. Mod. Phys. 30, 257 (1958). [21] R. de Diego, E. Garrido, D.V. Fedorov, A.S. Jensen, [1] E.E. Salpeter, Astrophys.J. 115, 326 (1952). Phys. Lett. B695, 324 (2011). [2] F. Hoyle, Astrophys.J. Suppl.1, 121 (1954). [22] R.deDiego,E.Garrido,D.V.Fedorov,andA.S.Jensen, [3] K.Nomoto,F.-K.Thielemann,S.Miyaji,Astron.Astro- Europhys. Lett. 90, 52001 (2010). phys.149, 239 (1985). [23] S. Ali and A.R. Bodmer, Nucl. Phys.80, 99 (1966). [4] K. Langanke, M. Wiescher, and F.-K. Thielemann, Z. [24] T. Yamada and P. Schuck, Eur. Phys. J. A 26, 185 Phys.A324, 147 (1986). (2005). [5] P. Descouvemont and D. Baye, Phys. Rev. C36, 54 [25] K. Nomoto, Astrophys.J. 253, 798 (1982). (1987).