First principles electromagnetic responses in medium-mass nuclei

TL;DR

This paper reviews recent advances in calculating electromagnetic responses in medium-mass nuclei using coupled-cluster theory, highlighting new computational methods and improved agreement with experimental data.

Contribution

It introduces a many-body formulation of the Lorentz integral transform method for electromagnetic responses and revisits key observables with enhanced coupled-cluster calculations.

Findings

Successful calculation of photoabsorption cross sections for $^{16, 22}$O and $^{40}$Ca.

Revised estimates of electric dipole polarizability and neutron-skin thickness in $^{48}$Ca.

Better agreement with experimental data for electromagnetic observables.

Abstract

We review the recent progress made in the computation of electromagnetic response functions in light and medium-mass nuclei using coupled-cluster theory. We show how a many-body formulation of the Lorentz integral transform method allows to calculate the photoabsorption cross sections of $^{16, 22}$O and $^{40}$Ca. Then, we discuss electromagnetic sum rules, with particular emphasis on the electric dipole polarizability, $\alpha_D$. By including triples corrections in coupled-cluster theory, we revisit $^{48}$Ca, for which, beside the electric dipole polarizability, we had previously investigated the neutron and proton radii, as well as the size of the neutron-skin thickness [1]. We show that correlations among these observables still hold, albeit a better agreement with experiment is obtained for $\alpha_D$ and the prediction of a small neutron-skin thickness is further corroborated.

Tables

Table 1. Table 1: List of results with NN + + 3N Hamiltonians  Hebeler2011 ; Ekstroem2015 for the charge radius and the electric dipole polarizability for 16 O and 40 Ca in the D/D approximation. For the notation of the potentials we follow Ref.  Hebeler2011 . Experimental values are taken from Ref.  angeli2013 (radius) and Ref.  Ahrens75 ; Miorelli2018 ; Birkhan17 (electric dipole polarizability).

16O
Interaction	$R_{\mathrm{ch}}$ [fm]	${\alpha }_D$ [fm3]
2.0/2.0(EM)	2.62	0.46
2.0/2.0(PWA)	2.74	0.54
1.8/2.0(EM)	2.60	0.44
2.2/2.0(EM)	2.63	0.48
2.8/2.0(EM)	2.67	0.52
NNLOsat	2.71	0.58
Experiment	2.6991(52) angeli2013	0.58(1) Ahrens75
		0.568(9) Miorelli2018
40Ca
2.0/2.0(EM)	3.35	1.67
2.0/2.0(PWA)	3.55	2.03
1.8/2.0(EM)	3.31	1.57
2.2/2.0(EM)	3.38	1.75
2.8/2.0(EM)	3.44	1.94
NNLOsat	3.48	2.08
Experiment	3.4776(19) angeli2013	2.23(3) Ahrens75
		1.87(3) Birkhan17

Table 2. Table 2: Effect of triples corrections for α D subscript 𝛼 𝐷 \alpha_{D} of 16 O in fm 3 using various NN + + 3N Hamiltonians  Hebeler2011 ; Ekstroem2015 . For the notation of the potentials we follow Ref.  Hebeler2011 . The experimental value is taken from Refs.  Ahrens75 ; Miorelli2018 .

16O
Interaction	D/D	T-1/D
2.0/2.0(EM)	0.46	0.42
2.0/2.0(PWA)	0.54	0.48
1.8/2.0(EM)	0.44	0.41
2.2/2.0(EM)	0.48	0.43
2.8/2.0(EM)	0.52	0.45
NNLOsat	0.58	0.50
Experiment		0.58(1) Ahrens75
		0.568(9) Miorelli2018

Table 3. Table 3: Effect of triples corrections in 48 Ca for R p subscript 𝑅 𝑝 R_{p} , R n subscript 𝑅 𝑛 R_{n} , R skin subscript 𝑅 skin R_{\rm skin} in fm and α D subscript 𝛼 𝐷 \alpha_{D} in fm 3 using various NN + + 3N Hamiltonians  Hebeler2011 ; Ekstroem2015 as in Table  2 . Calculations were performed for the harmonic-oscillator frequency of ℏ ​ ω = 16 Planck-constant-over-2-pi 𝜔 16 \hbar\omega=16 MeV.

48Ca, $R_p$
Interaction	D	T-1
2.0/2.0(EM)	3.271	3.296
2.0/2.0(PWA)	3.502	3.530
1.8/2.0(EM)	3.225	3.247
2.2/2.0(EM)	3.307	3.332
NNLOsat	3.431	3.453
48Ca, $R_n$
Interaction	D	T-1
2.0/2.0(EM)	3.423	3.445
2.0/2.0(PWA)	3.656	3.681
1.8/2.0(EM)	3.375	3.395
2.2/2.0(EM)	3.461	3.483
NNLOsat	3.571	3.587
48Ca, $R_{\mathrm{skin}}$
Interaction	D	T-1
2.0/2.0(EM)	0.151	0.149
2.0/2.0(PWA)	0.154	0.150
1.8/2.0(EM)	0.150	0.148
2.2/2.0(EM)	0.154	0.151
NNLOsat	0.139	0.134
48Ca, ${\alpha }_D$
Interaction	D/D	T-1/D
2.0/2.0(EM)	2.118	1.952
2.0/2.0(PWA)	2.693	2.414
1.8/2.0(EM)	1.979	1.860
2.2/2.0(EM)	2.233	2.021
NNLOsat	2.645	2.259

Equations

$$|{\Psi }_0⟩=\exp (\hat{T})|{\Phi }_0⟩,$$

$${\hat{T}}_1$$

$${\hat{T}}_2$$

$${\hat{T}}_3$$

$$\cdots$$

$${\overline{H}}_N|{\Phi }_0⟩=E_0|{\Phi }_0⟩,$$

$${\overline{H}}_N=\exp (-\hat{T}){\hat{H}}_N\exp (\hat{T}),$$

$$0$$

$$R(\omega )=\sum _n|⟨{\Psi }_0|\hat{\Theta }|{\Psi }_n⟩{|}^2\delta (E_n-E_0-\omega ).$$

$$L({\omega }_0,\Gamma )=\frac{\Gamma }{\pi }\int d\omega \frac{R(\omega )}{(\omega -{\omega }_0{)}^2+{\Gamma }^2}\text{mbox},$$

$$\sum _n|{\Psi }_n⟩⟨{\Psi }_n|=1\text{mbox},$$

$$L({\omega }_0,\Gamma )=\frac{\Gamma }{\pi }\times$$

$$⟨{\Psi }_0|{\hat{\Theta }}^{†}\frac{1}{\hat{H}-E_0-{\omega }_0+i\Gamma }\frac{1}{\hat{H}-E_0-{\omega }_0-i\Gamma }\hat{\Theta }|{\Psi }_0⟩$$

