1 Effective field theory for systems near unitarity

Nuclear physics at very low energies hosts a fascinating emergent phenomenon: out of the tremendously complicated dynamics of quarks and gluons, governed by the strong interaction (Quantum Chromodynamics, QCD) that is highly nonperturbative in this regime, ultimately arise strikingly simple features for systems of few nucleons. It was realized many decades ago [1,2,3,4] that the low-energy two-nucleon system can be parameterized with a formula that has become famous, in nuclear physics and beyond, as the effective range expansion (ERE):

$$\begin{aligned} k\cot \delta _0(k) = {-}\frac{1}{a} + \frac{r}{2} k^2 + \cdots . \end{aligned}$$
(1)

Here \(\delta _0(k)\) denotes the S-wave scattering phase shift for two particles (here, nucleons in a single fixed spin configuration) with relative momentum k. The leading parameter in this expansion, called “scattering length” and denoted by a, governs the nucleon-nucleon (\(N\!N\)) cross section at low energies and completely determines it in the limit where the relative momentum k of the nucleons goes to zero:

$$\begin{aligned} \sigma = 4\pi a^2 + {\mathcal {O}}(k^2) . \end{aligned}$$
(2)

Empirically, in the \({}^3S_1\) (“t”) and \({}^1S_0\) (“s”) \(N\!N\) spin channels the values are known to be \(a_t\simeq 5.4~\mathrm {fm}\) and \(a_s\simeq {-}23.7~\mathrm {fm}\), respectively. Compared to the typical range of the nuclear interaction, set by one-pion exchange providing the longest-range component as \(R\sim M_\pi ^{-1} \simeq 1.4~\mathrm {fm}\), these scattering lengths are unnaturally large, \(a_{s,t}\gg R\), qualitatively. Through Eq. (2) this implies that the nuclear force is particularly strong as the energy goes to zero. The fact that this is so can be understood as an accidental “fine tuning” of the QCD parameters [5,6,7,8,9] (the quark masses, in particular) to be close to a critical point where the scattering lengths diverge. This point is called the “unitarity (or unitary) limit,” and it is the heart of the emergent simplicity mentioned at the outset.

Systems near the unitarity limit exhibit universal features. As the two-body scattering length becomes large, the details of the underlying interaction largely cease to matter and to a very good approximation the behavior of the system is determined qualitatively by the fact that a is large, and quantitatively by how large exactly it is. This phenomenon places low-energy nuclear systems in a common universality class with other systems near unitarity, such as cold atomic gases (where the scattering length can be tuned experimentally via Feshbach resonances [10]), nuclei with a halo/cluster structure [11], or certain mesons which can be interpreted as hadronic molecules [12].

In the two-body sector, universality relates scattering parameters to shallow bound and virtual states. This is a consequence of Eq. (1) and the principle of analyticity: the ERE provides an expansion of the S matrix, so whenever poles at complex momenta—in particular bound and virtual states, which reside at purely imaginary momenta—fall within the radius of convergence of the expansion, they are described by the same parameters. Since, schematically, the S-wave S matrix is given by \(1 + \mathrm {i}t\) with

$$\begin{aligned} t(k) \sim \frac{1}{k\cot \delta _0(k) -\mathrm {i}k} , \end{aligned}$$
(3)

the pole condition is \(\cot \delta _0(k) = \mathrm {i}\) for some \(k = \mathrm {i}\kappa \). Keeping only the first term in the ERE this gives \(\kappa = 1/a\) and one sees that for positive (negative) a one has a bound (virtual) state in the system. Having a sufficiently large \(\left| a\right| \) ensures that indeed these momentum scales lie within the radius of convergence of the ERE, which for nucleons is determined by the position of the pion cut, \(M_\pi /2\). For the two-nucleon system at the physical point one has the deuteron as a shallow bound state (\(B_D\simeq 2.224~\mathrm {MeV}\), with the difference to \(1/(M_Na_t^2) \approx 1.41~\mathrm {MeV}\) being due to range corrections) in the \({}^3S_1\) channel, and a very shallow virtual state at \(B_{N\!N^*}\simeq 0.068~\mathrm {MeV}\) (with a relatively small range correction since \(\left| a_s\right| \) is so large). In the unitarity limit, \(a_{s,t}\rightarrow \infty \), both of these poles become zero-energy S-wave states.

A more striking universal behavior is encountered for three and more particles: in the unitarity limit there exists an infinite tower of three-body bound states, geometrically spaced (the binding energy of each subsequent state is given by a fixed factor times the previous level) and accumulating at zero energy, a phenomenon that has become famous as the Efimov effect [13]. At large but finite scattering length the spectrum is cut off in the infrared due to the existence of a two-body pole in the S matrix. It was shown in Refs. [14,15,16] that for physical values of the \(N\!N\) scattering lengths the triton can be interpreted as the single remaining bound state of such an Efimov tower. More recently it was established in a model-independent way [17] that a virtual state in the three-nucleon (3N) system, known to exist for a long time [18, 19], is in fact an S-matrix pole that would be an excited Efimov state instead if nature were just slightly closer to the unitarity limit. This confirms a relation previously observed in a separable potential model [20].

The phenomenon continues at the four-body level. At unitarity, each three-boson Efimov state (with binding energy \(B_3\)) is associated with two four-boson states [21]. One of these is almost five times as deeply bound as the trimer, \(B_4/B_3\simeq 4.611\), while the other resides just below the particle-trimer threshold, \(B_{4^*}/B_3\simeq 1.002\) [22]. Universality implies that if the \(N\!N\) scattering lengths were infinite, the ground state of \({}^{4}\hbox {He}\) would be located at 4.6 times the binding energy of the triton (neglecting Coulomb and other isospin breaking effects). In nature, the ground state is at \(B_\alpha /B_H\simeq 3.66\), and there exists a \(0^+\) resonance state just above the proton-triton threshold, i.e., one has \(B_{\alpha ^*}/B_H\simeq 1.05\), where \(B_H\simeq 7.72~\mathrm {MeV}\) is the \({}^3\mathrm {He}\) binding energy, taken as reference here to at least partially account for isospin breaking effects. The closeness of these ratios to the unitarity-limit values is a strong indication that nature may be perturbatively close to unitarity for systems of at least four nucleons.

In the following, this work discusses how to construct an effective field theory (EFT) that captures all phenomena mentioned above. EFTs are a powerful tool widely used in modern theoretical physics. In nuclear physics they enable the consistent construction of nuclear forces systematically connected to QCD by choosing a “theoretical resolution” at which effective interactions between degrees of freedom appropriate for the energy scales of interest are constructed. The richness of nuclear phenomena implies that there are a number of different EFTs relevant for nuclear physics, forming the “tower” of theories that gives rise to the name of the topical issue this work is contributed to. A recent review of EFTs that use nucleons and mesons as degrees of freedom can be found in Ref. [23]. Here the focus is on setting up an EFT that systematically expands light nuclei around the unitarity limit, expanding on previous work considering the unitarity expansion [24, 25] by considering charge radii of light nuclei in addition to binding energies. The following sections present the setup and implementation of the unitarity expansion. After a discussion of the main results, Sect. 5 will give a summary and outlook to address the question where the unitarity expansion resides within the tower, or landscape, of nuclear EFTs.

As a variant of what has become known as “Pionless EFT,” the unitarity expansion is defined in terms of a Lagrange density

$$\begin{aligned} {\mathcal {L}}&= N^\dagger \left( \mathrm {i}{{{\mathcal {D}}}}_0+\frac{\varvec{{{\mathcal {D}}}}^2}{2M_N}\right) N \nonumber \\&\quad + \sum \nolimits _{{\mathbf {i}}}C_{0,{\mathbf {i}}} \left( N^T P_{{\mathbf {i}}} N\right) ^\dagger \left( N^T P_{{\mathbf {i}}} N\right) \nonumber \\&\quad + D_0 \left( N^\dagger N\right) ^3 + \cdots , \end{aligned}$$
(4)

