Electromagnetic character of the competitive γγ/γ-decay from 137mBa

Second-order processes in physics is a research topic focusing attention from several fields worldwide including, for example, non-linear quantum electrodynamics with high-power lasers, neutrinoless double-β decay, and stimulated atomic two-photon transitions. For the electromagnetic nuclear interaction, the observation of the competitive double-γ decay from 137mBa has opened up the nuclear structure field for detailed investigation of second-order processes through the manifestation of off-diagonal nuclear polarisability. Here, we confirm this observation with an 8.7σ significance, and an improved value on the double-photon versus single-photon branching ratio as 2.62 × 10−6(30). Our results, however, contradict the conclusions from the original experiment, where the decay was interpreted to be dominated by a quadrupole-quadrupole component. Here, we find a substantial enhancement in the energy distribution consistent with a dominating octupole-dipole character and a rather small quadrupole-quadrupole component in the decay, hindered due to an evolution of the internal nuclear structure. The implied strongly hindered double-photon branching in 137mBa opens up the possibility of the double-photon branching as a feasible tool for nuclear-structure studies on off-diagonal polarisability in nuclei where this hindrance is not present.

P olarisability 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 investigation 1 . 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 obtained 2 from the photonuclear population of excited states, where the transition matrix elements of the wave functions correspond to the electric dipole transition, E1, between the ground state, I 0 , and an excited state, I n , with e the elementary unit charge and E n 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 offdiagonal nuclear polarisability can be defined analogously to Eq. (1) in terms of either electric and magnetic components, or components of different multipolarities as or, Due to the parity conserving properties of the strong force, these decays can only be observed between two different states, I i and I f . 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-Meyer 3 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 CaF 2 :Eu 2+ crystals 4 . 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 16 O 5,6 , 40 Ca 7 , and 90 Zr 7 . 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 systems 8,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 required 10 that the observation of the γγ/γ decay mode was announced 11 . The set-up used for that experiment consisted of five LaBr 3 :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) e 2 fm 4 /MeV and α E3M1 = 10.1(4.2) e 2 fm 4 /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 component 12 . The decay diagram of this process is shown in Fig. 1.
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-up 13,14 or highly charged radioactive ions 15 . Here we report on an experiment using the ELI Gamma Above Neutron Threshold (ELIGANT) detector system 16 (ELI-NP) facilities [18][19][20] in a configuration 21 similar to what was used in reference 11 . The experimental set-up was optimised for obtaining a clean signal over a wide angular range 21 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 137 Ba, contradicting the conclusions of the original experiment 11 . 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 sources 22 , or even in-beam experiments within reasonable beam times, to follow the evolution of the quadrupole-quadrupole strength.
While ELIGANT consists of both LaBr 3 :Ce and CeBr 3 detectors, the CeBr 3 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 toolkit 23 based on the GEANT4 framework 24 , and a 152 Eu source with an activity of 460 kilo Becquerel (kBq) and a 60 Co source with an activity of 60 kBq. For a comprehensive overview, see reference 21 . The γγ/γ-decay data on 137 Ba were collected using a Table 1 Experimental and calculated α coefficients and γγ/γ decay branching ratios. The Γ γγ /Γ γ decay branching ratio is shown both with unquenched (g eff s ¼ g bare s ) and quenched gyromagnetic spin factors (g eff s ¼ 0:6g bare s ). The latter limit was chosen based on the reproduction of individual reduced transition probabilities. Depending on the calculation the values of g eff s to best reproduce nuclear data are typically within this range. Thus, these limits should be representative of the uncertainties in the theoretical calculations, giving a range of~50% for both the α M2E2 and α E3M1 values for both models between the two extremes. The listed values closest to the measured branching are shown in bold font. The best fit for the decay branching ratio for the EDF + QPM calculations, not listed here, is obtained when choosing g eff Energy spectra. From the data set obtained with the 137 Cs 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, Δt 1,2 ≤ 655 ps. This condition was obtained from the coincident 444 kilo electron-volt (keV) and 245 keV γ rays from the 2 þ 2 ! 4 þ 1 ! 2 þ 1 decay chain in 152 Sm following the electron capture decay of 152 Eu. Corrections for detector efficiencies were done on an event-by-event basis 21 . A time difference of 20 ns ≤ Δt 1,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 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 E 1 + E 2 with these conditions imposed is shown in Fig. 3.
Branching. As experimental observable to evaluate the relative decay probability, we use the definition of the integrated differential branching ratio 11 , 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 E 1 = 180 keV and E 2 = 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. 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.  Fig. 3 Summed double-γ energy spectrum and data reduction. Black data points show the summed energy of two coincident photons detected in the CeBr 3 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 jE γ;1 À E γ;2 j<300 keV. c Final data with the additional condition that jE γ;1 À E γ;2 j<960 À ðE γ;1 þ E γ;2 Þ keV to remove the cosmic-ray induced background. The error bars represent the one standard-deviation statistical uncertainty.
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 dΓ γγ dω / ω 5 ω 0 5 for M2E2 and as dΓ γγ dω / ω 3 ω 0 7 for E3M1 with ω 0 ¼ 662 À ω. It is clear from these relations that the energy-sharing distributions are expected to have a maximum at E γ ¼ E 0 ¼ 331 keV for the M2E2 type transitions, while an asymmetric maximum is expected at E γ = 200 keV and E 0 ¼ 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 highenergy 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 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 .