$$=\frac{\Gamma }{\pi }⟨\tilde{\Psi }|\tilde{\Psi }⟩\text{mbox}.$$

$$(\hat{H}-z)|\tilde{\Psi }⟩=\hat{\Theta }|{\Psi }_0⟩,$$

$$R(\omega )={\omega }^{3/2}\exp \left(-\alpha \pi (Z-1)\sqrt{\frac{2\mu }{\omega }}\right)\sum _i^N{c}_i{e}^{-\frac{\omega }{\beta i}},$$

$$({\overline{H}}_N-z)|{\tilde{\Psi }}_R(z)⟩={\overline{\Theta }}_N|{\Psi }_0^R⟩,$$

$${\overline{\Theta }}_N=\exp (-\hat{T}){\hat{\Theta }}_N\exp (\hat{T}).$$

$$|{\tilde{\Psi }}_R(z)⟩$$

$$\hat{\mathcal{R}}(z)$$

$$+$$

$${\sigma }_{\gamma }(\omega )=4{\pi }^2\alpha \omega R(\omega ),$$

$$\hat{\Theta }=\sum _k^A\left({\mathbf{r}}_k-{\mathbf{R}}_{\mathrm{cm}}\right)\left(\frac{1+{\tau }_k^3}{2}\right).$$

$$S_n={\int }_{{\omega }_{th}}^{\infty }d\omega R(\omega ){\omega }^n,$$

$${\alpha }_D=2\alpha {\int }_{{\omega }_{th}}^{\infty }d\omega \frac{R(\omega )}{\omega }.$$

$${\alpha }_D=⟨{\Psi }_0|{\hat{\Theta }}^{†}\frac{1}{\hat{H}-E_0}\hat{\Theta }|{\Psi }_0⟩.$$

$$L({\omega }_0,\Gamma \to 0)=\int R(\omega )\delta (\omega -{\omega }_0)d\omega =R({\omega }_0),$$

$${\alpha }_D=2\alpha \int \frac{L({\omega }_0,\Gamma \to 0)}{{\omega }_0}d{\omega }_0.$$

$${\overline{\Theta }}_N$$

$$R_{\mathrm{skin}}=R_n-R_p.$$

$${\delta }_{\mathcal{O}}=\sqrt{({\delta }_{\mathcal{O}}^{\hslash \omega }{)}^2+({\delta }_{\mathcal{O}}^{\text{CC}}{)}^2},$$

$${\delta }_{\mathcal{O}}^{\hslash \omega }$$

Full text

††offprints:

1 Introduction

Electromagnetic probes are invaluable tools to study the nature of composite quantum-mechanical systems, such as nuclei. Due to the small value of the electromagnetic coupling constant $\alpha$, one can cleanly relate measured cross sections to properties of the composite system via perturbation theory, leading to a more complete understanding of the internal dynamics. Indeed, electromagnetic probes have historically enabled important discoveries regarding the nucleus and the strong dynamics governing its multifaceted properties. Most notably, the study of photonuclear reactions lead to the discovery of giant dipole resonances and to their interpretation in terms of collective modes GoT48 ; steinwedel1950 . A complete body of data has been collected over the past decades for stable nuclei, and some selected studies were performed even on unstable nuclei, leading, e.g., to the discovery of pygmy resonances, see, e.g., Ref. Bracco2019 and references therein. But where do we stand with the theory today?

One of the goals of modern nuclear theory is to be able to explain nuclear phenomena starting from protons and neutrons as degrees of freedom and by connecting their interactions to quantum chromodynamics via the use of chiral effective field theories Weinberg90 ; Epelbaum09 ; Machleidt11 ; Epelbaum12 or with other more traditional potentials Wiringa1995 . This research path is called “ab initio approach” in low-energy nuclear physics, see, e.g., Refs. Leidemann12 ; Bacca:2014tla ; Hebeler2015 . What is meant by that is that, for a given interaction Hamiltonian, the quantum-mechanical problem of protons and neutrons interacting with each other is either solved exactly, or within controlled approximations morten .

In the last years, ab initio computations in nuclear physics have advanced tremendously. While until one or two decades ago it was possible to only deal with very few nucleons, today one can reach even mass number $A\sim 100$ hagen2016b ; simonis2017 ; Sn ; gysbers2019 and above. The resulting growth of first principle calculations is well captured in Fig. 1, where we display the trend for ab initio calculations for the nuclear many-body problem as a function of $A$.

In early decades the progress was linear in $A$, mostly due to exponentially expensive algorithms implemented on machines with exponentially growing computing power. More recently, newly developed polynomial scaling algorithms, such as coupled-cluster theory hagen2014 , made it possible to exploit the exponential growth in computing power, pushing the boundaries of first principle calculations. However, one has to note that this fast growing progress regards mostly the computation of bulk properties such as binding energies and radii. Where does this ab initio approach stand with respect to the computation of electroweak reactions with nuclei? Or, in other words, are we able to describe the above mentioned collective modes in a microscopic way, starting from protons and neutrons interacting through realistic forces?

Only recently we have paved the road to a qualitative and quantitative description of electromagnetic reactions with nuclei from first principles. In this work, we will review the progress made in this respect thanks to the introduction of a new technique obtained by merging coupled-cluster theory hagen2014 with the Lorentz integral transform approach efros1994 ; Efl07 , that led to a novel computational tool to address electromagnetic reactions and related observables in an ab initio fashion Bacca:2013dma .

The paper is structured in the following way. We will first describe the computational tools in Section 2, then, after a brief explanation of the chiral interactions in Section 3, we will present results in Section 4. Finally, we will draw our conclusions and present an outlook for the future in Section 5.

2 Computational tools

2.1 Coupled-cluster theory

Coupled-cluster theory, as originally introduced by Coester and Kümmel coester1960 , is aimed at solving the Schrödinger equation for a many-body system. For a given Hamiltonian $\hat{H}$ describing the quantum system, one assumes that the correlated many-body wave function can be written with an exponential ansatz

[TABLE]

where $|\Phi_{0}\rangle$ is a Slater determinant of any kind. The operator $\hat{T}$, typically expanded in $n$-particle–$n$-hole ($np$–$nh$) excitations (or clusters) $\hat{T}=\hat{T}_{1}+\hat{T}_{2}+\dots+\hat{T}_{A}$, is responsible for introducing correlations. Using the formalism of second quantization the correlation operators can be written as

[TABLE]

where indices $i,j,k,\dots$ label occupied single–particle (hole) states in the reference Slater determinant, whereas the $a,b,c,\dots$ indicate unoccupied (particle) states.