where the notation has been taken over from Refs. [24, 26]. The degrees of freedom are nonrelativistic nucleon isospin doublets \(N=(p\;n)^T\), coupled to photon fields \(A_\mu \)via the covariant derivative \({\mathcal {D}}_\mu = \partial _\mu + \mathrm {i}eA_\mu (1+\tau _3)/2\), where e is the proton charge and \(\tau _a\) is used to label isospin Pauli matrices. Of the electromagnetic interactions, only the static Coulomb potential is relevant up to high orders (defined later), where corrections from transverse photons will eventually enter. The strong interaction is parameterized in Eq. (4) by the “low-energy constants (LECs)” \(C_{0,{\mathbf {i}}}\) and \(D_{0}\), defining contact (zero-range) interactions without derivatives between two and three nucleons, respectively. The \(P_{{\mathbf {i}}}\) denote projectors onto the \(N\!N\)S waves, \({\mathbf {i}}= {}^1S_0,{}^3S_1 \). Contact interactions with derivatives as well as higher-body forces are contained in the ellipses in Eq. (4), along with other interactions not shown explicitly here.

The ellipses in Eq. (4) represent a fundamental feature of an EFT, namely that the Lagrangian contains all possible terms which are allowed by the symmetries of the system at hand. For the EFT of nucleons considered here these symmetries are inherited from QCD as the underlying theory: each term in Eq. (4) is required to be invariant under Galilean boosts (plus systematic relativistic corrections), rotations, isospin, and other discrete symmetries. It is of course not arbitrary that Eq. (4) explicitly shows some terms but not others. In order to be predictive, each EFT comes with an organizational principle called “power counting,” which attributes the various terms to increasingly higher orders. A starting point for this organization is typically a naïve dimensional analysis (NDA): fields and derivatives acting on them are assigned their canonical dimensions, defining the exponent of a typical low-momentum scale Q. In order for each term to have overall dimension four, appropriate powers of the EFT breakdown scale \(M_{\text {hi}}\) are included in the prefactor. For the standard Pionless EFT expansion, \(M_{\text {hi}}\sim R^{{-}1}\sim M_\pi \), and this is kept for the construction of the unitarity expansion. However, while standard Pionless EFT assumes \(Q \sim 1/a_{s,t}\), Ref. [24] suggested to count these scales separately as \(\aleph \sim 1/a_{s,t}\) while assuming that

$$\begin{aligned} Q \sim Q_A = \sqrt{2M_NB_A/A} . \end{aligned}$$
(5)

This is a momentum scale associated with the binding energy per nucleon in an A-nucleon system, which for \(A=2\) coincides with the canonical definition of the binding momentum. With this assumption one obtains \(\aleph< Q < 1/R\) such that it is possible to set up a combined expansion in two parameters \(\aleph /Q\) and QR. Coulomb effects are perturbative for momenta of order \(Q_A\) as well and are naturally captured by the expansion if one takes into account that the Coulomb momentum scale \(k_C = \alpha M_N\), with the fine-structure constant \(\alpha \approx 1/137\), is naturally included in the \(\aleph \) scale [26].

For a calculation of few-body states it is convenient to switch from the Lagrangian formulation to standard quantum mechanics expressed in terms of potentials. In the two-body sector, it is possible to write

$$\begin{aligned} V_{2,{\mathbf {i}}}^{(0)} = C_{0,{\mathbf {i}}}^{(0)}|g\rangle \langle g| , \end{aligned}$$
(6)

where \(C_{0,{\mathbf {i}}}^{(0)}\) is the leading-order (LO) piece of the non-derivative contact LEC \(C_{0,{\mathbf {i}}}\) in Eq. (4), which has an expansion of the form

$$\begin{aligned} C_{0,{\mathbf {i}}} = C_{0,{\mathbf {i}}}^{(0)} + C_{0,{\mathbf {i}}}^{(1)} + \cdots . \end{aligned}$$
(7)

Apart from this, \(V_{2,{\mathbf {i}}}^{(0)}\) is defined in terms of a separable Gaussian regulator function, given by

$$\begin{aligned} \langle p|g\rangle = g(p^2) = \exp ({-}p^2/\varLambda ^2) \end{aligned}$$
(8)

in momentum space. This makes the zero-range theory well defined by regularizing the otherwise divergent interaction via the introduction of a cutoff scale \(\varLambda \). Both the value of \(\varLambda \) and the particular form of the regulator function are arbitrary and renormalization, discussed below, ensures that observables are independent of these choices. The separable form is however particularly convenient for the formalism explained in the following. It makes it possible to algebraically solve the Lippmann-Schwinger equations for the LO T matrices \(t_{\mathbf {i}}^{(0)}\),

$$\begin{aligned} t_{\mathbf {i}}^{(0)} = V_{2,{\mathbf {i}}}^{(0)} + V_{2,{\mathbf {i}}}^{(0)} G_0 t_{\mathbf {i}}^{(0)} , \end{aligned}$$
(9)

where

$$\begin{aligned} G_0(z) = \frac{1}{z - H_0} \end{aligned}$$
(10)

defines, for an arbitrary energy z, the two-body Green’s function in terms of the free (purely kinetic) Hamiltonian \(H_0\). The result is:

