Introduction

Polarisability is a fundamental concept in physics and chemistry defined from the principles of electromagnetic interaction. It describes how applied electric or magnetic fields induce an electric or magnetic dipole, or higher-order multipole, moment in the matter under investigation1. In nuclear physics, the simple concept of polarisability influences observables over a broad range of topics. For example, the static dipole polarisation of the shape of the ground and excited states in atomic nuclei is influenced by the coupling to high-energy collective modes like the giant dipole resonance (GDR) via virtual excitations. In this case, the nuclear static dipole polarisability, αd, is obtained2 from the photonuclear population of excited states,

$${\alpha }_{{\rm{d}};{\rm{E1}}}=2e\sum _{n}\frac{{\left|\langle {I}_{0}| | {\rm{E1}}| | {I}_{n}\rangle \right|}^{2}}{{E}_{n}-{E}_{0}},$$
(1)

where the transition matrix elements of the wave functions correspond to the electric dipole transition, E1, between the ground state, I0, and an excited state, In, with e the elementary unit charge and En the energy of the state.

By expanding the concept of polarisability beyond the scalar case, one can divide the polarisability tensor into separate components. Typically, these are either spatial components like the birefringence properties of crystals or electric and magnetic multipole components. Within the nuclear structure framework, this type of off-diagonal polarisabilities can appear in very weak second-order processes. In the electromagnetic case, the off-diagonal nuclear polarisability can be defined analogously to Eq. (1) in terms of either electric and magnetic components, or components of different multipolarities as

$${\alpha }_{{\rm{M2E2}}}=\sum _{n}\frac{\langle {I}_{{\rm{f}}}| | {\rm{E2}}| | {I}_{n}\rangle \langle {I}_{n}| | {\rm{M2}}| | {I}_{{\rm{i}}}\rangle }{{E}_{n}-\omega }$$
(2)

or,

$${\alpha }_{{\rm{E3M1}}}=\sum _{n}\frac{\langle {I}_{{\rm{f}}}| | {\rm{M1}}| | {I}_{n}\rangle \langle {I}_{n}| | {\rm{E3}}| | {I}_{{\rm{i}}}\rangle }{{E}_{n}-\omega }.$$
(3)

Due to the parity conserving properties of the strong force, these decays can only be observed between two different states, Ii and If. In the definition above, the denominator depends on the interference frequency, ω, of the emitted γ rays and is approximated to be half of the initial state energy. This type of second-order electromagnetic processes of atoms was discussed in the doctorate dissertation of Maria Göppert-Meyer3 where she estimated a probability for an atomic two-photon absorption process relative to the single-photon process to be approximately 10−7, later to be confirmed with the observation of this effect in CaF2:Eu2+ crystals4.

For many years, double-γ decay was only observed in exceptional cases where both the ground state and the initial state have a spin-parity Jπ = 0+ character for the doubly magic nuclei 16O5,6, 40Ca7, and 90Zr7. Here single γ-emission is blocked, and only conversion-electron decay and double-γ decay are allowed. In these experiments, the obtained information consists of correlations between energies and angles of these γ-rays, used to determine the decay probabilities of electric and magnetic dipoles. For a generalisation of this phenomenon and the possibility to use it as a spectroscopic tool for a more fundamental understanding of the underlying physics, large state-of-the-art high-purity germanium (HPGe) detector systems8,9 have been used to search for the competitive double-γ (γγ/γ) decay, where also the single γ decay is allowed. Even though unsuccessful in that respect, these experiments successfully measured an E5 transition with the branching of 1.12(9) × 10−7. It is only with instrumentation developments of detector materials that can provide both the energy and time resolution required10 that the observation of the γγ/γ decay mode was announced11. The set-up used for that experiment consisted of five LaBr3:Ce detectors arranged in a planar configuration with relative angles of 72∘ between the detectors, providing angular distribution data points at 72∘ and 144∘. Thus, the collaboration could announce a γγ/γ decay signal with 5.5 σ (standard deviations) statistical significance, near but above the typical discovery limit of 5 σ. From the two angular data points as well as the energy spectrum of the individual γ rays at 72∘ angle, the off-diagonal polarisabilities αM2E2 = 33.9(2.8) e2fm4/MeV and αE3M1 = 10.1(4.2) e2fm4/MeV were favoured. While the observation of the peak associated with γγ/γ decay was statistically clear, the nature of this decay was more uncertain, having the two dominating multipolarity combinations separated only by a small statistical difference, favouring the αM2E2 component12. The decay diagram of this process is shown in Fig. 1.

