# Analyzing rotational bands in odd-mass nuclei using effective field theory and Bayesian methods

^{1}Institute of Nuclear and Particle Physics and Department of Physics and Astronomy, Ohio University, Athens, OH, United States^{2}Department of Physics and Astronomy, King Saud University, Riyadh, Saudi Arabia^{3}Lawrence Livermore National Laboratory, Livermore, CA, United States

We recently developed an Effective Field Theory (EFT) for rotational bands in odd-mass nuclei. Here we use EFT expressions to perform a Bayesian analysis of data on the rotational energy levels of ^{99}Tc, ^{155,157}Gd, ^{159}Dy, ^{167,169}Er, ^{167,169}Tm, ^{183}W, ^{235}U and ^{239}Pu. The error model in our Bayesian analysis includes both experimental and EFT truncation uncertainties. It also accounts for the fact that low-energy constants (LECs) at even and odd orders are expected to have different sizes. We use Markov Chain Monte Carlo (MCMC) sampling to explore the joint posterior of the EFT and error-model parameters and show both the LECs and the breakdown scale can be reliably determined. We extract the LECs up to fourth order in the EFT and find that, provided we correctly account for EFT truncation errors in our likelihood, results for lower-order LECs are stable as we go to higher orders. LEC results are also stable with respect to the addition of higher-energy data. We extract the expansion parameter for all the nuclei listed above and find a clear correlation between the extracted and the expected value of the inverse breakdown scale, *W*, based on the single-particle and vibrational energy scales. However, the *W* that actually determines the convergence of the EFT expansion is markedly smaller than would be naively expected based on those scales.

## 1 Introduction

Rotational bands are ubiquitous in the spectra of medium-mass and heavy nuclei. As has been known for 70 years [1], they emerge in a description of the nucleus as a nearly rigid axially-symmetric rotor [2]. For even-even nuclei the simplest rotational bands consist of 0^{+}, 2^{+}, … states and their energies are described by an expansion in powers of *I*(*I* + 1), where *I* is the spin of the rotational state [3, 4]. This behavior has recently been obtained in *ab initio* calculations of the Be isotope chain [5–9] and ^{34}Mg [10].

Odd-mass neighbors of a rotor nucleus can then be understood as a fermion coupled to the rotor. The fermion dynamics is simpler in the intrinsic frame in which the nucleus is not rotating, but this frame is non-inertial, so solving the problem there induces a Coriolis force proportional to *I* appear, and produce staggering between adjacent levels. Which powers of *I* are present depends on the value of the quantum number, *K*, the projection of the fermion angular momentum on the rotor axis. For *K* = 1/2 bands the energy-level formula is:

where *A*_{K}, *E*_{K}, *A*_{1}, *B*_{1}, and *B*_{K} are parameters, related to rotor properties and single-particle matrix elements, that need to be either derived from a microscopic model or estimated from data.