The many-body Schrödinger equation for the ground state becomes

[TABLE]

employing the similarity transformed Hamiltonian

[TABLE]

where $\hat{H}_{N}$ is normal ordered with respect to the reference Slater determinant. The amplitudes of the $\hat{T}$ operator, such as $t_{i}^{a}$ , $t_{ij}^{ab}$, $t_{ijk}^{abc}$, etc., are found by solving the non-linear equations obtained by

[TABLE]

Here, $|\Phi_{i}^{a}\rangle$, $|\Phi_{ij}^{ab}\rangle$ and $|\Phi_{ijk}^{abc}\rangle$ are Slater determinants constructed as $1p$–$1h$, $2p$–$2h$, $3p$–$3h$, $\dots$ excitations on top of the reference state, respectively hagen2014 ; shavittbartlett2009 .

Coupled-cluster theory is exact when the expansion of the $\hat{T}$ operator is performed up to $Ap$–$Ah$ excitations. However, due to the exponential ansatz of Eq. (11), even when truncations are introduced, the result is very close to the exact one. For closed (sub-) shell nuclei the coupled-cluster method truncated at the $2p$–$2h$ level, at the so called coupled-cluster singles and doubles level – labeled with D in this work–, captures about $90\%$ of the full correlation energy, while adding triples excitations, about 99$\%$ of the correlation energy is accounted for bartlett2007 ; hagen2009b . The advantage of the method is that it scales polynomial with increasing system size (by system size we mean a measure of the number of particles and the size of the employed basis). For example, in the D approximation, the algorithm scales as $n_{o}^{2}n_{u}^{4}$, where $n_{o}$ and $n_{u}$ are the number of occupied and unoccupied orbitals, respectively. Coupled-cluster theory has been successfully applied to study properties of closed (sub-) shell nuclei and their neighbors (see for example Ref. hagen2014 for a review, and references therein).

2.2 Lorentz integral transform method

While enormous progress has been made in first-principles computations of ground-state properties of nuclei with increasing mass number $A$, advances in the calculation of electromagnetic reactions with nuclei have been slower, because of the additional challenges one has to face.

Electromagnetic cross sections are typically related to the nuclear response function, defined as

[TABLE]

Here, $|\Psi_{0}\rangle$ and $|\Psi_{n}\rangle$ denote ground and final state wave functions of the nucleus, and $E_{0}$ and $E_{n}$ are their respective energies, with $\omega=E_{n}-E_{0}$. The operator ${\hat{\Theta}}$ is a generic electromagnetic operator. One has to note that the $\sum_{n}$ indicates both the sum over discrete states and an integration over continuum Hamiltonian eigenstates. The calculation of the latter is the main bottleneck in the computation of response functions and thus of electromagnetic cross sections. In particular, the calculation of excited states in the continuum for medium-mass nuclei constitutes an open problem. At a given energy, the wave function is composed by many different channels, corresponding to all possible partitions into different fragments, which are difficult to calculate.

A method that allows to circumvent this issue is the Lorentz integral transform (LIT) efros1994 ; Efl07 , which reformulates the problem in such a way that the explicit knowledge of all $|\Psi_{n}\rangle$ in the continuum is not necessary. The Lorentz integral transform is defined as

[TABLE]

where $\omega_{0}$ and $\Gamma$ are parameters, with $\Gamma>0$. By substituting $R(\omega)$ in Eq. (7) with the expression from Eq. (6) and using the completeness relation of the Hamiltonian eigenstates,

[TABLE]

one obtains

[TABLE]

In this way, the LIT of the response function is basically proportional to the squared norm of the state $|\widetilde{\Psi}\rangle$. This state is found as solution of the Schrödinger-like equation

[TABLE]

where $z=E_{0}+\omega_{0}+i\Gamma$, for different values of the parameters $\omega_{0}$ and $\Gamma$. Because of the fact that ${L}(\omega_{0},\Gamma)$ is finite, the unique solution $|\widetilde{\Psi}\rangle$ of Eq. (10) has the same asymptotic boundary conditions as a bound state. Thus, only bound-state methods are required to solve this equation.

In the LIT approach one first computes ${L}(\omega_{0},\Gamma)$ in a direct way, without requiring the knowledge of $R(\omega)$. In a second step, the response function is obtained from a numerical inversion of ${L}(\omega_{0},\Gamma)$ efros1999 ; andreasi2005 . The typical inversion procedure is based on a least-squares fit. First, we make an ansatz for the shape of the response function, e.g., as

[TABLE]

where the exponential prefactor is a Gamow factor, assuming that the first channel is the one proton knock out, so that the remaining nucleus has $(Z-1)$ protons. Here $\mu\sim\frac{A-1}{A}m$ is the reduced mass with $m$ being the nucleon mass. The least-squares fit is optimizing the coefficients $c_{i}$ so that the LIT of Eq. (11) is coinciding with the calculated one. The coefficient $\beta$ is a non-linear fit parameter which is also varied in the fit procedure. Typically, one can change the ansatz in Eq. (11), vary $N$ and perform the inversion of LITs at different $\Gamma$ values to obtain an uncertainty of the inversion procedure.

For few-body systems, where a direct calculation of $|\Psi_{n}\rangle$ even in the continuum case is possible, it has been shown that the LIT method leads to an exact response function lapiana2000 ; golak2002 with the full final state interaction included. Because the relevant equation to solve, Eq. (10), is a bound-state equation, this method essentially circumvents the obstacle of the continuum calculation, the price to pay being that a very good precision in the calculation of $L(\omega_{0},\Gamma)$ is needed in order to obtain a stable inversion Efl07 .

The application of the LIT method used in conjunction with hyperspherical harmonics expansions to solve Eq. (10) allowed, e.g., to perform studies of the photodisintegration of the six- and seven-body nuclei bacca2002 ; BaB04 ; BaA04 .

2.3 Merging coupled-cluster theory with the Lorentz integral transform method

Predictive ab initio calculations of electromagnetic reactions have traditionally been limited to relatively light mass number, as discussed above. Medium-mass and heavy nuclei are typically studied with other theories, such as mean field based approaches and density functional theory Erler2011 ; Nakatsukasa2012 ; Piekarewicz2012 ; RocaMaza_Paar , which, despite being extremely useful, have a less direct connection to quantum chromodynamics.

To surpass previous limitations of the ab initio approach, we have merged the advantage of the LIT method of reducing the continuum problem to the solution of a bound-state equation with the mild computational scaling that characterizes coupled-cluster theory with increasing mass number. This led to the introduction of a new technique, which we call LIT-CC Bacca:2013dma and which essentially is a coupled-cluster formulation of the LIT method.

