Extracting free-space observables from trapped interacting clusters

TL;DR

This paper develops an improved formula for relating trapped two-particle energy spectra to free-space scattering phase shifts, enabling more accurate extraction of nuclear scattering information from finite trap calculations.

Contribution

It introduces a systematically improved formula for low-energy dynamics in finite traps, facilitating the extraction of nuclear scattering phase shifts from ab initio calculations.

Findings

Provides a new formula valid for finite traps

Enables extraction of nuclear phase shifts from many-body calculations

Establishes effective field theory as a useful framework

Abstract

The energy spectrum of two short-range interacting particles in a harmonic potential trap has previously been related to free-space scattering phase shifts. But the existing formula for this purpose is exact only in the limit of an infinitely shallow trap. Here we provide a systematically improved formula---describing the low-energy dynamics---that enables the use of finite traps. This paves the way for extracting nuclear scattering phase shifts from {\it ab initio} nuclear many-body structure calculations, a long-sought goal in nuclear physics. The derivation establishes effective field theory as a powerful framework for studying the connection between structure information of a trapped system (with two or more sub-clusters) and continuum physics in the fields of both nuclear and condensed-matter physics.

Figures

Tables

Table 1. Table 1: C ¯ i , j subscript ¯ 𝐶 𝑖 𝑗 \overline{C}_{i,j} for the s 1 / 2 subscript 𝑠 1 2 s_{1/2} channel.

	$i=0$	$1$	$2$
$j=0$	$-$0.3536	$-$0.4097	0.004206
$1$	0.7898	$-$0.1756	$-$0.4063
$2$	0.1743	0.2010
$3$	$-$0.01278	0.2949

Table 2. Table 2: C ¯ i , j subscript ¯ 𝐶 𝑖 𝑗 \overline{C}_{i,j} for the p 3 / 2 subscript 𝑝 3 2 p_{3/2} channel.

	$i=0$	$1$	$2$
$j=0$	0.02304	$-$1.870	0.5400
$1$	$-$0.5571	$-$0.8676	0.3118
$2$	0.5609	$-$0.2829
$3$	0.08738	$-$0.03911

Table 3. Table 3: C ¯ i , j subscript ¯ 𝐶 𝑖 𝑗 \overline{C}_{i,j} for the p 1 / 2 subscript 𝑝 1 2 p_{1/2} channel.

	$i=0$	$1$	$2$
$j=0$	0.1390	$-$2.084	0.7400
$1$	$-$0.4427	$-$1.090	0.5210
$2$	0.6225	$-$0.4064
$3$	0.1114	$-$0.08700

Equations

$$p^{2\ell +1}\cot {\delta }_{\ell }(E)$$

$$\sum _{i,j=0}^{\infty }{C}_{i,j}{\left(M_{{}_{\mathrm{R}}}\omega \right)}^{2i}{p}^{2j}$$

$$p^{2\ell +1}\cot {\delta }_{\ell }\left(E\right)$$

$${\mathcal{L}}_0$$

$${\mathcal{L}}_I$$

$$\Sigma (E)=\frac{{\mathcal{A}}_{\ell }}{\pi }{\int }_0^{T_{\Lambda }}d{T}_q\frac{(2M_{{}_{\mathrm{R}}}{T}_q{)}^{\ell +\frac{1}{2}}}{E-T_q+i{0}^{+}}$$

$$\displaystyle\quad\quad\ =-\mathscr{A}_{\ell}\big{[}ip^{2\ell+1}+\sum_{j=0}^{+\infty}L_{\ell,j}(\Lambda)p^{2j}\big{]}\;,$$

$${\mathcal{A}}_{\ell }\equiv \frac{g_{\ell }^2}{M_{{}_{\mathrm{R}}}^{2\ell -1}}\frac{2^{\ell -1}\ell {!}^2}{\pi (2\ell +1)!},L_{\ell ,j}(\Lambda )\equiv \frac{2{\Lambda }^{2\ell -2j+1}}{\pi \left(2\ell -2j+1\right)}.$$

$$D=\frac{1}{{\sigma }_{\ell }(E+{\Delta }_{\ell })-d_j^{(\ell )}{E}^j-\Sigma }=\frac{-{\mathcal{A}}_{\ell }^{-1}}{p^{2\ell +1}[\cot {\delta }_{\ell }-i]},$$

$$\text{with}{p}^{2\ell +1}\cot {\delta }_{\ell }=\sum _{j=0}^{\infty }{C}_{0,j}{p}^{2j},\text{and}$$

$$C_{0,j}=\frac{{\mathcal{A}}_{\ell }^{-1}}{(2M_{{}_{\mathrm{R}}}{)}^j}{\left\{-{\sigma }_{\ell }{\Delta }_{\ell },-{\sigma }_{\ell },d_2^{(\ell )},d_3^{(\ell )}\dots \right\}}_j-L_{\ell ,j}(\Lambda ).$$

$${\Sigma }_{\omega }(E)$$

$$=\frac{{\mathcal{A}}_{\ell }}{\pi }{\left(4M_{{}_{\mathrm{R}}}\omega \right)}^{\ell +\frac{1}{2}}\sum _{n=0}^{n_{\Lambda }}{f}_{\ell }\left(z_{{}_E},n\right),$$

$$f_{\ell }\left(z_{{}_E},n\right)$$

$$D_{\omega }=\frac{1}{{\sigma }_{\ell }(E+{\Delta }_{\ell })-d_j^{(\ell )}{E}^j-{\Sigma }_{\omega }(E)-d_{j,k}^{(\ell )}{E}^j(M_{{}_{\mathrm{R}}}\omega {)}^{2k}}=\frac{(-){\mathcal{A}}_{\ell }^{-1}}{p^{2\ell +1}\cot {\delta }_{\ell }+\frac{1}{{\mathcal{A}}_{\ell }}\left[{\Sigma }_{\omega }(E)-\mathcal{P}\Sigma (E)+d_{j,k}^{(\ell )}{E}^j(M_{{}_{\mathrm{R}}}\omega {)}^{2k}\right]}.$$

$$p^{2\ell +1}\cot {\delta }_{\ell }(E)+\frac{d_{j,k}^{(\ell )}}{{\mathcal{A}}_{\ell }}{E}^j(M_{{}_{\mathrm{R}}}\omega {)}^{2k}=\frac{\mathcal{P}\Sigma (E)-{\Sigma }_{\omega }(E)}{{\mathcal{A}}_{\ell }}.$$

$$\displaystyle-\frac{1}{\pi}\left(4M_{{}_{\mathrm{R}}}\omega\right)^{\ell+\frac{1}{2}}\bigg{[}\sum_{n=0}^{n_{\Lambda}}f_{\ell}\left(z_{{}_{E}},n\right)+\pi\sum_{j=0}^{\ell}z_{{}_{E}}^{j}L_{l,j}\!\!\left(\!\!\sqrt{\tfrac{\bar{T}_{\Lambda}}{2}}\right)\bigg{]}$$

$$\equiv -\frac{1}{\pi }{\left(4M_{{}_{\mathrm{R}}}\omega \right)}^{\ell +\frac{1}{2}}\sum _{n=0}^{(\mathcal{R})}{f}_{\ell }\left(z_{{}_E},n\right),$$

$$\sum _{n=0}^{(\mathcal{R})}{f}_{\ell }(z,n)=(-{)}^{\ell }\pi \frac{\Gamma \left(\frac{\ell }{2}+\frac{3}{4}-z\right)}{\Gamma \left(\frac{1}{4}-\frac{\ell }{2}-z\right)},$$

$$\left[-\frac{d^2}{d{\bar{r}}^2}+\frac{\ell (\ell +1)}{{\bar{r}}^2}+{\bar{r}}^2\right]u_{\ell }=2\bar{E}{u}_{\ell },$$

$$u_{\ell }$$

$$M\left(a,b,z\right)\approx \Gamma \left(b\right)\left[\frac{z^{a-b}{e}^z}{\Gamma \left(a\right)}+\frac{(-z{)}^{-a}}{\Gamma \left(b-a\right)}\right]\left[1+O(z^{-1})\right],$$