Fig. 1: Decay diagram from the 137Cs ground state to the 137Ba ground state.
figure 1

Illustration of the single-γ and the two types of double-γ decay, as fed by the β decay of 137Cs, including half lives of 137Cs and 137mBa. The energy of the 137Cs ground state (Qβ) is given relative to the 137Ba ground state. Here, M4 corresponds to the single-photon decay. The blue and pink decays show the lowest octupole-dipole and quadrupole–quadrupole components, respectively.

Full size image

Given the nature of this experiment to observe a longstanding prediction of a fundamental concepts in quantum mechanics and quantum electrodynamics, and the possibility to extract nuclear structure observables from this, it is highly desirable to independently confirm this observation. Some possibilities that have been under discussion to perform this independent confirmation is to either return to the HPGe approach with complex detector systems and event processing like the Advanced GAmma Tracking Array (AGATA) set-up13,14 or highly charged radioactive ions15. Here we report on an experiment using the ELI Gamma Above Neutron Threshold (ELIGANT) detector system16,17 at the Extreme Light Infrastructure-Nuclear Physics (ELI-NP) facilities18,19,20 in a configuration21 similar to what was used in reference11. The experimental set-up was optimised for obtaining a clean signal over a wide angular range21 based on the reported intensities of the decay mode. Here, we can confirm the existence of the competitive double-photon decay process in atomic nuclei with an 8.7σ significance. We, however, find a significant octupole-dipole, E3M1, matrix element product contribution to the double-γ decay mode of 137Ba, contradicting the conclusions of the original experiment11. From our calculations using the energy-density-functional (EDF) + quasiparticle-phonon model (QPM) and the Monte Carlo shell model (MCSM), we find that both models reproduce the octupole-dipole component consistently, but the nature and the strength of the quadrupole–quadrupole component, differ significantly. It is interesting to note that this suggests an additional hindrance, reducing the γγ/γ branching with almost an order of magnitude in the most extreme case of Table 1, compared to calculations without this hindrance. This most extreme case is also the case that best reproduces the αE3M1 polarisability. This opens for the possibility of a significant increase of the γγ/γ branching in nuclei in this region that do not exhibit this hindrance. In this case, experiments would be feasible also with more exotic sources22, or even in-beam experiments within reasonable beam times, to follow the evolution of the quadrupole–quadrupole strength.

Table 1 Experimental and calculated α coefficients and γγ/γ decay branching ratios.
Full size table

Results

Experimental set-up

The experiment was performed using eleven 3″ × 3″ CeBr3 detectors from ELIGANT, shown in Fig. 2a. While ELIGANT consists of both LaBr3:Ce and CeBr3 detectors, the CeBr3 detectors were chosen to remove any possible source of background contribution from the natural radioactivity in lanthanum. The detector configuration was a circle with an inner radius to the front-face of the scintillators of 40 cm. This distance was enough to separate true coincidences from multiple Compton scattering of single γ rays using the photon time-of-flight (TOF), see Fig. 2b. The relative angles between the eleven detectors were 32.7∘, with an opening angle, given by the lead shielding, of  ±3.4∘. This gave five independent γγ-correlation angles centred at: 32.7∘, 65.5∘, 98.2∘, 130.9∘, and 163.6∘. The detectors were separated with a minimum of approximately 15 cm of effective lead shielding between two neighbouring detectors to remove any contribution from single Compton scattering between detector pairs at low angles. The set-up was characterised both with an in-house toolkit23 based on the GEANT4 framework24, and a 152Eu source with an activity of 460 kilo Becquerel (kBq) and a 60Co source with an activity of 60 kBq. For a comprehensive overview, see reference21. The γγ/γ-decay data on 137Ba were collected using a 137Cs source with an activity of 336  kBq for 49.5 days of active data taking. The source had a thin circular active area with a diameter of 3 mm, encapsulated in the centre of a cylindrical polymethylmethacrylate capsule with a diameter of 25 and 3 mm thickness, see Fig. 2.

Fig. 2: Experimental set-up.
figure 2