Given that in coupled-cluster theory one introduces the exponential ansatz and then one works with similarity transformation, in the LIT-CC method Eq. (10) becomes simply

[TABLE]

where $|\Psi_{0}^{R}\rangle\equiv|\Phi_{0}\rangle$ is the right ground state, while $\overline{\Theta}_{N}$ is the similarity transformed normal-ordered electromagnetic operator

[TABLE]

The solution $|\tilde{\Psi}_{R}(z)\rangle$ of Eq. (12) is found as linear superposition of particle-hole excitations on top of the reference Slater determinant as

[TABLE]

where

[TABLE]

These operators build up particle–hole ($p$–$h$) excitations and analogous expressions exist for the bra states.

If the expansions in the $\hat{T}$ and the $\hat{{\cal R}}$ operators are performed up to $Ap$–$Ah$, then the theory is exact. However, in practical applications one has to truncate these expansions. The most common approximation is the singles and doubles scheme D, for both the ground state and the excited states. In this work we will investigate different approximation schemes in the ground and excited states. In order to keep the notation concise we therefore denote each scheme with a pair of labels (separated by a ‘$/$’ symbol), with the largest order of correlation included in the ground state on the left, and the largest order of correlation included in the excited states on the right. For example, when truncating $\hat{T}=\hat{T}_{1}+\hat{T}_{2}$ and $\hat{{\mathcal{R}}}=r_{0}+\hat{{\mathcal{R}}}_{1}+\hat{{\mathcal{R}}}_{2}$, we will denote the calculation with D/D.

3 Interactions from $\chi$EFT

In the last decades we have observed the emergence, and its increased application, of chiral effective field theories ($\chi$EFTs)  Weinberg90 ; Epelbaum09 ; Machleidt11 ; Epelbaum12 to systematically derive the interactions of nucleons among themselves and with external electroweak probes. This approach allows to maintain a deeper connection to the underlying fundamental theory of quantum chromodynamics (QCD). Effective Lagrangians, expressed in terms of nucleons and pions, are constructed so as to preserve all symmetries, in particular the chiral symmetry, characterizing QCD in the limit of vanishing quark masses. The expansion in powers of $(Q/\Lambda_{\chi})^{\nu}$ is based on a separation of scales, where $Q$ is the low momentum characterizing low energy nuclear physics and $\Lambda_{\chi}\sim 1$ GeV is the chiral-symmetry breaking scale. The coefficients of the expansion are called Low-Energy Constants (LECs). They encapsulate the unresolved short-range physics and are adjusted to experimental data. This then makes it possible to predict nuclear observables to any degree $\nu$ of desired accuracy, with an associated theoretical error roughly given by $(Q/\Lambda_{\chi})^{(\nu+1)}$.

In this approach, see, e.g., Refs. Epelbaum09 ; Epelbaum12 ; Machleidt11 , three-nucleon (3N) forces and higher-body forces arise naturally and consistently with two-nucleon (NN) interactions. As such, they play an important role in consistent calculations. In this paper, we will first show results obtained only with NN forces Entem03 and then also present results using NN$+$3N forces. On the one hand we use a set of NN$+$3N Hamiltonians starting from the next-to-next-to-next-to-leading order (N3LO) NN potential of Ref. Entem03 evolved to lower resolution scales using the similarity renormalization group bogner2007 . These low-momentum interactions are then supplemented with a non-local 3N force at next-to-next-to-leading order (N2LO), adjusting the 3N LECs $c_{D}$ and $c_{E}$ to reproduce the 3H binding energy and the 4He charge radius (for further details see Ref. Hebeler2011 ). On the other hand we use the NNLOsat NN$+$3N Hamiltonian Ekstroem2015 which was adjusted to reproduce few-body observables as well as binding energies and radii in selected nuclei up to mass number $A\approx 25$. For this latter case, it is to be noted that the NNLOsat potential has been fit to experimental data using the $\Lambda$-CCSD(T) approximation in a model space consisting of 15 oscillator shells and a frequency of 22 MeV. Here, we do not keep these model-space parameters fixed, but vary them in a truly ab initio spirit.

The advantage of using a variety of interactions is that, by doing so, we can give an estimate of systematic uncertainties of the employed Hamiltonians. It also allows us to study whether different observables are correlated, which in turn could be used to make predictions for relevant quantities. Recent examples include the neutron radius and dipole polarizability of 48Ca Hagen2016 , the $2^{+}$ excited state in 78Ni hagen2016b , and electromagnetic transitions in light nuclei calci2016 .

In all the calculations we will show, we use a model-space truncation in $N_{\rm max}$, the number of harmonic-oscillator shells. While for NN only we reach model-space sizes up to $N_{\rm max}=18$, the truncation for NN$+$3N calculations is $N_{\rm max}=14$, if not otherwise specified. Additionally, the 3N matrix elements are truncated in the sum of the three-particle energies with typically $E_{\rm 3max}\leq 16$. The coupled-cluster computations start from a Hartree-Fock reference state. For NN$+$3N calculations the 3N contributions are included in normal-ordered two-body approximation, discarding residual 3N forces. This approximation is valid for light and medium-mass nuclei as shown in Refs. Hagen:2007ew ; Roth:2011vt .

4 Results

In this section we will present results obtained with the LIT-CC method. We will first focus on the dipole response function and photodisintegration cross section in Subsection 4.1. Then, we will show our studies of sum rules concentrating on the electric dipole polarizability $\alpha_{D}$ in Subsection 4.2. Finally, we will address the role of triples corrections in Subsection 4.3 and revisit the correlations among $\alpha_{D}$ and nuclear radii for 48Ca in Subsection 4.4.

4.1 Dipole response functions

The first electromagnetic reaction observables for which we exploited the power of coupled-cluster theory to address medium-mass nuclei has been the photodisintegration cross section. In the unretarded dipole approximation valid at energies below the pion-production threshold, the photodisintegration cross section can be written as

[TABLE]

where $\omega$ is the excitation energy and $R(\omega)$ is the dipole response function, basically Eq. (6) where the electromagnetic operator is the translationally invariant dipole

[TABLE]

Here ${\bf r}_{k}$ and ${\bf R}_{\rm cm}$ are the coordinates of the $k$-th particle and the center-of-mass, respectively, while $(1+\tau^{3}_{k})/2$ defines the projection operator on the $Z$ protons, with $\tau^{3}_{k}$ being the third component of the $k$-th nucleon isospin.

