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

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 ³⁴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 $\vec{j}\cdot \vec{I}$, the dot product of the single-fermion angular momentum and the total angular momentum of the fermion-rotor system. When combined with other mechanisms, such as excitation of the fermion to higher-single particle states and the fermion disturbing the rotor, this induces a string of terms in the energy-level formula [11]. Odd powers of 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:

$\begin{array}{ll}\hfill E\left(I\right)& =A_{K}I\left(I+1\right)+E_K+A_1{\left(-1\right)}^{I+1/2}\left(I+1/2\right)\hfill \\ \hfill & +B_{1}I\left(I+1\right){\left(-1\right)}^{I+1/2}\left(I+\frac{1}{2}\right)+B_K{\left[I\left(I+1\right)\right]}^2\hfill \end{array}(1)$

where A_K, E_K, A₁, B₁, 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⁴LO) result for the energies, and the first corrections to it are $\mathcal{O}(E{Q}^4)$. The EFT analysis of Eq. 1 organizes it in terms of increasingly accurate predictions: the N^kLO energy-level formula has accuracy $\mathcal{O}(E{Q}^k)$. All short-distance/high-energy physical mechanisms that affect the energies up to that accuracy are subsumed into the parameters or low-energy constants (LECs) that multiply the 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 $\mathcal{O}(E{Q}^k)$ truncation uncertainty present at order N^kLO [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 ⁹⁹Tc, ^167,169Er, ^167,169Tm, ¹⁸³W, ²³⁵U and ²³⁹Pu as well as K = 3/2 bands in ^155,157Gd and ¹⁵⁹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, $\vec{j}$, as building blocks. The resulting Lagrangian corrects that of a rigid rotor with contributions arranged as a series in powers of a small expansion parameter, 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

$E_{\text{LO}}\left(I,K\right)=A_{\text{rot}}I\left(I+1\right)+E_K,(2)$

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 $\vec{j}\cdot \vec{I}$ form as the Coriolis force. This produces:

$\begin{array}{cc}\hfill E_{\text{NLO}}\left(I,K\right)& =A_{\text{rot}}I\left(I+1\right)+E_K\hfill \\ \hfill & +A_1{\left(-1\right)}^{I+1/2}\left(I+\frac{1}{2}\right){\delta }_{1/2}^K,\hfill \end{array}(3)$

where ${\delta }_{K^{\prime }}^K$ is the Kronecker delta. The LEC A₁ 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₁/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²LO) is

$\begin{array}{cc}\hfill E_{{\mathrm{N}}^2\mathrm{L}\mathrm{O}}\left(I,K\right)& =A_{K}I\left(I+1\right)+E_K\hfill \\ \hfill & +A_1{\left(-1\right)}^{I+1/2}\left(I+\frac{1}{2}\right){\delta }_{1/2}^K.\hfill \end{array}(4)$

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_rotW. 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³LO corrections to the energy of a rotational band are both ∼ I³ for I ≫ 1, but take a different form in the K = 1/2 and K = 3/2 bands:

$\begin{array}{cc}\hfill \mathrm{\Delta }{E}_{{\mathrm{N}}^3\mathrm{L}\mathrm{O}}\left(I,K\right)& =B_1{\left(-1\right)}^{I+1/2}\left(I+\frac{1}{2}\right)I\left(I+1\right){\delta }_{1/2}^K\hfill \\ \hfill & +A_3{\left(-1\right)}^{I+3/2}\left(I+\frac{1}{2}\right)\left(I-\frac{1}{2}\right)\left(I+\frac{3}{2}\right){\delta }_{3/2}^K.\hfill \end{array}(5)$

with B₁ and A₃ expected to be of order A₁W². Last, at N⁴LO we have the additional term:

$\mathrm{\Delta }{E}_{{\mathrm{N}}^4\mathrm{L}\mathrm{O}}\left(I,K\right)=B_K{\left[I\left(I+1\right)\right]}^2.(6)$

with B_K expected to be of order A_rotW³.

This pattern continues: at odd orders we add terms that correct the staggering term and have LECs of order A₁Wⁿ⁻¹, while the even-order terms provide the overall trend with I and have LECs of order A_rotWⁿ⁻¹. (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₁, ΔA, B₁, 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₃, and B_K, and a₁ = 0.) We then divide the nth-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:

$c_n=\frac{a_n}{A_{\text{rot}}{W}^{n-1}}.(7)$

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 ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$.

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, ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$.

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, ${\vec{y}}_{\mathit{exp}}$, and prior information on the model, P∗, as

