Global description of beta-minus decay in even-even nuclei with the axially-deformed Skyrme finite amplitude method
GGlobal description of beta-minus decay in even-even nuclei with the axially-deformedSkyrme finite amplitude method
M. T. Mustonen
1, 2, ∗ and J. Engel † Department of Physics and Astronomy, CB 3255,University of North Carolina, Chapel Hill, NC 27599-3255 Center for Theoretical Physics, Sloane Physics Laboratory,Yale University, New Haven, Connecticut 06502, USA (Dated: October 18, 2018)We use the finite amplitude method for computing charge-changing Skyrme-QRPA transitionstrengths in axially-deformed nuclei together with a modern Skyrme energy-density functional to fitseveral previously unconstrained parameters in the charge-changing time-odd part of the functional.With the modified functional we then calculate rates of beta-minus decay for all medium-mass andheavy even-even nuclei between the valley of stability and the neutron drip line. We fit the Skyrmeparameters to a limited set of beta-decay rates, a set of Gamow-Teller resonance energies, and a setof spin-dipole resonance energies, in both spherical and deformed nuclei. Comparison to availableexperimental beta-decay rates shows agreement at roughly the same level as in other global QRPAcalculations. We estimate the uncertainty in our rates all the way to the neutron drip line through aconstruction that extrapolates the errors of known beta-decay rates in nuclei with intermediate Q values to less stable isotopes with higher Q values. PACS numbers: 21.60.Jz, 23.40.HcKeywords: finite amplitude method, beta decay, neutron-rich nuclei
I. INTRODUCTION
Beta-decay rates are an important ingredient in simula-tions of the astrophysical r -process. Because parts of the r -process path are still not accessible to experiment, it isup to theoretical models to produce approximate rates formany relevant neutron-rich isotopes. Models could alsohelp resolve the issues raised by Ref. [1], which arguedthat the flux of antineutrinos from nuclear reactors doesnot agree with the standard model. Ref. [2] pointed outthat the discrepancy could be due to an overly simpletreatment of first-forbidden beta decay in fission products.Several schemes/methods for calculating beta-decayrates across almost the entire nuclear chart have beendevised. They include the proton-neutron quasiparticlerandom phase approximation (pnQRPA) with a residualGamow-Teller interaction [3], gross theory [4, 5], semi-gross theory [6], the extended Thomas-Fermi plus Struti-nsky integral (ETFSI) method [7], and artificial neuralnetworks [8], and, very recently, the relativistic sphericalpnQRPA [9]. Full beta-decay tables for neutron-rich iso-topes have been only published by M¨oller et al [4, 5] andMarketin et al [9].Many other authors have applied more sophisticatedand/or computationally intensive methods to smaller setsof nuclei, focusing on some of those important for the r process. Recent examples of such work include a deformedpnQRPA computation with the Bonn-CD interaction ofthe decay of neutron-rich isotopes with Z = 36 − ∗ [email protected] † [email protected] and Sr [12]; similar calculations with Gogny interactionin the N = 82 , ,
184 isotonic chains [13]; relativisticpnQRPA [14] for 20 ≤ Z ≤
50; and relativistic pnQRPAfor N ≈ ,
82 [15].Computational barriers have thus far prevented theproduction of a beta-decay table for the entire nuclearchart in a fully self-consistent Skyrme mean-field approachthat allows deformation. Recently, however, we reported[16] an implementation of the charge-changing finite am-plitude method, which sidesteps the QRPA eigenvalueproblem. We obtain beta-decay rates by directly com-puting the required sums and integrals over allowed finalstates of the response to charge-changing perturbations.We will soon make available a code called pnfam thatimplements the method.We could proceed by choosing an existing density func-tional, interpreting it as a density-dependent effectiveinteraction, and calculating beta-decay rates. If we wereinterested in, e.g., the effects of tensor terms discussedin Refs. [17] and [18], we could take them from already-parametrized functionals. Such a procedure would requiresome parameter fitting because pairing interactions andstrengths, especially those associated with isoscalar pair-ing, are not usually specified alongside Skyrme particle-hole effective interactions. But that approach is still toolimiting because not all Skyrme functionals can be consis-tently represented as effective interactions. In particular,the time-even and time-odd parts of the functional, whichare related if the functional is the mean-field expecta-tion value of a Hamiltonian, need not be related in moregeneral constructions.Our main goal here is to assess the ability of the SkyrmeQRPA with deformation to predict β − decay, and to useexisting data (decay rates and resonance energies) to a r X i v : . [ nu c l - t h ] O c t constrain the isoscalar-pairing strength and the othertime-odd coupling constants, which in the general energy-density functional (EDF) picture, are not fixed by fitsto masses, precisely because they are independent of thetime-even functional. In much of what follows, therefore,we will not assume that the functional results from mean-field theory with an interaction and so will have to fit asignificant number of parameters. After presenting ourmethods and assessment, we display the (summarized)results of a a full table of beta-minus rates, computed withthe pnFAM, for even-even nuclei in all medium-mass andheavy isotopic chains. We use a simple but apparentlyaccurate model to quantify and extrapolate theoreticaluncertainty.This article is organized as follows: Section II is a briefoverview of the theoretical background, and Section IIIdetails our computational approach and parameter-fittingprocedures. In Section IV, we assess the quality of our re-sults, comparing them to earlier work and to experimentaldata where available. Section V contains conclusions. II. THEORETICAL BACKGROUNDA. Finite amplitude method
The finite amplitude method (FAM), a formulation ofthe random-phase approximation that speeds the compu-tation of nuclear response functions, was introduced innuclear physics in Ref. [19]. It was later generalized tothe quasiparticle random phase approximation (QRPA)in Ref. [20] and to the relativistic QRPA in Ref. [21].In Ref. [16] we applied the method to charge-exchangetransitions, in particular allowed and first-forbidden betadecay.The FAM solves equations for the amplitude of thelinear response to a small but finite perturbation. As aresult, the method does not directly yield the poles andresidues of the response, which are the central objects inthe matrix version of the QRPA. But if the goal of thecomputation is to get transition strength functions in alarge model space the FAM can yield results in orders ofmagnitude less CPU time than matrix QRPA.The FAM offers another advantage for beta decay: theweighted sums or integrals of transition strength can be ex-pressed as contour integrals. This fact was first exploitedby Hinohara et al. [22], who evaluated the response at arelatively small number of complex frequencies to com-pute sum rules. In Ref. [16] we used the idea to evaluatethe more complicated beta-decay phase-space-weightedintegrals, which are not analytic and contain interferenceterms between first-forbidden operators. With the FAMwe can thus use typical supercomputer resources to calcu-late many observables in a large number of nuclei. We areable to extend systematic Skyrme parameter fitting frommean-field calculations to deformed QRPA calculations,at least in a preliminary way.
B. Model parameters and fitting targets
Our starting point in the particle-hole channel is ageneric Skyrme EDF: E = (cid:88) t =0 , t (cid:88) t = − t (cid:90) d r (cid:0) H even tt ( r ) + H odd tt ( r ) (cid:1) . (1)Here H even tt contains products of time-even local densi-ties only, with coefficients fixed by fitting to masses andperhaps a few other quantities, and H odd tt is given by H odd tt ( r ) ≡ C st [ ρ ] s tt + C ∆ st s tt · ∇ s tt + C jt j tt + C Tt s tt · T tt + C s ∇ jt s tt · ∇ × j tt + C Ft s tt · F tt + C ∇ st ( ∇ · s tt ) , (2)with the spin density s tt , the current density j tt , the spin-kinetic density T tt , and the tensor-kinetic density F tt defined e.g. in Ref. [23]. If one requires the functional to bethe mean-field expectation value of a Skyrme interaction,EDF coupling constants are completely determined bythe (fewer) parameters that specify the interaction, asdiscussed in Refs. [23, 24] and mentioned above. In thiswork, we adopt the view that the effective “interaction”comes from the energy-density functional rather than theother way around. Consequently we are free to fit all thetime-odd coupling constants without spoiling the mass fitsgenerated the time-even couplings. We do, however, adoptthe values obtained from the Skyrme parametrization asour starting point for the fits, unless we note otherwise.The subset { C s , C T , C F } of time-odd coupling con-stants maps directly to the parameters of Landau-Migdalinteraction for infinite homogeneous nuclear matter withtensor terms [25]. In that sense, these couplings are in-timately related to the bulk properties of the nuclearmatter. The constant C F is purely tensor in character,and it determines the tensor term in the Landau-Migdal in-teraction. The spin-density coupling constant C s stronglyaffects the Gamow-Teller strength distribution and can befit to the locations of Gamow-Teller resonances [24]. Thelast term C T maps to a linear combination of the Lan-dau parameters of both the tensor term and a term thatdepends on the scattering angle of the Laundau-Migdalquasiparticles.Two other parameters have a large effect in the QRPA:the strengths of the residual particle-particle (or pairing)interaction between protons and neutrons, V pp = (cid:16) V ˆΠ T =0 + V ˆΠ T =1 (cid:17) (cid:18) − α ρ ( r ) ρ c (cid:19) δ ( r ) , (3)where ρ c = 0 .
16 fm − is the saturation density of nu-clear matter and α ∈ [0 ,
1] controls the pairing density-dependence (throughout this work we use mixed pairing,i.e. α = 0 . T = 1) proton-neutron pair-ing mainly affects Fermi beta decay, which plays almostno role in heavy nuclei; we simply set its strength V tobe the average of the neutron-neutron and proton-protonpairing strengths (both fixed in the HFB part of the calcu-lation). The isoscalar ( T = 0) pairing is a different story;it has a strong effect on Gamow-Teller decay, and deter-mining a reasonable value for its strength is a commonissue in single and double beta decay computations. TheHFB mean field is independent of the T = 0 pairing termas long as explicit proton-neutron mixing is neglected.We are thus free to fit V in the QRPA.We include three types of observables in our fittingprocedure: Gamow-Teller resonance energies, spin-dipoleresonance energies, and total beta decay rates, all in bothspherical and deformed isotopes. Although we pick severalsets of target nuclei for fitting the decay half-lives, eachset includes a large range of mass values so that our fitscan be global. To assess the success of the fits and tocompare them to earlier work in very different models, wealso compute the following metrics for beta-decay tables,as laid out e.g. in [4, 5]: the residual of each computedlg t value (where t is the half-life and the logarithm is10-based), r = lg (cid:18) t th. t exp. (cid:19) , (4)the average of the residuals, M r = 1 n n (cid:88) i =1 r i , (5)(where we’ve used an index i to indicate that there is an r for every nucleus) and the standard deviation of theresiduals around the average, σ r = (cid:118)(cid:117)(cid:117)(cid:116) n n (cid:88) i =1 (cid:2) ( r i − M r (cid:3) . (6)Ref. [8] contains an excellent compilation of these quanti-ties. C. Uncertainty analysis
As theoretical approaches grow more sophisticated, theanalysis of theoretical uncertainty is growing in impor-tance. Here we attempt to provide reasonable estimatesfor the uncertainty in our predicted half-lives, particularlyin isotopes that are too short-lived to allow measurement.The standard prescription for assigning a theoretical(statistical) uncertainty ∆ O to a computed observable O is [26] ∆ O = (cid:115)(cid:88) ab ∂ O ∂x a (cid:12)(cid:12)(cid:12)(cid:12) x = x C ab ∂ O ∂x b (cid:12)(cid:12)(cid:12)(cid:12) x = x , (7)where C is the covariance matrix C = ( J T J ) − (8) and x = ( x , . . . , x N x ) are the N x parameters of the model.The partial derivatives of all the observables {O a } withrespect to all the parameters evaluated at the result ofthe fit x form the Jacobian J : J ab = ∂ O a ∂x b (cid:12)(cid:12)(cid:12)(cid:12) x = x . (9)When not analytically accessible, the needed partialderivatives can be estimated through finite central dif-ferences. (In general, the theoretical uncertainty relatedto the Jacobian must be supplemented by numerical andexperimental uncertainties. We have assumed that thetheoretical uncertainty is much larger than the other two,which we therefore neglect.)To use Eq. (7) to assign an uncertainty to every beta-decay rate in our table, we would need to evaluate theJacobian in Eq. (9) for (the logarithms of) all the half-lives t in our table. Unfortunately, the required 2 N x full decay-table computations are still not possible in a reasonableamount of computer time, even with the efficiency of theFAM. We therefore attempt to gauge the uncertainties andtheir Q dependence in a more na¨ıve way. We constructa simple few-parameter model for the uncertainties thatwe can then fit to the observed differences between ournumerical results and known experimental values. Theresulting approach is agnostic about how the decay isactually calculated; it treats the nuclear model as a blackbox that produces predictions for Q values and beta-decay rates. It does, however, make several assumptionsabout both the calculated and experimental strengthdistributions that are only approximately correct.The first assumption is that the final states that con-tribute significantly to a decay rate lie not too far fromthe ground state in a relatively small window of excita-tion energy, so that we can reasonably approximate them,again either in our calculation or in the real world, byone effective state with an effective Q value q eff (that isnot too different from the ground-state Q value): q eff = (cid:80) k C k f k ( q k + 1 , Z f ) q k (cid:80) k C k f k ( q k + 1 , Z f ) . (10)Here C k is the standard integrated shape function forthe transition to the final state k — for an allowed statedecay this is simply the transition strength, and for anon-unique forbidden decay it encapsulates several terms— and f k is the usual allowed phase-space integral. Thequantity q k is the Q value of the decay to the state k divided by m e c , and Z f is the charge of the daughternucleus.With these definitions, we can proceed to define aneffective shape factor C eff as t = κ (cid:80) k C k f ( q k + 1 , Z f ) = κC eff f ( q eff + 1 , Z f ) . (11)The C k depend on the q k for forbidden transitions, but thedependence is weak compared to that of the correspondingphase-space integral, and so we neglect it in our effectiveshape factor.All these definitions can be made independently for theexperimental strength distribution and the theoreticalone. The quantity we wish to understand is the ratio r of the theoretical and experimental lifetimes: r = lg t th t exp = lg C expeff C theff + lg f ( q expeff + 1 , Z f ) − lg f ( q theff + 1 , Z f ) , (12)where the meanings of q expeff and q theff and the correspondingquantities C expeff and C theff should be clear. We omit thenucleus index i in them for brevity. Because the phasespace grows quickly with decay energy, the effective Q values will usually be close to ground-state-to-ground-state Q value. We therefore expand both logarithms inEq. (12) of the phase-space factors about the theoreticalground-state-to-ground-state Q value q thg.s. to first orderin q :lg f ( q + 1 , Z f ) ≈ lg f ( q thg.s. + 1 , Z f )+ d lg f ( q + 1 , Z f ) dq (cid:12)(cid:12)(cid:12)(cid:12) q = q thg.s. · ( q − q thg.s. )= lg f ( q thg.s. + 1 , Z f )+ f (cid:48) ( q thg.s. + 1 , Z f ) f ( q thg.s. + 1 , Z f ) · q − q thg.s. ln 10 . (13)This approximation is best when the Q value is high,because the curvature of lg f ( q + 1 , Z f ) is small at high q .For Q values lower than about 2–3 MeV the first-orderapproximation is poor, so we exclude such data pointsfrom our analysis.Replacing the two logarithms in Eq. (12) by the first-order expressions, we find that several terms cancel in thedifference, so that r ≈ lg C expeff C theff + f (cid:48) ( q thg.s. + 1 , Z f ) f ( q thg.s. + 1 , Z f ) · q expeff − q theff ln 10 . (14)We now make another set of assumptions, this timeabout the distribution of the errors in the theoreticalvalues: First, we assume that the relative error in theeffective shape factor is normally distributed with a slightsystematic component. That is, we assume that for eachnucleus, we have lg C expeff C theff = c + δc, (15)where c is a constant, independent of the nucleus inwhich the decay occurs, that captures the systematic error and the set of all nucleus-dependent δc ’s is normallydistributed with standard deviation ∆ c , which is a mea-sure of statistical error that is independent of decay energy.Similarly, we assume that the errors in the theoreticaleffective Q value follow a normal distribution, so that q expeff − q theff ln 10 = q + δq, (16)where q is now a nucleus-independent systematic error inthe effective Q value and the set of δq ’s (again, one foreach nucleus) is normally distributed around zero with thestandard deviation ∆ q that represents statistical error,again independent of Q . (Both q and δq are expressedin units of m e c , and the factor 1 / ln 10 is absorbed into q and δq for convenience.) Finally, we assume that theerrors δq and δc are independent.With these assumptions, we then have r = c + f (cid:48) ( q thg.s. + 1 , Z f ) f ( q thg.s. + 1 , Z f ) · q + δr, (17)where δr is a random error, the set of which must benormally distributed at each q g.s.th with width∆ r ( q thg.s. ) = ∆ c + (cid:18) f (cid:48) ( q thg.s. + 1 , Z f ) f ( q thg.s. + 1 , Z f ) (cid:19) ∆ q . (18)We can now use Eqs. (17) and (18) to determine the valuesof c , q , ∆ c , and ∆ q . We obtain the first two through afit to Eq. (17), which expresses r i as a function of f (cid:48) /f (and thus of q thg.s. ), with the set of ratios r i given by ourcalculations (and experiment) and δr i set to zero. Finally,we insert the fit values of c and q into the square of Eq.(17), which then expresses δr i as a function of ( f (cid:48) /f ) ,and determine the values of ∆ c and ∆ q by requiringthat the line for ∆ r ( q thg.s. ) as a function of ( f (cid:48) /f ) thatexpresses Eq. (18) goes through the middle of the data,with as many points below ∆ r as above it. A secondleast squares fit of ∆ c and ∆ q in the right hand side ofEq. (18) to the set of δr i ’s accomplishes the task nicely.The explicitly Q -dependent ∆ r can then be extrapolatedto large Q values where there are no data.For both fits, in practice, we only include data pointswith Q ≥ | r i | > f (cid:48) /f a simpleratio of two polynomials, independent of Z f : f (cid:48) ( x, Z f ) f ( x, Z f ) ≈ x − x + 15 x − x + 15 x − . (19)Once the four parameters of the uncertainty model havebeen determined, the uncertainties of the theoretical pre-dictions are obtained the inverse relation of (4) t exp = t th · r (20)with r = r ( q thg.s. ) ± ∆ r ( q thg.s. ) . (21) -1 0 1 2 3 4 0 2 4 6 8 10 12 14 Q c o m p . - Q e x p . [ M e V ] Q exp. [MeV](a) SV-min0.626 � p o i n t s Q comp. -Q exp. [MeV]-2-1 0 1 2 3 4 0 2 4 6 8 10 12 14 Q c o m p . - Q e x p . [ M e V ] Q exp. [MeV](b) SkO'0.154 � p o i n t s Q comp. -Q exp. [MeV] Figure 1. The differences between computed and experimental Q values, with SkO’ (top) and with SV-min (bottom). Theinsets show the distribution of differences. The errors inour computed Q values follow a normal distribution with anaverage of 0 .
154 MeV and a standard deviation of 0 .
576 MeV,with no noticeable bias when moving to higher Q values. III. SKYRME FUNCTIONAL ANDCOMPUTATIONAL METHOD
The first step in our computational procedure is theconstruction of ground states in the doubly-even mothernucleus with hfbtho , a well-established HFB solver work-ing in a (transformed) harmonic-oscillator basis [36]. Wecut off the single-particle space at 60 MeV to avoid diver-gences from our zero-range pairing. For each nucleus, wesearch for a prolate, an oblate, and a spherical solution,and take the most bound of these to be the ground state.Because the beta-decay rates are very sensitive to the Q value of the decay, we look for a modern Skyrme func-tional that reproduces Q values well. Because we don’texplicitly treat odd nuclei, we use the prescription of Ref.[37] to approximate the Q value; we have checked theprescription against odd-A calculations in the equal fillingapproximation, and the two procedures generally agreeto within about 0.5 MeV. Of the several functionals weexamine, SkO’ [38] (with the strengths of proton-protonand neutron-neutron pairing fit to the experimental pair-ing gaps of ten isotopes picked in a wide mass range50 ≤ A ≤ Q values, produc-ing errors for ground-state-to-ground-state Q values thatare normally distributed, with an average systematic errorof 0 .
154 MeV and statistical error of 0 .
576 MeV. Figure 1compares the Q values produced by SkO’ with those ofthe next best functional, SV-min.To compute beta-decay rates and resonance energieswe use the code pnfam , an implementation of the charge-changing finite amplitude method presented in [16]. Builtto work together with hfbtho , pnfam allows us to com-pute properties of axial deformed nuclei, including bothallowed and first-forbidden beta decay.We obtain our most robust fits by fixing all but twoof the time-odd coupling constants of the functional atvalues implied by its interpretation as an interaction.(The mapping of the Skyrme coupling constants to theenergy-density functional coupling constants is explicitlydiscussed in Refs. [23, 24].) The exceptions are C s and C ∆ s . We set the latter to zero to avoid known finite-size instabilities [39] which, in the case of hfbtho and pnfam , manifest themselves as divergences in the iter-ative solution. That leaves C s , which, along with theisoscalar pairing strength V , we fit to a set of Gamow-Teller resonance energies, spin-dipole resonance energies,and beta-decay rates selected from a wide mass range withno particular region favored. We use the code pounders ,based on a derivative-free algorithm [40] designed for op-timizing computer-time-consuming penalty functions, toefficiently minimize the weighted the sum of the leastsquares, simultaneously fitting both parameters. We take C s to be independent of the density and V to have thesame density dependence as the proton and neutron pair-ing. The axial-vector coupling constant g A is known tobe quenched in nuclei, but the source and magnitude ofthe quenching is an open problem; a variety of very dif-ferent values for an effective g A have been used. For lackof a better prescription, we use the commonly adoptedquenched value g A = 1 . χ to the number of degrees of freedom, followingthe recommendation in Ref. [26]. We assume that thetheoretical error dominates the experimental error, andthus assign equal weight to observable of the same type.We select these weights based on how well the differentobservables are reproduced in an initial test fit, so that set GT resonances SD resonances beta-decay half-livesA Pb,
Sn, Ge,
Te, Zr, Ca none Ar, Cr, Ni, Zn, Kr,
Sr,
Ru,
Cd,
Sn,
BaB same as A none Ti, Zn, Sr,
Pd,
Te,
Sm,
Yb,
Pt,
Rn,
UC same as A none Ti, Ni, Sr,
Ru,
Te,
Nd,
Yb,
Pt,
Rn,
UD those of A and
Nd none Ti, Zn, Kr,
Cd,
Ce,
Gd,
PtE same as D Zr, Pb Ti, Zn, Kr,
Cd,
Ce,
Gd,
RnTable I. The sets of fitting targets used in this work. The beta-decay half-lives in set A range from 0 .
069 s (
Sr) to 1 .
84 s( Kr), in set B from 95 . Zn) to 45360 s (
Pt), in set C from 0 .
54 s (
Ru) to 9399 . Sr), and in set E from 0 .
046 s( Kr) to 444 s (
Rn). The nuclei selected for fitting the beta-decay half-lives in sets D and E all exhibit an excitation spectrumclearly associated with either a spherical or a well-deformed shape; set E only consists of open-shell nuclei. The experimentalGamow-Teller resonance energies are from Refs. [28–32], the spin-dipole resonance energies from Refs. [33, 34], and the half-lifesfrom Ref. [35]. each type of observable is approximately equally weightedin the actual fit. Typical fits then take 10—20 thousandCPU hours, and we use XSEDE supercomputers [41] tocarry them out.Following the fit, we proceed to compute the beta-decayrates of all even-even neutron-rich nuclei with 28 ≤ Z ≤ A ≥
50, all the way to the neutron drip line, omittingjust a few very stable isotopes for which the Q value isnegative in our HFB calculations. IV. RESULTS AND DISCUSSIONA. Fit results
To assess how sensitive our fit is to the set of targetbeta-decay rates and resonance energies, we repeat theprocess with four sets of rates, summarized in Table I.Each of these sets spans a large range of masses. Set Acontains beta-decay isotopes with relatively short half-lives only, set B relatively long half-lives only, and set Ca wide range of half-lives. (Short half-lives should be lesssensitive to details in nuclear structure). Set D containsonly nuclei that are known to be rather rigid, with anexcitation spectrum characteristic of a spherical or a well-deformed nucleus. (The QRPA, which is based on a singlemean field, should work best in rigid nuclei). In set E, weinclude only open-shell rigid nuclei (for which isoscalarpairing should be most effective), swapping out
Pt for
Rn, and including two spin-dipole resonances. Figure 2shows the quality of the fits for set E both with computed(1D) and experimental (1E) Q values (see Table III fordefinitions of the number-letter combinations). The twoprocedures yield very similar results. The comparisonin Fig. 3, which shows the results of the two-parameterfits (along with those of more-parameter fits discussedshortly) when the resulting functionals are applied to theset of all measured even-even β − -decay half-lives, showsthat all these two-parameter fits (1A through 1E) yieldthe same level of predictive accuracy.Can we do better by including some of the other time-odd coupling constants in our fit? To find out, we refit the four-parameter set { V , C s , C T , C F } , which determinesthe Landau parameters { g (cid:48) , g (cid:48) , h (cid:48) } and thus allows usto incorporate infinite-nuclear-matter stability conditions[25] as constraints. The parameter C F introduces a time-odd tensor term that is not present in the two-parameterfits. The result, however, improves the description ofthe fitting targets only marginally, as the points labeled3A and 3B in Fig. 2 show, and actually worsens theagreement with half-life measurements overall (as Fig. 3shows). The situation gets even worse when we use theresults of this fit as a starting point to fit three more time-odd coupling constants, { C j , C ∇ j , C ∇ s } (fit 4). Then thebeta-decay rates to which we fit are reproduced better,but the agreement with all measured rates deteriorates.Improvement in the fitting targets accompanying a de-terioration in overall agreement with data is a symptomof overfitting. To better understand why this happens,we evaluate the Jacobian matrix (9) at the parametervalues produced by the two-parameter fit 1E. The Jaco-bian appears in Table II, with the values of the couplingconstants in natural units following the prescription ofRef. [42] and the natural scale of isoscalar pairing takento be the strength of isovector pairing. A clear columnstructure appears in both the resonance energies and thehalf-lives, signaling that the members of each individualset move largely in unison when the parameters are varied.Thus there are essentially just two meaningful degrees offreedom that we can expect to fix with this experimentaldata: V and C s .To see this in more detail, we carry out a singular valuedecomposition of the Jacobian. The largest singular value,122 .
53, corresponds to a vector pointing nearly in the di-rection of C s in parameter space, and the second largest,10 .
85, to a vector pointing nearly in the direction of V .The third largest value, 1 . C T , with many other directions mixed in. The charge-changing data we have available — Gamow-Teller andspin-dipole resonances, and half-lives — are not enoughto reasonably constrain more than the parameters V and C s in our initial fit.Figure 3, besides containing the results of our fits, -8-6-4-2 0 2 4 6 8 Ca Ge Zr Sn Te Nd Pb E c o m p u t e d - E e x p . GT resonances -4-2 0 2 4 Zr Pb E c o m p u t e d - E e x p . SD resonances 10 -2 -1 Ti Zn Kr Cd Ce Gd Rn t c o m p . / t e x p . half-lives 1E1D3A3B45 Figure 2. Reproduction of target data in those of our fits that use the target set E. The four-parameter fits 3A and 3B yieldalmost the same results as the two-parameter fits 1D and 1E, with only a small decrease in the penalty function. contains results from other work: Refs. [3, 6, 8] and [9].Of all the these computations, the one by Homma et al.[3] seems to best reproduce the known β − half-lives, eventhough it neglects non-unique first-forbidden decay anduses simple separable interactions. As Figure 4 shows,in our computation the non-unique 1 − contribution isquite important (even dominant) in many experimentallyinaccessible nuclei, so it is far from clear how the variouscalculations will fare with data in the future. In anyevent, the most striking fact is that all the computationsmanage to existing data at roughly the same level ofprecision. It may not be possible to do much betterwithout moving beyond Skyrme QRPA, at least whileusing a global parameter set as we have done here. B. Extrapolation to neutron-rich isotopes
Figure 4 displays the relative contribution to the decayrate from each multipole. Except in the immediate vicinityof the valley of stability, the changes appear quite gradualas a function of Z and N . In nuclei with large Q values,the details of single-particle structure are less importantthan in isotopes for which transitions to only a few low-energy states are possible.Fig. 4 also demonstrates the importance of going be-yond the allowed approximation. In many heavy nuclei,the computed rates are dominated by the first-forbiddenchannel. Towards the drip line, both allowed and forbid-den channels are important for all masses. The figure alsoshows that the non-unique 1 − channel is usually the mostimportant of the forbidden multipoles. Thus any quench-ing of the (unique) 2 − channel and anti-quenching of the (non-unique) 0 − channel from meson-exchange currents[43] would not have a significant impact on our overallresults. In the 1 − channel the contributions of severaldifferent operators makes the effects of quenching hardto estimate.In Figure 5 we compare our half-lives to those of M¨oller et al. [5] in all medium and heavy even-even isotopes. Ourhalf-lives tend to be longer then those of Ref. [5] close tothe valley of stability in light nuclei and somewhat shorterin heavy nuclei (with significant forbidden contributions).Approaching the neutron drip line, the two computationsyield similar results up to a constant offset in those ofRef. [5] in even-even nuclei. All models can expect to dobetter near the drip line, where a significant fraction ofthe total β -decay strength can be below threshold.Because our na¨ıve model for uncertainties is based onseveral assumptions that are only approximately corrector cannot easily be verified, we check its predictions wherethere are enough data to do so. Figure 6 shows the ratiosof our half-lives to those of experiment together withthe uncertainty model’s mean value and one- and two-standard-deviation bands, all as a function of ground-state Q value. We can discount the model at very low Q butit appears to work well above Q ≈ Q is generally large.A recent RIBF measurement [44] of 110 neutron-richisotopes, 40 of them previously unknown, allows us totest the reliability of our predictions and especially our O d O /dC s d O /dV d O /dC F d O /dC T d O /dC ∇ s d O /dC ∆ s d O /dC j d O /dC ∇ j Pb E GTR Sn E GTR Ge E GTR Te E GTR Zr E GTR Ca E GTR Pb E SDR Zr E SDR Ti lg t Zn lg t Kr lg t Cd lg t Ce lg t Gd lg t Pt lg t t values are hence dimensionless and those of the resonance energies are in the units of MeV.fit starting point target set Q values fitted parameters1A SkO’ A comp. V = − . C s = 128 . V = − . C s = 133 . V = − . C s = 126 . V = − . C s = 129 . V = − . C s = 99 . V = − . C s = 132 . V = − . C s = 144 . C T = − . C F = − . V = − . C s = 120 . C T = − . C F = − . C j = 54 . C ∇ j = − . C ∇ s = − .
55 SkO’ E comp. V = − . C s = 146 . C j = − . C ∆ s , which is set to zero everywhere to avoid finite-size instabilities. The units of V and C s are MeV fm , and the units of the other coupling constants are MeV fm . model for theoretical uncertainties. Because the dataare so recent, we did not include them in any of ourfits, and hence we are effectively using older data topredict the results of these new measurements. We have28 even-even nuclei with which to compare rates; forhalf of these there are earlier data in the ENSDF set(Figure 7). Our predictions agree with experiment towithin our theoretical uncertainty (though our error barsmay be a bit too large here). Our uncertainty model thusappears to be reasonable. V. CONCLUSIONS
We have explored the ability of the axially-deformedSkyrme QRPA to provide a global description of beta-decay rates in even-even neutron-rich nuclei. With ex-perimental rates and charge-exchange resonance energiesas fitting targets we have found that among time-oddcouplings, only those multiplying the isoscalar-pairingand the spin-density parts of the functional are well con-strained; attempts to fit more than these two constantslead to overfitting. The tensor contributions to the energy-density functional are, in particular, not well constrainedby this data. To get more accurate Skyrme-QRPA pre-dictions, one can resort to local fits, i.e. A -dependentcouplings. The recent work of Ref. [9], for example, at- -1.5-1-0.5 0 0.5 1 1.5 1A 1B 1C 1E 2 3A 3B 4 5 Mö Ho Na Co Ma l g ( t c o m p . / t e x p . ) (a) t ≤ l g ( t c o m p . / t e x p . ) (b) t ≤
100 s-1.5-1-0.5 0 0.5 1 1.5 1A 1B 1C 1E 2 3A 3B 4 5 Mö Ho Na Co Ma l g ( t c o m p . / t e x p . ) (c) t ≤
10 s-1.5-1-0.5 0 0.5 1 1.5 1A 1B 1C 1E 2 3A 3B 4 5 Mö Ho Na Co Ma l g ( t c o m p . / t e x p . ) (d) t ≤ Figure 3. Comparison of the mean and standard deviationof the lg t values in our fits with those of previous work. Thelabels for our fits correspond to those of Table I. The resultsfrom prior work are contained in [5] (M¨o), [3] (Ho), [6] (Na),[8] (Co), and [9] (Ma). Only even-even isotopes are considered. Figure 4. The contributions of different allowed and first-forbidden multipoles to the total computed beta-decay rates.Only even-even nuclei are considered. Z Nlg (t
FRDM /t our ) -4-2 0 2 4 Figure 5. Comparison of our computed half-lives in neutron-rich nuclei with those of Ref. [5]. taches a sensible A dependence to the strength of isoscalarpairing.The level of agreement between our calculations anddata throughout the isotopic chart is similar to that pro-duced by other recent computations, in spite of our con-sistent inclusion of deformation, tensor terms in the func-tional, etc. It could be difficult to do much better withoutan account of multiphonon effects, as in the work, e.g., ofRef. [45].The most glaring shortcoming of our work here is the re-striction to even-even nuclei. An extension of the FAM toodd-mass nuclei will be the subject of a future publication[46]. For the moment, we make our results for the 1387even-even neutron-rich nuclei, with crudely estimated the-oretical uncertainties, available as supplementary materialto this article. -4-3-2-1 0 1 2 3 4 0 2 4 6 8 10 12 14 r Q (MeV) uncertaintySkO'bias
Figure 6. The fit 3A with the mean (green line) and one- andtwo-standard-deviation bands from our uncertainty model. -3 -2 -1 S r S r M o Z r M o M o M o R u R u R u P d P d C d S n h a l f - l i f e ( s ) FAM predictionRIBF measurement
Figure 7. Our predictions for the 14 half-lives of neutron-richeven-even nuclei measured only recently [44] and not includedin the ENSDF data set. All the measured half-lives fall withinour one-sigma error bars, suggesting that our uncertaintyestimates are too pessimistic in this particular region.
ACKNOWLEDGMENTS
We thank W. Nazarewicz, J. Dobaczewski, and S. Wildfor useful discussions. Support for this work was pro-vided through the Scientific Discovery through AdvancedComputing (SciDAC) program funded by US Departmentof Energy, Office of Science, Advanced Scientific Com-puting Research and Nuclear Physics, under ContractNo. DE-SC0008641, ER41896. This work used the Ex-treme Science and Engineering Discovery Environment(XSEDE), which is supported by National Science Foun-dation grant number ACI-1053575. [1] G. Mention, M. Fechner, T. Lasserre, T. A. Mueller,D. Lhuillier, M. Cribier, and A. Letourneau, Phys. Rev. D , 073006 (2011). [2] A. C. Hayes, J. L. Friar, G. T. Garvey, G. Jungman, andG. Jonkmans, Phys. Rev. Lett. , 202501 (2014).[3] H. Homma, E. Bender, M. Hirsch, K. Muto, H. V.Klapdor-Kleingrothaus, and T. Oda, Phys. Rev. C ,2972 (1996).[4] P. M¨oller, J. R. Nix, and K. L. Kratz, Atomic Data andNuclear Data Tables , 131 (1997).[5] P. M¨oller, B. Pfeiffer, and K.-L. Kratz, Phys. Rev. C ,055802 (2003).[6] H. Nakata, T. Tachibana, and M. Yamada, Nucl. Phys.A , 521 (1997).[7] I. N. Borzov and S. Goriely, Phys. Rev. C , 035501(2000).[8] N. J. Costiris, E. Mavrommatis, K. A. Gernoth, andJ. W. Clark, Phys. Rev. C , 044332 (2009).[9] T. Marketin, L. Huther, and G. Mart´ınez-Pinedo, (2015),arXiv:1507.07442.[10] D.-L. Fang, B. A. Brown, and T. Suzuki, Phys. Rev. C , 024314 (2013).[11] D. Ni and Z. Ren, Phys. Rev. C , 064320 (2014).[12] D. Ni and Z. Ren, J. Phys. G: Nucl. Part. Phys. ,125102 (2014).[13] M. Martini, S. P´eru, and S. Goriely, Phys. Rev. C ,044306 (2014).[14] Z. M. Niu, Y. F. Niu, H. Z. Liang, W. H. Long, T. Nikˇsi´c,D. Vretenar, and J. Meng, Phys. Lett. B , 172 (2013).[15] T. Nikˇsi´c, T. Marketin, D. Vretenar, N. Paar, and P. Ring,Phys. Rev. C , 014308 (2005).[16] M. T. Mustonen, T. Shafer, Z. Zenginerler, and J. Engel,Phys. Rev. C , 024308 (2014).[17] F. Minato and C. L. Bai, Phys. Rev. Lett. , 122501(2013).[18] V. De Donno, G. Co’, M. Anguiano, and A. M. Lallena,Phys. Rev. C , 024326 (2014).[19] T. Nakatsukasa, T. Inakura, and K. Yabana, Phys. Rev.C , 024318 (2007).[20] P. Avogadro and T. Nakatsukasa, Phys. Rev. C , 014314(2011).[21] T. Nikˇsi´c, N. Kralj, T. Tutiˇs, D. Vretenar, and P. Ring,Phys. Rev. C , 044327 (2013).[22] N. Hinohara, M. Kortelainen, W. Nazarewicz, andE. Olsen, Phys. Rev. C , 044323 (2015).[23] E. Perli´nska, S. G. Rohozi´nski, J. Dobaczewski, andW. Nazarewicz, Phys. Rev. C , 014316 (2004).[24] M. Bender, J. Dobaczewski, J. Engel, and W. Nazarewicz,Phys. Rev. C , 054322 (2002).[25] S. O. B¨ackman, O. Sj¨oberg, and A. D. Jackson, Nucl.Phys. A , 10 (1979).[26] J. Dobaczewski, W. Nazarewicz, and P. G. Reinhard, J.Phys. G: Nucl. Part. Phys. , 074001 (2014).[27] H. Primakoff and S. P. Rosen, Rep. Prog. Phys. , 121(1959).[28] H. Akimune, I. Daito, Y. Fujita, M. Fujiwara, M. B. Green-field, M. N. Harakeh, T. Inomata, J. J¨anecke, K. Katori,S. Nakayama, H. Sakai, Y. Sakemi, M. Tanaka, andM. Yosoi, Phys. Rev. C , 604 (1995).[29] B. D. Anderson, T. Chittrakarn, A. R. Baldwin, C. Lebo,R. Madey, P. C. Tandy, J. W. Watson, B. A. Brown, and C. C. Foster, Phys. Rev. C , 1161 (1985).[30] R. Madey, B. S. Flanders, B. D. Anderson, A. R. Baldwin,J. W. Watson, S. M. Austin, C. C. Foster, H. V. Klapdor,and K. Grotz, Phys. Rev. C , 540 (1989).[31] K. Pham, J. J¨anecke, D. A. Roberts, M. N. Harakeh,G. P. A. Berg, S. Chang, J. Liu, E. J. Stephenson, B. F.Davis, H. Akimune, and M. Fujiwara, Phys. Rev. C ,526 (1995).[32] T. Wakasa, H. Sakai, H. Okamura, H. Otsu, S. Fujita,S. Ishida, N. Sakamoto, T. Uesaka, Y. Satou, M. B. Green-field, and K. Hatanaka, Phys. Rev. C , 2909 (1997).[33] K. Yako, H. Sagawa, and H. Sakai, Phys. Rev. C ,051303(R) (2006).[34] T. Wakasa, M. Okamoto, M. Dozono, K. Hatanaka,M. Ichimura, S. Kuroita, Y. Maeda, H. Miyasako, T. Noro,T. Saito, Y. Sakemi, T. Yabe, and K. Yako, Phys. Rev.C , 064606 (2012).[35] “Evaluated nuclear structure data file (ENSDF),”http://nndc.bnl.gov/ensdf (2014).[36] M. V. Stoitsov, N. Schunck, M. Kortelainen, N. Michel,H. Nam, E. Olsen, J. Sarich, and S. Wild, Comp. Phys.Comm. , 1592 (2013).[37] J. Engel, M. Bender, J. Dobaczewski, W. Nazarewicz,and R. Surman, Phys. Rev. C , 014302 (1999).[38] P. G. Reinhard, D. J. Dean, W. Nazarewicz,J. Dobaczewski, J. A. Maruhn, and M. R. Strayer, Phys.Rev. C , 014316 (1999).[39] N. Schunck, J. Dobaczewski, J. McDonnell, J. Mor´e,W. Nazarewicz, J. Sarich, and M. V. Stoitsov, Phys.Rev. C , 024316 (2010).[40] T. Munson, J. Sarich, S. Wild, S. Benson, and L. C.McInnes, TAO 2.0 Users Manual , Tech. Rep. ANL/MCS-TM-322 (Mathematics and Computer Science Division,Argonne National Laboratory, 2012).[41] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither,A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D.Peterson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr,Computing in Science & Engineering , 62 (2014).[42] M. Kortelainen, R. J. Furnstahl, W. Nazarewicz, andM. V. Stoitsov, Phys. Rev. C , 011304(R) (2010).[43] E. K. Warburton, Phys. Rev. C , 233 (1991).[44] G. Lorusso, S. Nishimura, Z. Y. Xu, A. Jungclaus,Y. Shimizu, G. S. Simpson, P.-A. S¨oderstr¨om, H. Watan-abe, F. Browne, P. Doornenbal, G. Gey, H. S. Jung,B. Meyer, T. Sumikama, J. Taprogge, Z. Vajta, J. Wu,H. Baba, G. Benzoni, K. Y. Chae, F. C. L. Crespi,N. Fukuda, R. Gernh¨auser, N. Inabe, T. Isobe, T. Ka-jino, D. Kameda, G. D. Kim, Y.-K. Kim, I. Kojouharov,F. G. Kondev, T. Kubo, N. Kurz, Y. K. Kwon, G. J.Lane, Z. Li, A. Montaner-Piz´a, K. Moschner, F. Naqvi,M. Niikura, H. Nishibata, A. Odahara, R. Orlandi, Z. Pa-tel, Z. Podoly´ak, H. Sakurai, H. Schaffner, P. Schury,S. Shibagaki, K. Steiger, H. Suzuki, H. Takeda, A. Wendt,A. Yagi, and K. Yoshinaga, Phys. Rev. Lett. , 192501(2015).[45] E. Litvinova, B. A. Brown, D. L. Fang, T. Marketin, andR. G. T. Zegers, Phys. Lett. B730