$$c_1\frac{\Gamma \left(\ell +\frac{3}{2}\right)}{\Gamma \left(\frac{\ell }{2}+\frac{3}{4}-\frac{\bar{E}}{2}\right)}+c_2\frac{\Gamma \left(-\ell +\frac{1}{2}\right)}{\Gamma \left(-\frac{\ell }{2}+\frac{1}{4}-\frac{\bar{E}}{2}\right)}=0.$$

$$c_1{\bar{r}}_c^{\ell +1}M$$

$$+$$

$$\displaystyle\!\!\!\!\!\frac{\left(2\bar{r}_{c}\right)^{2\ell+1}}{(-)^{\ell}\mathscr{N}_{\ell}}\frac{\Gamma\!\left(\frac{\ell}{2}+\frac{3}{4}-\frac{\bar{E}}{2}\right)}{\Gamma\!\left(\frac{1}{4}-\frac{\ell}{2}-\frac{\bar{E}}{2}\right)}\!=\!\frac{M\!\left(\frac{1}{4}\!-\!\frac{\ell}{2}\!-\!\frac{\bar{E}}{2},\frac{1}{2}\!-\!\ell,\bar{r}_{c}^{2}\right)}{M\!\left(\frac{\ell}{2}\!+\!\frac{3}{4}\!-\!\frac{\bar{E}}{2},\ell\!+\!\frac{3}{2},\bar{r}_{c}^{2}\right)}.\$$

$$\frac{{\left(p{r}_c\right)}^{2\ell +1}}{(-){\mathcal{N}}_{\ell }}\frac{y_{\ell }(p{r}_c)}{j_{\ell }p{r}_c)}\stackrel{?}{=}\frac{M\left(\frac{1}{4}-\frac{\ell }{2}-\frac{\bar{E}}{2},\frac{1}{2}-\ell ,{\bar{r}}_c^2\right)}{M\left(\frac{\ell }{2}+\frac{3}{4}-\frac{\bar{E}}{2},\ell +\frac{3}{2},{\bar{r}}_c^2\right)}.$$

$$1+\frac{(2\ell +1){\left(p{r}_c\right)}^2}{4\ell (\ell +1)-3}+\frac{(\ell +3)(2\ell +1){\left(p{r}_c\right)}^4}{(2\ell -3)(2\ell +3{)}^2(2\ell +5)}+\cdots$$

$$\underset{\omega \to 0}{\lim }M$$

$$=\underset{\omega \to 0}{\lim }M\left(-\frac{E}{2\omega },-\ell +\frac{1}{2},r_c^2{M}_{{}_{\mathrm{R}}}\omega \right)$$

$$={}_0{F}_1\left(;-\ell +\frac{1}{2};-\frac{1}{2}{M}_{{}_{\mathrm{R}}}E{r}_c^2\right),$$

$$\underset{\omega \to 0}{\lim }M$$

$$={}_0{F}_1\left(;\ell +\frac{3}{2};-\frac{1}{2}{M}_{{}_{\mathrm{R}}}E{r}_c^2\right),$$