The photodisintegration cross section of Eq. (16) has been calculated with the LIT method for a variety of light nuclei, from deuteron lapiana2000 to 7Li BaA04 , where the most extensive studies have been performed on 4He PRL_4He . In these cases the Schrödinger-like equation was solved using exact few-body techniques. Before the new LIT-CC method is used to study medium-mass nuclei, it is useful to benchmark it against the above mentioned exact calculations. 4He presents itself as an interesting case study, because it is a closed-shell nucleus, where the LIT-CC method can be most easily applied, and it has been substantially studied both from the theoretical and the experimental point of view. From the theory side, it is particularly instructive to test the D/D approximation against exact calculations performed with effective interaction hyperspherical harmonics Barnea2000 ; Barnea2001 . In particular, we have performed such a comparison by using a two-body interaction derived in chiral effective field theory at N3LO Entem03 .

Results of this comparison are shown in Fig. 2, where we present the response function obtained after the inversion of the LIT. In the inversion, we impose the response function to be zero before the threshold energy $\omega_{\rm th}$ which is the difference between the binding energy of 4He and 3H, with $\gamma$ + 4He $\rightarrow$ 3H + $p$ being the first open reaction channel in the photodisintegration process. Because, with the N3LO potential the binding energies of 4He and 3H are not correctly reproduced due to the missing 3N forces, a simple way to correct for that is to shift the curves from the theoretical threshold $\omega_{\rm th}=17.54$ MeV to the experimental one of 19.82 MeV. In this way, one can focus on the shape of the cross section in the comparison to data. Three-nucleon forces were included, e.g., in Ref. PRL_4He , and their effect is to correct the threshold energy and decrease the peak height by a few percent.

The coupled-cluster D/D result is shown by the thick curve. The thickness of the curve is obtained from inverting the LIT with $\Gamma=10$ MeV and $\Gamma=20$ MeV and varying $N$ in Eq. (11). The EIHH results are instead represented by the thin curve. In this case, by inverting LITs with $\Gamma=10$ and 20 MeV, the two results overlap exactly. Interestingly, the D/D response function is close to the EIHH result, proving to be a very good approximation. Only small deviations for energies between about $\omega=30$ and $50$ MeV are seen. These are though much smaller than the uncertainties of most of the experimental data available for the dipole response function.

In Fig. 2, we show a selected set of data for comparison, see Ref. Bacca:2014tla for a more extensive discussion. The Arkatov et al. Arkatov data from the 70s cover the broadest energy range. They are in relatively good agreement with more recent data by Raut et al. Raut:2012zz ; Tornow:2012zz , where exclusive experiments were carried out measuring ($\gamma,p$) and ($\gamma,n$) at the High Intensity Gamma-Ray source. We summed up the two channels and show just the points in the energy range below the three-body break-up. These first two sets of data obtained from traditional photoabsorption experiments, agree rather well with dipole response function data from Nakayama et al. Nakayama obtained via the study of the 4He(7Li,7Be) reaction. A completely different trend is instead shown by the Shima et al. data from Ref. Shima:2005ix , obtained with a quasi-monoenergetic photon beam and a time projection chamber, where a simultaneous measurement of both the 4He$(\gamma,n)^{3}$He and 4He$(\gamma,p)^{3}$H reactions has been performed. Obviously the uncertainties in the experimental data (and the disagreement of about a factor of 2 of the Shima et al. data), are much larger than the slight differences we observed between the EIHH and D/D coupled-cluster calculation. We impute this latter difference to the missing triples and quadruple correlations and conclude that they can be safely neglected as their effect is small.

After we successfully benchmarked the new method with exact hyperspherical harmonics on 4He, we exploit the mild computational scaling of the LIT-CC method with increasing mass number and investigate medium-mass nuclei, for the first time with an ab initio approach.

In Fig. 3, we show the dipole response function of 16O calculated with the LIT-CC method at the D/D approximation level, using the same NN interaction derived in chiral effective field theory at N3LO Entem03 . The curve is shown starting from the experimental threshold $\omega_{th}=12.1$ MeV. The band thickness is obtained by inverting the LIT with width $\Gamma=10$ MeV and by varying the number of basis functions employed in the inversion. Similar results are obtained by inverting the LIT at $\Gamma=20$ MeV. We compare the theoretical results with experimental data by Ahrens et al. Ahrens75 , who measured the total photoabsorption cross section $\sigma_{\gamma}$ on an oxygen target with natural abundance ($99.762\%$ 16O) with an attenuation method. We also compare to a more recent evaluation by Ishkhanov et al. ishkhanov2004 . We observe that the theoretical result is smeared compared to data, but overall the total dipole strength is correctly reproduced, as well as the bulk of the strength is in the right energy range. It is to note though that with this two-body interaction and with the D/D approximation, we do not see the structures at around 10 MeV, which are found in experiment.

Another advantage of coupled-cluster theory is that one can also compute neutron-rich nuclei. In particular, in the oxygen isotopic chain, a few of its neutron-rich isotopes have been studied at rare isotope beam facilities. The comparison of stable and unstable nuclei can provide key information about nuclear forces at the extremes of matter. Thus, it is interesting to apply these ab initio methods to exotic nuclei. We will focus on 22O, which is a closed sub-shell nucleus, for which the LIT-CC method can be applied.

In Fig. 4 we show the photoabsorption cross section of 22O computed in the D/D approximation PRC2014 with the chiral N3LO NN force Entem03 . The width of the curve is obtained by performing several inversions, for $\Gamma=20,10$ and 5 MeV. We compare our results to experimental data taken at GSI by Leistenschneider et al. leistenschneider2001 . These are obtained from a Coulomb excitation experiment and are turned into the equivalent of a photoabsorption cross section. The experimental data show a small peak at low energy. This structure is often named pygmy dipole resonance and was experimentally observed in neutron-rich nuclei Bracco2019 . Interestingly, our first principle calculation also presents a low-energy peak. The curve is shifted to start from the experimental threshold energy, as done for the other nuclei shown here. We remind the reader that the implemented two-body interaction was tuned only on two-nucleon data. Despite that, we see the emergence of a pronounced substructure at low energy as showed by the data. At higher energies, the D/D results are larger than the data. This is expected because, while the experiment measured a semi-inclusive cross section the theoretical calculation is for an inclusive cross section, where proton emission channels are included.

The heaviest nucleus for which we calculated a response function with the LIT-CC method is 40Ca. In Fig. 5 we show the results in the D/D approximation obtained using the chiral N3LO NN force Entem03 used for all other nuclei so far. The width of the curve in Fig. 5 is obtained by inverting LITs with different $\Gamma$ parameters. We compare it to the measured cross section by Ahrens et al. Ahrens75 , where a natural calcium target was used with a photoabsorption attenuation method. The data show a very pronounced peak, referred to as the giant dipole resonance and located at an excitation energy of around 20 MeV. This structure is well reproduced by the LIT-CC results.