$$\begin{aligned}&t_{\mathbf {i}}^{(0)}(z;{\mathbf {k}},{\mathbf {k}}') = \langle {\mathbf {k}}|t_{\mathbf {i}}^{(0)}|{\mathbf {k}}'\rangle = g(k^2) \tau _{\mathbf {i}}(z) g(k'^2) , \end{aligned}$$
(11)
$$\begin{aligned}&\tau _{\mathbf {i}}(z) = \left[ 1/{C_{0,{\mathbf {i}}}^{(0)}} - \langle g|G_0|g\rangle \right] ^{-1} . \end{aligned}$$
(12)

The regulator ensures that \(\langle g|G_0|g\rangle \) is finite. At the on-shell point, \(E = k^2/M_N\) and \({\mathbf {k}}={\mathbf {k}}'\), this solution can be matched directly to the ERE, yielding

$$\begin{aligned} C_{0,{\mathbf {i}}}^{(0)} \rightarrow C_{0,{\mathbf {i}}}^{(0)}(\varLambda ) = \frac{1}{2\pi ^2M_N}\left( \frac{1}{a_{\mathbf {i}}} -\theta _0\varLambda \right) ^{{-}1} . \end{aligned}$$
(13)

With this running coupling appearing in Eq. (6), the two-body sector of the theory is renormalized. The number \(\theta _0\) in general depends on the choice of regulator; for the Gaussian form used here one finds \(\theta _0 = 1/\sqrt{2\pi }\). In the unitarity limit, \(1/a_{\mathbf {i}}= 0\), such that the leading-order two-body interaction is parameter free:

$$\begin{aligned} C_{0,{\mathbf {i}}}^{(0)}(\varLambda ) = \frac{{-}1}{2\pi ^2M_N}\frac{1}{\theta _0\varLambda } . \end{aligned}$$
(14)

Perturbative higher orders are defined by formally expanding the full T matrices \(t_{\mathbf {i}}\) as

$$\begin{aligned} t_{\mathbf {i}}= t_{\mathbf {i}}^{(0)} + t_{\mathbf {i}}^{(1)} + t_{\mathbf {i}}^{(2)} + \cdots , \end{aligned}$$
(15)

where \(t_{\mathbf {i}}^{(0)}\) is defined by Eq. (12). The corrections \(t_{\mathbf {i}}^{(n)}\) for \(n>0\) can conveniently be obtained by solving similar integral equations [25, 27]. For the unitarity expansion, corrections from the finite scattering length enter at NLO via\(C_{0,{\mathbf {i}}}^{(1)}\), yielding a separable potential \(V_{2,{\mathbf {i}}}^{(1)}\) with the same form as Eq. (6). For \(t_{\mathbf {i}}^{(1)}\), this gives rise to

$$\begin{aligned} t_{\mathbf {i}}^{(1)} = V_{2,{\mathbf {i}}}^{(1)} + V_{2,{\mathbf {i}}}^{(1)} G_0 t_{\mathbf {i}}^{(0)} + V_{2,{\mathbf {i}}}^{(0)} G_0 t_{\mathbf {i}}^{(1)} , \end{aligned}$$
(16)

which, just like the LO equation, can be solved algebraically (see Ref. [25] for explicit details). From this procedure one obtains

$$\begin{aligned} C_{0,{\mathbf {i}}}^{(1)}(\varLambda ) = \frac{{-}2\pi ^2M_N}{a_{\mathbf {i}}} \left[ C_{0,{\mathbf {i}}}^{(0)}(\varLambda )\right] ^2 . \end{aligned}$$
(17)

Range corrections enter at NLO together with \(C_{0,{\mathbf {i}}}^{(1)}\). They are generated by contact interactions involving quadratic derivatives acting on the nucleon fields, included in the ellipses in Eq. (4). The corresponding potential can be written in momentum space as \(C^{(1)}_{2,{\mathbf {i}}} g(k^2) \left( k^2+k'^2\right) g(k'^2)\). By virtue of this still being a separable interaction, the corresponding version of Eq. (16) with

$$\begin{aligned}&\langle \mathbf {k}|V^{(1)}_{2,{\mathbf {i}}}|\mathbf {k}'\rangle \nonumber \\&\quad = C_{0,{\mathbf {i}}}^{(1)} g(k^2) g(k'^2) + C^{(1)}_{2,{\mathbf {i}}} g(k^2) \left( k^2+k'^2\right) g(k'^2) \end{aligned}$$
(18)

can still be solved algebraically. Matching the result to the ERE (1) up to the quadratic term gives

$$\begin{aligned} C_{2,{\mathbf {i}}}^{(1)}(\varLambda )= & {} {\pi ^2M_N} \left( \frac{r}{2} - \frac{1}{\theta _2\varLambda }\right) \left[ C_{0,{\mathbf {i}}}^{(0)}(\varLambda )\right] ^2 , \end{aligned}$$
(19a)
$$\begin{aligned} C_{0,{\mathbf {i}}}^{(1)}(\varLambda )= & {} {4\pi ^2M_N} {\theta _2\varLambda ^3} C_{0,{\mathbf {i}}}^{(0)}(\varLambda ) C_{2,{\mathbf {i}}}^{(1)}(\varLambda ) , \end{aligned}$$
(19b)

with \(\theta _2 = \theta _0/4\) for the Gaussian regulator used here. Going to higher orders is straightforward, proceeding in the same way via integral equations that can be solved algebraically, recursively using the solutions of previous orders [25, 27]. At second order, the T-matrix correction is obtained from

$$\begin{aligned} t_{\mathbf {i}}^{(2)} = V_{2,{\mathbf {i}}}^{(2)} + V_{2,{\mathbf {i}}}^{(2)} G_0 t_{\mathbf {i}}^{(0)} + V_{2,{\mathbf {i}}}^{(1)} G_0 t_{\mathbf {i}}^{(1)} + V_{2,{\mathbf {i}}}^{(0)} G_0 t_{\mathbf {i}}^{(2)} . \end{aligned}$$
(20)

The treatment of Coulomb contributions (which are neglected in this work), in particular the matching between perturbative and nonperturbative Coulomb regimes within Pionless EFT, is discussed in Refs. [25, 26].

Leading order for \(A>2\) is not complete with only the \(C_{0,{\mathbf {i}}}^{(0)}\) interactions. It is a distinct feature of Pionless EFT, intimately related to the Efimov effect [14,15,16], that a three-nucleon interaction enters at LO. Naïvely it would be expected to contribute only much later in the power counting because the larger number of fields, according to NDA, implies more inverse powers of \(M_{\text {hi}}\) in the prefactor. Analogously to the two-body interactions, the potential induced by the term involving \(D_0\) in Eq. (4) can be written in a separable form,

$$\begin{aligned} V_3^{(0)} = D_0^{(0)} \, |{}^3\mathrm {H} \rangle |\xi \rangle \langle \xi |\langle {}^3\mathrm {H} | \end{aligned}$$
(21)

at LO, where \(|{}^3\mathrm {H} \rangle \) projects onto a three-nucleon state and the regulator \(|\xi \rangle \) is defined, for Jacobi momenta \({\mathbf {u}}_1 = \frac{1}{2}({\mathbf {k}}_1-{\mathbf {k}}_2)\) and \({\mathbf {u}}_2 = \frac{2}{3}[{\mathbf {k}}_3-\frac{1}{2}({\mathbf {k}}_1+{\mathbf {k}}_2)]\), as

$$\begin{aligned} \langle {\mathbf {u}}_1 {\mathbf {u}}_2|\xi \rangle = g\big (u_1^2+\frac{3}{4}u_2^2\big ). \end{aligned}$$
(22)

The \({\mathbf {k}}_i\) label the individual nucleon momenta. An NLO correction \(V_3^{(1)}\) has the same form as Eq. (21), but involves the LEC \(D_0^{(1)}\). Both \(D_0^{(0)}\) and \(D_0^{(1)}\) are determined by the triton binding energy and then enter in other calculations of \(A\ge 3\) observables. These are described in the following.

2 Faddeev and Faddeev-Yakubovsky equations

This section gives an overview of the three- and four-body formalism, implementing a unified framework to solve, respectively, the Faddeev and Faddeev-Yakubovsky equations used to obtain the results presented in the following sections. The main aspects are explained in broad strokes, referring the reader to the references given for more background. Developments needed to calculate charge radii along with perturbative corrections are, however, elaborated on further, with key results explained in the main text and additional details provided in Appendices A and B.

The basis for a description of the three-nucleon system are Jacobi momenta

$$\begin{aligned} {\mathbf {u}}_1= & {} \frac{1}{2}({\mathbf {k}}_1-{\mathbf {k}}_2) , \end{aligned}$$
(23a)
$$\begin{aligned} {\mathbf {u}}_2= & {} \frac{2}{3}\left[ {\mathbf {k}}_3-\frac{1}{2}({\mathbf {k}}_1+{\mathbf {k}}_2)\right] , \end{aligned}$$
(23b)

where the \({\mathbf {k}}_i\) label the individual nucleon momenta, conjugate to position vectors \({\mathbf {x}}_i\). Projecting these momenta onto partial waves yields states \(|u_1u_2;s\rangle \), where

$$\begin{aligned} |s\rangle = | \left( {l_2}{\left( {\left( {l_1}{s_1}\right) \!{j_1}}{\tfrac{1}{2}}\right) \!{s_2}}\right) \!{J}; \left( {t_1}{\tfrac{1}{2}}\right) \!{T} \rangle \end{aligned}$$
(24)

collects angular momentum, spin, and isospin quantum numbers. They are coupled such that \(\left( {l_1}{s_1}\right) \!{j_1}\) and \(t_1\) describe the two-nucleon subsystem, whereas \(l_2\) denotes the orbital angular momentum associated with the Jacobi momentum \(u_2\) and \(s_2\) is an intermediate quantum number. For the trinucleon bound states, the total spin and isospin are \(J=T=1/2\). These states are determined by solving the Faddeev equation [28]

$$\begin{aligned} |\psi ^{(0)}\rangle= & {} G_0\,t^{(0)}\,P |\psi ^{(0)}\rangle \nonumber \\&+ \frac{1}{3}(G_0 + G_0t^{(0)}\,G_0)V_3^{(0)}(1+P)|\psi ^{(0)}\rangle , \end{aligned}$$
(25)

where \(|\psi ^{(0)}\rangle = |\psi _{(12)3}^{(0)}\rangle \) is one of three equivalent two-body Faddeev components. As already done in the discussion of the two-body sector, explicit superscripts “(0)” are used to denote leading-order quantities. Alternatively, one can incorporate the three-body interaction \(V_3^{(0)}\) by writing [29]

$$\begin{aligned} |{{\tilde{\psi }}}^{(0)}\rangle= & {} G_0\,t^{(0)}\,P |{{\tilde{\psi }}}^{(0)}\rangle + G_0\,t^{(0)}\,|\psi _3^{(0)}\rangle , \end{aligned}$$
(26a)
$$\begin{aligned} |\psi _3^{(0)}\rangle= & {} G_0\,t_3^{(0)}\,(1+P) |{{\tilde{\psi }}}^{(0)}\rangle , \end{aligned}$$
(26b)

where \(|\psi _3^{(0)}\rangle \) is an auxiliary amplitude, and

$$\begin{aligned} t_3^{(0)} = V_3^{(0)} + V_3^{(0)}\,G_0\,t_3^{(0)} . \end{aligned}$$
(27)

The tilde is used to distinguish \(|{{\tilde{\psi }}}^{(0)}\rangle \) from \(|\psi ^{(0)}\rangle \). In either form of the Faddeev equations, \(G_0\) denotes the free three-body Green’s function and \(P = P_{12}P_{23} + P_{13}P_{23}\) generates the non-explicit components through permutations. \(t^{(0)}\) collectively denotes the two-body T-matrices \(t_{\mathbf {i}}^{(0)}\). Note that \(|\psi _3^{(0)}\rangle \) can be eliminated by inserting Eq. (26b) into Eq. (26a), yielding an equation of the form

$$\begin{aligned} |\psi ^{(0)}\rangle = K^{(0)}|\psi ^{(0)}\rangle \end{aligned}$$
(28)

with

$$\begin{aligned} K^{(0)} = G_0\,t^{(0)}\,P + G_0\,t^{(0)}\,G_0\,t_3^{(0)}\,(1+P) , \end{aligned}$$
(29)

and alternatively a similar kernel can be obtained from Eq. (25) in terms of \(V_3^{(0)}\). Either form of the Faddeev equations is solved by representing it within the space of states \(|u_1u_2;s\rangle \), discretizing the momenta \(u_{1,2}\) on a quadrature mesh. The binding energy is determined by varying the energy E, entering as an argument to both \(G_0\) and \(t^{(0)}\) until the kernel K has a unit eigenvalue. At that energy, \(E = {-}B_0\), \(|\psi ^{(0)}\rangle \) can then be determined as the corresponding eigenstate and finally one constructs the full wavefunction as

$$\begin{aligned} |\varPsi ^{(0)}\rangle = (1 + P) |\psi ^{(0)}\rangle + |\psi _3^{(0)}\rangle . \end{aligned}$$
(30)

If one starts from Eq. (25), there is no explicit three-body Faddeev component and one simply has \(|\varPsi ^{(0)}\rangle = (1 + P) |\psi ^{(0)}\rangle \).

Fundamentally, the permutation operator P leads to a coupling of different partial waves (see Ref. [30] for an excellent pedagogical discussion of both this and the Faddeev equations in general), and for the construction of the full wavefunction (30) it is important to include higher partial waves: the proper antisymmetry of \(|\varPsi ^{(0)}\rangle \) is only recovered as more and more states are included, which means that in principle all observables calculated from \(|\varPsi ^{(0)}\rangle \) have to be checked for convergence.

The full wavefunctions \(|\varPsi ^{(0)}\rangle \) are used to calculate both perturbative shifts for the binding energy,

$$\begin{aligned} B_1 = \langle \varPsi ^{(0)}|V^{(1)}|\varPsi ^{(0)}\rangle , \end{aligned}$$
(31)

as well as the radius at LO (see Sect. 4). Note that \(|\varPsi ^{(0)}\rangle \) is assumed here to be properly normalized,

$$\begin{aligned} \langle \varPsi ^{(0)}|\varPsi ^{(0)}\rangle = 1 . \end{aligned}$$
(32)

In calculating the matrix elements in Eqs. (31) and (32) it is advantageous to exploit the antisymmetry of \(|\varPsi ^{(0)}\rangle \) as much as possible because that will speed up convergence of results with respect to the number of partial-wave channels. For example, the two-body part of the full potential \(V^{(1)}\) can be expressed through permutations in terms of only the potential between the pair of nucleons 1 and 2, and it holds that \((1+P)(1+P) = 3(1+P)\).

Note furthermore that Eq. (26), and similarly Eq. (25), can be significantly simplified by exploiting the fact that the two- and three-body interactions are separable and act only within S waves. As a result, it suffices to work with merely two coupled equations for the triton [29], and using the procedure described in the previous paragraph it is furthermore possible to eliminate all intermediate higher partial-wave components in the calculation of \(B_1\). These simplifications are used in the practical implementation to fit the three-body LECs \(D_0^{(0)}\) and \(D_0^{(1)}\).

Describing the four-nucleon system requires an additional Jacobi momentum

$$\begin{aligned} {\mathbf {u}}_3 = \frac{3}{4}\left[ {\mathbf {k}}_4-\frac{1}{3}({\mathbf {k}}_1+{\mathbf {k}}_2+{\mathbf {k}}_3)\right] , \end{aligned}$$
(33)

as well as an alternative set of momenta \(({\mathbf {v}}_1,{\mathbf {v}}_2,{\mathbf {v}}_3)\), describing a 2+2 cluster setup, i.e., \({\mathbf {v}}_1={\mathbf {u}}_1\), \({\mathbf {v}}_3\) denotes the relative momentum in the (34) system, and \({\mathbf {v}}_2\) is defined as the relative momentum between the (12) and (34) subsystems.

Including the remaining quantum numbers, this leads to channel states

$$\begin{aligned} |a\rangle= & {} |\! \left( {l_2}{\left( {\left( {l_1}{s_1}\right) \!{j_1}}{\tfrac{1}{2}}\right) \!{s_2}}\right) \!{j_2}, \left( {l_3}{\tfrac{1}{2}}\right) \!{j_3}, \left( {j_2}{j_3}\right) \!{J}; \left( {\left( {t_1}{\tfrac{1}{2}}\right) \!{t_2}}{\tfrac{1}{2}}\right) \!{T} \rangle ,\nonumber \\ \end{aligned}$$
(34a)
$$\begin{aligned} |b\rangle= & {} |\! \left( {\lambda _2}{\left( {\lambda _1}{\sigma _1}\right) \!{\iota _1}}\right) \!{\iota _2}, \left( {\lambda _3}{\sigma _3}\right) \!{\iota _3}, \left( {\iota _2}{\iota _3}\right) \!{J}; \left( {\tau _2}{\tau _3}\right) \!{T} \rangle . \nonumber \\ \end{aligned}$$
(34b)

The \(|a\rangle \) are a straightforward extension of three-nucleon states (24), including the angular momentum \(l_3\) associated with \({\mathbf {u}}_3\) as well as spin and isospin \(\tfrac{1}{2}\) for the fourth nucleon. For the b states, \((\lambda _1,\sigma _1,\tau _1)\) and \((\lambda _3,\sigma _3,\tau _3)\) are two-body quantum numbers for the (12) and (34) subsystems, respectively, whereas \(\lambda _{1,2,3}\) are the angular momenta associated with \(v_{1,2,3}\). For \({}^4\mathrm {He}\) the total spin and isospin are \(J=T=0\).

Following Refs. [29, 31], the Faddeev-Yakubovsky equations can be written as

$$\begin{aligned} |\psi _A^{(0)}\rangle= & {} G_0 t^{(0)} P \big [(1 - P_{34})|\psi _A^{(0)}\rangle + |\psi _B^{(0)}\rangle \big ] \nonumber \\&+ \frac{1}{3}(1 + G_0 t^{(0)}) G_0 V_3^{(0)} |\varPsi ^{(0)}\rangle \end{aligned}$$
(35a)
$$\begin{aligned} |\psi _B^{(0)}\rangle= & {} G_0 t^{(0)} {\tilde{P}} \big [(1 - P_{34})|\psi _A^{(0)}\rangle + |\psi _B^{(0)}\rangle \big ] , \end{aligned}$$
(35b)

corresponding to the decomposition

$$\begin{aligned} |\varPsi ^{(0)}\rangle= & {} (1 + P)\Big [(1 - P_{34} - P_{34}P)|\psi _A^{(0)}\rangle \nonumber \\&+ (1 + {\tilde{P}})|\psi _B^{(0)}\rangle \Big ] \end{aligned}$$
(36)

of the full four-body wavefunction. The two distinct Faddeev-Yakubovsky components \(|\psi _A\rangle \) and \(|\psi _B\rangle \) correspond to, respectively, 3+1 and 2+2 cluster configurations of the four-body system, with the former naturally expressed in terms of the Jacobi momenta \(u_i\) and states \(|a\rangle \), and the latter in terms of \(v_i\) and \(|b\rangle \). Note that the same notation is used here as for the three-body case, and \(G_0\) in Eqs. (35a, b) now represents the free four-body Green’s function. In addition to the operator P already encountered in the three-body system, Eqs. (35a, b) include the further permutations \(P_{34}\) and \({\tilde{P}} = P_{13}P_{24}\) to ensure proper antisymmetry. In fact, the overall symmetry is determined by the sign in front of \(P_{34}\): to study a bosonic system, one would use \((1+P_{34})\) acting on \(|\psi _A^{(0)}\rangle \) to construct fully symmetric states.

The structure of Eqs. (35a, b) can be made clearer by rewriting them in a generic matrix form,

$$\begin{aligned} \left( \mathbb {1}- {\hat{K}}^{(0)}\right) |\varvec{\psi }^{(0)}\rangle = 0 , \end{aligned}$$
(37)

with \(|\varvec{\psi }^{(0)}\rangle = (|\psi _A^{(0)}\rangle , |\psi _B^{(0)}\rangle )^T\) and the kernel

$$\begin{aligned} {\hat{K}}^{(0)} = G_0 t^{(0)} {\hat{P}} + \frac{1}{3}(G_0 + G_0 t^{(0)} G_0) V_3^{(0)} {\hat{P}}_3 . \end{aligned}$$
(38)

In Eq. (38), \(G_0\) and \(t^{(0)}\) are understood to be diagonal matrices, and the permutation operators are collected in

$$\begin{aligned} {\hat{P}}= & {} \mathrm {diag}(P, {\tilde{P}}) \otimes \begin{pmatrix} (1-P_{34}) &{} 1 \\ (1-P_{34}) &{} 1 \end{pmatrix} , \end{aligned}$$
(39a)
$$\begin{aligned} {\hat{P}}_{3}= & {} \begin{pmatrix} (1+P)(1-P_{34}-P_{34}P) &{} (1+P)(1+{\tilde{P}}) \\ 0 &{} 0 \end{pmatrix} ,\nonumber \\ \end{aligned}$$
(39b)

where “\(\otimes \)” represents matrix multiplication. From this form the structural analogy to the three-body Faddeev case becomes obvious.

Just like for the Faddeev equations, Eqs. (35a, b) are solved by projecting onto states \(|u_1u_2u_3;a\rangle \), \(|v_1v_2v_3;b\rangle \) [31], discretizing all momenta on a grid, and looking for a unit eigenvalue of the resulting kernel matrix as a function of the energy. However, the set of coupled equations does not naturally truncate even if all interactions are pure S wave. This means that already for a determination of the binding energies it is necessary to truncate the sums in Eq. (34) (by choosing all total angular momenta \(j_i\) and \(\iota _i\) less than some \(j_{\text {max}}\)) and study the numerical convergence of results as \(j_{\text {max}}\) is increased.

3 Binding energies of light nuclei

By construction, the unitarity expansion renders the deuteron a zero-energy bound state at leading order. Since the expansion is set up in powers of the inverse scattering lengths, it corresponds to the zero-range binding momentum \(\kappa _t = 1/a_t\) in the \({}^3S_1\) channel. As demonstrated explicitly in Ref. [25] by using the perturbative formalism discussed in Sect. 1, this implies that the deuteron remains at zero energy at NLO and only moves to \(1/(M_Na_t^2)\) in an N\(^2\)LO calculation. This is so for both the pure expansion in \(1/a_{t}\), neglecting range corrections, but interestingly also for the full unitarity expansion that includes, via\(C_{2,{\mathbf {i}}}^{(1)}\), range corrections starting at NLO. This is so because the unitarity LO shifts all corrections that mix ERE parameters to a higher order compared to where they would be with a finite scattering length at LO. Overall, the dominant source of uncertainty for the deuteron energy comes from the \(1/(Q_2 a_t)\) expansion, which still amounts to a 50% effect at \(\text {N}^{2}\text {LO} \). Conservatively taking the experimental binding energy as reference value for the uncertainty estimate yields \(B_d^{\text {N}^{2}\text {LO}} = 1.41 \pm 1.12~\mathrm {MeV}\).

The triton, being the “anchor point” of the expansion that determines the value of three-body parameter \(D_0^{(0)}\), stays fixed at the physical binding energy at each order. With the finite physical scattering lengths entering through \(C_{0,{\mathbf {i}}}^{(1)}\) at NLO, the three-body LEC \(D_0^{(1)}\) compensates the shift in the triton energy to keep it in place. This leaves the \(^{3}\hbox {He}\) binding energy as a nontrivial prediction. While at LO by construction the trinucleon bound states are degenerate,Footnote 1 finite-scattering-length corrections together with Coulomb effects (specifically, one-photon exchange) produce a triton-helion splitting \((B_{T}-B_{H})^{\text {NLO}} \simeq (0.92 \pm 0.18)~\mathrm {MeV}\) at NLO [24], in good agreement with the experimental value \((B_T-B_H)^{\mathrm{exp}}\simeq 0.764~\mathrm {MeV}\). Details regarding the perturbative treatment of Coulomb effects are discussed in Refs. [25, 26].

The \({}^4\mathrm {He}\) nucleus provides a more serious test of the unitarity expansion. Since \(Q_A R \sim 0.8\) for \({}^4\mathrm {He}\), it is the standard Pionless EFT part of the unitarity scheme which naïvely one might doubt to work, while the pure unitarity expansion, \(\aleph /Q_A\), should indeed work better with increasing \(Q_A\). Figure 1 shows the \({}^4\mathrm {He}\) binding energy as a function of the momentum cutoff \(\varLambda \). The observed convergent behavior as \(\varLambda \) increases indicates that the EFT calculation is properly renormalized, as established originally in Refs. [29, 32, 33]. Results for the standard pionless LO, given by the (blue) diamonds in Fig. 1, are consistent with this earlier work.

While any \(\varLambda \) above the breakdown scale (of order \(M_\pi \)) is a valid cutoff choice in principle, polynomials in \(1/\varLambda \) are fitted to the points in Fig. 1 to quantitatively assess the convergence and conveniently extrapolate \(\varLambda \rightarrow \infty \). This procedure gives \(B_\alpha = 39(12)~\mathrm {MeV}\) in the unitarity limit, with the uncertainty estimated as \({\mathcal {O}}(r_{s,t}/a_{s,t}) \simeq 30\%\) based on the expectation that range effects \({}^3S_1\) are the dominant correction. Including the finite-scattering lengths as NLO corrections gives the (green) circles in Fig. 1, very close to the standard pionless LO, indicating that the \(1/(Q_4a_{s,t})\) expansion appears to converge remarkably well up to this order, and indeed the extrapolated result 30(9) \(\mathrm {MeV}\) comes out very near the standard pionless LO value of \(31(9)~\mathrm {MeV}\).Footnote 2

Fig. 1
figure 1

\({}^4\mathrm {He}\) binding energy as function of the Gaussian cutoff parameter \(\varLambda \). The (blue) diamonds and (green) squares show, respectively, the results for standard Pionless EFT and the unitarity expansion at LO. Inclusion of first-order corrections in \(1/a_{s,t}\) (i.e., an incomplete NLO that neglects range and Coulomb effects) gives the (green) circles. The closeness to the standard leading order demonstrates how well this part of the expansion convergences. Large symbols on the right edge indicate results for an extrapolation \(\varLambda \rightarrow \infty \) (see text)

Full size image

It should be stressed that an NLO including only finite-a corrections is incomplete: in the full unitarity expansion, range corrections and Coulomb effects enter at the same time. While the latter are expected to be small given that \({}^4\mathrm {He}\) is rather deeply bound, it turns out that the inclusion of range effects actually has a profound consequence: in Ref. [35] it is shown that a four-body interaction is required to renormalize the universal four-boson system once range corrections are included at NLO, and universality implies that this conclusion carries over directly to \({}^4\mathrm {He}\) in Pionless EFT. The implication is that a four-nucleon input datum is required to fix the unknown four-body parameter at NLO, and this is most naturally taken to be \({}^4\mathrm {He}\) ground-state energy. Other properties, such as the ground-state radius (discussed below) or the position of the \(0^+\) excited state, will remain predictions at NLO  – unless it turns out that additional many-body forces are required for a renormalized NLO calculation of these observables. While it may seem unlikely and is certainly not to be expected based on NDA, such a possibility cannot a priori be excluded: each calculation needs to be carefully checked if it is properly renormalized.

The rapid convergence of the pure unitarity expansion persists off the physical point. Figure 2 shows the correlation between 3N and 4N binding energies, known as the Tjon line [36]. Its existence is explained by the three-body parameter largely governing the physics of the system [33]. It is seen that the result starting from unitarity is shifted very close to having the exact scattering lengths at LO over a significant range or triton energies. This observation provides further evidence that the unitarity expansion converges well and that the results found at the physical point are not merely accidental.

The unitarity expansion for ground-state energies up to \(A=4\) is summarized in Table 1. Observables fixed as input data at a given order are shown as underlined text. This in particular includes the \({}^3\mathrm {He}\) binding energy at \(\text {N}^{2}\text {LO} \)  since according to Refs. [25, 26] an isospin-breaking three-body force is required once perturbation theory mixes Coulomb effects and range corrections.Footnote 3 Taking into account the findings of Ref. [35], the \({}^4\mathrm {He}\) binding energy should be underlined for a complete NLO as well.

Fig. 2
figure 2

Tjon line: correlation between the \({}^4\mathrm {He}\) and \({}^3\mathrm {H}\) binding energies. (Blue) dotted curve: standard pionless LO result; (green) dashed upper curve: unitarity limit at LO. Additional points nearly on top of the blue curve: inverse scattering lengths added in first-order perturbation theory. Star: experimental point

Full size image
Table 1 Unitarity expansion convergence pattern. Underlined values indicate energies which are used as input values to determine three-body LECs. An asterisk superscript indicates an incomplete NLO result which only includes the finite-scattering length but no contributions from effective ranges or electromagnetic interactions
Full size table

4 Charge radii and form factors

It has so far been established that the unitarity expansion describes well the ground-state energies of light nuclei. While certainly impressive given how simple the LO of the expansion is, it is still merely a first step towards showing that the scheme comprehensively captures the properties of light nuclei.

Further insight can be gained by considering charge radii of the \(A=3,4\) systems as well. Within the momentum-space framework employed in this work, this is achieved by calculating charge form factors

$$\begin{aligned} F_C(q^2) = \langle \varPsi |{\hat{\rho }}({\mathbf {q}})|\varPsi \rangle \end{aligned}$$
(40)

for the \({}^3\mathrm {H}\) and \({}^4\mathrm {He}\) ground states, from which one obtains point charge radii as

$$\begin{aligned} \langle r_0^2\rangle = {-}\dfrac{1}{6}\dfrac{\mathrm {d}}{\mathrm {d}(q^2)} F_C(q^2)\Big |_{q^2=0} \ \ \text {,}\ \ \langle r_0\rangle \equiv \sqrt{\langle r_0^2\rangle } . \end{aligned}$$
(41)

In Eq. (40), the total charge operator \({\hat{\rho }} \equiv J_0\), i.e., the zero component of the electric current \(J_\mu \), is given by the sum of the individual nucleon contributions,

$$\begin{aligned} {\hat{\rho }} = \sum _{i=1}^A {\hat{\rho }}_i . \end{aligned}$$
(42)

The (anti-)symmetry of the wavefunction makes it possible to replace this sum by \(A{\hat{\rho }}_i\) for any fixed i. A particularly convenient choice for three nucleons is \(i=3\) because it holds that

$$\begin{aligned} {\mathbf {x}}_3 = {\mathbf {R}}^{(3)} + \frac{2}{3}{\mathbf {r}}_2 , \end{aligned}$$
(43)

where \({\mathbf {r}}_2\) is the relative distance conjugate to \({\mathbf {u}}_2\) and \({\mathbf {R}}^{(3)}\) is the overall center-of-mass coordinate. With this choice, the momentum-space expression for the current operator involves a momentum transfer only onto the Jacobi momentum \({\mathbf {u}}_2\). Likewise, for four nucleons a good choice is \(i=4\) because with

$$\begin{aligned} {\mathbf {x}}_4 = {\mathbf {R}}^{(4)} + \frac{3}{4}{\mathbf {r}}_3 \end{aligned}$$
(44)

one obtains a momentum transfer only onto \({\mathbf {u}}_3\).

To use \({\hat{\rho }}\) within the Faddeev-Yakubovsky framework, it is necessary to represent it within the appropriate partial-wave basis. Since \({\hat{\rho }}\) does not depend on spin, two-body matrix elements of \({\hat{\rho }}\) are given by

$$\begin{aligned}&\langle u;\left( {l}{s}\right) \!{j} m|{\hat{\rho }}({\mathbf {q}})|u';\left( {l'}{s'}\right) \!{j}'m'\rangle \nonumber \\&\quad = \delta _{jj'}\delta _{mm'} \langle u;l|\;\!\!|{\hat{\rho }}({\mathbf {q}})|\;\!\!|u';l\rangle , \end{aligned}$$
(45)

where the reduced matrix element on the right-hand side is given by

$$\begin{aligned}&\langle u;l m|{\hat{\rho }}({\mathbf {q}})|u';l'm'\rangle \nonumber \\&\quad = \delta _{ll'}\delta _{mm'} \frac{1}{2}\sum _{k} \sqrt{\left( {\begin{array}{c}2l\\ 2k\end{array}}\right) } C_{k0,(l-k)0}^{l0} \nonumber \\&\quad \quad \times \int _{-1}^1 \mathrm {d}x\,P_k(x) \frac{\delta \big (u'-\iota (u,q,x)\big )}{u'^2} \frac{u^{l-k}\big ({-}\tfrac{1}{2}q\big )^{k}}{\iota (u,q,x)^{l}} , \end{aligned}$$
(46)

without the \(\delta _{mm'}\), where

$$\begin{aligned} \iota (p,q,x) = \sqrt{ p^2-pqx+q^2/4} . \end{aligned}$$
(47)

A detailed derivation of Eq. (46) is provided in Appendix A.

Embedding \({\hat{\rho }}\) into the three- and four-body bases merely leads to additional Dirac and Kronecker deltas, as well as to kinematic prefactors multiplying the momentum transfer q which can be read off from Eqs. (43) and (44):

$$\begin{aligned}&\langle u_1u_2;s|{\hat{\rho }}_3({\mathbf {q}})|u_1'u_2';s'\rangle \nonumber \\&\quad = \delta _{l_1l_1'} \delta _{l_1l_1'} \delta _{s_1s_1'} \delta _{t_1t_1'} \delta _{s_2s_2'} \delta _{JJ'} \delta _{TT'}\nonumber \\&\qquad \times \frac{\delta (u_1-u_1')}{u_1^2} \langle u_2;l_2|\;\!\!|{\hat{\rho }}(\tfrac{4}{3}{\mathbf {q}})|\;\!\!|u_2';l_2'\rangle , \end{aligned}$$
(48)
$$\begin{aligned}&\langle u_1u_2u_3;a|{\hat{\rho }}_4({\mathbf {q}})|u_1'u_2'u_3';a'\rangle \nonumber \\&\quad = \delta _{l_1l_1'} \delta _{s_1s_1'} \delta _{j_1j_1'} \delta _{t_1t_1'} \delta _{s_2s_2'} \delta _{j_2j_2'} \delta _{t_2t_2'} \delta _{JJ'} \delta _{TT'}\nonumber \\&\qquad \frac{\delta (u_1-u_1')}{u_1^2} \frac{\delta (u_2-u_2')}{u_2^2} \langle u_3;l_3|\;\!\!|{\hat{\rho }}(\tfrac{3}{2}{\mathbf {q}})|\;\!\!|u_3';l_3'\rangle . \end{aligned}$$
(49)

At leading order in the unitarity expansion, the form factor is given by

$$\begin{aligned} F_C^{(0)}(q^2) = \langle \varPsi ^{(0)}|{\hat{\rho }}({\mathbf {q}})|\varPsi ^{(0)}\rangle , \end{aligned}$$
(50)

and analogously Eq. (41), with added superscripts “(0),” yields the LO point charge radii. At NLO, the correction to the form-factor isFootnote 4

$$\begin{aligned} F_C^{(1)}(q^2) = 2\,\langle \varPsi ^{(1)}|{\hat{\rho }}({\mathbf {q}})|\varPsi ^{(0)}\rangle , \end{aligned}$$
(51)

where the factor 2 follows from symmetry. Perturbatively expanding Eq. (41) gives

$$\begin{aligned} \langle r_0\rangle ^{(1)} = \frac{1}{2}\frac{\langle r_0^2\rangle ^{(1)}}{\langle r_0\rangle ^{(0)}} , \end{aligned}$$
(52)

with \(\langle r_0^2\rangle ^{(1)}\) calculated from the slope of \(F_C^{(1)}(q^2)\) at \(q^2=0\).

Evaluating the perturbative radius shifts defined above requires the NLO correction to the wavefunctions that enter in Eq. (51). Is is possible to obtain these for three- and four-body systems from inhomogeneous versions of the Faddeev and Faddeev-Yakubovsky equations, respectively. As is derived in Appendix B, for three particles one has \(|\varPsi _1\rangle = (1+P)|\psi _1\rangle \), with the NLO Faddeev component \(|\psi _1\rangle \) defined as a solution ofFootnote 5

$$\begin{aligned}&\left[ 1 - G_0 t^{(0)} P - (G_0 + G_0t^{(0)}\,G_0)V_3^{(0)} \right] |\psi _1\rangle \nonumber \\&\quad = (G_0 + G_0 t^{(0)} G_0)\left[ V_2^{(1)}(1+P) + V_3^{(1)}+ B_1 \right] |\psi _0\rangle .\nonumber \\ \end{aligned}$$
(53)

In Eq. (53), \(G_0\) and \(t^{(0)}\) are understood to be evaluated at the LO binding energy, \(E = {-}B_0\). Similarly, for four nucleons NLO Faddeev-Yakubovsky equations can be written as

$$\begin{aligned} \left( \mathbb {1}- {\hat{K}}^{(0)}\right) |\varvec{\psi }^{(1)}\rangle = {\hat{K}}^{(1)} |\varvec{\psi }^{(0)}\rangle , \end{aligned}$$
(54)

with the kernel \({\hat{K}}^{(0)}\) as defined in Eq. (38) and

$$\begin{aligned} {\hat{K}}^{(1)}= & {} B_1 (G_0 + G_0t^{(0)}G_0) + G_0 t^{(1)} {\hat{P}} \nonumber \\&+ G_0 t^{(1)} G_0 V_3^{(0)} {\hat{P}}_3 + (G_0 + G_0t^{(0)}G_0) V_3^{(1)} {\hat{P}}_3 .\nonumber \\ \end{aligned}$$
(55)

Note that as explained in Appendix B special care has to be taken when solving Eqs. (53) and (54) to account for the fact that the operators on the left-hand sides are singular at \(E={-}B_0\). From the components \(|\varvec{\psi }^{(1)}\rangle = (|\psi _A^{(1)}\rangle , |\psi _B^{(1)}\rangle )^T\) the full NLO correction \(|\varPsi ^{(1)}\rangle \) is obtained analogously to Eq. (36). For practical calculations Eq. (54) is simplified based on the fact that all interactions are chosen to be separable. This can be achieved with the same factorization as used in Ref. [29] at LO.

Fig. 3
figure 3

Point-charge radii for \({}^3\mathrm {H}\) and \({}^4\mathrm {He}\) as function of the Gaussian cutoff parameter \(\varLambda \). The (purple) diamonds and (cyan) pentagons show, respectively, \({}^3\mathrm {H}\) results for standard Pionless EFT and the unitarity expansion at LO. The (cyan) hexagons are obtained by perturbatively including \(1/*\) correction on top of the unitarity LO. For \({}^4\mathrm {He}\), results are shown in the upper part, with symbols as in Fig. 1. Large symbols on the right edge indicate results for an extrapolation \(\varLambda \rightarrow \infty \) (see text)

Full size image

As for the \({}^4\mathrm {He}\) energy discussed in the previous section, the focus here is on the 1/a part of the unitarity expansion while the inclusion of range corrections is postponed to future work. Results for the ground-state radii of both \({}^3\mathrm {H}\) and \({}^4\mathrm {He}\) are shown in Fig. 3 as a function of the UV cutoff \(\varLambda \). Convergence as \(\varLambda \) increases is evident from the plot, and just like it was done for the binding energies, polynomials in \(1/\varLambda \) are fitted to the data points in order to extrapolate \(\varLambda \rightarrow \infty \). The horizontal lines in Fig. 3 show the experimental values of the point charge radii, which, following Ref. [38], are defined as

$$\begin{aligned} \langle r_0^2\rangle _{{}^3\mathrm {H}} = \langle r^2\rangle _{{}^3\mathrm {H}} - \langle r^2\rangle _p - 2\langle r^2\rangle _n \end{aligned}$$
(56)

for the triton, and

$$\begin{aligned} \langle r_0^2\rangle _{{}^4\mathrm {He}} = \langle r^2\rangle _{{}^4\mathrm {He}} - 2\langle r^2\rangle _p - 2\langle r^2\rangle _n \end{aligned}$$
(57)

for \({}^4\mathrm {He}\). That is, contributions from the root-mean-square radii of the individual nucleons are subtracted from the experimental nuclear charge radii.

Using experimental values from Ref. [39] for the quantities appearing on the right-hand sides of the above definitions one obtains, with error bars negligible compared to those of the present theoretical calculation, \(\langle r_0^2\rangle _{{}^3\mathrm {H}}^{\text {exp}} = 1.59~\mathrm {fm}\) and \(\langle r_0^2\rangle _{{}^4\mathrm {He}}^{\text {exp}} = 1.72~\mathrm {fm}\).

The lower part of Fig. 3 shows results for the triton. In the limit \(\varLambda \rightarrow \infty \), indicated as points on the right border of the plot, the \({}^3\mathrm {He}\) radius comes out as \(\langle r_0\rangle _{{}^3\mathrm {H}}^{(0)} = 1.15(35)~\mathrm {fm}\) for the standard pionless LO, and \(\langle r_0\rangle _{{}^3\mathrm {H}}^{(0)} = 1.04(31)~\mathrm {fm}\) at unitarity. Perturbative corrections shift the unitarity LO result more than half way towards the value obtained for physical scattering lengths at leading order,

$$\begin{aligned} \langle r_0\rangle _{{}^3\mathrm {H}}^{(0)} + \langle r_0\rangle _{{}^3\mathrm {H}}^{(1)} = 1.10(33)~\mathrm {fm}, \end{aligned}$$
(58)

indicating that the unitarity expansion works well for this observable. This is in line with the results of Ref. [40], where good convergence is found for a perturbative expansion of \({}^3\mathrm {H}\) and \({}^3\mathrm {He}\) radii around an SU(4) symmetric leading order (of which the unitarity limit is a special case). The result obtained for standard pionless LO is furthermore in excellent agreement with the calculation of Ref. [38], at unitarity the radius satisfies well the universal relation [38, 40, 41]

$$\begin{aligned} M_NB_{{}^3\mathrm {H}} \langle r_0^2\rangle _{{}^3\mathrm {H}} = (1 + s_0^2)/9 \approx 0.224 . \end{aligned}$$
(59)

As done for binding energies it is assumed here that the \(Q_AR\) part dominates the overall expansion, yielding a \(30\%\) uncertainty both at LO as well as NLO. Indeed, from Ref. [38] it is known that range corrections contribute significantly to the triton radius in Pionless EFT and shift the result close to the experimental value once they are included. This uncertainty assignment places the experimental value outside the error band of the unitarity LO result. Since it is purely based on omitted range corrections, this is not actually a reason for concern and merely indicates that to be yet more conservative one should consider adding the uncertainties from the \(Q_AR\) and \(\aleph /Q_A\) expansions coherently.

Results for \({}^4\mathrm {He}\) radius, shown in the upper part of Fig. 3, look equally good. In fact, consistent with what is found for the binding energy, the result obtained from a standard pionless LO calculation, \(\langle r_0\rangle _{{}^4\mathrm {He}}^{(0)} = 1.69(51)~\mathrm {fm}\), comes out surprisingly close to the experimental value. At unitarity the radius comes out smaller, \(\langle r_0\rangle _{{}^4\mathrm {He}}^{(0)} = 1.49(45)\)\(\mathrm {fm}\), consistent with the observed overbinding in the unitarity limit. Inclusion of perturbative 1/a corrections shifts this value to

$$\begin{aligned} \langle r_0\rangle _{{}^4\mathrm {He}}^{(0)} + \langle r_0\rangle _{{}^4\mathrm {He}}^{(1)} = 1.73(52)~\mathrm {fm}, \end{aligned}$$
(60)

in excellent agreement with the standard pionless LO result and therefore providing yet more evidence for the good convergence properties of the unitarity expansion. In particular, the fact that convergence appears to be somewhat faster for \({}^4\mathrm {He}\) than for the triton is in good agreement with \(\aleph /Q_4 < \aleph /Q_3\) obtained from Eq. (5), and it therefore reinforces confidence in this estimate.

5 Summary and perspectives

Superficially, the unitarity expansion may seem like merely a minor departure from standard Pionless EFT. It is rather well known that Pionless EFT, unlike Chiral EFT, is the ideal EFT to describe few-nucleon systems at low energies since its expansion explicitly embraces implications from the scattering lengths being large, basing its power counting explicitly on this fact. Chiral EFT is limited at low energies by its simultaneous expansion in both momenta and around the chiral limit, with \(M_\pi \ne 0\) parametrizing the distance from it. This combination yields a power counting for \(Q\sim M_\pi \) which does not easily capture the physics of the regime \(Q\ll M_\pi \). Notably, one-pion exchange only contributes to the \(N\!N\) scattering lengths through loop effects.

However, the unitarity expansion does in fact constitute a significant paradigm shift in the EFT-based description of light nuclei: it goes as far as saying that the details of the two-body sector, represented by the experimental values of the scattering lengths, do not actually matter much to describe the gross properties of light nuclei. Instead, it fully embraces universality and uses the three-body sector as anchor point, constructing a leading order with just a single parameter and an exact manifestation of the Efimov effect. In this work it has been shown that the 1/a expansion of the unitarity scheme works well not only for binding energies of up to four nucleon systems, but that similarly good convergence is obtained for the \({}^3\mathrm {H}\) and \({}^4\mathrm {He}\) point charge radii as well. This finding solidifies the picture drawn in the introduction of this work, placing few-nucleon systems in a universal regime perturbatively close to the unitarity limit.

It is an important next step to include range corrections and Coulomb effects, thus considering the full unitarity scheme that pairs the expansions in \(\aleph /Q_A\) and \(Q_A/M_{\text {hi}}\)\(=Q_AR\). So far, this has been investigated only for the \({}^3\mathrm {H}\)-\({}^3\mathrm {He}\) energy splitting, where by construction range corrections cancel at NLO  [24]. An isospin-breaking three-nucleon force is required once range corrections mix with Coulomb contributions at \(\text {N}^{2}\text {LO} \) [25]. Range corrections are known to significantly contribute to the triton point charge radius [38]. For the \({}^4\mathrm {He}\) radius, the closeness to the experimental point already without range corrections found in this work leaves little room (about \(0.01~\mathrm {fm}\) in the central value) for a significant shift at full NLO, suggesting that several of the contributions neglected here might cancel. From Ref. [35] it is known that full NLO requires a four-nucleon force to be included along with range corrections, which is most conveniently fit to reproduce the \({}^4\mathrm {He}\) energy at the experimental point at NLO. It is conceivable that this fixing of the energy will maintain a good reproduction of the radius, just like it is observed for the standard Pionless LO result. Coulomb contributions should also be included in a complete NLO calculation, along with isospin breaking in the \({}^1S_0\)NN scattering lengths. This is expected to be a small effect for a bound state as deep as \({}^4\mathrm {He}\), but it is interesting to note that fitting the four-nucleon force to exactly reproduce the \({}^4\mathrm {He}\) binding energy at NLO will inevitably absorb isospin-breaking contributions as well. Clearly a careful overall consideration of the LEC fitting procedure is called for in light of this to avoid possible overfitting of individual parameters. Bayesian methods stand ready as a powerful tool to address this [42, 43].

Apart from such more technical issues, it is an exciting question how far into the nuclear chart the unitarity expansion can reach and what exactly its place is in the tower of nuclear EFTs. The observation that bosonic systems at unitarity exhibit saturation for large numbers of particles [44] and recent calculations of nuclear matter using interactions guided by unitarity [45] provide reason to be optimistic that universality, and in particular discrete scale invariance [46], is able to inform more than just few-nucleon calculations. On the other hand, Refs. [47, 48] indicate that few-nucleon systems beyond \(A>4\) may not be bound in the unitarity limit. To further assess this situation one should investigate whether these states can be found as resonance (or virtual state) poles at unitarity, and if so, if these poles are perturbatively close to the situation in the real world.

In the bosonic sector, the promotion of many-body forces to lower orders than where they would be expected according to NDA is a fascinating consequence of universality [35], but it does impose practical limitations on many-body calculations. Beyond four nucleons the influence of Fermi statistics is expected to become important, which will most likely limit the promotion of A-nucleon forces with \(A>4\). However, at the same time one should wonder how much this might constrain the usefulness of universality in general, and at which point the Fermi momentum becomes a relevant scale for the description of nuclei. A calculation of, for example, n-\(\alpha \) scattering within the unitarity expansion will be an important next test of the framework and help assess its exact place in the tower of nuclear EFTs.