a Coincident γ rays could originate either from true double-γ decay events illustrated with red cones, or from multiple Compton scattering between detectors illustrated with blue cones. The location and geometry of the source are also shown, with the active area marked in black. b Multiple Compton scattering events were rejected by the time difference (Δt) between the γ-ray interactions, shown in the blue histogram. The time condition for prompt γ-rays is shown as red dashed lines and verified with a 152Eu source.

Full size image

Energy spectra

From the data set obtained with the 137Cs source, a (γ1, γ2) coincidence matrix was constructed where the γ rays were considered coincident if the time difference between them were less than one standard deviation from the prompt time distribution, Δt1,2 ≤ 655 ps. This condition was obtained from the coincident 444 kilo electron-volt (keV) and 245 keV γ rays from the \({2}_{2}^{+}\to {4}_{1}^{+}\to {2}_{1}^{+}\) decay chain in 152Sm following the electron capture decay of 152Eu. Corrections for detector efficiencies were done on an event-by-event basis21. A time difference of 20 ns ≤ Δt1,2 ≤ 820 ns was used to estimate the uncorrelated background events with two detected γ rays and subtracted after applying an appropriate scaling factor. To remove the background contribution from electron–positron pairs produced by cosmic rays a multiplicity-two condition was assigned together with an additional energy condition that \(\left|{E}_{1}-{E}_{2}\right|<960-\left({E}_{1}+{E}_{2}\right)\) keV. The full data set, as well as the different angular groups, were used to construct the summed energy spectra. The peaks were fitted assuming a quadratic background both with a Gaussian distribution as well as GEANT4 simulated data. Both fitting methods gave consistent results. The full summed spectrum of E1 + E2 with these conditions imposed is shown in Fig. 3.

Fig. 3: Summed double-γ energy spectrum and data reduction.
figure 3

Black data points show the summed energy of two coincident photons detected in the CeBr3 detectors for events with a multiplicity of two. Grey data points show the sum energy spectrum when the multiplicity is larger than two, which mainly correspond to the background induced by cosmic ray showers. We also show the fit to the data of a quadratic background as a dashed red line and the fit of the background plus a Gaussian peak as a solid red line. a Raw data before any conditions. b Reduced data with the condition that the energy difference between two γ rays is \(|{E}_{\gamma ,1}-{E}_{\gamma ,2}|<300\) keV. c Final data with the additional condition that \(|{E}_{\gamma ,1}-{E}_{\gamma ,2}|<960-({E}_{\gamma ,1}+{E}_{\gamma ,2})\) keV to remove the cosmic-ray induced background. The error bars represent the one standard-deviation statistical uncertainty.

Full size image

Branching

As experimental observable to evaluate the relative decay probability, we use the definition of the integrated differential branching ratio11,

$$\delta ({E}_{1},{E}_{2},{\theta }_{1,2})=\frac{{(4\pi )}^{2}}{{\Gamma }_{\gamma }}\mathop{\int}\nolimits_{{E}_{1}}^{{E}_{2}}{\rm{d}}\omega {\left.\frac{{\rm{d}}{\Gamma }_{\gamma \gamma }^{5}}{{\rm{d}}\omega {\rm{d}}\Omega {\rm{d}}\Omega ^{\prime} }\right|}_{{\theta }_{1,2}}.$$
(4)

In this definition, Γγ is the total single-gamma decay width, proportional to the size of the single-gamma peak. Given an angle, θ1,2, the differential decay is integrated over the frequency of the γ ray, ω. The frequency is proportional to the energy, and the integration limits are taken as the edges of the energy bin of interest. In the experimental spectrum, a natural low-energy limit comes from the low-energy threshold of the detectors around 120 keV. However, to reduce the contamination from the 511 keV γ-rays originating from electron–positron annihilation, the integration limits E1 = 180 keV and E2 = 331 keV were chosen. The upper limit was chosen as half of the total energy as we are not able to distinguish any relative ordering of the γ rays. This procedure was performed for all combinations of θ1,2 and δ was evaluated as a function of angle. The results of this evaluation are shown in Fig. 4.

Fig. 4: Angular distribution.
figure 4

The angular correlation of the two photons emitted in the double-γ decay from this work and reference11, compared with the expected angular distributions of pure M2E2 and E3M1 decay. The error bars represent the one standard-deviation statistical uncertainty.

Full size image

