Energies and radii of light nuclei around unitarity

Light nuclei fall within a regime of universal physics governed by the fact that the two-nucleon scattering lengths are large compared to the typical nuclear interaction range set by one-pion exchange. This places nuclear physics near the so-called unitarity limit in which the scattering lengths are exactly infinite. Effective field theory provides a powerful theoretical framework to capture this separation of scales in a systematic way. It is shown here that the nuclear force can be constructed as a perturbative expansion around the unitarity limit and that this expansion has good convergence properties for both the binding energies of $A=3,4$ nuclei as well as for the radii of these states.


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 lowenergy two-nucleon system can be parametrized with a formula that has become famous, in nuclear physics and beyond, as the effective range expansion (ERE): Here δ 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 (NN ) cross section at low energies and completely determines it in the limit where the relative momentum k of the nucleons goes to zero: Empirically, in the 3 S 1 ("t") and 1 S 0 ("s") NN spin channels the values are known to be a t 5.4 fm and a s −23.7 fm, respectively. Compared to the typical range of the nuclear interaction, set by one-pion exchange providing the longest-range component as R ∼ M −1 π 1.4 fm, these scattering lengths are unnaturally large, a s,t R. 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]), or certain mesons which can be interpreted as hadronic molecules [11].
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 momentafall 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 + it with the pole condition is cot δ 0 (k) = i for some k = iκ. Keeping only the first term in the ERE this gives κ = 1/a and one sees that for positive (negative) a one has a bound (virtual) state in the system. Having a sufficiently large |a| ensures that indeed these momentum scales lie within arXiv:1910.12627v1 [nucl-th] 28 Oct 2019 the radius of convergence of the ERE, which for nucleons is determined by the position of the pion cut, M π /2. For the two-nucleon system at the physical point one has the deuteron as a shallow bound state (B D 2.224 MeV, with the difference to 1/(M N a 2 t ) ≈ 1.41 MeV being due to range corrections) in the 3 S 1 channel, and a very shallow virtual state at B NN * 0.068 MeV (with a relatively small range correction since |a s | is so large). In the unitarity limit, a s,t → ∞, both of these poles become zeroenergy 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 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 Efimov effect [12]. 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. [13,14,15] that for physical values of the NN 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 [16] that a virtual state in the three-nucleon (3N ) system, known to exist for a long time [17,18], 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 [19].
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 [20]. One of these is almost five times as deeply bound as the trimer, B 4 /B 3 4.611, while the other resides just below the particle-trimer threshold, B 4 * /B 3 1.002 [21]. Universality implies that if the NN scattering lengths were infinite, the ground state of 4 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 α /B H 3.66, and there exists a 0 + resonance state just above the proton-triton threshold, i.e., one has B α * /B H 1.05, where B H 7.72 MeV is the 3 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 con-tributed to. A recent review of EFTs that use nucleons and mesons as degrees of freedom can be found in Ref. [22]. 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 [23,24] by considering charge radii of light nuclei in addition to binding energies. After the setup and implementation of the expansion as well as results are presented in the following sections, Sec. 5 will come back to 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 where the notation has been taken over from Refs. [25,23]. The degrees of freedom are nonrelativistic nucleon isospin doublets N = (p n) T , coupled to photon fields where e is the proton charge and τ 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 parametrized in Eq. (4) by the "low-energy constants (LECs)" C 0,i and D 0 , defining contact (zero-range) interactions without derivatives between two and three nucleons, respectively. The P i denote projectors onto the NN S waves, i = 1 S 0 , 3 S 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 hi are included in the prefactor. For the standard Pionless EFT expansion, M hi ∼ R −1 ∼ M π , and this is kept for the construction of the unitarity expansion. However, while standard Pionless EFT assumes Q ∼ 1/a s,t , Ref. [23] suggested to count these scales separately as ℵ ∼ 1/a s,t while assuming that 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 ℵ < Q < 1/R such that it is possible to set up a combined expansion in two parameters ℵ/Q and QR. Coulomb effects are perturbative for momenta of order Q A as well and naturally captured by the expansion if one takes into account that the Coulomb momentum scale k C = αM N with the finestructure constant α ≈ 1/137, is naturally included in the ℵ scale [25]. 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 possible to write where C (0) 0,i is the leading-order (LO) piece of the nonderivative contact LEC C 0,i in Eq. (4), which has an expansion of the form Apart from this, V 2,i is defined in terms of a separable Gaussian regulator function, given by in momentum space. This makes the zero-range theory well defined by regularizing the otherwise divergent interaction via the introduction of a cutoff scale Λ. Both the value of Λ 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 where 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: The regulator ensures that g|G 0 |g is finite. At the onshell point, E = k 2 /M N and k = k , this solution can be matched directly to the ERE, yielding With this running coupling appearing in Eq. (6), the twobody sector of the theory is renormalized. The number θ 0 in general depends on the choice of regulator; for the Gaussian form used here one finds θ 0 = 1/ √ 2π. In the unitarity limit, 1/a i = 0, such that the leading-order twobody interaction is parameter free: Perturbative higher orders are defined by formally expanding the full T matrices t i as where t (0) i is defined by Eq. (12). The corrections t (n) i for n > 0 can conveniently be obtained by solving similar integral equations [26,24]. For the unitarity expansion, corrections from the finite scattering length enter at NLO via C (1) 0,i , yielding a separable potential V (1) 2,i with the same form as Eq. (6). For t (1) i , this gives rise to which, just like the LO equation, can be solved algebraically (see Ref. [24] for explicit details). From this procedure one obtains Range corrections enter at NLO together with C 0,i . 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,i g(k 2 ) k 2 + k 2 g(k 2 ). By virtue of this still being a separable interaction, the corresponding version of Eq. (16) with can still be solved algebraically. Matching the result to the ERE (1) up to the quadratic term gives with θ 2 = θ 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 [26,24]. At second order, the T -matrix correction is obtained from For a perturbative treatment of Coulomb contributions, which are neglected in this work, see Refs. [25,24]. Leading order is however not complete with only the C (0) 0,i interactions. It is a distinct feature of Pionless EFT, intimately related to the Efimov effect [13,14,15], 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 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, at LO, where | 3 H projects onto a J = T = 1 /2 threenucleon state and the regulator |ξ is defined, for Jacobi momenta The k i label the individual nucleon momenta. An NLO correction V has the same form as Eq. (21), but involves the LEC D

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 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 a further, with key results explained 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 where the k i label the individual nucleon momenta, conjugate to position vectors x i . Projecting these momenta onto partial waves yields states |u 1 u 2 ; s , where collects angular momentum, spin, and isospin quantum numbers. They are coupled such that (l 1 s 1 )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 the Faddeev equation [27] where |ψ (0) = |ψ (0) (12)3 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 where |ψ is an auxiliary amplitude, and 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 with and alternatively a similar kernel can be obtained from Eq. (25) in terms of V 3 . Either form of the Faddeev equations is solved by representing it within the space of states |u 1 u 2 ; s , 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 , |ψ (0) can then be determined as the corresponding eigenstate and finally one constructs the full wavefunction as If one starts from Eq. (25), there is no explicit three-body Faddeev component and one simply has |Ψ (0) = (1 + P )|ψ (0) . Fundamentally, the permutation operator P leads to a coupling of different partial waves (see Ref. [29] 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 |Ψ (0) is only recovered as more and more states are included, which means that in principle all observables calculated from |Ψ (0) have to be checked for convergence.
The full wavefunctions |Ψ (0) are used to calculate both perturbative shifts for the binding energy, as well as the radius at leading order (see Sec. 4). Note that |Ψ (0) is assumed here to be properly normalized, In calculating the matrix elements in Eqs. (31) and (32) it is advantageous to exploit the antisymmetry of |Ψ (0) 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 Eqs. (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 [28], and using the procedure described in the previous paragraph it is furthermore possible to eliminate all intermediate higher partial-wave components int the calculation of B 1 . These simplifications are used in the practical implementation to fit the three-body LECs D Describing the four-nucleon system requires an additional Jacobi momentum as well as an alternative set of momenta (v 1 , v 2 , v 3 ), describing a 2+2 cluster setup, i.e., v 1 = u 1 , v 3 denotes the relative momentum in the (34) system, and v 2 is defined as the relative momentum between the (12) and (34) subsystems.
Including the remaining quantum numbers, this leads to channel states The |a are a straightforward extension of three-nucleon states (24), including the angular momentum l 3 associated with u 3 as well as spin and isospin 1 2 for the fourth nucleon. For the b states, (λ 1 , σ 1 , τ 1 ) and (λ 3 , σ 3 , τ 3 ) are two-body quantum numbers for the (12) and (34) subsystems, respectively, whereas λ 1,2,3 are the angular momenta associated with v 1,2,3 . For 4 He the total spin and isospin are J = T = 0.
Following Refs. [28,30], the Faddeev-Yakubovsky equations can be written as corresponding to the decomposition of the full four-body wavefunction. The two distinct Faddeev-Yakubovsky components |ψ A and |ψ B 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 , and the latter in terms of v i and |b . Note that the same notation is used here as for the three-body case, and G 0 in Eqs. (35) now represents the free four-body Green's function. In addition to the operator P already encountered in the three-body system, Eqs. (35) include the further permutations P 34 andP = 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 |ψ A to construct fully symmetric states. The structure of Eqs. (35) can be made clearer by rewriting them in a generic matrix form, with In Eq. (38), G 0 and t (0) are understood to be diagonal matrices, and the permutation operators are collected in From this form, the structural analogy to the three-body Faddeev case becomes obvious. Just like for the Faddeev equations, Eqs. (35) are solved by projecting onto states |u 1 u 2 u 3 ; a , |v 1 v 2 v 3 ; b [30], 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 Eqs. (34) (by choosing all total angular momenta j i and ι i less than some j max ) and study the numerical convergence of results as j max is increased. 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 κ t = 1/a t in the 3 S 1 channel. As demonstrated explicitly in Ref. [24] by using the perturbative formalism discussed in Sec. 1, this implies that the deuteron remains at zero energy at NLO and only moves to 1/(M N a 2 t ) 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,i , 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  [23]. Details regarding the perturbative treatment of Coulomb effects are discussed in Refs. [25,24].
The 4 He nucleus provides a more serious test of the unitarity expansion. Since Q A R ∼ 0.8 for 4 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, ℵ/Q A , should indeed work better with increasing Q A . Figure 1 shows the 4 He binding energy as a function of the momentum cutoff Λ. The observed convergent behavior as Λ increases indicates that the EFT calculation is properly renormalized, as established originally in Refs. [28,31,32]. Results for the standard pionless LO, given by the (blue) diamonds in Fig. 1, are consistent with this earlier work.
While any Λ above the breakdown scale (of order M π ) is a valid cutoff choice in principle, polynomials in 1/Λ are fitted to the points in Fig. 1 to quantitatively assess the convergence and conveniently extrapolate Λ → ∞. This procedure gives B α = 39(12) MeV in the unitarity limit, with the uncertainty estimated as O(r s,t /a s,t ) 30% based on the expectation that range effects are the dominant correction. Including the finite-scattering lengths as NLO corrections gives the (green) circles in Fig. 1 up to this order, and indeed the extrapolated result 30 (9) MeV comes out very near the standard pionless LO value of 31 (9) MeV. 1 It should be stressed that an NLO including only finitea 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 He is rather deeply bound, it turns out that the inclusion of range effects actually has a profound consequence: in Ref. [34] 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 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 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 to be 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 [35]. Its existence is explained by the three-body parameter largely governing the physics of the system [32]. 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 3. Observables fixed 1 Note that the NLO data points at finite cutoffs shown in Fig. 1 differ slightly from the results shown in Refs. [23,33]. This is due to a small error that has been fixed in the numerical implementation. Incidentally, the NLO result in the limit Λ → ∞ is almost unaffected by this, with the value merely changing from the 29.5 MeV reported in Ref. [23] to 30.2 MeV. The difference is negligible compared to the 9 MeV uncertainty estimated at this order. Therefore, all conclusions of the previous work remain unchanged.
as input data at a given order are shown as underlined text. This in particular included the 3 He binding energy at N 2 LO since according to Refs. [25,24] an isospin-breaking three-body is required once perturbation theory mixes Coulomb effects and range corrections. 2 Taking into account the findings of Ref. [34], the 4 He binding energy should be underlined for a complete NLO as well.

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 for the 3 H and 4 He ground states, from which one obtains point charge radii as In Eq. (40), the total charge operatorρ ≡ J 0 , i.e., the zero component of the electric current J µ , is given by the sum of the individual nucleon contributions, The (anti-)symmetry of the wavefunction makes it possible to replace this sum by Aρ i for any fixed i. A particularly convenient choice for three nucleons is i = 3 because it holds that where r 2 is the relative distance conjugate to u 2 and 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 u 2 . Likewise, for four nucleons a good choice is i = 4 because with one obtains a momentum transfer only onto u 3 .
To useρ within the Faddeev-Yakubovsky framework, it is necessary to represent it within the appropriate partialwave basis. Sinceρ does not depend on spin, two-body matrix elements ofρ are given by u; (ls)jm|ρ(q)|u ; (l s )j m = δ jj δ mm u; l||ρ(q)||u ; l , where the reduced matrix element on the right-hand side is given by without the δ mm , where A detailed derivation of Eq. (46) is provided in Appendix A. Embeddingρ 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):  At leading order in the unitarity expansion, the form factor is given by and analogously Eq. (41), with added superscripts "(0)," yields the LO point charge radii. At NLO, the correction to the form-factor is 3 where the factor 2 follows from symmetry. Perturbatively expanding Eq. (41) gives r 0 (1) = 1 2 with r 2 0 (1) calculated from the slope of F (1) 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 threeand 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 |Ψ 1 = (1+P )|ψ 1 , with the NLO Faddeev component 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 with the kernelK (0) as defined in Eq. (38) and  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. [28] at leading order.
As for the 4 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 H and 4 He are shown in Fig. 3 as a function of the UV cutoff Λ. Convergence as Λ increases is evident from the plot, and just like it was done for the binding energies polynomials in 1/Λ are fitted to the data points in order to extrapolate Λ → ∞. The horizontal lines in Fig. 3 show the experimental values of the point charge radii, which, following Ref. [37], are defined as for the triton, and for 4 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. [38] for the quantities appearing on the left-hand sides of the above definitions one obtains, with error bars negligible compared to those of the present theoretical calculation, r 2 0 exp 3 H = 1.59 fm and r 2 0 exp 4 He = 1.72 fm. The lower part of Fig. 3 shows results for the triton. In the limit Λ → ∞, indicated as points on the right border of 4 Note that Eq. (53) makes explicit use of the fact that the three-body force considered here is symmetric under all permutations, such that V indicating that the unitarity expansion works well for this observable. This is in line with the results of Ref. [39], where good convergence is found for a perturbative expansion of 3 H and 3 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. [37], at unitarity the radius satisfies well the universal relation [40,37,39] As done for binding energies it is assumed here that the Q A R part dominates the overall expansion, yielding a 30% uncertainty both at LO as well as NLO. Indeed, from Ref. [37] 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 experiment 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 A R and ℵ/Q A expansions coherently.
Results for 4 He radius, shown in the upper part of 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 He than for the triton is in good agreement with ℵ/Q 4 < ℵ/Q 3 obtained from Eq. (5), and it therefore reinforces confidence in this estimate.

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 π = 0 parametrizing the distance from it. This combination yields a power counting for Q ∼ M π which does not easily capture the physics of the regime Q M π . Notably, one-pion exchange only contributes to the NN 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 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 obtain for the 3 H and 4 He point charge radii as well. This finding solidifies the picture drawn in the introduction of this work, placing fewnucleon 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 ℵ/Q A and Q A /M hi = Q A R. So far, this has been investigated only for the 3 H-3 He energy splitting, where by construction range corrections cancel at NLO [23]. An isospin-breaking threenucleon force is required once range corrections mix with Coulomb contributions at N 2 LO [24]. Range corrections are known to significantly contribute to the triton point charge radius [37]. For the 4 He radius, the closeness to the experimental point already without range corrections found in this work leaves little room for a significant shift at full NLO. From Ref. [34] it is known that such a calculation will require a four-nucleon force to be included, which is most conveniently fit to reproduce the 4 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 1 S 0 N N scattering lengths. This is expected to be a small effect for a bound state as deep as 4 He, but it is interesting to note that fitting the four-nucleon force to exactly reproduce the 4 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 [41,42].
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 [43] and recent calculations of nuclear matter using interactions guided by unitarity [44] provide reason to be optimistic that universality, and in particular discrete scale invariance [45], is able to inform more than just fewnucleon calculations. On the other hand, Refs. [46,47] 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 [34], 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-α 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. A Partial-wave decomposition of charge operators A generic one-body charge operator between two-body states |u; lm can be written as u; lm|ρ(q)|u ; l m = d 3 p d 3 p u; lm|p p|ρ(q)|p p |u ; l m The dependence on the angles of q can be isolated by using a procedure analogous to the one described in Ref. [29] for the permutation operators appearing in the Faddeev equations. It holds that and furthermore wherek = 2k + 1.
In these expressions, Y LM l1,l2 is used to denote two coupled spherical harmonics. The resulting product can be reduced: At this point, it is possible to perform the integral over p in Eq. (61), yielding: (65) Sinceρ is a scalar operator that cannot connect l = l, only f 2 = m 2 = 0 can contribute and Y * f2m2 (q) reduced to a factor 1/ √ 4π. This leads to a cascade of further simplifications: C l m lm,00 = δ ll δ mm (66a) 0 l l λ 1 λ 2 k = (−1) l +λ 1 +λ 2 δ ll δ λ 2 k / lk (66b) Putting everything together leads to Eq. (46) in the main text. The next step is embedding the current into the threenucleon system. To that end one decouples the states |s defined in Eq. (24) to isolate the u 2 part: Taking matrix elements, it is possible to exploit thatρ is diagonal, so one simply gets a number of Kronecker deltas from the reduction of the Clebsch-Gordan coefficients. In particular, For the isospin part, one finds such that overall one arrives at Eq. (48). Analogously, for the four-nucleon system, the first step is decoupling the u 3 part from the states (34a). Omitting the intermediate quantum numbers for the (123) subsystem as well as the isospin part, one finds leading to Eq. (49).

B Perturbative expansion of few-body bound states
This appendix discusses methods to obtain perturbative corrections for few-body wavefunctions. In the main text of this work, these corrections are used to calculate firstorder shifts for three-and four-nucleon charge radii within the unitarity expansion, but having access to wavefunctions generally enables a variety of further calculations. For example, second-order corrections to binding energies can be obtained in a way that is numerically much simpler than the procedure used in Ref. [24] to extract such second-order shifts from off-shell T -matrix corrections. An abstract operator notation is used throughout this section, and conventions for sub-and superscripts differ from the usage in the main text in order to simplify the notation. After a general discussion that does not assume any fixed number of particles, concrete first-order equations for three-and four-body states are derived in Sec. B.2. While explicit results are given for first-order perturbation theory, it is clear from the following discussion that higher-order equations can be derived analogously in a recursive fashion, much like it is done in Ref. [37] for a calculation of the triton charge radius up to next-to-next-to leading order. The derivations presented here are built on the concept of using inhomogeneous equations to calculate perturbative corrections, introduced in Ref . [37] specifically for three-body vertex functions, and previously in Ref. [26] for scattering calculations.

B.1 Generic discussion
In principle, one can base a perturbative corrections for bound-state wavefunctions on the corresponding expansion of Lippmann-Schwinger equation following the discussion in Sec. 1. Using the fact that as the energy approaches a bound-state at E = −B, the T matrix has a pole: In this expression, the vertex function |B is related to the wavefunction |Ψ via In order to develop a formalism that is connected to the bound-state Faddeev formalism, it is however more instructive to start directly from the Schrödinger equation. Expanding (74c) gives the well-known first-order equation In principle, an explicit solution is given by where S denotes the whole spectrum and the sum is generally a sum over discrete states plus an integral over the continuous spectrum (note that for simplicity it is assumed here that there are no degeneracies in the spectrum). Using Eq. (76) however requires a complete diagonalization of the leading-order Hamiltonian, which is not convenient for calculations based on Faddeev-and Faddeev-Yakubovsky equations that by construction only give access to specific individual states. Moreover, for the systems considered in this work almost all contributions in Eq. (76) would come from discretized continuum of scattering states, which is numerically challenging. It is therefore interesting and relevant to look for alternative ways of solving for |Ψ 1 . As a first step, one can rewrite Eq. (75) as and recognize from this that Eq. (76) involves the full leading-order Green's function in spectral representation, with the bound-state |Ψ 0 subtracted. This subtraction is crucial because implies the operator on the right-hand side of Eq. (77) is singular. It is the term −B 1 |Ψ 0 on the right-hand side that generates this subtraction. With this insight it is possible to derive methods to deal with the problem. Using the definition of the free Green's function, it is possible to further rewrite Eq. (77) as with z fixed at −B 0 . For any z = −B 0 , the kernel on the left-hand side of Eq. (80) is regular and can therefore be solved, after discretization, as a linear system of equations. Knowing from Eq. (76) that the solution is actually well defined in the limit z → −B 0 , numerical extrapolation can be used to obtain |Ψ 1 . An alternative procedure is to consider a potential with the leading-order bound state removed. This can be achieved using the replacement [48,49,50] in Eq. (80), with λ a large positive constant that moves the bound state far away from the low-energy spectrum we are interested in. This procedure is valid because the orthogonality condition for the first-order correction, Ψ 1 |Ψ 0 = 0, implies that the modified equation is equivalent to original one. Numerically, it gives excellent agreement with the extrapolation method for a two-body test case.

B.2 Faddeev and Faddeev-Yakubovsky decomposition
For a three-body state one can write, in analogy to Eq. (30), |Ψ 1 = (1 + P )|ψ 1 and insert this decomposition into Eq. (77). Three-body interactions are neglected here to keep the discussion as transparent as possible. Application of the remaining two-body interaction can be simplified by exploiting the symmetry of the full wavefunction, which implies that whereṼ is the potential acting only on the specific pair of particles used to define the Faddeev component |ψ . This gives: Since (1+P ) commutes with the free Hamiltonian H 0 , one can bring an application of (1 + P ) to the left of each term in this equation. While (1 + P ) can in general have zero eigenvalues, components that are in ker(1 + P ) are clearly irrelevant for the description of the physical bound state and it is therefore possible to proceed with a simplified equation: where G 0 (−B 0 ) −1 = −B 0 − H 0 . Acting with G 0 (−B 0 ) on both sides of this one can perform a partial inversion by using the Lippmann-Schwinger equation for the full two-body Green's function G associated with V . This can then be eliminated in favor of the corresponding t matrix by using for t =Ṽ +Ṽ G 0 t. Overall one arrives at: . This is an inhomogeneous integral equation which involves exactly the same kernel as the Faddeev equation at leading order, whereas the terms on the left-hand side are straightforward to calculate from known quantities. As a final step one calculates |Ψ 1 from |ψ 1 in exactly the same way |Ψ 0 is calculated from |ψ 0 . By keeping three-body forces in the derivation one analogously obtains Eq. (53) given in the main text.
It is useful to note that the same result can alternatively be derived directly from a perturbative expansion of the Faddeev equation setting |ψ = |ψ 0 + |ψ 1 + · · · . To proceed it is important to carefully regard the energy arguments of the operators, noting that B = B 0 + B 1 + · · · : The term involving the energy derivative in Eq. (89b) looks peculiar, but, assuming the potential to be energy independent, it can be shown by differentiating the Lippmann-Schwinger equation that d dz t 0 (z) = −t 0 (z)G 0 (z) 2 t 0 (z) .
Using this and making use of the leading-order Faddeev equation, to simplify some terms, one obtains [1 − G 0 t 0 P ] |ψ 1 = B 1 (G 0 + G 0 t 0 G 0 )|ψ 0 + G 0 t 1 P |ψ 0 (92) All energy arguments are −B 0 and have therefore been omitted as well. Using again Eq. (91) as well as the Lippmann-Schwinger equation for t 1 , Eq. (16) in the main text, it possible to show that G 0 t 1 P |ψ 0 = (G 0 + G 0 tG 0 )Ṽ 1 (1 + P )|ψ 0 , so that indeed Eq. (92) is equivalent to Eq. (87) The previous result is both reassuring and very useful because it implies that in order to derive perturbative corrections for four-body states one can simply start from the Faddeev-Yakubovsky equations, avoiding a tedious detour through the Schrödinger equation.
A convenient starting point for this calculation is generic matrix form of the Faddeev-Yakubovsky equations, analogous to Eq. (37): with |ψ = (|ψ A , |ψ B ) T and where three-body forces are again neglected for simplicity. As before, one expands all quantities in this equation, carefully keeping track of the energy arguments. This gives All energy arguments in Eqs. (97) are fixed now at −B 0 and have been omitted. Again the term involving t 1 can be rewritten in terms ofṼ 1 as With three-body forces included one obtains the slightly more complicatedK 1 shown in Eq. (55) in the main text.
The fact that within the Faddeev(-Yakubovsky) formalism one has to calculate B 1 with an increasing number of partial-wave components included to check convergence of the numerical calculation renders the extrapolation method discussed in Sec. B.1 unstable beyond the two-body sector. At the same time, the projection method which modifies the potential to remove the leading-order bound state from the low-energy spectrum becomes impractical because for an A-body state Eq. (81) involves an A-body potential, which is expensive to handle computationally (at leading order) for A > 3. Fortunately, it is possible to employ the projection method directly at the level of Faddeev(-Yakubovsky) components and make the replacement 5 in Eq. (87), where K 0 = G 0 t 0 P , and analogously for Eq. (96). While this ensures that the kernel becomes regular, the solution of the modified equation, denoted by |Ψ 1 , will not in general directly provide the correct solution |Ψ 1 (obtained from |ψ 1 by adding the appropriate permutations). In particular, |Ψ 1 may not be orthogonal to |Ψ 0 . It is however possible to simply project out the undesired component by setting which gives the desired the solution. The correctness of the result can be checked by evaluating which follows from Eq. (75), since B 0 is known from the leading-order calculation.