$\begin{array}{cc}\hfill \mathrm{pr}\left({\mathbf{a}}_k,W,{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}}|{\vec{y}}_{\mathit{exp}},P\ast \right)& =\mathrm{pr}\left({\mathbf{a}}_k|W,{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},{\vec{y}}_{\mathit{exp}},P\ast \right)\hfill \\ \hfill & \times \mathrm{pr}\left(W|{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},{\vec{y}}_{\mathit{exp}},P\ast \right)\hfill \\ \hfill & \times \mathrm{pr}\left({\bar{c}}_{\mathit{even}}|{\bar{c}}_{\mathit{odd}},{\vec{y}}_{\mathit{exp}},P\ast \right)\hfill \\ \hfill & \times \mathrm{pr}\left({\bar{c}}_{\mathit{odd}}|{\vec{y}}_{\mathit{exp}},P\ast \right).\hfill \end{array}(8)$

Marginalization of this posterior distribution over W, ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$ yields the posterior distribution for a_k. Other marginalizations can be carried out to obtain posteriors for W, ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$.

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 $\mathcal{O}$ as

$\mathrm{pr}\left(\mathcal{O}|{\vec{y}}_{\mathit{exp}},P\ast \right)=\int d\vec{\theta }\delta \left(\mathcal{O}-\mathcal{O}\left(\vec{\theta }\right)\right)\mathrm{pr}\left(\vec{\theta }|{\vec{y}}_{\mathit{exp}},P\ast \right)(9)$

where $\vec{\theta }$ represents the LECs and the error model parameters. Calculating the observable at each point in the parameter space $\vec{\theta }$ and then integrating over the parameters $\vec{\theta }$ allows one to carefully account for correlations between the parameters.

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

$\begin{array}{cc}\hfill \mathrm{pr}\left({\mathbf{a}}_k,W,{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}}|{\vec{y}}_{\mathit{exp}},P\ast \right)& =\mathrm{pr}\left({\vec{y}}_{\mathit{exp}}|{\mathbf{a}}_k,W,{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},P\ast \right)\hfill \\ \hfill & \times \mathrm{pr}\left({\mathbf{a}}_k|W,{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},P\ast \right)\hfill \\ \hfill & \times \mathrm{pr}\left(W|{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},P\ast \right)\hfill \\ \hfill & \times \mathrm{pr}\left({\bar{c}}_{\mathit{even}}|P\ast \right)\mathrm{pr}\left({\bar{c}}_{\mathit{odd}}|P\ast \right)\hfill \\ \hfill & \times \frac{1}{\mathrm{pr}\left({\vec{y}}_{\mathit{exp}}|P\ast \right)}.\hfill \end{array}(10)$

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

1. $\mathrm{pr}({\vec{y}}_{\mathit{exp}}|{\mathbf{a}}_k,W,{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},P\ast )$ is the likelihood of the experimental data given specific values of both the LECs that appear in the energy formula at order k and the parameters in our error model.

2. $\mathrm{pr}({\mathbf{a}}_k|W,{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},P\ast )$ is the prior distribution of the LECs given the parameters encoding the systematic expansion of the EFT.

3. $\mathrm{pr}(W|{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},P\ast )$ is the prior distribution of the inverse breakdown scale given the characteristic sizes of even and odd natural coefficients.

4. $\mathrm{pr}({\bar{c}}_{\mathit{even}}|P\ast )$ and $\mathrm{pr}({\bar{c}}_{\mathit{odd}}|P\ast )$ are the prior distributions of the even and odd characteristic sizes (In Eq. 10 we assume an uncorrelated prior on ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$.).

5. $\mathrm{pr}({\vec{y}}_{\mathit{exp}}|P\ast )$ is the evidence, which we drop in what follows as it does not depend on the parameters we are interested in extracting and functions only as a normalization constant.

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

$\begin{array}{c}E\left(I\right)=A_{\text{rot}}I\left(I+1\right)\left\{1+\sum _{n=odd}^k{c}_n{W}^{n-1}{\left(-1\right)}^{I+1/2}\left(I+\frac{1}{2}\right){\left[I\left(I+1\right)\right]}^{\left(n-3\right)/2}+\sum _{n=even}^k{c}_n{W}^{n-1}{\left[I\left(I+1\right)\right]}^{\left(n-2\right)/2}\right\}.\hfill \end{array}(11)$