This data can be directly fitted to the generalised polarisability functions of Eq. (6) discussed in the methods section, using only αM2E2 and αE3M1 as free parameters. Other components like αE2M2 or αM3E1 could, in principle, also contribute. However, the general polarisability functions are linearly dependent in the exchange of terms, weighted by the coefficients given by the Wigner 6j symbols, and this experiment is not sensitive to this ordering. These additional components are, furthermore, expected to be small. Thus, we restrict the discussion to the αM2E2 and αE3M1 polarisabilities from here on.

Energy-sharing distributions

The angular distributions themselves are not enough to completely distinguish between the contribution from the different polarisabilities. When calculating the goodness-of-fit (χ2), two local minima corresponding to either a large αM2E2 component or a large αE3M1 component appear. Instead, it is necessary to study the energy-sharing distributions between the two individual γ-rays. From Eqs. (6) and (7) in the methods section it is clear that the energy dependence of the decay for the two different cases follows \(\frac{{\rm{d}}{\Gamma }_{\gamma \gamma }}{{\rm{d}}\omega }\propto {\omega }^{5}\omega {^{\prime} }^{5}\) for M2E2 and as \(\frac{{\rm{d}}{\Gamma }_{\gamma \gamma }}{{\rm{d}}\omega }\propto {\omega }^{3}\omega {^{\prime} }^{7}\) for E3M1 with \(\omega ^{\prime} =662-\omega\). It is clear from these relations that the energy-sharing distributions are expected to have a maximum at \({E}_{\gamma }=E^{\prime} =331\) keV for the M2E2 type transitions, while an asymmetric maximum is expected at Eγ = 200 keV and \(E^{\prime} =442\) keV for the E3M1 type transitions.

For this purpose, δ from Eq. (4) was evaluated in separate slices of 30 keV energy difference between the low- and high-energy limit of Eγ. Figure 5a shows the results of these evaluations. A χ2 value was then calculated based on (4) for different values of αM2E2 and αE3M1 simultaneously using the energy-integrated angular data points and the angle-summed energy data points as

$${\chi }^{2}=\sum_{{{{\theta_{1,2}^{i}}\atop {{E}_{1}=181}}\atop {{E}_{2}=331}}}\frac{{\delta }^{2}({E}_{1},{E}_{2},{\theta }_{1,2}^{i})}{{\sigma }_{\delta }^{2}({E}_{1},{E}_{2},{\theta }_{1,2}^{i})}+\sum_{{{{E}_{1}={E}^{i}_{{\rm{low}}}}\atop {{E}_{2}={E}^{i}_{{\rm{high}}}}}}\frac{{\delta }^{2}({E}_{1},{E}_{2},\Sigma {\theta }_{1,2})}{{\sigma }_{\delta }^{2}({E}_{1},{E}_{2},\Sigma {\theta }_{1,2})},$$
(5)

where σδ is the statistical uncertainty in each data point, including both the signal and the subtracted background. The systematic uncertainty mainly originates from uncertainties in the intrinsic and geometric efficiencies of the set-up and is expected to be on the order of a few %, much smaller than the statistical uncertainties, and have been neglected in this expression. When including the data from Walz et al., only the energy distribution of the 72∘ data was included in the second part of Eq. (5), and the lower energy summation limit for the 144∘ data point was set to 206 keV in the first part of Eq. (5). The resulting χ2 surface is shown in Fig. 5b. As seen here, the χ2 analysis from this data favours a large αE3M1 component, in contradiction with both the experimental interpretation and theoretical conclusions reported in ref. 11.

Fig. 5: Multipole nature of the γγ/γ decay.
figure 5

a Energy-sharing distribution for the two photons in the double-γ decay compared with the expected energy distributions of pure M2E2 and E3M1 decay, as well as the two-step E5 + M1 decay where the measured intensity8 has been subtracted from the 300 keV data point. The data points correspond to the sum of the differential branching ratio defined in Eq. (4) over all the available angles. b Two-dimensional goodness-of-fit, χ2, plot for the two α parameters with the experimental data. The contours are separated by one standard deviation. The data from the present work is shown as a black point, data from reference11 is shown as a red point, and a fit with the two data sets combines is shown as a purple point. The error bars represent the one standard-deviation statistical uncertainty, except the error bars in Eγ,low that represent the width of the energy bin.

Full size image

Discussion