As mentioned in the Introduction, first interpretations of such resonances, were given in terms of collective models GoT48 ; steinwedel1950 . Now, with the advent of novel computational techniques such as the LIT-CC method we are able to show the emergence of these collective modes, both for the giant as well as for the pygmy resonance, from first principles.

Clearly, more work needs to be done to better assess the uncertainties included in the theoretical framework, as well as to include 3N forces. Nevertheless, the results shown in this Section constitute the first successful attempt to describe dipole response functions with ab initio methods. The main difficulty of this approach is that one needs to be able to compute the LIT very precisely, in order to be able to stably invert it. One limiting factor to achieve a sub-percentage convergence in the LIT calculation, is the availability of large model spaces, i.e., of computer memory and of matrix elements of the starting interaction. In all the above shown computations, the precision we could reach employing model spaces of $N_{\rm max}=18$ was of the order of one to two percent, and that is the main reason why we obtain quite thick bands in the inversion procedure.

4.2 Electric dipole polarizability

If one is able to calculate the dipole response function, one can then easily compute its existing sum rules as

[TABLE]

where $n$ is an integer number. However, in several instances, it is easier to directly compute some selected sum rules than the whole response function, as the former may be computed as the expectation value of an operator on the ground state.

A particularly interesting sum rule of the dipole response function is the electric dipole polarizability $\alpha_{D}$, which is defined as

[TABLE]

It is basically an inverse energy weighted sum rule of the dipole response function. As written in Eq. (19), it is evident that $\alpha_{D}$ contains the information on $R(\omega)$ at all energies $\omega$, including those in the continuum. A calculation of $\alpha_{D}$ would then require to be able to solve the many-body scattering problem in the continuum or to use alternative approaches as the LIT method. Indeed, for all the nuclei discussed in Subsection 4.1 it is easy to compute $\alpha_{D}$ by simply performing the integral in Eq. (19).

On the other hand, using the completeness relations it is possible to rewrite $\alpha_{D}$ as

[TABLE]

So, if one is able to deal with the operator $\hat{\Theta}^{\dagger}\frac{1}{\hat{H}-E_{0}}\hat{\Theta}$, the calculation of $\alpha_{D}$ reduces to an expectation value of the ground state, where no excited state is involved at all. However, that operator is not easy to write down, so practically what one does is to include completeness of the eigenstates (or a bound-state representation of them) and use an eigenrepresentation of the Hamiltonian. At this point, one can then see a relation between this equation with the second line of Eq. (9). In fact, because in the limit $\Gamma\rightarrow 0$ the Lorentzian kernel becomes a delta function

[TABLE]

$\alpha_{D}$ can be computed from the LIT as

[TABLE]

Equation (21) is nothing else than a discretized response function and Eq. (22) is just an integral of that with an inverse energy weight. In Ref. Miorelli2016 we showed that this method is equivalent to computing first the response function with the LIT approach and then directly integrate the response function. The advantage of using Eq. (22) is that it does not require to invert $L(\omega_{0},\Gamma)$ and have, consequently, an additional numerical error. For this reason, if one is interested just in the electric dipole polarizability, it is preferable to use this method. In a slightly different spirit, one can understand that it is more complicated to compute the whole response function $R(\omega)$ because knowing the response functions means knowing all of its existing sum rules, not just one.

The $L(\omega_{0},\Gamma)$ in the LIT-CC approach is always computed using the Lanczos algorithm for non-symmetric matrices PRC2014 and in Ref. Miorelli2016 we showed that $\alpha_{D}$ can be computed as continued fraction of the Lanczos coefficients Lanczos50 . With this technology at hand, we can now explore the dependence of the polarizability on the employed nuclear interaction. For this purpose, we will also use the potentials and normal-ordered 3N forces introduced in Section 3. For the notation of the potentials we follow Ref. Hebeler2011 .

In Table 1 we show results obtained within the D/D approximation for the nuclear charge radius $R_{\rm ch}$ and for $\alpha_{D}$ for 16O and 40Ca. Calculations were performed with a model-space size of $N_{\rm max}=14$, for which we show the central value obtained with the optimal harmonic-oscillator frequency. The uncertainty associated with these D/D calculations were estimated to be of about 1$\%$ for the radius and 2$\%$ for the polarizability Miorelli2016 .

One can readily see that, as expected, $\alpha_{D}$ and $R_{\rm ch}$ are strongly correlated. At the bottom of the table we also show the experimental results from Angeli and Marinova angeli2013 for the radius and from Ahrens et al. Ahrens75 for the polarizability. When integrating the experimental cross section from Ahrens et al. Ahrens75 , a lower value for the polarizability was obtained in Refs. Miorelli2018 ; Birkhan17 . The agreement with data is particularly good for the interaction NNLOsat, obviously so for the radius of 16O, since it was fit to reproduce the experimental value, while the polarizability and both observables for 40Ca are predictions of this potential. As shown in Ref. Miorelli2016 , if one used two-body forces only, one would still see a strong correlation, but generally an underestimation of the experimental value for both the charge radius and the polarizability.

We remind the reader that, while in Ref. Miorelli2016 we were able to compute the full $L(\omega_{0},\Gamma)$ and invert it for 4He and 16O with the NNLOsat potential, the inversion procedure was found to be much more unstable than for the calculations with two-body forces only shown in Subsection 4.1, primarily due to the fact that we are limited in the model-space size to about $N_{\rm max}=14$ (as opposed to $N_{\rm max}=18$ for two-body forces alone), which translate in a larger uncertainty in the $L(\omega_{0},\Gamma)$ and consequently on the inversion. Thus, from now on, we will continue our discussion on the electric dipole polarizability only.

4.3 Adding coupled-cluster triples corrections

Coupled-cluster theory is a systematically improvable method. So far we showed results obtained within the coupled-cluster singles and doubles approximation and have assessed uncertainties due to the model-space expansion and residual dependence on the harmonic-oscillator frequency, and by comparing to exact results in light nuclei Hagen2016 ; Ekstroem2015 . In order to improve the results and more rigorously assess the uncertainty associated to the coupled-cluster expansion, the next natural step is to include triples corrections. We include leading-order $3p$–$3h$ excitations using the so called CCSDT-1 iterative triples approach lee1984 which is a good approximation to the full triples. In the ground state, it typically accounts for about $99\%$ of the correlation energy and includes the leading-order contribution $\left(\hat{H}_{N}\hat{T}_{2}\right)_{C}$ (here the index $C$ denotes connected terms shavittbartlett2009 ) to the $\hat{T}_{3}$ amplitudes with an energy denominator given by the Hartree-Fock single-particle energies, while all $\hat{T}_{3}$ contributions to the $\hat{T}_{1}$ and $\hat{T}_{2}$ amplitudes are fully included. We also solve for the corresponding left ground state in the CCSDT-1 following Ref. watts1995 . Here, we denote this approximation with T-1, consistently with Ref. Miorelli2018 . Analogously, for the equation of motion and excited states we will add leading-order $3p$–$3h$ excitations in the T-1 approach watts1995 ; jansen2016 . When T-1 is used both in the ground state and excited states, we will denote this scheme as T-1/T-1.