Over the years a number of models have had success describing this pattern from underlying density functional theory [12–15] or shell-model [15–18] dynamics. The models also predict specific values for the coefficients that appear in Eq. 1. In Ref. [19] we took a different approach, organizing formula (1) as an effective field theory (EFT) expansion in powers of the small parameter, *Q*. For values of *I* appreciably larger than one the expansion parameter should be modified to *Q* = *I*/*I*_{{br}}, with *I*_{{br}} the spin of the nuclear state at which dynamical effects associated with single-particle and/or vibrational degrees of freedom cause the polynomial expansion in powers of *I* to break down. To simplify our later presentation we notate the inverse of the breakdown scale as *W* ≡ 1/*I*_{{br}}. We then have *Q*=*IW*. This description of rotational bands in odd-mass nuclei builds on the successful EFT developed for even-even nuclei in Refs. [3, 4]. Other efforts to develop an EFT for these rotational bands can be found in Refs. [20, 21].

In the odd-mass rotor EFT, Eq. 1 is the next-to-next-to-next-to-next-to leading order (N^{4}LO) result for the energies, and the first corrections to it are ^{k}LO energy-level formula has accuracy *I*-dependent terms in Eq. 1. In Ref. [19] we determined these LECs by fitting the lowest levels in the different rotational bands we analyzed. However, this runs the risk of fine-tuning the values of the LECs to those levels, and it does not provide uncertainty estimates for them. Better parameter estimation would use all the data available on a particular band, and account for the ^{k}LO [31, 32].

Bayesian methods for EFT parameter estimation do just that [32–35]. Reference [34] showed that the effect of neglected terms in the EFT expansion could be included in the error model by modifying the likelihood so that the covariance matrix that appears there includes both experimental uncertainties and EFT truncation errors. More recently, Ref. [35] showed that MCMC sampling of that likelihood enabled the simultaneous determination of the LECs and the parameters of the error model, i.e., the value of *W* and the typical size of the “order one” dimensionless coefficients that appear in the EFT expansion.

In this work we apply the EFT parameter estimation technology developed in Refs. [32–35] to the problem of rotational bands in odd-mass nuclei. We consider *K* = 1/2 bands in ^{99}Tc, ^{167,169}Er, ^{167,169}Tm, ^{183}W, ^{235}U and ^{239}Pu as well as *K* = 3/2 bands in ^{155,157}Gd and ^{159}Dy. Section 2 summarizes the elements of the EFT that are relevant for this paper. Section 3 then develops the Bayesian statistical model we use to analyze data on rotational bands. We first write down the likelihood that includes both experimental and theory uncertainties, and then explain how we use known information on the expected size of the LECs and the expansion parameter to set priors. A novel feature of this work, compared to earlier Bayesian EFT parameter-estimation studies, is that our statistical model incorporates the possibility that the LECs at even and odd orders have different typical sizes. This reflects the physics of odd-order LECs that are associated with matrix elements of the fermion spin, while even-order LECs contain a combination of effects from the rotor and the fermion. Section 4 contains details of our Markov Chain Monte Carlo sampler, and then Section 5 presents the results for LECs and the inverse breakdown scale, *W*, that we obtain from sampling the Bayesian posterior. We conclude in Section 6. All the results and figures generated from this work can be reproduced using publicly available Jupyter notebooks [36].

## 2 Rotational EFT background

Here we summarize the results of the EFT for rotational bands in odd-mass nuclei that was developed up to fourth order in the angular velocity of the system in Ref. [19]. This theory constructs the Lagrangian of the particle-rotor system using its angular velocity and the angular momentum of the unpaired fermion, *Q* = *W*∗*I*, according to a power-counting scheme that counts powers of the system’s angular velocity. Naively, we expect *W* to be of order *E*_{rot}/*E*_{high}, where *E*_{rot} is the energy scale at which rotational excitation take place and *E*_{high} is the scale of high-energy physics not explicitly taken into account by the EFT. At leading order (LO), the energy of a rotational band on top of a bandhead with spin *K* is

where *I* is the spin of the rotational state (or, equivalently, the total angular momentum of the fermion-rotor system), and *A*_{rot} and *E*_{K} are LECs that must be fitted to experimental data. *A*_{rot} is determined by the moment of inertia of the even-even nucleus (the rotor) to which the unpaired fermion is coupled.

At next-to-leading order (NLO) rotational bands with *K* = 1/2 are affected by a term that takes the same

where *A*_{1} is expected to be of order *A*_{rot} times a sum of matrix elements involving the fermion’s total angular momentum operator (for details see Ref. [19]). From previous studies we see that *A*_{1}/*A*_{rot} < 1. This correction, sometimes called the signature term, causes staggering between adjacent states in *K* = 1/2 bands.

The energy of a rotational band at next-to-next-to-leading order (N^{2}LO) is

The term proportional to *A*_{K} combines the LO term proportional to *A*_{rot} and corrections entering at this order with the same spin dependence. From our power counting we expect the shift Δ*A* = *A*_{rot} − *A*_{K} to be of order *A*_{rot}*W*. In contrast to *A*_{rot}, *A*_{K} is band dependent and so should be fitted to data on the rotational band of interest.

The N^{3}LO corrections to the energy of a rotational band are both ∼ *I*^{3} for *I* ≫ 1, but take a different form in the *K* = 1/2 and *K* = 3/2 bands:

with *B*_{1} and *A*_{3} expected to be of order *A*_{1}*W*^{2}. Last, at N^{4}LO we have the additional term:

with *B*_{K} expected to be of order *A*_{rot}*W*^{3}.

This pattern continues: at odd orders we add terms that correct the staggering term and have LECs of order *A*_{1}*W*^{n−1}, while the even-order terms provide the overall trend with *I* and have LECs of order *A*_{rot}*W*^{n−1}. (In both cases *n* is the order of our expansion.) This difference in the expected sizes of odd and even LECs comes from the physics. Odd-order LECs are associated with operators in the effective Lagrangian that couple rotor and fermionic degrees of freedom, while even-order LECs encode both rotor-fermion interactions and effects coming from the non-rigidity of the rotor itself.

In what follows we denote the LECs *A*_{1}, Δ*A*, *B*_{1}, and *B*_{K} generically as {*a*_{n}: *n* = 1, … , *k*}≡**a**_{k}, where *k* is the order of the EFT calculation. (In the case of *K* = 3/2 bands the set is Δ*A*, *A*_{3}, and *B*_{K}, and *a*_{1} = 0.) We then divide the *n*th-order LEC, *a*_{n}, by the reference scale and the power of the inverse breakdown scale assigned to it by the EFT power counting, i.e., construct:

We expect these coefficients *c*_{n} to be of order one, i.e., they should be natural coefficients. However, because sets of odd and even natural coefficients seem to have different sizes we will assume the even and odd *c*_{n}’s are drawn from two different distributions with different characteristic sizes that we denote by

## 3 Building the Bayesian model

### 3.1 Building the posterior

Our goal in this analysis is to use the information on the expected size of LECs to stablize the extraction of their values as we add more levels to the analysis, or as we use energy-level formulae computed at different EFT orders. At the same time, we want to estimate the inverse breakdown scale, *W*, of the theory, as well as the characteristic sizes for even and odd coefficients,

We want to obtain the posterior distribution for all the LECs that appear at order *k*, a set we collectively denote by **a**_{k}. Here we will obtain the joint posterior pdf of **a**_{k}, the inverse breakdown scale, *W*, and the characteristic sizes. To do this we follow the successful endeavor by the BUQEYE collaboration in Refs. [33–35], and write the posterior, given experimental data, *P*∗, as

Marginalization of this posterior distribution over *W*, **a**_{k}. Other marginalizations can be carried out to obtain posteriors for *W*,

This joint posterior distribution tells us the probability of the LECs and the error model parameters given experimental data. We could use this posterior distribution to get other quantities or observables, such as the energy of a particular rotational level, which depend on the LECs or the error model parameters. These are now represented by distributions and not single numbers. Their distributions are called posterior predictive distributions (PPD). We write the PPD of an observable

where

Using Bayes’ theorem, we can express the posterior of Eq. 8 as

The terms on the right-hand side of Eq. 10 have the following interpretations:

1. *k* and the parameters in our error model.

2.

3.

4.

5.

### 3.2 Building the likelihood

We now build the likelihood function accounting for the expected error between the experimental and theoretical values, for data on *K* = 1/2 rotational bands. The corresponding likelihood for *K* = 3/2 bands is built analogously. Following [34] we start by writing our observable (the energy of a particular rotational level) at order *k* as

We choose the leading-order energy for each level, *A*_{rot}*I*(*I* + 1), to be the reference scale *E*_{ref} for the observable. The dimensionless coefficients *c*_{n} (see Eq. 7) are assumed to be *k* + 1 is even, and anticorrelated for adjacent levels if *k* + 1 is odd. To account for this correlation or anticorrelation we write the theory covariance matrix as the outer product of a vector representing the theory error, *m* energy levels that enter the likelihood. We also account for experimental errors by writing the covariance matrix as

where we take

where *m* is the number of levels included in the likelihood estimation.

We note that since the theory error is the outer product of the theory error with itself, the theory covariance Σ_{th} is singular. Including the experimental error solves this singularity problem for the covariance Σ. However, Σ can still become ill-conditioned for higher values of *W* if the experimental errors are too small; numerical issues then arise when we try to invert the covariance matrix.

Including more terms in the estimate for the theoretical error produces a steeper peak in the likelihood function, see Figure 1, which, in turn, restricts the values sampled for *W* to a narrower region. Because it precludes the sampler exploring large values of *W*, this inclusion of more omitted terms in the model of the theoretical error solves the numerical problem of ill-conditioned matrices and gives a more accurate extraction of the LECs and the error-model parameters.

**FIGURE 1**. Comparing the log of the likelihood when accounting for different number of omitted terms, *p*, in the theory error. Apart from *W*, the parameters that enter the likelihood were chosen to be the median parameters after we had sampled the posterior distribution for ^{169}Er.

In what follows we estimate the theory error including omitted terms up to a certain cutoff order *k*_{max}. Our theory error estimate for the level with spin *I* is then

where the *l* and *I*-dependence of the *l*th term is chosen to match that in Eq. 11, and is denoted here by *P*_{l}(*I*), a polynomial of power *l*. We arrange the contributions to the theory error, Eq. 14 as the *p* columns of a *m* × *p* matrix *σ*_{th}, where *p* = *k*_{max} − *k* is the number of omitted terms. Each column in this matrix then corresponds to the theory-error structure, while each row corresponds to a different energy level. To obtain Σ_{th} we then again take the outer product of *σ*_{th} with itself, i.e., we construct an outer product in our *m*-dimensional data space, while also taking an inner product in order space. This results in the theory error associated with different orders being added in quadrature, while maintaining the correlation structure of the theory error across the data space.

### 3.3 Building the priors

The prior distributions for an order-*n* LEC is taken to be a Gaussian with mean zero and standard deviation

encoding the EFT expectations for the sizes of the LECs arising from the power counting described in Section 2. The standard deviation in Eq. 15 allows the possibility for even and odd LECs to have different typical sizes. Combining the Gaussian priors for the LECs yields

The LEC *E*_{K} is just an energy shift and its size is not determined by the EFT power counting. We set the prior on it to be Gaussian with mean zero and a standard deviation, *A*_{rot} by means of Eq. 3.

We choose not to impose any expectations regarding the size of the expansion parameter in the prior for *W* and so take it to be flat between two limits:

Limiting *W* from above restricts the sampler from going to high values of *W*, as they make the covariance matrix ill-conditioned and harder to invert. For all cases we check that the posterior for *W* is confined to values well below *W*_{cut}.

The priors on the characteristic sizes *χ*^{2} distributions

where the cutoff *χ*^{2} distribution, given by

is shown for different values of *ν* and *τ* in Figure 2. We stress that we chose identical priors for

**FIGURE 2**. Prior distribution of the size of the dimensionless natural coefficients,

The scaled-inverse-*χ*^{2} favors small values of *K* = 3/2 bands where we expect smaller odd-order corrections to the leading-order energy than in *K* = 1/2 bands.

## 4 Running the sampler

To sample the posterior distribution in Eq. 10 we use the Python ensemble sampling toolkit for affine-invariant MCMC (emcee) [37]. We run the sampler for each nucleus at a certain EFT order using the *m* rotational levels from the bandhead up to some *I*_{max} and accounting for *p* omitted terms in the theory error. We use 64 walkers to sample the posterior distribution for an initial 10,000 steps. We then continue running the sampler with 3,000 step increments. After every 3,000 steps we calculate the autocorrelation time, *τ*_{α}, where *α* indexes an LEC or an error-model parameter. We declare the sampler to be converged if the sampler meets two criteria. First, the number of steps has to be more than 50 times the highest *τ*_{α}. Second, the change in any of the *τ*_{α}’s has to be less than 2% from its value after the last 3,000 step increment.

To get the posterior distributions we discard 2 ×max (*τ*_{α}) steps from the beginning of the chain (*burn-in*) and 0.5 ×min (*τ*_{α}) steps in between steps we accept (*thinning*).

A sample corner plot of the marginalized distributions of the LECs and the error-model parameters *W*, ^{167}Er is shown in Figure 3. This figure clearly shows that the posterior distributions for all parameters are fully converged. For this particular case we set *W*_{cut} = 0.16 and *W* and the characteristic sizes prevent the covariance matrix from being ill-conditioned. We also ran the sampler for ^{167}Er at different values of *W*_{cut} and

**FIGURE 3**. Corner plot for the marginalized distributions of the LECs and the error-model parameters at N^{4}LO for ^{167}Er including all adopted rotational levels (*I*_{max} =16.5) and accounting for six omitted terms in the theory error. The inset in the top right corner shows the correlations between posterior parameters. The order of the parameters on the corner plot is the same on the correlations plot. (Here *E*_{K} and all the EFT LECs are expressed in keV. The error-model parameters are dimensionless.)

For some cases, namely ^{99}Tc and ^{183}W, the posterior distribution of *W* was initially at the upper limit of the prior. We then ran into numerical problems when increasing *W*_{cut} trying to encompass the entire posterior. This problem was solved by decreasing the number of levels included in the analysis, i.e., decreasing *I*_{max}. It was then possible to increase *W*_{cut} without encountering problems with degenerate matrices. This means that for ^{99}Tc we were only able to extract the LECs and *W* at *I*_{max} = 11.5. We note that this is beyond the breakdown scale for this particular nucleus and therefore we believe that the extraction of the LECs and the error model parameters is not as reliable as for the other nuclei considered in this work. For ^{183}W we needed to remove two levels from the upper end of the data set for the sampler to be numerically stable.

In Figure 3 we see clear correlations between *E*_{K}, *A*, and *B* and also between *A*_{1} and *B*_{1}. (Here we have dropped the subscript *K* on *A* and *B*; it is to be understood that all LECs are band dependent.) The correlation coefficients given in the inset in the top-right corner of the figure make the block-diagonal structure of the covariance matrix clear. To a good approximation the correlation matrix can be decomposed into a correlation matrix for even-order LECs, one for odd-order LECs, and one for the error-model parameters.

We note that, as expected, *I* are smaller than those carrying even powers of *I*. This size difference is connected to different physics correcting the effective Lagrangian at even and odd orders.

To see which of the parameters has the narrowest distribution and therefore places the strongest constraint on the posterior distribution, we did a Singular Value Decomposition (SVD) of the Hessian matrix. We found that the eigenvector with the highest eigenvalue, i.e., the parameter combination with smallest absolute error, is made up mostly of the highest-order LEC. This is unsurprising, since that LEC, *B*, is markedly smaller than the others (we note that its relative error is actually larger than that on, e.g., *A*_{1}).

We initially found a peculiar correlation between LECs in some cases where the rotational band was built on the ground state of the nucleus we were looking at. There we found the eigenvector with the highest eigenvalue was a very particular linear combination that involved all the LECs. We ultimately traced this correlation to the fact that the ground state experimental error had been set to zero, and so the combination of LECs that entered the formula for the ground-state energy was very well constrained (theory error is also very small there). This problem was solved by adding a small experimental error to the ground state. We chose it to be equal to the error that the NNDC quotes on the energy of the first excited state.

## 5 Results

In this section we show results for our Bayesian analysis of the rotational energy levels in ^{99}Tc, ^{155,157}Gd, ^{159}Dy, ^{167,169}Er, ^{167,169}Tm, ^{183}W, ^{235}U and ^{239}Pu. The experimental data are taken from the National Nuclear Data Center (NNDC) [22–30]. Except for the cases of ^{99}Tc and ^{183}W noted above, we included all levels in a certain rotational band according to the adopted level determination in the NNDC.

### 5.1 Stable LEC extraction across EFT orders and additional data

In this subsection we show that lower-order LECs extracted for the selected rotational bands are stable across EFT orders and with the addition of high-energy data, provided that we account for enough omitted terms when treating the theory error. For ^{169}Er, ^{167}Er, ^{169}Tm, and ^{239}Pu including omitted terms up to *k*_{max} = 10, i.e., accounting for six omitted terms at N^{4}LO, was enough to stabilize the extraction of the LECs.

As an example, we show the stability of the extracted LEC, *A*_{1}, across number of levels included at different EFT orders in Figure 4. In this figure, *I*_{max} is the spin of the highest-energy level included in a particular analysis. The central values of the resulting posteriors are consistent with each other within 68% credible intervals, shown as error bars in the figure. Adding more levels to the analysis narrows the posteriors for the LECs up to a certain *I*_{max}, after which the widths of these distributions saturate. Figure 4 also demonstrates striking agreement between the distributions obtained at low and high EFT orders: they are almost identical as long as omitted terms up to the same *k*_{max} are accounted for in both analyses.

**FIGURE 4**. Posteriors for *A*_{1} describing ^{169}Er as a function of *I*_{max} at different EFT orders. The solid line connects the median values and the error bands encompass the 16th and 84th percentiles of the marginalized distribution.

The importance of including more than one omitted term in the theory error estimate is evident in Figure 5. The top and bottom panels of the figure show the way that posteriors for *B*_{1} and *B* evolve as *I*_{max} increases. This is done using three error models that include different numbers of omitted terms. These results show that including more omitted terms in the model of the theory error removes the drifting and staggering of the central values.

**FIGURE 5**. Posteriors for *B*_{1} and *B* describing ^{239}Pu a function of *I*_{max} for different values of *k*_{max}. The solid line connects the median values and the error bands encompass the 16th and 84th percentiles in the marginalized distribution.

For both cases the distributions at *k*_{max} = 10 agree within errors as we go higher in *I*_{max}. The narrowing of the distribution as we go higher in *I*_{max} is clearly seen in those two figures. In addition to having less data, the broadening of the error bands at low *I*_{max} comes from the fact that including less levels in the analysis leads to highly correlated LECs. This allows the numerically larger errors on the lower-order LECs to contribute to the errors on the higher-order LECs, thereby enhancing them.

In Figure 6 we show the decrease in the correlations between the LECs as *I*_{max} increases. The high correlation between the LECs at low *I*_{max} occurs because these analyses do not include enough data to constrain all LECs independently. Furthermore, the high correlation between the LECs at low *I*_{max} also results in an unreliable extraction of the inverse breakdown scale *W*. This comes from the fact that at low energies the theory truncation error is very small compared to the experimental error. Indeed, adding more terms to our EFT error model (i.e., increasing *k*_{max}) leads to higher correlation between the LECs at low *I*_{max}. Thus, the number of levels required to reliably extract *W* increases with increasing *k*_{max}.

**FIGURE 6**. Correlations between LECs and error-model parameters as a function of *I*_{max}, resulting from the analysis on the lowest *K* =1/2 rotational band in ^{239}Pu at N^{4}LO with *k*_{max} =10.

Starting instead at the low-*I* end of the data: when we progressively remove the lowest-energy levels from the data set *D* used to construct the likelihood we rapidly lose the ability to reliably extract the LECs. Figure 7 shows that the distribution for *E*_{K} starts narrow and broadens as we remove levels from below. When we remove the six lowest energy levels the distribution of *E*_{K} is exactly the same as the prior distribution: the likelihood is making no contribution to the *E*_{K} posterior.

**FIGURE 7**. The distribution of *E*_{K} for ^{169}Er at N^{4}LO and *k*_{max} =10 as we successively remove the lowest energy levels from the data set *D*. The solid blue line connects the median values and the error bands encompass uncertainties between on the 16th and 84th percentiles of the samples in the marginalized distribution. The solid black lines show to size of the standard deviation set with the Gaussian prior on *E*_{K}.

The previous results were nearly the same for all cases considered in this work. However, even for *k*_{max} = 10, staggering and shifting of the LECs remains sizable for the *K* = 1/2 bands in ^{183}W, ^{167}Tm and ^{235}U. In ^{183}W and ^{167}Tm, these effects could be attributed to large expansion parameters, as they translate to large omitted contributions to the energies of the rotational levels. In ^{235}U, the fermionic matrix elements could be larger than naively expected, causing the systematic expansion of the EFT to be questionable as discussed in Ref. [19]. For ^{167}Tm we needed to go to *k*_{max} = 12 to get stable results, and for ^{235}U and ^{183}W we needed to go to *k*_{max} = 18.

For *K* = 3/2 bands, we were able to extract stable LECs from the ^{159}Dy analysis by setting *I*_{max} = 15.5. This extraction required us to consider omitted terms up to *k*_{max} = 16. This is because the spin at which the EFT breaks in this nucleus is *I*_{br} ≈ 15.5. (This, then, is the third case in which we do not use all the NNDC energy-level data available on a particular band.) ^{157}Gd is stable across orders and *I*_{max} and we get stable results at *k*_{max} = 12, while ^{155}Gd exhibits shifting and staggering due to a larger inverse breakdown scale, *W* ≈ 0.07, and we needed to go to *k*_{max} = 18 to get stable results.

The values of the LECs and the error-model parameters at N^{4}LO for the nuclei considered in this work are given in Tables 1, 2 respectively.

**TABLE 1**. The median value of the LECs at N^{4}LO compared with the standard deviation of their Gaussian priors with zero mean [see Eq. 15]. We see that, for nearly all cases, the LECs fall within this standard deviation. The uncertainties encompass the 16th and 84th percentiles of the samples in the marginalized distributions. *K* =3/2 rotational bands do not have a parameter *A*_{1} and the parameters (*B*_{1}, *A*_{3}) refer to *K* =1/2 and *K* =3/2 bands respectively. All the numbers are in units of keV.

**TABLE 2**. The median value of the error model parameters at N^{4}LO and the estimated expansion parameters based on rotational and single particle energy scales. The uncertainties encompass the 16th and 84th percentiles of the samples in the marginalized distributions.

### 5.2 Prior sensitivity

In addition to using the scaled-inverse-*χ*^{2} distribution as a prior for *σ* = 7 and *σ* = 3 respectively for all cases. These truncated Gaussian priors allow for smaller values of the characteristic sizes. But the standard deviations were chosen to still allow values for *χ*^{2} priors with *ν* = 1 and *τ*^{2} = 1.

The change in prior for *W* differ by less than 15%, and were similarly consistent.

The strongest dependence on the prior is that exhibited by the posteriors for

The changes in the posteriors of *W* on one hand and ^{157}Gd, where the difference is about 40%.

For all results that follow we used the scaled-inverse-*χ*^{2} distribution with *ν* = 1 and *τ*^{2} = 1 as the prior for both

### 5.3 Posterior predictive distributions

Figure 8 shows the PPDs (in blue) of the energy residuals as a function of the spin *I* for two cases considered in this work. These distributions are calculated using Eq. 9. In each figure, translucent blue lines connect energy residuals resulting from different LECs sets sampled from the posterior distribution in Eq. 8. The solid black line represents the median of the PPD, and the dashed lines encompass the region between the 16th and 84th percentiles. Meanwhile, the dark and light red bands show the truncation error and the experimental error added in quadrature at 68% and 95% credible levels respectively. To calculate the truncation error, we consider a theory error that accounts for *p* omitted terms. The omitted terms are combined in quadrature, just as they are in the likelihood defined in Section 3. This calculation was done using Eq. 9 i.e., by calculating the theory error at each point in the sample space and then marginalizing over the error parameters. The dependence of the size of the theory error on the prior on

**FIGURE 8**. The posterior predictive distribution for energy-level residuals at N^{4}LO and *k*_{max} =10 for ^{169}Er and *I*_{max} =17.5 (top panel) and at *k*_{max} =18 for ^{235}U and *I*_{max} =23.5 (bottom panel). The dark and light red bands show the truncation error plus the experimental error at 68% and 95% credible levels respectively. The lighter blue lines connect the energy residuals calculated from the distribution of the LECs. The solid black line represents the median of the distribution and the dashed lines indicate the 16th and 84th percentiles. The correlation shown on the plot is the highest correlation between any LEC and any error-model parameter. *I*_{br} was determined from the distribution of *W*. The dashed purple line shows the lower limit of *I*_{br}. The inset on the plot shows the residuals on the first five levels with an altered *y*-axis scale.

The correlation coefficient written in the legend in Figure 8 is the largest between any LEC and any error-model parameter for the shown analysis. When this value is small, the truncation error and the propagated LEC error could in principle be added together in quadrature.

In viewing Figure 8 it is important to remember that the truncation error on the energy residuals is highly correlated across levels. This comes from the high correlation between levels when building the correlation matrix that goes into the likelihood. This correlation also flows into a correlation between levels in the PPD of the energies. A correlation plot between two energy levels, like the ones in Figure 9, gives a 2D cut of this multi-dimensional correlation.

**FIGURE 9**. A 2D cut of the posterior predictive distribution at N^{4}LO and *k*_{max} =10 for ^{169}Er and *I*_{max} =17.5 (top panel) and at *k*_{max} =18 for ^{235}U and *I*_{max} =23.5 (bottom panel). The blue dots show the energies calculated from the distribution of the LECs. The black cross shows the experimental value and the black lines and black ellipse shows the corresponding experimental uncertainty. The remaining ellipses and lines show the truncation error and the experimental error added in quadrature. (All the ellipses are centered at the experimental value.) The orange, red and green account for 1, two and six omitted terms in the theory error respectively (In the top panel the red ellipse is completely covered by the green ellipse.)

In both panels we see the importance of accounting for more than one omitted term in the theory error. This is clearly shown in the reverse in the direction of the correlation from a negative to a positive correlation when going from the orange ellipse to the red ellipse. The orange ellipse is obtained when we account for only one omitted term, while the red ellipse includes the effect of two omitted terms. After accounting for six omitted terms the green ellipse is obtained and the 68% ellipse in principle expands. This is more clearly seen when we go to high-energy levels plotted in the lower panel in Figure 9. Note also that for lower-energy levels the correlation is smaller since the experimental error dominates over the truncation error, and we assumed that the experimental errors are not correlated across energy levels.

### 5.4 Model checking

In Figures 10, 11 we compare the marginalized posterior distributions of the LECs, on the *y*-axis, with their expected sizes from the EFT power counting, on the *x*-axis. Since we also extract the theory error parameters from the sampler and they are highly correlated among themselves, we calculate the expected size from the distributions of the error model parameters using Eq. 9. We notice that the error on the distribution of the LECs is very small compared to the error on the expected sizes that comes from the distribution of the theory error parameters.

**FIGURE 10**. The size of the NLO LEC, *A*_{1} (top panel) and the N^{3}LO LEC, *B*_{1} for *K* =1/2 bands and *A*_{3} for *K* =3/2 bands (bottom panel), on the *y*-axis, compared to its expected size from the EFT power counting, on the *x*-axis. Error bands on the LEC distribution are small and can not been seen on the plot. The error bands on the *x*-axis encompass the 16th and 84th percentiles. Different nuclei are labeled in the legend of the plot. The black dashed line has slope =1 and is plotted to facilitate comparison of prior expectation and results from the posterior. The yellow colored symbols are results for rotational bands with bandheads *K* =3/2, all the others are *K* =1/2 bands. *K* =3/2 rotational bands do not have a parameter *A*_{1} and we do not have them in the top panel. ^{99}Tc and ^{155}Gd are outliers and we exclude them from the plots (LECs values for these nuclei can be found in Table 1).

**FIGURE 11**. The size of the N^{2}LO LEC, Δ*A* (top panel) and the N^{4}LO LEC, *B* (bottom panel), on the *y*-axis, compared to its expected size from the EFT power counting, on the *x*-axis. Error bands on the LEC distribution are small and can not been seen on the plot. The error bands on the *x*-axis encompass uncertainties between the 16th and 84th percentiles. Different nuclei are labeled in the legend of the plot. The black dashed line has slope =1 and is plotted to facilitate comparison of prior expectation with results from the posterior. The yellow colored symbols are results for rotational bands with bandheads *K* =3/2, all the others are *K* =1/2 bands. ^{99}Tc and ^{183}W are outliers and we exclude them from the both plots. We also exclude ^{155}Gd from the bottom plot only (LEC values for these nuclei can be found in Table 1).

As these graphs are model-checking graphs, and since the estimates of LEC sizes plotted on the *x*-axis are meant as order-of-magnitude estimates, we do not expect perfect linear correlations. Nevertheless, the top panel in Figure 11 shows that, for all *K* = 1/2 bands considered, the LEC Δ*A* agrees with its expected size within error bands. This result is surprisingly better than expected. In contrast, the size of Δ*A* for *K* = 3/2 bands is larger than expected, especially for ^{155}Gd (see yellow symbols in Figure 11). There are two factors that could contribute to this. First, the *K* = 3/2 bands have larger fermionic matrix elements. This could hinder the systematic expansion of the EFT. Second, the *K* = 3/2 bands have relatively larger expansion parameters, see Figure 12.

**FIGURE 12**. The extracted inverse breakdown scale *W* from the marginalized posterior distribution obtained by sampling compared to its naively expected size. The expected size is taken to be the maximum of *E*_{rot}/*E*_{sp} and *E*_{rot}/*E*_{vib}. The dashed black line shows the best linear fit and its parameters are printed on the plot. The yellow colored symbols are results for rotational bands with bandheads *K* =3/2, all the others are *K* =1/2 bands. ^{99}Tc is an outlier and we exclude it from the plot (its values can be seen in Table 2).

The same discussion applies to the results in the remaining panels in Figures 10, 11, where we see good agreement between the LECs and their expected sizes for *K* = 1/2 bands. The disagreement with power-counting estimates for *K* = 3/2 bands at N^{3,4}LO is less of a concern than the one at N^{2}LO seen in the top panel of Figure 11, since these higher-order LECs are smaller than their expected sizes. This doesn’t undermine the convergence of the EFT expansion.

We also note here that the scale of the *x*-axis is prior dependent and could change by more than 50% in some nuclei, depending on the choice of prior on ^{157}Gd changing the prior on *A*_{3} was then equal to the expected size (i.e., the point for ^{157}Gd then falls exactly on the line in the bottom panel of Figure 10). This did not happen when the truncated normal is chosen as a prior for the analysis in ^{155}Gd and ^{159}Dy; this may occur because there is strong N^{5}LO energy-level staggering present in the data for these nuclei.

The size of *K* = 3/2 bands due to large fermionic matrix elements, this tension becomes clear. One example of this is seen in Figure 11 for *K* = 3/2 bands. There Δ*A* is large and favors large values of *B* is small and favors smaller values of

### 5.5 Higher than expected break-down scale

In Figure 12 we see a clear correlation between the extracted values of *W* and those that are expected based on each nucleus’ single-particle and vibrational energy scales, *E*_{sp} and *E*_{vib}. The expected *W* is the larger of *E*_{rot}/*E*_{sp} and *E*_{rot}/*E*_{vib}, while the extracted *W* comes from sampling the posterior in Eq. 10. This extracted *W* is what actually determines the convergence of the EFT expansion. It is markedly smaller than would be naively expected. The break-down scale of the theory is thus higher than naively expected: our rotational EFT works to much higher *I* than energy-scale arguments would suggest. This could occur because coupling between the higher rotational states explicitly included in the EFT and the high-energy states that are not explicitly included in our EFT is hindered by the large difference in angular momentum between them.

## 6 Conclusion

We performed a Bayesian analysis to extract the LECs and inverse breakdown scale *W* describing the rotational energy levels of diverse odd-mass nuclei within a recently developed EFT. This analysis corroborates the EFT organization for energy-level formulae which results from the assumed power-counting scheme: the extracted LECs of order *k* scale as *W*^{k−1}, i.e., according to EFT expectations. While our analysis reached this conclusion for both *K* = 1/2 and *K* = 3/2 rotational bands, the sizes of the LECs describing the latter exhibit larger deviations from their expected values than those describing the former. We attribute this behavior to the size of fermionic matrix elements, assumed to be of order one while organizing energy-level formulae. Since these matrix elements involve the angular momentum of the fermion, *K*. For the *K* = 3/2 bands studied in this work, however, this discrepancy does not destroy the systematic improvement of calculated energies up to N^{4}LO, as the sizes of extracted LECs are smaller that expected.

In order to ensure that the extracted values are independent of the EFT order and number of energy levels entering the analysis, we employed a theory error beyond the first-omitted-term approximation, considering omitted terms in the expansion for the energy of rotational levels up to order *k*_{max}. As we increased the number of omitted terms considered in the theory error, the corresponding log likelihood exhibited steeper and steeper peaks. Therefore, the ‘widths’ of the sampled posteriors decrease as *k*_{max} increases. Considering up to fourteen omitted terms at N^{4}LO enabled a stable extraction of the LECs and breakdown scale describing the levels of interest. The shapes of posteriors for low-order LECs extracted at this order and those extracted using lower-order energy formulae are, for all practical purposes, identical. On the other hand, the shapes of the posteriors depend strongly on the number of levels informing the model, narrowing as more levels are included. Nevertheless, the 68% credible intervals of these posteriors possess significant overlap, facilitating reliable LEC extraction.

In addition to the posteriors for the LECs and the inverse breakdown scales, our analysis yielded distributions for the characteristic sizes of even and odd *c*_{n}’s,

Although the distributions of

These considerations mean that our extractions of the LECs and the theory error parameters in the EFT of rotational bands in odd-mass nuclei are robust under the choice of prior. The formalism presented here also gives robust results for LECs across orders and as more data is added to the analysis. We conclude that a Bayesian framework that incorporates theory errors in the likelihood offers significant advantages for LEC extraction in EFTs. This methodology has already been used for the extraction of LECs in the NN potential from phase shifts [34] and to constrain parameters of the three-nucleon force [35]. But it is a very general approach which should improve the parameter estimation for LECs in any EFT.

## Data availability statement

The datasets presented in this study can be found in the online repository: https://github.com/inamlah/brb.

## Author contributions

This work has been conceptualized by DP and IA. Both IA and EACP have contributed to writing the code. IA, EACP and DP have all contributed to analyzing the results and writing the manuscript.

## Funding

This work was supported by the US Department of Energy, contract DE-FG02-93ER40756 (IKA, DRP) by the National Science Foundation CSSI program Award OAC-2004601 (DRP), and under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 (EACP). IKA acknowledges the support of King Saud University and the Ministry of Education in Saudi Arabia.

## Acknowledgments

We thank Dick Furnstahl and Jordan Melendez for useful discussions. We are also grateful for Daniel Odell’s significant conceptual and computational assistance.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor CF declared a past co-authorship with the author DRP.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2022.901954/full#supplementary-material

## References

1. Bohr A. On the quantization of angular momenta in heavy nuclei. *Phys Rev* (1951) 81:134–8. doi:10.1103/PhysRev.81.134

2. Rowe DJ. *Nuclear collective motion: Models and theory*. Singapore: World Scientific (2010). Chapter 6.

3. Papenbrock T. Effective theory for deformed nuclei. *Nucl Phys A* (2011) 852:36–60. [arXiv:1011.5026 [nucl-th]]. doi:10.1016/j.nuclphysa.2010.12.013

4. Coello Pérez EA, Papenbrock T. Effective theory for the nonrigid rotor in an electromagnetic field: Toward accurate and precise calculations ofE2transitions in deformed nuclei. *Phys Rev C* (2015) 92(no. 1):014323. [arXiv:1502.04405 [nucl-th]]. doi:10.1103/PhysRevC.92.014323

5. Caprio MA, Maris P, Vary JP. Emergence of rotational bands in *ab initio* no-core configuration interaction calculations of light nuclei. *Phys Lett B* (2013) 719:179–84. [arXiv:1301.0956 [nucl-th]]. doi:10.1016/j.physletb.2012.12.064

6. Maris P, Caprio MA, Vary JP. Erratum: Emergence of rotational bands in *ab initio* no-core configuration interaction calculations of the Be isotopes. *Phys Rev Cphys Rev C* (2015) 91(no.2):014310. 029902. [erratum:[arXiv:1409.0881 [nucl-th]]. 10.1103/physrevc.99.029902. doi:10.1103/PhysRevC.91.014310

7. Jansen GR, Signoracci A, Hagen G, Navrátil P. Open s d-shell nuclei from first principles. *Phys Rev C* (2016) 94(no.1):011301. [arXiv:1511.00757 [nucl-th]]. doi:10.1103/PhysRevC.94.011301

8. Caprio MA, Fasano PJ, Maris P, McCoy AE, Vary JP. Probing *ab initio* emergence of nuclear rotation. *Eur Phys J A* (2020) 56(no.4):120. [nucl-th]]. doi:10.1140/epja/s10050-020-00112-0

9. McCoy AE, Caprio MA, Dytrych T, Fasano PJ. Emergent sp(3, R) dynamical symmetry in the nuclear many-body system from an *ab initio* description. *Phys Rev Lett* (2020) 125(no.10):102505. [nucl-th]]. doi:10.1103/physrevlett.125.102505

10. Hagen G, Novario SJ, Sun ZH, Papenbrock T, Jansen GR, Lietz JG, et al. Angular-momentum projection in coupled-cluster theory: Structure of Mg34. *Phys Rev C* 105(6). [arXiv:2201.07298 [nucl-th]]. doi:10.1103/physrevc.105.064311

11. Bohr A, Mottelson B. *Nuclear structure, volume II: Deformations*. Singapore: World Scientific (1998).

12. Dudek J, Nazarewicz W, Szymanśki J. Independent quasiparticle analysis of rotational bands in^{156}Er. *Phys Scr* (1981) 24:309–11. doi:10.1088/0031-8949/24/1b/029

13. Cwiok S, Nazarewicz W, Dudek J, Szymanski Z. Analysis of the backbending effect inYb166, Yb168, andYb170within the Hartree-Fock-Bogolyubov cranking method. *Phys Rev C* (1980) 21:448–52. doi:10.1103/PhysRevC.21.448

14. Afanasjev AV, Abdurazakov O. Pairing and rotational properties of actinides and superheavy nuclei in covariant density functional theory. *Phys Rev C* (2013) 88(no.1):014320. [arXiv:1307.4131 [nucl-th]]. doi:10.1103/PhysRevC.88.014320

15. Zhang ZH, Huang M, Afanasjev AV. Rotational excitations in rare-earth nuclei: A comparative study within three cranking models with different mean fields and treatments of pairing correlations. *Phys Rev C* (2020) 101(no.5):054303. arXiv:2003.07902 [nucl-th]]. doi:10.1103/PhysRevC.101.054303

16. Inglis DR. Particle derivation of nuclear rotation properties associated with a surface wave. *Phys Rev* (1954) 96:1059–65. doi:10.1103/PhysRev.96.1059

17. Velazquez V, Hirsch JG, Sun Y, Guidry MW. Backbending in Dy isotopes within the projected shell model. *Nucl Phys A* (1999) 653:355–71. [arXiv:nucl-th/9901044 [nucl-th]]. doi:10.1016/S0375-9474(99)00238-9

18. Liu SX, Zeng JY, Yu L. Particle-number-conserving treatment for the backbending in Yb isotopes. *Nucl Phys A* (2004) 735:77–85. doi:10.1016/j.nuclphysa.2004.02.007

19. Alnamlah IK, Pérez EAC, Phillips DR. Effective field theory approach to rotational bands in odd-mass nuclei. *Phys Rev C* (2021) 104(no.6):064311. [arXiv:2011.01083 [nucl-th]]. doi:10.1103/PhysRevC.104.064311

20. Papenbrock T, Weidenmüller HA. Effective field theory for deformed odd-mass nuclei. *Phys Rev C* (2020) 102(no.4):044324. [nucl-th]. doi:10.1103/physrevc.102.044324

21. Chen QB, Kaiser N, Meißner UG, Meng J. Effective field theory for triaxially deformed odd-mass nuclei. [arXiv:2003.04065 [nucl-th]].

22. Baglin CM. Nuclear data sheets for A = 169. *Nucl Data Sheets* (2008) 109(no.9):2033–256. doi:10.1016/j.nds.2008.08.001

23. Baglin CM. Nuclear data sheets for A = 167. *Nucl Data Sheets* (2000) 90(3):431–644. doi:10.1006/ndsh.2000.0012

24. Browne E, Tuli JK. Nuclear data sheets for A = 239. *Nucl Data Sheets* (2014) 122:293–376. doi:10.1016/j.nds.2014.11.003

25. Browne E, Tuli JK. Nuclear data sheets for A = 235. *Nucl Data Sheets* (2014) 122:205–92. doi:10.1016/j.nds.2014.11.002

26. Reich CW. Nuclear data sheets for A = 159. *Nucl Data Sheets* (2012) 113(no.1):157–363. doi:10.1016/j.nds.2012.01.002

27. Nica N. Nuclear data sheets for A=155. *Nucl Data Sheets* (2019) 160:1–404. doi:10.1016/j.nds.2019.100523

28. Nica N. Nuclear data sheets for A = 157. *Nucl Data Sheets* (2016) 132:1–256. doi:10.1016/j.nds.2016.01.001

29. Browne E, Tuli JK. Nuclear data sheets for A = 99. *Nucl Data Sheets* (2017) 145:25–340. doi:10.1016/j.nds.2017.09.002

30. Baglin CM. Nuclear data sheets for A = 183. *Nucl Data Sheets* (2016) 134:149–430. doi:10.1016/j.nds.2016.04.002

31. Furnstahl RJ, Phillips DR, Wesolowski S. A recipe for EFT uncertainty quantification in nuclear physics. *J Phys G: Nucl Part Phys* (2015) 42(no.3):034028. [nucl-th]]. doi:10.1088/0954-3899/42/3/034028

32. Schindler MR, Phillips DR. Bayesian methods for parameter estimation in effective field theories. *Ann Phys (N Y)* (2009) 324:682–708. [arXiv:0808.3643 [hep-ph]]. doi:10.1016/j.aop.2008.09.003

33. Wesolowski S, Klco N, Furnstahl RJ, Phillips DR, Thapaliya A. Bayesian parameter estimation for effective field theories. *J Phys G: Nucl Part Phys* (2016) 43(7):074001. [nucl-th]]. doi:10.1088/0954-3899/43/7/074001

34. Wesolowski S, Furnstahl RJ, Melendez JA, Phillips DR. Exploring Bayesian parameter estimation for chiral effective field theory using nucleon–nucleon phase shifts. *J Phys G: Nucl Part Phys* (2019) 46(no.4):045102. [arXiv:1808.08211 [nucl-th]]. doi:10.1088/1361-6471/aaf5fc

35. Wesolowski S, Svensson I, Ekström A, Forssén C, Furnstahl RJ, Melendez JA, et al. Rigorous constraints on three-nucleon forces in chiral effective field theory from fast and accurate calculations of few-body observables. *Phys Rev C* (2021) 104(no.6):064001. arXiv:2104.04441 [nucl-th]]. doi:10.1103/PhysRevC.104.064001

36. Alnamlah IK, Coello Pérez EA, Phillips DR. Bayesian rotational bands. *GitHub repository* (2022). Available at: http://github.com/inamlah/brb.

Keywords: EFT, bayesian analysis, rotational bands, collective models, nuclear structure

Citation: Alnamlah IK, Coello Pérez EA and Phillips DR (2022) Analyzing rotational bands in odd-mass nuclei using effective field theory and Bayesian methods. *Front. Phys.* 10:901954. doi: 10.3389/fphy.2022.901954

Received: 22 March 2022; Accepted: 05 July 2022;

Published: 22 August 2022.

Edited by:

Christian Forssén, Chalmers University of Technology, SwedenReviewed by:

Norbert Kaiser, Technical University of Munich, GermanyChristopher Koerber, Ruhr University Bochum, Germany

Copyright © 2022 Alnamlah, Coello Pérez and Phillips. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: I. K. Alnamlah, ialnamlah@ksu.edu.sa