To understand these results, we performed theoretical calculations of the polarisation functions from Eq. (2), using the QPM25 approach. The application of the QPM in the case of odd-mass spherical nuclei is discussed in detail in reference26. In particular, the nuclear structure of 137Ba was studied within the framework of this model in references27,28 and in reference11. In the work presented here, the calculations were built on the EDF theory coupled with the QPM29 to obtain magnetic and electric spectral distributions. The model parameters of the EDF + QPM approach are firmly determined from nuclear structure data or derived fully microscopically30,31,32. The theoretical results are shown in Table 1 and agree with the data in terms of absolute branching strength, Γγγ/Γγ. In addition, the EDF + QPM used here predicts a significantly larger αE3M1 than the value reported in ref. 11 from the QPM, close to our experimental observations. However, the relative magnitude of the αM2E2 and αE3M1 coefficients obtained from the EDF + QPM theory, as well as from reference11, are different than the experimental results obtained in this work. In particular, the present measurement indicates that the αM2E2 coefficient is significantly smaller than previously reported, and at this level of complexity EDF + QPM it is not able to account for the apparent discrepancy with the experimental data.

To understand the origin of this discrepancy, the properties of the dominant, low-lying, states were investigated from another perspective using the state-of-the-art nuclear MCSM33,34. These calculations were used to extract information from the three lowest-energy Jπ = 7/2+ states, the five lowest-energy Jπ = 5/2+ states, as well as the ground Jπ = 3/2+ and isomeric Jπ = 11/2− states. The neutron component in the MCSM wave function of the isomeric Jπ = 11/2− state is dominated by a single neutron hole in νh11/2. The Jπ = 7/2+ state is different, however, with most of the neutron hole occupation is in νd3/2, coupled to a 2+ state of six valence protons. The νg7/2 orbital itself is almost full. This is in contrast with the EDF + QPM results where the 2+ ⊗ νd3/2 contribution is 38.7% and the νg7/2 single-particle component is 51.3%. Thus, the odd-neutron contribution to the M2 transition rate in the MCSM would require a highly hindered transition between νh11/2 and νd3/2, or by utilising a minor νg7/2 vacancy. This, gives rise to a strongly hindered M2 transition within the MCSM, with a reduced transition probability, B(M2) = 13.5 × 10−3 μ2fm2, three orders of magnitudes less than predicted by the EDF + QPM model where B(M2) = 14.9 μ2fm2. This can explain the observed suppression of αE2M2. It is interesting to note that with increasing excitation energy, the MCSM predict a smooth change in orbital occupation from νd3/2 to νd5/2, constructively adding to the M2 transition strength for all the calculated 7/2+ transitions in contrast to the EDF + QPM where all higher-lying states act destructively. Table 2 lists the contributing low-lying matrix elements discussed here.

Table 2 Calculated matrix elements.
Full size table

Regarding the αE3M1 component of the decay, the main components obtained from the EDF + QPM calculations are from the coupling of the single-particle mode with the surface vibrations of the even-even core. As a consequence, due to the exchange of the collective \({3}_{1}^{-}\) octupole phonon, we obtain a rather strong E3 transition, consistent with our experimental observations. For these states, the EDF + QPM and the MCSM give a consistent picture with a constructive addition to the strength for each successive state among the first three excited states with the main difference that in the EDF + QPM, the main contribution comes from the \(5/{2}_{1}^{+}\) state while the MCSM predicts that the \(5/{2}_{2,3}^{+}\) states are dominating.

Methods

Experimental set-up

The set-up consisted of eleven 3” × 3” CeBr3 detectors coupled with Hamamatsu R6233 photomultiplier tubes and built-in voltage dividers. The high voltage for the photomultiplier tubes was provided by a CAEN SY4527 power supply. The signals were read out using one CAEN V1730 digitiser operating with a 14-bit resolution at a 500 MS/s sampling rate and a dynamic range of 0.5 Vpp, running PSD firmware. The digitisers were controlled using the Multi Instance Data Acquisition System software and triggered individually. Each event consisted of the energy, the time-stamp, and the the digitised voltage pulse from the detector. The sub-nanosecond time information was obtained from the value of the time-stamp corrected by a digital interpolation of the sampling points in the recorded pulse, at half of the maximum value of the pulse and interpolated using a quadratic polynomial.

Polarisation functions