When including $\hat{T}_{3}$ contribution in the similarity transformation of a normal-ordered one-body operator $\hat{\Theta}_{N}$ one has

[TABLE]

where $\overline{\Theta}^{D}_{N}$ is the similarity-transformed operator in the D approximation. In Ref. Miorelli2018 we have shown that $\overline{\Theta}_{N}$ can safely be approximated with $\overline{\Theta}^{D}_{N}$, dramatically simplifying the calculations. Thus, we will only report on results obtained in this scheme.

At this point, it is interesting to compare calculations for the electric dipole polarizability with triples with those only including singles and doubles. Furthermore, because in the computation of $\alpha_{D}$ we need to specify an approximation scheme for the ground state and for the excited states, we also explore the case where T-1 triples are included only in the ground state, and denote this scheme with T-1/D. Obviously, this calculation is less computationally intensive than the T-1/T-1 case.

In Fig. 6, we show $\alpha_{D}$ for 4He and 16O with the D/D, T-1/D and T-1/T-1 schemes. We employ the chiral two-body force at N3LO from Entem03 for 4He because we want to benchmark the various approximation schemes against virtually exact EIHH, which cannot easily employ the NNLOsat Hamiltonian due to the non-locality of the 3N forces. For 16O instead we use the NNLOsat Hamiltonian Ekstroem2015 , since our goal is to compare with experimental data from Ahrens et al. Ahrens75 .

For 4He we see that the T-1/D calculation agrees best with the hyperspherical harmonics result shown by the black line, surprisingly better than the T-1/T-1 approach. The width of the bands reflects the residual harmonic-oscillator dependence for the largest model space $N_{\rm max}=14$ and no cut has been done on the $3p$–$3h$ configurations. Overall, the effect of $3p$–$3h$ excitations in 4He is small, at the level of 1 $\%$, which is comparable to the uncertainty obtained from the model-space variation and dependence on the harmonic-oscillator frequency.

For 16O, while the D/D value is obtained at $N_{\rm max}=12$, when adding $3p$–$3h$ configurations we used a model space of $N_{\rm max}=12$ and $E^{F}_{\rm 3max}=14$ for T-1/T-1 and T-1/D. Here, $E^{F}_{\rm 3max}$ is an energy cut on the allowed $3p$–$3h$ excitations, defined as $|N_{a}-N_{F}|+|N_{b}-N_{F}|+|N_{c}-N_{F}|\leq E^{F}_{\rm 3max}$ and $|N_{i}-N_{F}|+|N_{j}-N_{F}|+|N_{k}-N_{F}|\leq E^{F}_{\rm 3max}$ with a given harmonic-oscillator shell $N_{p}$ and the harmonic-oscillator shell at the Fermi surface $N_{F}$. The bands are obtained by assigning a $2\%$ uncertainty, accounting for the combined uncertainty from the $E^{F}_{\rm 3max}$ cut and the residual harmonic-oscillator dependence. Correlations arising from $3p$–$3h$ excitations reduce the size of $\alpha_{D}$, in this case both for T-1/T-1 and T-1/D, by $8\%$ and $10\%$, respectively. Interestingly, T-1/T-1 and T-1/D are very close to each other showing that including triples corrections into the ground state is more important than including them in the excited states. The few percent difference between the T-1/D and T-1/T-1 results can be taken as an estimate of neglected higher-order correlations. In comparison with the experimental data from Ahrens75 shown by the grey band, we see that the addition of triples leads to a further deviation of $\alpha_{D}$ with respect to the experimental data, which agreed better in the D/D approximation for this interaction. While the comparison to experiment has to be seen as a judgment on the Hamiltonian itself, one also has to keep in mind that the experimental value has been extracted from sum rules from photoabsorption data that may be prone to larger systematic uncertainties than those quoted, because it is difficult to estimate the role of multipoles beyond the dipole.

In Table 2 we include new results for 16O obtained within the T-1/D scheme for the electric dipole polarizability, computed with all the Hamiltonians as in Table 1. Interestingly, we find that with triples corrections our results compare better with the value of 0.4959 fm3 obtained by Raimondi and Barbieri Raimondi with the self consistent Green’s function method for the same NNLOsat interaction. This might be coincidental as explained in Ref. Miorelli2018 , since coupled-cluster theory in the simple D/S approximation also yields a similar value, namely 0.503 fm3. Overall, we find that the effect of coupled-cluster triples corrections is lower for the softer interactions and higher for harder interactions, reaching about 15$\%$ for the NNLOsat case. For all the interactions, we get values that slightly under-estimate the experimental data, which as mentioned before, could potentially also include the effect of higher multipoles.

4.4 Revisiting correlations in 48Ca

Recently, the first ab initio computation of the dipole polarizability and the neutron-skin thickness of 48Ca were provided Hagen2016 . 48Ca is a neutron-rich isotope of calcium, with 20 protons and 28 neutrons. Despite the fact that it is one of the candidates for neutrino-less double beta decay, its half life is so long, that it can be considered as a stable nucleus and studied in fixed target experiments. In particular, an important observable that is of great interest is the so-called neutron-skin thickness, which is defined as the difference of the root mean square distribution of neutrons, $R_{n}$, and of protons, $R_{p}$, in the nucleus as

[TABLE]

The neutron-skin thickness basically measures the extent of the neutron distribution in a nucleus. While nuclear charge (thus proton) distributions are traditionally easily measured, e.g., with elastic electron scattering, the distribution of neutrons is much more elusive and difficult to access directly. Presently, the cleanest way to measure this quantity is expected to be parity violating electron scattering, which is mediated by the weak force through the exchange of a $Z^{0}$ boson. Due to the fact that the proton-weak current is approximately zero, this probe is very sensitive to the location of the neutrons. There are plans to measure neutron skins at Jefferson Laboratory with the CREX CREX and PREX experiments PREX , as well as in Mainz with the MREX n-skin project @MESA MREX .

In Ref. Hagen2016 a study of the correlations of $\alpha_{D}$ with $R_{n}$ and $R_{\rm skin}$ was presented, that provided theoretical predictions for all the three above mentioned quantities targeted by modern experiments. The electric dipole polarizability of 48Ca was indeed recently measured with $(p,p^{\prime})$ reactions in Osaka and was found to be in rather good agreement with coupled-cluster calculations Birkhan17 , albeit somewhat smaller.