$$J_{\alpha }(z)=\frac{{\left(\frac{z}{2}\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}{}_0{F}_1\left(;\alpha +1;-\frac{z^2}{4}\right).$$

$$j_{\ell }(z)=\sqrt{\frac{\pi }{2z}}{J}_{\ell +\frac{1}{2}}(z),y_{\ell }(z)=(-{)}^{\ell +1}\sqrt{\frac{\pi }{2z}}{J}_{-\ell -\frac{1}{2}}(z),$$

$$\frac{{}_0{F}_1\left(;-\ell +\frac{1}{2};-\frac{1}{2}{M}_{{}_{\mathrm{R}}}E{r}_c^2\right)}{{}_0{F}_1\left(;\ell +\frac{3}{2};-\frac{1}{2}{M}_{{}_{\mathrm{R}}}E{r}_c^2\right)}$$

$$\displaystyle\frac{(-1)^{\ell+1}}{\mathscr{N}_{\ell}}\Bigl{[}$$

Full text

Introduction

Nuclear experiments at low energy can not manipulate many-body systems to the extent possible in condensed-matter or cold-atom experiments. However, with progress in many-body methods Barbieri and Carbone (2017); Barrett et al. (2013); Carlson et al. (2015); Lee (2009); Stroberg et al. (2019) and increasing computing power (and quantum computers Preskill (2018)), one can start manipulating nuclear systems computationally. Here, I show how trapping two clusters at low energy in a harmonic potential well tells us about their free-space scattering through a formula connecting low-energy phase shifts with the confined spectrum. In this approach the trap compacts the system and reduces the required degrees of freedom enough to allow controlled ab initio calculations, as will be demonstrated elsewhere. (See e.g., Nollett et al. (2007); Navrátil et al. (2016); Elhatisari et al. (2015); Shirokov et al. (2018) for other ab initio approaches of computing light-nucleus scatterings.)

A formula for particles in a harmonic-potential trap was derived in Ref. Busch et al. (1998), and later generalized to include the full energy dependence of the phase-shift (besides the scattering length term in Busch et al. (1998)) and for partial waves beyond s-wave Blume and Greene (2002); Block and Holthaus (2002); Bolda et al. (2002); Idziaszek and Calarco (2006); Stetcu et al. (2007); Suzuki et al. (2009); Stetcu et al. (2010); Rotureau et al. (2010); Blume (2012). The result for angular momentum $\ell$ Suzuki et al. (2009); Stetcu et al. (2010); Luu et al. (2010) (called the BERW formula here) is

[TABLE]

This holds at the eigenenergies $E\equiv$ ${p^{2}}/2M_{{}_{\mathrm{R}}}$—with the center-of-mass (CM) energy subtracted—in a trap where each particle experiences a potential $\omega^{2}\bm{r}^{2}/2$ times its mass; $M_{{}_{\mathrm{R}}}$ is the reduced mass and $\delta_{\ell}$ the phase shift. Equation (1) is analogous to the Luscher formula Luscher (1991); Beane et al. (2004) that is widely applied in Lattice Quantum Chromodynamcs (for a system on a space-time torus).

Refs. Luu et al. (2010); Rotureau et al. (2012) have used Eq. (1) to extract nuclear scattering from ab initio spectrum calculations.

However, away from the infinitely-shallow-trap limit (i.e., for $\omega\neq 0$), Eq. (1) does not capture the external potential’s modifications to the interaction at short distances. To illustrate the impact on extracting phase shifts, I use a two-body potential model Ali et al. (1985) designed for describing neutron-$\alpha$ scattering [see the supplemental materials (SM) for details]. Figure 1(a) shows $3/2{}^{-}$ p-wave phase shifts extracted using Eq. (1) at the eigenenergies with $\omega=3,4,6,9,16$ MeV (typical values applicable in ab initio calculations): they fail to align on a smooth curve and systematically deviate from the exact curve Luu et al. (2010).

Here, I remedy the BERW formula using pionless effective field theory (EFT) van Kolck (1999); Bedaque and van Kolck (2002); Hammer et al. (2017), which enables low-energy dynamics to be studied without specifying the details of the short-distance physics (e.g., potential or cluster structure and excitation). This EFT was used to re-derive and generalize the Luscher formula Beane et al. (2004); Briceno et al. (2013). The improved formula for a harmonic trap is

[TABLE]

The constants $C_{i,j}$ depend implicitly on $\ell$ but are independent of $\omega$ and $p$; they are dimensionful and scale as proper powers of a high-momentum scale $M_{H}$ (as dictated by, e.g., the cluster excitations), unless there is fine tuning. When $\sqrt{M_{{}_{\mathrm{R}}}\omega}$ and $p$ are smaller than $M_{H}$, the series sum converges and thus can be truncated with a controlled error. However, outside the convergence domain, where the details of the finite-range interaction and its interplay with the trap potential are being probed, the series expansion becomes infeasible. It is also worth noting that in the case of physically zero-range interactions (i.e., $M_{H}\rightarrow\infty$), the $C_{i\neq 0,j}$ terms, which capture the trap-induced modifications, would vanish, and Eq. (2) becomes equivalent to Eq. (1).

To infer the phase shifts from Eq. (2) given the eigenenergies, the $C_{i\neq 0,j}$ terms must be simultaneously calibrated with the $C_{0,j}$. The latter determine the free-space phase shifts via the effective range expansion (ERE) van Kolck (1999); Hammer et al. (2017). Knowing the full potential in the n-$\alpha$ model, one can fix $C_{i,j}$ (see the SM) and generate Fig. 1(b): the inset shows that the phase shifts extracted from Eq. (1) for a given $\omega$ sit on a curve parameterized by a generalized ERE, in which the $j^{\textrm{th}}$-order coefficient is given by $\sum_{i=0}C_{i,j}(M_{{}_{\mathrm{R}}}\omega)^{2i}$. After subtracting the trap-induced modifications, the extracted phase shifts agree with the “True” curve.

The essence of Eq. (2), that the trap-induced modifications can be parameterized using a Taylor expansion in the $\omega^{2}$ and $p^{2}$ variables, can be seen in Fig. 1(c). The symbols show the differences between the right-hand-side (RHS) and left-hand-side (LHS) in Eq. (1)—scaled by $M_{H}^{-3}$—at the $\omega$-dependent eigenenergies, while the solid lines plot the summation of the $C_{i\neq 0,j}$ terms in Eq. (2)’s LHS with $i\leq 2$ and $0\leq j\leq 3$. Indeed, they interpolate those symbols. Of course, outside the convergence domain, the series expansion would fail, as shown in Fig. 1(d). The differences between Eq. (2)’s two sides, based on two series truncations ($j\leq 3$ for “Up to $E^{3}$” and $j\leq 2$ for “Up to $E^{2}$”, and $i\leq 2$ in both cases), are plotted against a large range of eigenenergies. The truncation errors behave as the leading terms left out of the summation in the low-energy region, but then increase to $100\%$ when the symbols reach the shaded region indicating the range of Eq. (2)’s RHS; this also suggests the breakdown scale for $E$ is between $20$ and $40$ MeV. (See the SM for more details on series convergence.)

To extract nuclear phase shifts (or $C_{i,j}$s) from ab initio spectra, Eq. (2) will play a crucial role, because ab initio calculations, developed to computing compact nuclei, have uncontrolled errors when $\omega\rightarrow 0$. To illustrate Eq. (2), two models are used in the SM: a hard-sphere potential model is solved exactly, while the n-$\alpha$ model is studied numerically. The rest of the paper is devoted to the derivation of Eq. (2), emphasizing a new set of interaction vertices between the external potential (or background field) and trapped particles, and renormalization.

Derivation through EFT

I start by constructing an EFT Lagrangian for two spin-[math] particles—for simplicity—in the $\ell^{\textrm{th}}$ partial wave, with a harmonic potential coupled to each particle. The framework is valid at low energies, where the details of the short-distance physics and its interplay with the trap are not resolved. I follow the conventions of Ref. Zhang et al. (2018). Let $c(x)$ and $n(x)$ be particle fields with masses $M_{\mathrm{c}}$ and $M_{\mathrm{n}}$ ($c^{\ast}$ and $n^{\ast}$ are the complex conjugations), while $\phi_{m_{\ell}}$ is the so-called dimer field Kaplan et al. (1998a); Beane and Savage (2001); Bedaque and van Kolck (2002); Bedaque et al. (2003); Braaten and Hammer (2006); Briceno et al. (2013); Hammer et al. (2017) with spin $l$, projection $m_{\ell}$, and mass $M_{nc}=M_{\mathrm{n}}+M_{\mathrm{c}}$ ($\phi^{\dagger m_{\ell}}\equiv\left(\phi_{m_{\ell}}\right)^{\ast}$). The dimer $\phi$ couples to $n$-$c$ and represents the compound system. The background field $\mathscr{B}(x)$ is $m\omega^{2}\bm{x}^{2}/2$ in the Lab frame with $m$ as a reference mass. The $\phi$ propagator—and the related self-energy corrections due to $n$-$c$ multiple scattering—will be the central piece in the derivation: in free space it is directly related to the $n$-$c$ scattering $T$-matrix, while in the trap its poles give the system’s spectrum.

The Lagrangian is $\mathscr{L}_{0}+\mathscr{L}_{I}$, where

[TABLE]

The building blocks of $\mathscr{L}_{0,I}$ are invariant under Galilean transformations, including rotation, translation, and boost. (A relevant discussion on Galilean invariance in EFT can be found e.g., in Ref. Braaten (2015).) In both lagrangians for $\psi=n$, $c$, or $\phi$, the $\psi^{\ast}[i\tilde{\partial}_{t}+\bm{\partial}^{2}/(2M_{\psi})]\psi$ structures with $i\tilde{\partial}_{t}\equiv i\partial_{t}-M_{\psi}\mathscr{B}(x)/m$ are $\psi$’s internal energies (i.e., total energies with kinetic and external potential energies subtracted), and therefore Galilean invariant.

The $g_{\ell}$ coupling in $\mathscr{L}_{I}$ uses $n$-$c$’s relative velocity $\bm{V}$, while $V^{\otimes\ell}$ denotes a rank-$\ell$ operator composed of $\ell$ copies of $\bm{V}$ normalized such that when $m_{\ell}=+\ell$, $[V^{\otimes\ell}]_{m_{\ell}}=[(V^{+1})^{\ell}]^{\ast}$ with $V^{+1}\equiv-(V^{x}+iV^{y})/\sqrt{2}$. This term means $\phi^{\dagger m_{\ell}}$ is coupled to a $n$-$c$ configuration having $\ell$ and $m_{\ell}$ as its relative angular quantum numbers. The $m_{\ell}$ indices are implicitly summed up so that this term is a scalar. (In general, the spin and vector indices need to be properly contracted to form scalars.) In addition, both $\bm{V}$ and spins are invariant under translation and boost, and thus the $g_{\ell}$ coupling preserves Galilean invariance. Note that repeated indices in the Lagrangian (and in Eqs. (6), (8), and (9)) are implicitly summed with specified ranges.

It should be mentioned that the interactions in $\mathscr{L}_{0}$ and $\mathscr{L}_{I}$ with the external potential turned off follow closely previous works using a dimer-field approach Kaplan et al. (1998a); Beane and Savage (2001); Bedaque and van Kolck (2002); Bedaque et al. (2003); Braaten and Hammer (2006); Briceno et al. (2013); Hammer et al. (2017): $\sigma_{\ell}$ ($=\pm 1$), $\Delta_{\ell}$, $g_{\ell}$, and $d_{j}^{(\ell)}$ together reproduce the ERE (see Eq. (6) and Briceno et al. (2013); Beane and Savage (2001)). This approach is equivalent Kaplan et al. (1998a); Bedaque and van Kolck (2002); Braaten and Hammer (2006) to the EFTs without dimer fields (see further discussion in the “further comments” section).

The $d^{(\ell)}_{j\geq 0,k\geq 1}$ terms are also Galilean invariant, considering the external potential $\mathscr{B}$ is a scalar field. However, their specific structures are severely constrained by a unique property of a harmonic potential: the CM of a multi-particle system is decoupled from its internal dynamics Caprio et al. (2012). (It can be understood based on that the external force on the multi-particle’s CM depends only on CM’s location in the harmonic potential well, i.e., not affected by any other degrees of freedom.) In these couplings with ${M_{{}_{\mathrm{R}}}^{2}}\bm{\partial}^{2}\!\mathscr{B}/(3m)=M_{{}_{\mathrm{R}}}^{2}\omega^{2}$, $\bm{\partial}^{2}$ ensures that they only shift the system’s energy by $\bm{r}$-independent but $\omega^{2}$-dependent functions so that the CM behaves as a free particle in traps, i.e., decoupled from internal dynamics.

Besides powers of $\bm{\partial}^{2}\!\mathscr{B}$, the other possible scalar objects built of $\mathscr{B}$ include (1) $\mathscr{B}^{2}$, $\mathscr{B}^{3}$, …(2) $(\bm{\partial}\mathscr{B})^{4}$, $(\bm{\partial}\mathscr{B})^{6}$, …[$(\bm{\partial}\mathscr{B})^{2}$ is proportional to $\mathscr{B}$], and (3) products of (1) and (2). (Derivatives higher than $\bm{\partial}^{2}$ applied on $\mathscr{B}$ would give zero, and therefore are not relevant here.) They would induce external potentials with powers of $\bm{r}^{2}$ higher than 1. Copies of $\bm{\partial}\mathscr{B}$ can also be used to construct tensor (vector) objects, which create anisotropic external potentials that need to be coupled to particles’ momenta or spins. As the result, the new scalar and tensor objects and their products would break the CM–internal-dynamics decoupling if they are present in any interaction terms with particles.

In principle, $\mathscr{B}$ can be coupled to the $\phi^{\ast}nc$ operators (e.g., the $g_{\ell}$ term), which again must take the form of $\left(\bm{\partial}^{2}\!\mathscr{B}\right)^{1,2,\ldots}$. However, these terms can be eliminated by rescaling the $\phi$ field by $1+\#(M_{{}_{\mathrm{R}}}\omega)^{2}+...$ Coleman et al. (1969). Since the rescaling-induced terms are already present as $d^{(\ell)}_{j,k}$ couplings in $\mathscr{L}_{I}$, the trap modification to $g_{\ell}$ is not included.

Lastly in the free space, defining energy relative to the $n$-$c$ threshold sets $\Delta_{c}=\Delta_{n}=0$. Both are modified by $\mathscr{B}$ through “polarization” effects as $\Delta_{\ell}$ by $d^{(\ell)}_{j=0,k}$ couplings, but they only affect the energy-references in traps and for simplicity not shown here.

To compute the propagator of the dimer $\phi$, its self-energy correction due to $n$-$c$ multiple scattering needs to be included. A cut-off on momentum is applied to regularize loops in free space, while in traps the cut-off is applied on the virtual excitation energy van Kolck (1999). (However, for fine-tuned systems other schemes would be preferred, e.g., power divergence subtraction Kaplan et al. (1998b).) Within time-independent perturbation theory Zhang et al. (2018), the one-loop self-energy bubble diagram in free space is $(2\pi)^{3}\delta(\bm{P}-\bm{P}^{\prime})\delta^{m_{\ell}}_{m_{\ell}^{\prime}}\Sigma(E_{{}_{L}},\bm{P})\equiv\langle\phi_{\bm{P}^{\prime}}^{m_{\ell}^{\prime}}|H_{g_{\ell}}\left(E_{{}_{L}}-H_{0}+i0^{+}\right)^{-1}H_{g_{\ell}}|\phi_{\bm{P}}^{m_{\ell}}\rangle$. $H_{0}$ and $H_{g_{\ell}}$ are the Hamiltonians derived from $\mathscr{L}_{0}$ and the $g_{\ell}$ term in $\mathscr{L}_{I}$, respectively Zhang et al. (2018). Both states are plane waves, with $\bm{P}$, $\bm{P}^{\prime}$, $E_{{}_{L}}$, $m_{\ell}$, and $m_{\ell}^{\prime}$ as $\phi$’s momenta and energy in the Lab frame, and its spin projections. (The realtionships between Feynman diagrams and the matrix elements defined here and below can be found in Ref. Zhang et al. (2018).) One then obtain

[TABLE]

$p\equiv\sqrt{2M_{{}_{\mathrm{R}}}\left(E+i0^{+}\right)}$, $T_{q}\equiv{\bm{q}^{2}}/(2M_{{}_{\mathrm{R}}})$, $\Lambda$ is the cut-off on $|\bm{q}|$, and $T_{\Lambda}\equiv\Lambda^{2}/(2M_{{}_{\mathrm{R}}})$. $E\equiv E_{{}_{L}}-{\bm{P}^{2}}/(2M_{nc})$ is the energy in the CM frame. Note $L_{\ell,j>\ell}(\Lambda)\rightarrow 0$ as $\Lambda\rightarrow\infty$.

The fully dressed free-space $\phi$ propagator, which is defined through $(2\pi)^{3}\delta(\bm{P}-\bm{P}^{\prime})\delta^{m_{\ell}}_{m^{\prime}_{\ell}}D(E_{{}_{L}},\bm{P})\equiv\langle\phi_{\bm{P}^{\prime}}^{m_{\ell}^{\prime}}|\left[E_{{}_{L}}-(H_{0}+H_{I})+i0^{+}\right]^{-1}|\phi_{\bm{P}}^{m_{\ell}}\rangle$, with $H_{I}$ from $\mathscr{L}_{I}$, can be computed by summing the self-energy-insertion diagrams due to $\Sigma$ and the $d_{j}^{(\ell)}$ vertices, yielding

[TABLE]

$D$ is related to $\delta_{\ell}$ through the scattering $T$-matrix, which is computed by multiplying $D$ with two $g_{\ell}$-vertices Zhang et al. (2018). The range of the index in $d_{j}^{(\ell)}$ in the implicit sum is fixed in $\mathscr{L}_{I}$, and in the $C_{0,j}$ definition $\{\ldots\}_{j}$ is the $j^{\textrm{th}}$ component of the list and $j$ is not summed.

Now let us turn to the trapped system. Based on $\mathscr{L}_{0}$, one can expand $n$, $c$ and $\phi$ fields using their corresponding harmonic-oscillator wave functions Stetcu et al. (2010). Again note that the $g_{\ell}$ coupling only picks up the $n$-$c$ configuration whose total angular momentum and projection equal those of the CM motion (i.e., $\phi$) and whose relative angular momentum and projection equal the $\phi$’s spin and projection ($\ell$ and $m_{\ell}$). Thus the matrix element between $\phi$’s eigenstates in a trap for defining its self-energy becomes $\delta^{{}_{\bm{N}_{\phi}^{\prime}}}_{{}_{\bm{N}_{\phi}}}\delta^{m_{\ell}^{\prime}}_{m_{\ell}}\Sigma_{\omega}(E)\equiv\langle\phi_{\bm{N}_{\phi}^{\prime}}^{m_{\ell}^{\prime}}|H_{g_{\ell}}\left(E_{{}_{L}}-H_{0}\right)^{-1}H_{g_{\ell}}|\phi_{\bm{N}_{\phi}}^{m_{\ell}}\rangle$ (note the absence of $i0^{+}$ in the Green’s function), with

[TABLE]

Here $z_{{}_{E}}\equiv E/(2\omega)$ and the relative energy $E\equiv E_{{}_{L}}-E^{(\phi)}_{\bm{N}_{\phi}}$, with $E^{(\phi)}_{\bm{N}_{\phi}}=(2N_{\phi}+\ell_{\phi}+\tfrac{3}{2})\omega$ as the CM’s energy. If $\Delta_{c}$ and $\Delta_{n}$ receive trap-dependent “polarization” corrections, these corrections also need to be subtracted in defining $E$. In the derivation, a unitary transformation between $n$ and $c$ single-particle and CM/relative motion eigenmodes has been used.

Summing over the quantum numbers associated with the intermediate state’s CM motion gives rise to the $\delta^{{}_{\bm{N}_{\phi}^{\prime}}}_{{}_{\bm{N}_{\phi}}}$ factor in defining $\Sigma_{\omega}$, since the CM’s decoupling property is preserved and thus so is $\bm{N}_{\phi}$. For the relative dynamics, $\bar{R}^{(r)}_{{}_{n,\ell}}$ is part of the eigenmode function $R^{(r)}_{{}_{\bm{N}_{r}}}$ Stetcu et al. (2010): $R^{(r)}_{{}_{\bm{N}_{r}}}(\bm{r})\equiv\bar{R}^{(r)}_{{}_{n,\ell}}(r)\,r^{\ell}Y_{\ell m_{\ell}}(\hat{r})$. $\bm{N}_{r}$ has $n,\ell$ for its radial excitation and angular momentum, and $E^{(r)}_{n,\ell}=(2n+\ell+\tfrac{3}{2})\omega$. A cut-off on $n$ is used to regularize the theory in a trap, which is in parallel with the regularization used in Eq. (5).

The $\phi$ propagator in the trap, defined as $D_{\omega}(E)\delta_{{}_{\bm{N}_{\phi}}}^{{}_{\bm{N}_{\phi}^{\prime}}}\delta^{m_{\ell}^{\prime}}_{m_{\ell}}$, can be computed by summing up all self-energy insertion diagrams, including insertions of $\Sigma_{\omega}$ and those of the $d_{j}^{(\ell)}$ and $d_{j,k}^{(\ell)}$ vertices. One get

[TABLE]

In the 2nd step, the principal value of the free-space self-energy $\mathscr{P}\Sigma$ is added and subtracted. Thus, the quantization condition can be derived by setting the denominator in Eq. (8) to zero:

[TABLE]

There exists a special relation between $\Lambda$ (or $\bar{T}_{\Lambda}\equiv T_{\Lambda}/\omega$) and $n_{\Lambda}$ such that the divergences in $\Sigma$ and $\Sigma_{\omega}$ cancel in Eq. (9), and thus $d^{(\ell)}_{j,k}$ are finite. For s-wave, $\bar{T}_{\Lambda}=2n_{\Lambda}(1+O(n_{\Lambda}^{-1}))$, but for p-wave the $n_{\Lambda}^{-1}$-order term needs to be specified: ${\bar{T}_{\Lambda}}=2n_{\Lambda}\left(1+\frac{7}{4}n_{\Lambda}^{-1}+O(n_{\Lambda}^{-2})\right)$; for d-wave another higher order term needs to be specified: $\bar{T}_{\Lambda}=2n_{\Lambda}\left(1+\frac{9}{4}n_{\Lambda}^{-1}-\frac{37}{32}n_{\Lambda}^{-2}+O(n_{\Lambda}^{-3})\right)$; for even larger $\ell$, more terms need to be specified accordingly. Details on the renormalization can be found in the SM. However, the above $n_{\Lambda}$-$\Lambda$ relations should be considered as a specific scheme; any alternative ones would need to ensure that the divergences in Eq. (9)’s right side can be absorbed by the $d^{(l)}_{j,k}$ terms in the left side so that phase shifts are cut-off independent and the CM-decoupling property is not violated.

The right side of Eq. (9) in this scheme then becomes

[TABLE]

with “$(\mathscr{R})$” labeling the renormalized series sum with $n_{\Lambda}\rightarrow+\infty$. To finish the derivation, this identity is needed:

[TABLE]

which holds in the entire complex $z$ plane (both sides have the same poles and residues, see the proof in the SM). By redefining $d^{(\ell)}_{j,i}\equiv\mathscr{A}_{\ell}(2M_{{}_{\mathrm{R}}})^{j}C_{i,j}$ in Eq. (9) and applying Eq. (11) in Eq. (10), Eq. (9) gives Eq. (2).

Further comments

It is worth comparing $D(E)$ in Eq. (6) and $D_{\omega}(E)$ in Eq. (8) in the complex $E$ plane. $1/D(E)$ has a branch cut—known as the unitary cut—on the positive real axis due to the $-ip^{2\ell+1}$ term, which changes into a series of poles—called “unitary” poles below—for $1/D_{\omega}(E)$ (from the term $[\Sigma_{\omega}(E)-\mathscr{P}\Sigma(E)]/\mathscr{A}_{\ell}$). Both non-analyticities are directly connected to unitarity and thus independent of framework, power counting, and fine tuning.

However, fine tuning and power counting do impact the behavior of the ERE function van Kolck (1999): in a natural case, $C_{0,j}\sim M_{{}_{H}}^{2\ell+1-2j}$; in a $1^{\textrm{st}}$ fine-tuned case, $C_{0,0}$ is enhanced; and in a $2^{\textrm{nd}}$ fine-tuned case, the function has low-energy poles. Note Ref. van Kolck (1999) uses EFTs without a dimer field; and the Lagrangian is the same for the three cases, except power countings. One can add the couplings between $\mathscr{B}(x)$ and particles to the Lagrangian, again by multiplying the short-distance-interaction terms with powers of $\bm{\partial}^{2}\mathscr{B}$, e.g., for s-wave $\mathscr{L}_{I}=[\tilde{d}_{0}+\sum_{k\geq 1}\tilde{d}_{0,k}(M_{{}_{\mathrm{R}}}^{2}\bm{\partial}^{2}\!\mathscr{B}/3m)^{k}](c\,n)^{\ast}(c\,n)+\dots$.

For each of these cases, the $T$-Matrix in a trap can be computed in the same way as the free-space one van Kolck (1999) but with the bare couplings substituted by the corresponding modified ones—e.g., $\tilde{d}_{0}\rightarrow\tilde{d}_{0}+\sum_{k\geq 1}\tilde{d}_{0,k}(M_{{}_{\mathrm{R}}}\omega)^{2k}$—and the unitary cut by the “unitary” poles. Since the EFT calculations reproduce the free-space $T$-matrix using ERE parameters $C_{i=0,j}$, the $T$-matrix in a trap can be parameterized in the same way but with $\omega$-dependent ERE parameters. The relation between these ERE parameters and the bare couplings is nonlinear, but the former’s $\omega$-dependence could be expanded in terms of $\omega^{2}$. (This expansion must be examined with care, if its convergence radius is much smaller than the naive estimate based on $M_{H}$, e.g., due to fine-tuning of the trap’s modification to the interaction at short distance.) Finally, by identifying the poles of the trap $T$-Matrix, one then reproduces Eq. (2) for the natural and the $1^{\textrm{st}}$ fine-tuned case; for the $2^{\textrm{nd}}$ fine-tuned case, a Laurent expansion of $p^{2\ell+1}\cot\delta_{\ell}$ was derived van Kolck (1999), so the same expansion should be used in Eq. (2)’s left side with the parameters carrying $\omega^{2}$ corrections. In other words in my approach with a dimer field, resumming of $d_{j}^{(\ell)}$ and $d_{j,k}^{(\ell)}$ terms is needed.

Summary

I have applied pionless EFT to two short-range interacting particles in an external harmonic trap to derive a systematically improved BERW formula that is exact even at finite $\omega$. It is valid when the infrared scale of the trap ($\sqrt{M_{{}_{\mathrm{R}}}\omega}$) and the relative momentum ($p$) are both smaller than the high momentum scale set by the dynamics. This provides a firm foundation for implementing a Luscher-formula-like approach to connect nuclear scattering and ab initio structure calculations. The derivation involved new coupling terms between the background field and particles, which lead to the improvements of the original BERW formula. Moreover, a careful analysis of renormalization shows a non-trivial relation between the cut-off $\Lambda$ on relative momentum in free space and cut-off $n_{\Lambda}$ on the number of radial excitation in a trap. The renormalization procedure is further confirmed by Eq. (11)’s proof. Both aspects are instructive for deducing connections between a trapped system (with two or more clusters) and free-space scattering/reactions for both nuclear and cold atom physics Blume (2012). It should also be interesting to apply this framework to study exotic atoms111Here the long-range interaction is the attractive Coulomb force. The so-called Deser–Trueman formula Combescure et al. (2007) relates exotic atom’s energy levels to the short-distance interaction’s scattering length. and quantum dots Combescure et al. (2007).

Acknowledgment

I would like to thank Dick Furnstahl, Chan Gwak, Jason Holt, David Kaplan, Ubirajara Van Kolck, Petr Navratil, Daniel Phillips, Martin Savage, Ragnar Stroberg, and Chieh-Jen Yang for helpful discussions, and Yuri Kovchegov for pointing out how to use contour integration to prove Eq. (11). I also thank Dick Furnstahl and Jordan Melendez for careful proofreading of the manuscript. The work was supported by the National Science Foundation under Grant No. PHY–1614460 and the NUCLEI SciDAC Collaboration under US Department of Energy MSU subcontract RC107839-OSU, the US Department of Energy under contract DE-FG02-97ER-41014, and the US Institute for Nuclear Theory.

I An exactly solvable case: hard-sphere potential

I demonstrate here that if the interaction is short-ranged and has the form of a hard sphere, the parameters in the improved BERW formula (i.e., Eq. (2) in the main text) can be found analytically. This model was studied in Ref. Block and Holthaus (2002) by using parabolic cylinder functions for s-wave channel. The hard-sphere potential is defined as $V_{s}(r)=+\infty$ if $r\leq r_{c}$ and [math] otherwise ($r\equiv|\bm{r}|$ with $\bm{r}$ the relative displacement between the two particles). In addition, each particle experiences an external harmonic potential. Because the CM motion is factorized, one can just focus on the relative motion; the corresponding external potential is $M_{{}_{\mathrm{R}}}\omega^{2}\bm{r}^{2}/2$, with $M_{{}_{\mathrm{R}}}$ the reduced mass.

Let us define $\bar{r}\equiv r/b$ and $\bar{r}_{c}\equiv r_{c}/b$ with $b\equiv 1/\sqrt{M_{{}_{\mathrm{R}}}\omega}$, and $\bar{E}\equiv E/\omega$. When $\bar{r}>\bar{r}_{c}$, the Schrödinger equation in the $\ell^{\textrm{th}}$ partial wave becomes

[TABLE]

with the radial wave function defined as $u_{\ell}(\bar{r})/r$. Thus, the wave function at $\bar{r}>\bar{r}_{c}$ is a linear combination of two independent solutions to the harmonic oscillator Schrödinger equation Suzuki et al. (2009):

[TABLE]

where $M(a,b,z)$ is the Kummer function Abramowitz and Stegun (1972).

When $\bar{r}\rightarrow+\infty$, $u_{\ell}$ should go to zero to guarantee that the wave function is normalizable. Since at large $z$ (DLMF, , Eq. 13.7.2)

[TABLE]

the following is required:

[TABLE]

Meanwhile, $u_{\ell}\left(\bar{r}_{c}\right)=0$ implies

[TABLE]

Equations (S4)–(S5) have nontrivial solutions only if the corresponding determinant is zero, which gives the quantization condition:

[TABLE]

Here, $\mathscr{N}_{\ell}\equiv(2\ell+1)!!(2\ell-1)!!$ [$(-1)!!\equiv 1$].

Note that the left side of Eq. (S6) is the right side of the main text’s Eq. (1) multiplied by a $(-){r_{c}^{2\ell+1}}/{\mathscr{N}_{\ell}}$ factor. Meanwhile, the phase shift due to $V_{s}(r)$ in the $\ell^{\textrm{th}}$ partial wave is $\tan\delta_{\ell}={j_{\ell}(pr_{c})}/{y_{\ell}(pr_{c})}$ Joachain (1975). Thus the discrepancy of the BERW formula is identified as the difference between the left side of the main text’s Eq. (1) multiplied by $-{r_{c}^{2\ell+1}}/{\mathscr{N}_{\ell}}$ and the right side of Eq. (S6), i.e.,

[TABLE]

The left side can be expanded in terms of $\left(pr_{c}\right)^{2}$:

[TABLE]

Now according to Ref. (DLMF, , Eq.(13.2.2)), $M(a,b,z)$ can be expanded as $1+\frac{a}{b}z+\frac{a(a+1)}{b(b+1)}\frac{z^{2}}{2!}+\dots$, indicating that Eq. (S7)’s right side, a function of $E$ and $r_{c}$, can be approximated by a double expansion in terms of powers of $\bar{r}_{c}^{2}$ and $2\bar{E}\bar{r}_{c}^{2}=(pr_{c})^{2}$ when both are small (note that the coefficient denominators in the expansion of the $M$ functions in Eq. (S7) are independent of $\bar{E}$ and $\bar{r}_{c}$). This also suggests the difference between left and right sides in Eq. (S7) can be expanded in terms of $\bar{r}_{c}^{2}$ and $2\bar{E}\bar{r}_{c}^{2}$.

It can be shown that with $\omega\to 0$ and $E$ and $r_{c}$ fixed, the right side of Eq. (S7) approaches its left side. In this limit, $\bar{E}\to\infty$,

[TABLE]

where ${}_{0}F_{1}$ is a Confluent Hypergeometric Limit Function. Reference Petkovšek et al. (1996) suggests

[TABLE]

Since

[TABLE]

one then obtain

[TABLE]

and thus prove that the difference between the two sides of Eq. (S7) disappears as $\omega\rightarrow 0$. This requires that the difference, if expanded in terms of $\bar{r}_{c}^{2}$ and $2\bar{E}\bar{r}_{c}^{2}$, must have positive powers of $\bar{r}_{c}^{2}$ ($\propto\omega$).

Moreover, both sides of Eq. (S7) are unchanged when $\bar{E}\to-\bar{E}$ and $\bar{r}_{c}^{2}\to-\bar{r}_{c}^{2}$ (i.e., $\omega\to-\omega$), because $M(a,b,z)=e^{z}\,M(b-a,b,-z)$(DLMF, , Eq.(13.2.39)), indicating that the powers of $\bar{r}_{c}^{2}$ are positive and even (the powers of $2\bar{E}\bar{r}_{c}^{2}$ are non-negative integers). Therefore, the following expansion is expected, with the first two terms explicitly given:

[TABLE]

Meanwhile, the ERE for $V_{s}$ can be derived easily from Eq. (S8): $\left({pr_{c}}\right)^{2\ell+1}\cot\delta_{\ell}\left(p\right)/\mathscr{N}_{\ell}=\sum_{j=0}^{\infty}\overline{C}_{i=0,j}\left({pr_{c}}\right)^{2j}$. If a high-momentum scale $M_{{}_{H}}=r_{c}^{-1}$ is identified, Eq. (S14) can be transformed to the improved BERW formula after redefining $C_{i,j}=\mathscr{N}_{\ell}M_{{}_{H}}^{2\ell+1-4i-2j}\overline{C}_{i,j}$.

II The 2nd model for numerical testing

In contrast to the analytical study in the last section, here a simple square-well model Ali et al. (1985) is used to numerically demonstrate that the discrepancy of the BERW formula can be expanded in terms of powers of $\omega^{2}$ and $p^{2}$, which is the essence of the improved BERW formula. The potential was constructed to qualitatively describe neutron-$\alpha$ scattering in the s- and p-waves Ali et al. (1985). The quantum numbers for the considered channels in $J^{\pi}$ notation are ${\frac{1}{2}}^{+}$, ${\frac{1}{2}}^{-}$ and ${\frac{3}{2}}^{-}$. The potential $V_{s}(r)=V_{0}(1+\beta\bm{L}\cdot\bm{\sigma})$ when $r<r_{c}$ and [math] when $r>r_{c}$, with $V_{0}=-33\,$MeV, $r_{c}=2.55\,$fm, $\beta=0.103$ Ali et al. (1985). $\bm{L}\cdot\bm{\sigma}$ is the spin-orbit coupling, which generates differences between the $p_{\frac{3}{2}}$ and $p_{\frac{1}{2}}$ channels.

In this exercise, the potential is treated as the exact underlying physics for the two particles. The phase shifts, calculated by solving the corresponding Schrödinger equations in the continuum, are shown in Fig. S1. These phase shifts are considered to be “exact” ones. (The ${\frac{3}{2}}^{-}$ channel’s exact phase shift is also shown in Fig. 1 in the main text.) Meanwhile, the spectra for the two particles in various harmonic potential traps can also be precisely computed. The goal is to test the original and the improved BERW formulas using these exact phase shifts and energy spectra.

In Fig. S2, the discrete symbols are the differences between the two sides in the main text’s Eq. (1) divided by $M_{H}^{2\ell+1}$ at eigenenergies associated with the trap frequency $\omega=0.1$, $\sqrt{0.1}$, $1$, $\sqrt{10}$ and $10$ MeV. Here, I take $M_{H}=200$ MeV ($\sim 1\;\mathrm{fm}^{-1}$), as motivated by the value of $r_{c}$. One can see that the leading-order difference does scale as $\omega^{2}$. Also, in the s-wave channel there is a deep bound state in free space—unphysical for the n-$\alpha$ system—and thus a distinct negative eigenenergy in all the traps. However, for the other two channels no bound state exists in free space.

To see how well the improved BERW formula works, one need to know the values of $C_{i,j}$ corresponding to this particular potential. First, another set of eigenenergies at extremely small trap frequencies (both on the order of $10^{-6}$ MeV) is computed. Second, they are used as inputs to fit $C_{i,j}$ values based on the improved BERW formula. A least-squares fit uses as the objective function the sum of the squares of the differences between the two sides in the improved BERW formula, as calculated at those small eigenenergies. The small values of $\omega$ and eigenenergies help separate the impact of $C_{i,j}$ at different orders, and enables a precise fit (it amounts to computing derivatives at $E=0$ and $\omega=0$ with points separated by tiny distances). The series in the improved BERW formula is nonetheless truncated to $i\leq 2$ and $j\leq 3$. To make the parameters dimensionless, one can rescale them by appropriate powers of $M_{H}$, $C_{i,j}\equiv\overline{C}_{i,j}M_{H}^{2\ell+1-4i-2j}$. The best-fit values for $\overline{C}_{i,j}$ are shown in Tables 1, 2, and 3.

A few higher-order $\bar{C}_{i,j=2}$ values are not shown in those tables, because over-fitting Furnstahl et al. (2015) in this simple exercise leads to values significantly larger than $1$. However, these contributions are very small in the plots shown if they on the order of $1$ (their natural size), and therefore they are set to zero in generating the plots.

Having determined the $C_{i,j}$, the $C_{i\neq 0,j}$ terms in the improved BERW formula are used in Fig. S2 to interpolate between the discrete symbols, i.e., the difference between the right and left sides of the BERW formula as computed at discrete eigenenergies up to $\approx 20$ MeV. Because the fitting of $C_{i,j}$ is carried for $E$ and $\omega$ near zero ($\sim 10^{-6}$ MeV), the agreement with the interpolating curves and discrete symbols demonstrates that the discrepancy of the BERW formula can indeed be expanded in terms of powers of $\omega^{2}$ and $p^{2}$, as implied by the improved BERW formula.

However, the truncation on the series sum in the improved BERW formula inevitably leads to truncation error. Fig. S3 plots such errors, i.e., the difference between the formula’s two sides with the $j$ index summed up to either $3$ (for the “Up to $E^{3}$” results) or $2$ (for the “Up to $E^{2}$” results) and with $i$ index summed up to $2$ for both cases. The shaded region shows the rough range of the formula’s RHS with $E$ above 5 MeV. In the low energy region, the error does behave like the leading terms left out of the series sum, but it increases with $E$ and eventually reaches the shaded region. The latter signals the break down of the series expansion, as the error reaches $100\%$. This also suggests the break down scale for $E$ is above $20$ MeV and perhaps below 40 MeV in all three channels, which are consistent with my choice of $M_{H}=200$ MeV. Note the behavior of the error above the break down scale could be more complicated than the simple power law in the low-energy region. In fact in all these channels, the errors oscillate and change their signs in the higher-energy range not shown in the plots.

III Details on the renormalization in the EFT derivation

This section discusses the particular renormalization scheme mentioned in the main text, in which the divergences in in $\Sigma$ and $\Sigma_{\omega}$ cancel in Eq. (9) in the main text. This scheme gives a specific relation between $\Lambda$ and $n_{\Lambda}$. To derive this relation, the $n_{\Lambda}$-dependence of $\Sigma_{\omega}$ needs to be studied.

Two formulas are useful for understanding $f_{\ell}(z_{E},n)$ at large $n$ and $\Sigma_{\omega}$ at large $n_{\Lambda}$. The first is (DLMF, , Eq. 5.11.13)

[TABLE]

Here $a$ and $b$ are real or complex constants, and $G_{k}(a,b)$ as a function of $a$ and $b$ is related to the generalized Bernoulli polynomials (DLMF, , Eq. 5.11.17) Temme (1996). The second is the Euler–-Maclaurin formula (DLMF, , Eq. 2.10.1), stating that for a smooth $f(x)$, its series sum can be approximated using an asymptotic expansion:

[TABLE]

where $B_{n}$ is a Bernoulli number. Only the $n_{\Lambda}$ dependent terms are shown. Like $\Sigma(E)$, the $0^{\textrm{th}}$ to $\ell^{\textrm{th}}$ derivatives of $\Sigma_{\omega}(E)$ diverge. So for s-waves, only $\Sigma_{\omega}(0)$ is considered:

[TABLE]

Thus $\bar{T}_{\Lambda}=2n_{\Lambda}(1+O(n_{\Lambda}^{-1}))$ ($\bar{T}_{\Lambda}\equiv T_{\Lambda}/\omega$) so that the divergence can be absorbed by $\sqrt{2\bar{T}_{\Lambda}}$ in the main text’s Eq. (10). The $O(n_{\Lambda}^{-1})$ term in the $n_{\Lambda}$-$\Lambda$ relation is not relevant when $n_{\Lambda}\to\infty$. For p-waves, the $0^{\textrm{t}h}$ and first derivatives are divergent:

[TABLE]

Thus, ${\bar{T}_{\Lambda}}=2n_{\Lambda}\left(1+\frac{7}{4}n_{\Lambda}^{-1}+O(n_{\Lambda}^{-2})\right)$, to have these divergences canceled by those in $\mathscr{P}\Sigma(E)$ in the main text’s Eq. (10). The $\frac{7}{4}n_{\Lambda}^{-1}$ piece in the expression must be specified, but higher-order terms are not needed. However for d-waves, another order higher needs to be specified: $\bar{T}_{\Lambda}=2n_{\Lambda}\left(1+\frac{9}{4}n_{\Lambda}^{-1}-\frac{37}{32}n_{\Lambda}^{-2}+O(n_{\Lambda}^{-3})\right)$; for even larger $\ell$, more terms need to be specified accordingly.

This requirement is tied to the fact that for a specific $E$-derivative, there is a tower of divergences with different degrees in $\Sigma_{\omega}$ (e.g., Eq. (S18)), in stark contrast with the same derivative of $\Sigma$, where the power of $\Lambda$ is fixed by the the dimensions.

IV A proof of Eq. (11) in the main text

Let us redefine $\tilde{z}\equiv z-\ell/2-3/4$, and

[TABLE]

The main text’s Eq. (11) then becomes

[TABLE]

The renormalization $\mathscr{R}$ is defined by Eq. (10) in the main text and the relationship between $T_{\Lambda}$ and $n_{\Lambda}$ that guarantees the cancellation of the divergences on that equation’s left side. The proof starts with integrating $g_{\ell}(\tilde{u})/(\tilde{u}-\tilde{z})$ in the complex $\tilde{u}$ plane over a large contour around the origin and crossing between the two singularities at $\tilde{u}=n_{\Lambda}$ and $n_{\Lambda}+1$ on the positive real axis. The contour $\overline{C_{\epsilon}}+C_{\epsilon}$ is plotted in Fig. S4. After rearranging terms one get

[TABLE]

The left side comes from the residue of the $1/(\tilde{u}-\tilde{z})$ pole in the contour integration while the series sum on the right side is from the residues of $g_{\ell}(\tilde{u})$’s poles (all on the positive real axis) in the same integration. Comparing Eq. (S23) to Eq. (S22) suggests that the contour integration in Eq. (S23) should cancel the series’s divergence as the $\bar{T}_{\Lambda}$ terms—i.e., those from $\mathscr{P}\Sigma(E)$—do in the main text’s Eq. (10). This is the focus of the following proof.

As preparation, it is important to understand the behavior of $g_{\ell}(\tilde{u})$ in two different regions (define $\theta_{\tilde{u}}\equiv\arg(\tilde{u})$ and $R$ as the radius of $\overline{C_{\epsilon}}$): $\pi\geq|\theta_{\tilde{u}}|\gg R^{-1}$ and $|\theta_{\tilde{u}}|\lesssim R^{-1}$. Let us look at the region close to the positive real axis first. Considering the presence of $g_{\ell}(\tilde{u})$’s poles, it is easier to use the following re-expression based on (DLMF, , Eq. 5.5.3):

[TABLE]

This expression moves the poles from the $\Gamma$ function to the $\cot$ function. When $R\rightarrow+\infty$,

[TABLE]

Here $\theta_{0}$ is a small positive number. Therefore, $\cot(\pi\tilde{u})\sim-i\operatorname{sgn}(\text{Im\,}\tilde{u})$ when $\pi-\theta_{0}>|\theta_{\tilde{u}}|\gg R^{-1}$ [the error scales as $\exp(-2\pi|\text{Im\,}(\tilde{u})|)$], but in the $2^{\textrm{nd}}$ region, $|\theta_{\tilde{u}}|\lesssim R^{-1}$, and its behavior is qualitatively different.

Since the $\Gamma$ function ratio in Eq. (S24) is the same as that ratio in $\tilde{f}_{\ell}(\tilde{z},n)$ with $n\rightarrow\tilde{u}$, by applying the asymptotic expansion from the main text’s Eq. (11) on the right side of Eq. (S24), I get, when $R\rightarrow+\infty$ such that $R\gg|\tilde{z}|$ and when $\pi-\theta_{0}>|\theta_{\tilde{u}}|\gg R^{-1}$,

[TABLE]

Here $\tilde{f}_{\ell}^{\mathrm{E}}(\tilde{z},\tilde{u})$ is the asymptotic expansion of $\tilde{f}_{\ell}(\tilde{z},\tilde{u})$ in terms of $1/\tilde{u}$ using the main text’s Eq. (11). It turns out the above expansion also holds when $\theta_{u}\rightarrow\pi$ or $-\pi$. However, using Eq. (S24) to understand $g_{\tilde{u}}$ close to the negative real axis is awkward, as it involves the cancellation of poles from the $\cot$ and the $\Gamma$ functions. Instead, Eq. (11) in the main text can be applied to analyze the original form of $g_{\ell}(\tilde{u})$ (see Eq. (S21)):

[TABLE]

Here $G_{k}(\ell+\frac{3}{2},1)=(-)^{k}G_{k}(0,-\ell-\frac{1}{2})$ is used, which can be inferred from its definition (DLMF, , Eq. 5.11.17) Temme (1996) and properties of the generalized Bernoulli polynomials. It is worth pointing out that the $\operatorname{sgn}(\text{Im\,}\tilde{u})$ factor moves the asymptotic series’s branch cut on the negative real axis due to the $\sqrt{\tilde{u}}$222As usual (DLMF, , Eq. 1.9.7), the $\sqrt{\tilde{u}}$ branch cut in the complex $\tilde{u}$ plane is on the negative real axis. Meanwhile $\Sigma(E)$ (and the related scattering $T$-matrix), as calculated in the main text’s Eq. (5), in the complex $E$ plane has a branch cut separating physical and unphysical sheets on the positive real axis Goldberger and Watson (1964), because it involves $-ip=\sqrt{-(2M_{{}_{\mathrm{R}}}E+i0^{+})}$. factor to the positive real axis, ensuring that the expansion series is analytic around the negative real axis. Therefore, Eq. (S26) holds when $\pi\geq|\theta_{\tilde{u}}|\gg R^{-1}$, but not in the $|\theta_{\tilde{u}}|\lesssim R^{-1}$ region. This suggests splitting the contour integration into two major pieces:

[TABLE]

In $\overline{C_{\epsilon}}+C_{\epsilon}-S$, $S$ is the point where $C_{\epsilon}$ crosses the real axis. Thus the contour means to integrate infinitely close to the upper and lower real axis, because of the corresponding integrand’s branch cut on the positive real axis. In the $2^{\textrm{nd}}$ piece, the first term should be integrated over $\overline{C_{\epsilon}}+C_{\epsilon}$. However, since the integrand is continuous across the point $S$ as long as $\epsilon\neq 0$, changing the contour to $\overline{C_{\epsilon}}+C_{\epsilon}-S$ does not affect the results.

Integrating the $1^{\textrm{st}}$ piece over $\overline{C_{\epsilon}}$ with $\epsilon\rightarrow 0^{+}$ gives

[TABLE]

Here $\tilde{f}_{\ell}^{\text{ET}}(\tilde{z},\tilde{u})$ truncates the expansion in $\tilde{f}_{\ell}^{\text{E}}(\tilde{z},\tilde{u})$ by only keeping terms with powers of $\tilde{u}$ from $\ell-1/2$ to $-1/2$ (the neglected terms’ integration vanishes no slower than $O(1/\sqrt{R})$ as $R\rightarrow\infty$). In the derivation, (1) the $-\pi\leq\theta_{\tilde{u}}<0$ branch has been rotated by $+2\pi$, which eliminated the $\operatorname{sgn}(\text{Im\,}\tilde{u})$ factor because of the $\sqrt{\tilde{u}}$ factor in the integrand; (2) a transformation, $\tilde{u}=v^{2}$, was used; and (3) the fact that $v\tilde{f}_{\ell}^{\text{E,T}}(\tilde{z}v^{2})$ is analytic in the upper complex $v$ plane and an even function on the real axis was used. Note that integrating the $1^{\textrm{st}}$ piece over $C_{\epsilon}-S$ gives [math] with $\epsilon\rightarrow 0^{+}$.

For the $2^{\textrm{nd}}$ piece, the $\overline{C_{\epsilon}}+C_{\epsilon}-S$ contour can be deformed to $C_{L}+C_{R-}+C_{\epsilon}+C_{R+}-S$ without crossing any singularities, with the left and right segments connecting at $\infty$. The integrand on the two sets of contours (including in the area enclosed by them) is 0 up to at most a $\sim e^{-R\#}$ correction (with $\#$ as a positive number), except on the segments with $|\theta_{\tilde{u}}|\lesssim 1/R$ along $C_{R\pm}$ and $C_{\epsilon}$. Therefore, $C_{L}$ can be safely ignored. Let us focus on $C_{R-}+C_{\epsilon}+C_{R+}-S$. Since

[TABLE]

the $2^{\textrm{nd}}$ piece is

[TABLE]

Note that $\cot(\pi\tilde{u})\sim 1/[\pi(\tilde{u}-n_{\Lambda})]$ when $\tilde{u}\rightarrow n_{\Lambda}$, its integration over $C_{\epsilon}-S$ with $\epsilon\rightarrow 0^{+}$ gives $(-)\tilde{f}_{\ell}^{\text{E}}(\tilde{z},n_{\Lambda})/2$.

Integrating over $C_{R\pm}$ with $\epsilon\rightarrow 0^{+}$ gives

[TABLE]

Because the integration is dominated by the $\Delta\ll n_{\Lambda}$ region, the two $\tilde{f}_{\ell}^{\mathrm{E}}$ are expanded in terms of Taylor series for the $2^{\textrm{nd}}$ argument at $n_{\Lambda}$ in the above derivation. The resulted integrations are proportional to Riemann $\zeta(2j)$ at even arguments, which are related to Bernoulli numbers (DLMF, , Eq. 25.5.1, 25.6.2, 24.2.2). Adding all the contributions, one can see that the $2^{\textrm{nd}}$ term on the right side of Eq. (S23) can be approximated by

[TABLE]

with the leading error scaling as $n_{\Lambda}^{-\frac{1}{2}}$ [see the discussion of $\tilde{f}_{\ell}^{\text{ET}}$. In comparison, the $e^{-n_{\Lambda}\#}$ error ($\#$ is a positive constant) is much smaller with large $n_{\Lambda}$]. The above expression is exactly the same as the divergent $n_{\Lambda}$-dependent pieces in $\Sigma_{\omega}$’s series sum derived using Eq. (S16) in the previous section, which completes the proof.