We choose the leading-order energy for each level, A_rotI(I + 1), to be the reference scale E_ref for the observable. The dimensionless coefficients c_n (see Eq. 7) are assumed to be $\mathcal{O}(1)$. The theory error ${\vec{\sigma }}_{th}$ at any order is due to terms omitted from the summations in Eq. 11. Its most significant contribution comes from the first omitted term in the EFT expansion. Accounting only for this term yields an estimate for the theory error that is fully correlated across levels if 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, ${\mathrm{\Sigma }}_{th}\equiv {\vec{\sigma }}_{th}\otimes {\vec{\sigma }}_{th}$. The vector ${\vec{\sigma }}_{th}$ contains the value of the first omitted term for each of the m energy levels that enter the likelihood. We also account for experimental errors by writing the covariance matrix as

$\mathrm{\Sigma }={\mathrm{\Sigma }}_{th}+{\mathrm{\Sigma }}_{\mathit{exp}}(12)$

where we take ${({\mathrm{\Sigma }}_{\mathit{exp}})}_{ij}\equiv {({\vec{\sigma }}_{\mathit{exp}})}_i^2{\delta }_{ij}$. The likelihood function is then

$\begin{array}{c}\hfill \mathrm{pr}\left({\vec{y}}_{\mathit{exp}}|{\mathbf{a}}_k,W,{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},P\ast \right)=\sqrt{\frac{1}{{\left(2\pi \right)}^m|\mathrm{\Sigma }|}}\exp \left(-\frac{1}{2}{\vec{r}}^T{\mathrm{\Sigma }}^{-1}\vec{r}\right),\end{array}(13)$

where $\vec{r}\equiv {\vec{y}}_{\mathit{exp}}-{\vec{y}}_{th}$ is the residual between the central experimental energy for a level and the theory result from Eq. 11 and 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

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 ¹⁶⁹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

${\sigma }_{th}\left(I\right)=A_{\text{rot}}\sum _{l=k+1}^{k_{\mathit{max}}}{\bar{c}}_{\mathit{even,odd}}{W}^{l-1}{P}_l\left(I\right),(14)$

where the $\bar{c}$ that is used here is ${\bar{c}}_{\mathit{even}}$ for even values of l and ${\bar{c}}_{\mathit{odd}}$ otherwise. The I-dependence of the lth 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

${\sigma }_n=\left\{\begin{array}{cc}{A}_{\text{rot}}{\bar{c}}_{\mathit{even}}{W}^{n-1}\hfill & \hfill \mathrm{i}\mathrm{f}n\mathrm{i}\mathrm{s}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n};\hfill \\ A_{\text{rot}}{\bar{c}}_{\mathit{odd}}{W}^{n-1}\hfill & \hfill \mathrm{i}\mathrm{f}n\mathrm{i}\mathrm{s}\mathrm{o}\mathrm{d}\mathrm{d},\hfill \end{array}\right.(15)$

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

$\begin{array}{c}\hfill \mathrm{pr}\left({\mathbf{a}}_k|W,{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},P\ast \right)=\frac{1}{\bar{E}\sqrt{2\pi }}\mathit{exp}\left(-\frac{E_k^2}{2{\bar{E}}^2}\right)\prod _{n=1}^k\frac{1}{{\sigma }_n\sqrt{2\pi }}\exp \left(-\frac{a_n^2}{2{\sigma }_n^2}\right).\end{array}(16)$

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, $\bar{E}$, that is wide enough to capture its value. The value for $\bar{E}$ is determined from the energy of the bandhead and 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:

$\mathrm{pr}\left(W|{\bar{c}}_{\mathit{even}},{\bar{c}}_{\mathit{odd}},P\ast \right)\propto \left\{\begin{array}{cc}1\hfill & W\in \left(0,W_{\mathit{cut}}\right)\hfill \\ 0\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{array}\right.(17)$

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 ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$, are taken to be identical scaled-inverse-χ² distributions

$\mathrm{pr}\left({\bar{c}}_l^2|P\ast \right)\propto \left\{\begin{array}{cc}{\chi }^{-2}\left(\nu =1,{\tau }^2=1\right)\hfill & {\bar{c}}_l^2\in \left(0,{\bar{c}}_{\mathit{cut}}^2\right)\hfill \\ 0\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\hfill \end{array}\right.(18)$

where the cutoff ${\bar{c}}_{\mathit{cut}}$ prevents numerical issues inverting the covariance matrix. The scaled-inverse-χ² distribution, given by

${\chi }^{-2}\left(x;\nu ,{\tau }^2\right)=\frac{{\left({\tau }^2\nu /2\right)}^{\nu /2}}{\mathrm{\Gamma }\left(\nu /2\right)}\frac{\exp \left[-\frac{\nu {\tau }^2}{2x}\right]}{x^{1+\nu /2}},(19)$

is shown for different values of ν and τ in Figure 2. We stress that we chose identical priors for ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$ even though we expect the former to be larger than the latter based on previous analyses of data on rotational bands [19]. We did not want to bias our analysis by imposing this hierarchy on the prior, instead anticipating that it will emerge naturally in the posteriors for those parameters.

FIGURE 2

FIGURE 2. Prior distribution of the size of the dimensionless natural coefficients, $\bar{c}$.

The scaled-inverse-χ² favors small values of ${\bar{c}}^2$ and has long tails. This allows the sampler to explore higher values of ${\bar{c}}^2$. The sharp decrease in this distribution for very small values of ${\bar{c}}^2$ could be a problem for cases where ${\bar{c}}_{\mathit{odd}}$ is much smaller than one. This is a concern in some 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, ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$, for the case of ¹⁶⁷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 ${\bar{c}}_{\mathit{cut}}=22$ for both ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$. As explained in Section 3.3, the cutoffs on W and the characteristic sizes prevent the covariance matrix from being ill-conditioned. We also ran the sampler for ¹⁶⁷Er at different values of W_cut and ${\bar{c}}_{\mathit{cut}}$ and found that different choices of these hyperparameters do not result in a significant change in the posterior distributions.

FIGURE 3

FIGURE 3. Corner plot for the marginalized distributions of the LECs and the error-model parameters at N⁴LO for ¹⁶⁷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 ⁹⁹Tc and ¹⁸³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 ⁹⁹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 ¹⁸³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₁ and B₁. (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, ${\bar{c}}_{\mathit{odd}}$ is smaller than ${\bar{c}}_{\mathit{even}}$. Corrections to the energy levels carrying odd powers of 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₁).

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 ⁹⁹Tc, ^155,157Gd, ¹⁵⁹Dy, ^167,169Er, ^167,169Tm, ¹⁸³W, ²³⁵U and ²³⁹Pu. The experimental data are taken from the National Nuclear Data Center (NNDC) [22–30]. Except for the cases of ⁹⁹Tc and ¹⁸³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 ¹⁶⁹Er, ¹⁶⁷Er, ¹⁶⁹Tm, and ²³⁹Pu including omitted terms up to k_max = 10, i.e., accounting for six omitted terms at N⁴LO, was enough to stabilize the extraction of the LECs.

As an example, we show the stability of the extracted LEC, A₁, 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

FIGURE 4. Posteriors for A₁ describing ¹⁶⁹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₁ 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

FIGURE 5. Posteriors for B₁ and B describing ²³⁹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

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 ²³⁹Pu at N⁴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

FIGURE 7. The distribution of E_K for ¹⁶⁹Er at N⁴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 ¹⁸³W, ¹⁶⁷Tm and ²³⁵U. In ¹⁸³W and ¹⁶⁷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 ²³⁵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 ¹⁶⁷Tm we needed to go to k_max = 12 to get stable results, and for ²³⁵U and ¹⁸³W we needed to go to k_max = 18.

For K = 3/2 bands, we were able to extract stable LECs from the ¹⁵⁹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.) ¹⁵⁷Gd is stable across orders and I_max and we get stable results at k_max = 12, while ¹⁵⁵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⁴LO for the nuclei considered in this work are given in Tables 1, 2 respectively.

TABLE 1

TABLE 1. The median value of the LECs at N⁴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₁ and the parameters (B₁, A₃) refer to K =1/2 and K =3/2 bands respectively. All the numbers are in units of keV.

TABLE 2

TABLE 2. The median value of the error model parameters at N⁴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-χ² distribution as a prior for ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$ we tried truncated Gaussians with mean zero and standard deviations σ = 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 ${\bar{c}}_{\mathit{odd}}$ and ${\bar{c}}_{\mathit{even}}$ larger than those resulting from scaled-inverse-χ² priors with ν = 1 and τ² = 1.

The change in prior for ${\bar{c}}_{\mathit{odd}}$ and ${\bar{c}}_{\mathit{even}}$ does not significantly change the posteriors for the LECs: the corresponding central values differ by less than 1%, and are consistent with each other within the 68% credible intervals. Central values of the posteriors for W differ by less than 15%, and were similarly consistent.

The strongest dependence on the prior is that exhibited by the posteriors for ${\bar{c}}_{\mathit{odd}}$ and ${\bar{c}}_{\mathit{even}}$: the central values differ in some cases by more than 50%. However, even these values are consistent with each other within 68% credible intervals, since the posteriors for the characteristic sizes are broad.

The changes in the posteriors of W on one hand and ${\bar{c}}_{\mathit{odd}}$ & ${\bar{c}}_{\mathit{even}}$ on the other are anticorrelated. We only care about combinations of them to set the size of the theory error and the expected size of the LECs. Thus, the dependence of the theory error and the expected size of the LECs on the prior for the characteristic sizes is less profound. The difference in sizes of the theory error resulting from the chosen priors is less than 20% for all cases except ¹⁵⁷Gd, where the difference is about 40%.

For all results that follow we used the scaled-inverse-χ² distribution with ν = 1 and τ² = 1 as the prior for both ${\bar{c}}_{\mathit{odd}}$ and ${\bar{c}}_{\mathit{even}}$, in keeping with the naturalness assumption.

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 ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$ is small in these cases: the theory error changes by about 10% when the prior is changed.

FIGURE 8

FIGURE 8. The posterior predictive distribution for energy-level residuals at N⁴LO and k_max =10 for ¹⁶⁹Er and I_max =17.5 (top panel) and at k_max =18 for ²³⁵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

FIGURE 9. A 2D cut of the posterior predictive distribution at N⁴LO and k_max =10 for ¹⁶⁹Er and I_max =17.5 (top panel) and at k_max =18 for ²³⁵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

FIGURE 10. The size of the NLO LEC, A₁ (top panel) and the N³LO LEC, B₁ for K =1/2 bands and A₃ 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₁ and we do not have them in the top panel. ⁹⁹Tc and ¹⁵⁵Gd are outliers and we exclude them from the plots (LECs values for these nuclei can be found in Table 1).

FIGURE 11

FIGURE 11. The size of the N²LO LEC, ΔA (top panel) and the N⁴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. ⁹⁹Tc and ¹⁸³W are outliers and we exclude them from the both plots. We also exclude ¹⁵⁵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 ¹⁵⁵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

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. ⁹⁹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,4LO is less of a concern than the one at N²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 ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$. For ¹⁵⁷Gd changing the prior on ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$ to a truncated normal allowed for smaller values of ${\bar{c}}_{\mathit{odd}}$ and A₃ was then equal to the expected size (i.e., the point for ¹⁵⁷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 ¹⁵⁵Gd and ¹⁵⁹Dy; this may occur because there is strong N⁵LO energy-level staggering present in the data for these nuclei.

The size of ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$ is constrained by both the sizes of the LECs and the size of the theory error. In a good systematic expansion the tension between those factors on setting the size of the $\bar{c}$’s would be small and one number apiece would suffice to represent the even and odd order corrections. However when the systematic expansion is hindered, as in the case for 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 ${\bar{c}}_{\mathit{even}}$, however, B is small and favors smaller values of ${\bar{c}}_{\mathit{even}}$. The eventual result is a compromise. This tension may be exacerbated by the truncation error also providing information on the size of the $\bar{c}$’s.

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, $\vec{j}$, we cannot exclude the possibility that the systematic behavior of the EFT is hindered in bands build on top of single-particle orbitals with larger values of 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⁴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⁴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, ${\bar{c}}_{\mathit{even}}$ and ${\bar{c}}_{\mathit{odd}}$. The values of ${\bar{c}}_{\mathit{odd}}$ are typically smaller than those for ${\bar{c}}_{\mathit{even}}$, in agreement with results from previous studies where the LECs were fitted to the smallest possible data sets. The difference of the characteristic sizes of even and odd LECs has its origin in the physics behind the corresponding contributions to the effective Lagrangian: while odd-order contributions correct the particle-rotor interaction, even-order contributions include terms that depend exclusively on the rotor degrees of freedom, thus correcting the physics of the core to which the particle is coupled. Here this conclusion was reached solely on the basis of experimental data; we assumed equal priors for both characteristic sizes.

Although the distributions of ${\bar{c}}_{\mathit{odd}}$ and ${\bar{c}}_{\mathit{even}}$ change depending on the choice of the their priors, that does not significantly change the distributions of the LECs. Altering the priors also does not have a large effect on the size of the theory error, which changes by less than 20% for nearly all cases.

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

37. Foreman-Mackey D, Hogg DW, Lang D, Goodman J. EMCEE: the MCMC hammer. Publ Astron Soc Pac (2013) 125:306–12. arXiv:1202.3665 [astro-ph.IM]]. doi:10.1086/670067

CrossRef Full Text | Google Scholar