To obtain the nuclear polarisabilites, \({\alpha }_{S^{\prime} L^{\prime} SL}\), from the differential decay probability we follow the theoretical treatment in refs. 11,35. Here the differential decay probability can be expressed in terms of generalised polarisation functions, \({P}_{J}^{\prime}(S,L,S^{\prime} ,L^{\prime} )\), and Legendre polynomials, \({P}_{l}(\cos \theta )\), as

$$\frac{{{\rm{d}}}^{5}{\Gamma }_{\gamma \gamma }}{{\rm{d}}\omega \Omega \Omega ^{\prime} }=\frac{\omega \omega ^{\prime} }{96\pi }\sum {P}_{J}^{\prime}({S}_{1}^{\prime}{L}_{1}^{\prime}{S}_{1}{L}_{1}){P}_{J}^{\prime}({S}_{2}^{\prime}{L}_{2}^{\prime}{S}_{2}{L}_{2})\sum {a}_{l}^{\xi }{P}_{l}(\cos \theta ),$$
(6)

where the generalised polarisation functions are defined as

$${P}_{J}^{\prime}(S^{\prime} L^{\prime} SL)= \, {(-1)}^{S+S^{\prime} }2\pi {(-1)}^{{I}_{i}+{I}_{f}}{\omega }^{L}\omega {^{\prime} }^{L^{\prime} }\cdot \\ \sqrt{\frac{L + 1}{L}}\sqrt{\frac{L^{\prime} + 1}{L^{\prime} }}\frac{\sqrt{2L + 1}\sqrt{2L^{\prime} + 1}}{(2L+1)!!(2L^{\prime} +1)!!}\cdot \\ \left(\left\{\begin{array}{lll}L&L^{\prime} &J\\ {I}_{f}&{I}_{i}&I\end{array}\right\}{\alpha }_{S^{\prime} L^{\prime} SL}+{(-1)}^{S+S^{\prime} }\left\{\begin{array}{lll}L^{\prime} &L&J\\ {I}_{f}&{I}_{i}&I\end{array}\right\}{\alpha }_{SLS^{\prime} L^{\prime} }\right).$$
(7)

The sums in Eq. (6) run over all the permutations of electric, S = 0, and magnetic, S = 1, combinations with multipolarity, L, allowed in the decay, and over all Legendre polynomials with non-zero coefficients. The general polarisability functions in Eq. (7) consist of a linear combination of the off-diagonal polarisabilities of the nucleus weighted by coefficients determined by the corresponding angular momentum algebra of the decay.

The quasiparticle-phonon model

The QPM Hamiltonian includes mean field, pairing interaction and separable multipole and spin-multipole interactions25. The mean field for protons and neutrons is defined as a Woods-Saxon potential with parameter sets derived self-consistently from a fully microscopic Hartree-Fock-Bogoljubov (HFB) calculations described in29,36. The method assures a good description of nuclear ground-state properties by enforcing that measured separation energies and nuclear radii are reproduced as close as possible36. The pairing and residual interaction parameters are fitted to reproduce the odd-even mass differences of neighbouring nuclei as well as the experimental values of the excitation energies and reduced transition probabilities of low-lying collective and non-collective states in the even-even core nucleus25. Of particular importance in these studies is the determination of the isovector spin-dipole coupling constant which is extracted from comparison to data from ref. 31 and fully self-consistent quasiparticle random phase approximation (QRPA) calculations using the microscopic EDF of ref. 30. Single-particle (s.p.) energies of the lowest-lying excited states in 137Ba are fine-tuned to experimental values to achieve the highest accuracy in the description of the experimental data. We point out that the s.p. energies problem is not a matter of the interaction parameters but originate in the quasiparticle spectrum, which indicates the necessity to go beyond the static mean-field formalism37,38.

In the QPM the wave functions of the excited states of an even-odd nucleus are constructed from a combination of quasiparticles originating from the single-particle orbitals and excitation phonons that are constructed from the excited states in the neighbouring even-even core nucleus :

$${\Psi }_{\nu }(JM)={C}_{J}^{\nu }\left\{{\alpha }_{JM}^{+}+\sum _{\lambda \mu i}{D}_{j}^{\lambda i}(J\nu ){[{\alpha }_{jm}^{+}{Q}_{\lambda \mu i}^{+}]} \, {JM}\right\}{\Psi }_{0}$$
(8)