Discussion
To understand these results, we performed theoretical calculations of the polarisation functions from Eq. (2), using the QPM 25 approach. The application of the QPM in the case of odd-mass spherical nuclei is discussed in detail in reference 26 . In particular, the nuclear structure of 137 Ba was studied within the framework of this model in references 27,28 and in reference 11 . In the work presented here, the calculations were built on the EDF theory coupled with the QPM 29 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 microscopically [30][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 reference 11 , 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  Fig. 5 Multipole nature of the γγ/γ decay. 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 intensity 8 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 reference 11 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. 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 MCSM 33,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 νh 11/2 . The J π = 7/2 + state is different, however, with most of the neutron hole occupation is in νd 3/2 , coupled to a 2 + state of six valence protons. The νg 7/2 orbital itself is almost full. This is in contrast with the EDF + QPM results where the 2 + ⊗ νd 3/2 contribution is 38.7% and the νg 7/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 νh 11/2 and νd 3/2 , or by utilising a minor νg 7/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 μ 2 fm 2 , three orders of magnitudes less than predicted by the EDF + QPM model where B(M2) = 14.9 μ 2 fm 2 . 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 νd 3/2 to νd 5/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.
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" CeBr 3 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 V pp , 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, α S 0 L 0 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 0 J ðS; L; S 0 ; L 0 Þ, and Legendre polynomials, P l ðcos θÞ, as where the generalised polarisation functions are defined as 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 interactions 25 . 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 in 29,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 possible 36 . 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 nucleus 25 . 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 137 Ba 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 formalism 37,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 singleparticle orbitals and excitation phonons that are constructed from the excited states in the neighbouring even-even core nucleus : The notation α þ jm is the quasiparticle creation operator with shell quantum numbers j ≡ [(n, l, j)] and projection m; Q þ λμ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 and D λi j ðJνÞ are the quasiparticle and 'quasiparticle ⊗ phonon' amplitudes for the ν state. The components ½α þ jm Q þ λμ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 E x = 20 MeV. The calculations of the α-coefficients of the double-γ decay probability of 137 Ba 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 E x = 10 MeV.
In the case of the E1 transitions, we have used effective charges e eff p ¼ ðN=AÞe (for protons) and e eff n ¼ ÀðZ=AÞe (for neutrons) to separate the centre of mass motion and 'bare' values for E2 and E3 transitions e p = e (for protons) and e n = 0 (for neutrons), where e is the electron charge. Following previous QPM calculations 32 , the magnetic transitions are calculated with a quenched effective spin-magnetic factor g eff s . The influence of the g eff s 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 bare s . The theoretical observations indicate that the values g eff s ¼ 0:6 À 0:7g bare s which are in agreement with our previous findings 27, 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, jΨ N b i, are obtained as a superposition of spin (I) and parity (π) projected Slater determinant basis states, jϕ n i, where N b is the number of basis states, P Iπ MK is the spin-parity projection operator, and the f N b n;K 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 1g 9/2 , 1g 7/2 , 2d 5/2 , 2d 3/2 , and 3s 1/2 even-parity orbitals, as well as the 1h 11/2 , 2f 7/2 , and 3p 3/2 oddparity orbitals. The two-body matrix elements were obtained from the JUN45 and SNBG3 data sets 39,40 , and the V MU interaction 41 . To obtain the transition matrix elements effective proton and neutron charges e p = 1.25 and e n = 0.75, and gyromagnetic factors g ℓ,p = 1, g ℓ,n = 0, g s,p = 5.586, and g s,n = − 3.826 was used. The calculations followed the procedure for the tin isotope chain closely 42 . Said reference and references within contains a detailed description of the procedure.

Data availability
Raw data were obtained at the Extreme Light Infrastructure -Nuclear Physics facility, Romania. All the data used to support the findings of this study are available from the authors upon reasonable request. The final data points can be obtained from https://doi. org/10.17632/skhmjshxdj.

Code availability
Sorting codes were developed at the Extreme Light Infrastructure -Nuclear Physics facility, Romania. All the codes for the experimental data used in this study are available from the authors upon reasonable request.