The radii and polarizability calculations in Refs. Hagen2016 were performed in the D/D scheme. Here, we revisit our results using the recently developed $3p$–$3h$ technology Miorelli2018 . In Ref. Miorelli2018 we already computed $\alpha_{D}$ for 48Ca in the T-1/D scheme for a couple of interactions and found that its values is reduced, improving the agreement with the experimental measurement by Birkhan et al. Birkhan17 . Here, we supplement our results with a few more Hamiltonians and present complete results for radii and the dipole polarizability in Table 3 (with an increased number of decimal digits compared to the other tables in this work for clarity of the following discussion). On the one hand, we find that triples corrections only mildly affect the calculations of $R_{p}$ and $R_{n}$, leading to an increase in both quantities of less than $1\%$. Triples effects are slightly larger on $R_{p}$ than on $R_{n}$, and as a consequence $R_{\rm skin}$ calculated in the T-1/D scheme is smaller than $R_{\rm skin}$ calculated in the D/D approximation. On the other hand, triples correlations affect $\alpha_{D}$ quite visibly, as already pointed out in Ref. Miorelli2018 . In particular, for the harder interactions, their effect is of about $15\%$.

In Fig. 7 we plot these results, showing the proton radius as a function of the skin radius, the neutron-distribution radius and the electric dipole polarizability. The squares correspond to coupled-cluster calculations in the T-1/D scheme using four different parameterization of the $\chi$EFT interactions Hebeler2011 , while the circles correspond to the results obtained with the NNLOsat Ekstroem2015 potential. Each theory point is plotted with error bars that include both the residual $\hbar\omega$-dependence, as well as an estimate of the coupled-cluster truncation error. Below we briefly explain the recipe we use to estimate uncertainties.

The uncertainty $\delta_{\cal O}$ for any observable ${\cal O}$ ($R_{p},R_{n},R_{\rm skin}$ and $\alpha_{D}$) is computed as the quadrature of

[TABLE]

where $\delta^{\hbar\omega}_{\cal O}$ is an estimate of the uncertainty related to the residual harmonic oscillator $\hbar\omega$-dependence, while $\delta^{\text{CC}}_{\cal O}$ is an estimate of the coupled-cluster truncation uncertainty. For the latter, in case of the D/D calculations we take the difference from the T-1/D and D/D results, while in case of the T-1/D calculations, given that we do not have any higher order coupled-cluster calculation available, we take half of the above mentioned difference. More formally stated, the two uncertainties summed in quadrature are taken to be

[TABLE]

for the D/D calculations, and

[TABLE]

for the T-1/D calculations. The values are computed for the maximum available model space size of $N_{\mathrm{m}ax}=14$, using two neighboring frequencies around the optimal value, namely $\hbar\omega_{1}=16$ MeV and $\hbar\omega_{2}=12$ MeV. The uncertainties of $R_{\rm skin}$ are obtained taking the correlation between $R_{p}$ and $R_{n}$ into account using standard covariance theory. Clearly, this is a rough uncertainty estimate, but we find it sensible for the purpose of updating Ref. Hagen2016 .

From Fig. 7, it is clear that the strong correlation between $R_{p}$, $R_{n}$ and $\alpha_{D}$, which was observed in Hagen2016 , is confirmed, as well as we reinforce the fact that $R_{\rm skin}$ is almost constant with respect to the employed interaction. As done in Ref. Hagen2016 , one can exploit the correlation among observables with the fact that the proton distribution radius is known experimentally, and from the intersection between the correlation band and the horizontal line, one can draw constraints on the neutron-skin thickness and the electric dipole polarizability. In Fig. 7 we do not explicitly show these constraints, but rather highlight the linear correlation bands and the experimental constraints with their corresponding uncertainty. With respect to Ref. Hagen2016 , revisiting the calculations with leading triples, we obtain that the constraints on the dipole polarizability move from $2.19\leq\alpha_{D}\leq 2.60$ fm3 to $1.92\leq\alpha_{D}\leq 2.38$ fm3, thus getting closer to the experimental value of 2.07(22) fm3 Birkhan17 , while the constraints for the neutron-skin thickness go from $0.12\leq R_{\rm skin}\leq 0.15$ fm to $0.13\leq R_{\rm skin}\leq 0.16$ fm, thus not varying much. Despite the fact that $R_{\rm skin}$ becomes smaller with triples for a given model space parameter set, the width of the neutron-skin thickness constraint, as well as the central value of $R_{\rm skin}$ are slightly increased with respect to Hagen2016 . This is due to a combination of facts: the use of symmetric correlation bands; the use of five instead of six Hagen2016 different interactions, and the use of optimal harmonic-oscillator frequencies, where the convergence in $N_{max}$ is faster. The difference is anyway small and the updated 48Ca study further establishes that the ab initio prediction of the neutron-skin thickness is much smaller than the values obtained from selected density functional theory calculations, shown in Fig. 7 by the diamonds. Thus, it is worth to reiterate that it is important to have a clean, as model-independent as possible, experimental determination of $R_{\rm skin}$.

5 Conclusions

In this paper we review the recent progress made in the computation of electromagnetic response functions and related sum rules using a coupled-cluster theory formulation of the Lorentz integral transform method. We present photoabsorption cross sections of 16,22O and 40Ca, showing that we obtain a reasonable description of the experimental data already with two-body forces. We also review our calculations of the electric dipole polarizability, for which we use NN$+$3N Hamiltonians. We also present new results for the polarizability of 16O and 48Ca that include leading triples coupled-cluster correlations. For the 48Ca case, we revisit our previous studies of correlations among the electric dipole polarizability, the neutron-, proton- and skin-radius. We show that correlations still hold and allow to draw improved predictions for the value of the polarizability, which is in better agreement with experiment, as well as for the neutron-skin thickness, for which we corroborate the earlier finding that ab initio theory predicts a smaller neutron skin than density functional theory.

Acknowledgments.– This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through the Collaborative Research Center [The Low-Energy Frontier of the Standard Model (SFB 1044)]; by the Cluster of Excellence “Precision Physics, Fundamental Interactions, and Structure of Matter” (PRISMA+ EXC 2118/1) funded by DFG within the German Excellence Strategy (Project ID 39083149) and by the Office of Nuclear Physics, U.S. Department of Energy, under grants desc0018223 (NUCLEI SciDAC-4 collaboration) and by the Field Work Proposal ERKBP72 at Oak Ridge National Laboratory (ORNL). Computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. The new calculations presented in this work were also performed on “Mogon II” at Johannes Gutenberg-Universität in Mainz.