The notation \({\alpha }_{jm}^{+}\) is the quasiparticle creation operator with shell quantum numbers j ≡ [(n, l, j)] and projection m; \({Q}_{\lambda \mu i}^{+}\) denotes the phonon creation operator with the angular momentum λ, projection m and QRPA root number i; Ψ0 is the ground state of the neighboring even-even nucleus and ν stands for the number within a sequence of states of given angular momentum Jπ and projection M. The coefficients \({C}_{J}^{\nu }\) and \({D}_{j}^{\lambda i}(J\nu )\) are the quasiparticle and ’quasiparticle ⊗ phonon’ amplitudes for the ν state. The components \({[{\alpha }_{jm}^{+}{Q}_{\lambda \mu i}^{+}]}_{JM}\) of the wave function (8) may violate the Pauli principle. The exact commutation relations between quasiparticle and phonon operators are used to solve this problem. The properties of the phonons are determined by solving QRPA equations from refs. 25,26. The model basis includes one-phonon states with spin and parity Jπ = 1±, 2±, 3±, 4±, 5± and excitation energies up to Ex = 20 MeV. The calculations of the α-coefficients of the double-γ decay probability of 137Ba include all low-energy excited states with spin and parity Jπ = 1/2±, 3/2±, 5/2±, 7/2±, 9/2± and excitation energies up to Ex = 10 MeV.

In the case of the E1 transitions, we have used effective charges \({e}_{{\rm{p}}}^{{\rm{eff}}}=(N/A)e\) (for protons) and \({e}_{{\rm{n}}}^{{\rm{eff}}}=-(Z/A)e\) (for neutrons) to separate the centre of mass motion and ’bare’ values for E2 and E3 transitions ep = e (for protons) and en = 0 (for neutrons), where e is the electron charge. Following previous QPM calculations32, the magnetic transitions are calculated with a quenched effective spin-magnetic factor \({g}_{{\rm{s}}}^{{\rm{eff}}}\). The influence of the \({g}_{{\rm{s}}}^{{\rm{eff}}}\) parameter on the experimental observables related to electromagnetic transitions of lowest-lying states and double-γ decay probability coefficients was investigated by carrying out EDF + QPM calculations for several choices of this parameter between 0.6 and 1 of the value of the ‘bare’ spin-magnetic moment, \({g}_{{\rm{s}}}^{{\rm{bare}}}\). The theoretical observations indicate that the values \({g}_{{\rm{s}}}^{{\rm{eff}}}=0.6-0.7{g}_{{\rm{s}}}^{{\rm{bare}}}\) which are in agreement with our previous findings27,28,32 reproduce quite well the experimental data on M1 and M2 transition strengths and the angular distribution of the two photons of the double-γ decay.

Monte Carlo shell model

In the MCSM, the approximated wave functions, \(|{\Psi }_{{N}_{b}}\rangle\), are obtained as a superposition of spin (I) and parity (π) projected Slater determinant basis states, \(|{\phi }_{n}\rangle\),

$$\left|{\Psi }_{{N}_{b}}\right\rangle =\mathop{\sum }\limits_{n=1}^{{N}_{b}}\mathop{\sum }\limits_{K=-I}^{I}{f}_{n,K}^{{N}_{b}}{P}_{MK}^{I\pi }\left|{\phi }_{n}\right\rangle ,$$
(9)

where Nb is the number of basis states, \({P}_{MK}^{I\pi }\) is the spin-parity projection operator, and the \({f}_{n,K}^{{N}_{b}}\) coefficients are obtained from diagonalising the Hamiltonian. The set of basis states are selected by Monte Carlo methods and iteratively refined to minimise the ground-state energy. The model space for these calculations included the 1g9/2, 1g7/2, 2d5/2, 2d3/2, and 3s1/2 even-parity orbitals, as well as the 1h11/2, 2f7/2, and 3p3/2 odd-parity orbitals. The two-body matrix elements were obtained from the JUN45 and SNBG3 data sets39,40, and the VMU interaction41. To obtain the transition matrix elements effective proton and neutron charges ep = 1.25 and en = 0.75, and gyromagnetic factors gℓ,p = 1, gℓ,n = 0, gs,p = 5.586, and gs,n = − 3.826 was used. The calculations followed the procedure for the tin isotope chain closely42. Said reference and references within contains a detailed description of the procedure.