Universal physics of bound states of a few charged particles

TL;DR

This paper explores the universal properties of shallow bound states of a few charged particles, revealing unique scaling laws influenced by Coulomb interactions, and emphasizing the importance of Coulomb-modified scattering parameters.

Contribution

It introduces universal scaling laws for charged particle bound states that depend on Coulomb-modified scattering length and effective range, extending understanding beyond neutral systems.

Findings

Bound states obey unique Coulomb-influenced scaling laws.

Both Coulomb-modified scattering length and effective range are essential for accurate descriptions.

States with large spatial extent exhibit universality independent of short-range potential shape.

Abstract

We study few-body bound states of charged particles subject to attractive zero-range/short-range plus repulsive Coulomb interparticle forces. The characteristic length scales of the system at zero energy are set by the Coulomb length scale $D$ and the Coulomb-modified effective range $r_{\mathrm{eff}}$. We study shallow bound states of charged particles with $D\gg r_{\mathrm{eff}}$ and show that these systems obey universal scaling laws different from neutral particles. An accurate description of these states requires both the Coulomb-modified scattering length and the effective range unless the Coulomb interaction is very weak ($D\to \infty$). Our findings are relevant for bound states whose spatial extent is significantly larger than the range of the attractive potential. These states enjoy universality -- their character is independent of the shape of the short-range potential.

Figures

Equations

$$-\frac{{\hslash }^2}{2\mu }\frac{{\partial }^{2}u}{\partial r^2}+\left[V_C(r)+V_S(r)\right]u(r)=-\frac{{\hslash }^2{\kappa }^2}{2\mu }u(r),$$

$$u(r)=\mathcal{N}\times \{\begin{array}{ll}\sin \left(\frac{r}{R}\sqrt{g-{\kappa }^2{R}^2}\right)& \text{if }r\le R\\ W_{-\eta ,1/2}(2\kappa r)\frac{\sin \left(\sqrt{g-{\kappa }^2{R}^2}\right)}{W_{-\eta ,1/2}(2\kappa R)}& \text{if }r>R,\end{array}$$

$$\sqrt{g}\cot (\sqrt{g})=(2\eta \psi (\eta )+2\eta \ln (2\kappa R)+4\eta \gamma +1)\kappa R,$$

$$a_C=D{\left[2\ln \left(\frac{2R}{D}{e}^{2\gamma }\right)-\sqrt{g}\cot (\sqrt{g})\frac{D}{R}\right]}^{-1},$$

$$D=\frac{1}{\kappa \eta }=\frac{{\hslash }^2}{k_e\mu Q_1{Q}_2}$$

$$2\eta \psi (\eta )-2\eta \ln (\eta )+1=-\frac{1}{\kappa a_C},$$

$$\kappa \simeq \frac{1}{a_C}\left[1-\frac{2a_C}{D}\ln \left(\frac{a_C{e}^{\gamma }}{D}\right)\right].$$

$${\kappa }^2\simeq \frac{6}{D{a}_C}+\frac{18}{5a_C^2},$$

$$\frac{\sqrt{⟨r^2{⟩}_0}}{D}=\eta \sqrt{\frac{{\int }_0^{\infty }{W}_{-\eta ,1/2}^2(2x)x^2\mathrm{d}x}{{\int }_0^{\infty }{W}_{-\eta ,1/2}^2(2x)\mathrm{d}x}},$$

$$2\eta \psi (\eta )-2\eta \ln (\eta )+1=-\frac{1}{\kappa }\left(\frac{1}{a_C}+\frac{r_{\mathrm{eff}}{\kappa }^2}{2}\right).$$

$${\kappa }^2\simeq \frac{6}{a_{C}D(1-3r_{\mathrm{eff}}/D)},$$

$$V_S^G=g_G{e}^{-r^2/(2R_G^2)},$$

$$\frac{⟨O⟩}{⟨O{⟩}_0}\simeq 1+\frac{r_{\mathrm{eff}}}{2D}\frac{W_{-\eta ,1/2}^2(0)}{{\int }_0^{\infty }{W}_{-\eta ,1/2}^2(2x)\mathrm{d}x},$$

$$\cot ({\delta }_0)=-\frac{1}{2\pi P_g}\left[\frac{D}{a_C}+\frac{k^2{D}^2}{6}\left(1-3\frac{r_{\mathrm{eff}}}{D}\right)\right]+....$$

$$\left[-\frac{{\hslash }^2}{2m_3}\frac{{\partial }^2}{\partial {{\mathbf{r}}_3}^2}+V_S^{lh}({\mathbf{r}}_1-{\mathbf{r}}_3)+V_S^{lh}({\mathbf{r}}_2-{\mathbf{r}}_3)\right]{\Psi }^{lh}=ϵ{\Psi }^{lh},$$

$$\left[-\frac{{\hslash }^2}{m_1}\frac{{\partial }^2}{\partial {\text{mathbfcal}R}^2}+ϵ(\text{mathbfcal}R)+V_C(\text{mathbfcal}R)+V_S^{hh}(\text{mathbfcal}R)\right]{\Phi }^{hh}=E{\Phi }^{hh},$$

$$\left(-\frac{{\hslash }^2}{m_1}\frac{{\partial }^2}{\partial \text{mathbfcal}{R}^2}-\frac{{\hslash }^{2}A}{m_1\text{mathbfcal}{R}^2}+k\frac{Q_1^2}{|\text{mathbfcal}R|}\right)\Phi =E\Phi ,$$

$$N\le \frac{2(R-b)}{D}+A\ln \left(\frac{AD}{2R}\right),$$

$$V_{3b}=g_{3b}{e}^{-(r_{12}^2+r_{23}^2+r_{13}^2)/(16R_G^2)}$$

$${\Psi }_{GEM}=\sum _{i=1}^M{a}_i\mathcal{S}{e}^{-{\alpha }_i{x}^2-{\beta }_i{y}^2},$$

$${\Psi }_{SVM}=\sum _{i=1}^{M_S}{a}_i^S{e}^{-{\alpha }_i^S{x}^2-{\beta }_i^S{y}^2-{\gamma }_i^S\mathbf{x}\cdot \mathbf{y}},$$

$${\Phi }_{m_1,m_2}=\sum _{i,l_1,l_2}{f}_{l_1,l_2}^i(x,y)Y_{l_1{m}_1}({\theta }_y,{\varphi }_y)Y_{l_2{m}_2}({\theta }_x,{\varphi }_x),$$

$$H_3{P}_l(\mathbf{x}\cdot \mathbf{y})=\sum _{l^{\prime }}{F}_{l^{\prime }}(x,y)P_{l^{\prime }}(\mathbf{x}\cdot \mathbf{y}),$$

$$\varphi =\sum f_l^i(x,y)P_l(\mathbf{x}\cdot \mathbf{y}).$$

$$f(r_{\text{eff}})=a_1{r}_{\text{eff}}+c_1{r}_{\text{eff}}^2+b_1$$

$$\tilde{f}(r_{\text{eff}})=a_2{r}_{\text{eff}}+c_2{r}_{\text{eff}}^2+d_2{r}_{\text{eff}}^3+b_2,$$

Full text

††journal: Physics Letters B

1 Introduction

Shallow bound states of two neutral particles with zero angular momentum live in a classically forbidden region and retain almost no information about binding interactions [1]. As a consequence, any short-range attractive potential, $V_{S}$, can be used to model these states as long as it fixes a few relevant parameters (e.g., the scattering length, effective range) to their physical values. A celebrated $V_{S}$ is a zero-range potential tuned to reproduce the scattering length [2, 3, 4]. It provides a powerful starting point for studying universal bound states (i.e., independent of the shape of $V_{S}$) in nuclear and atomic physics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12].

In this Letter we consider particles that interact via $V_{S}+V_{C}$, where $V_{C}$ is a repulsive Coulomb potential. Potentials $V_{S}+V_{C}$ are typical for cluster models of nuclei [5, 13, 14], e.g., in 17F between 16O and a proton [15, 16]. Furthermore, they provide an effective description of interactions between charged quasi-particles, e.g., between dressed electrons in crystals [17, 18, 19, 20, 21, 22]. We focus on $V_{S}$ of zero range, and explore it as a possible starting point for understanding realistic charged systems.

The main finding of our study is that there is a new family of universal few-body bound states for charged systems interacting via a potential $V_{C}+V_{S}$. Their properties are fully determined by the Sommerfeld parameter, the Coulomb-modified scattering length and effective range, as well as the three-body parameter. Note that, in contrast to neutral particles, the zero-range approximation to $V_{S}$ is not guaranteed to be useful for realistic shallow bound states: The Coulomb barrier makes the spatial extent of the wave function finite [23], forcing particles to explore the landscape of the short-range binding potential. For shallow two-body bound states, we show that finite-range corrections to the energy can be accounted for by an effective range parameter. For weakly-bound three- and four-body systems, we study these corrections numerically using the Gaussian Expansion Method [24] and Stochastic Variational Method with Gaussians [25, 26, 27, 28]. Details of the numerical methods are given in A.

2 Two-Body System

We consider two particles whose relative motion is described by the radial Schrödinger equation

[TABLE]

where $\mu$ is the reduced mass, $E_{2}=-\frac{\hbar^{2}\kappa^{2}}{2\mu}$ with $\kappa>0$ is the two-body energy, $V_{C}=k_{e}\frac{Q_{1}Q_{2}}{r}$ is the Coulomb potential energy ($Q_{1}$, $Q_{2}$ are the particle charges, $k_{e}$ is Coulomb’s constant), and $V_{S}$ is a binding potential of range $R$. We consider only zero angular momenta since our focus is on the bound states for $R\to 0$ (later referred to as the zero-range or universal limit). Moreover, we are free to choose any shape of $V_{S}$, which is irrelevant as long as the limit $R\to 0$ is well-defined for neutral particles interacting via $V_{S}$ (cf. [3, 29]). For simplicity, we assume that $V_{C}+V_{S}$ is a square well for $r\leq R$, i.e., $V_{S}(r)=-\frac{\hbar^{2}g}{2\mu R^{2}}-V_{C}(r)$, and $V_{C}$ otherwise. The dimensionless parameter $g>0$ sets the interaction strength. The wave function $u$ for this potential reads

[TABLE]

where $\cal{N}$ is a normalization constant that ensures $\int_{0}^{\infty}u^{2}\mathrm{d}r=1$, $\eta\equiv\frac{k_{e}\mu Q_{1}Q_{2}}{\hbar^{2}\kappa}$ is the Sommerfeld parameter, and $W$ is the Whittaker $W$-function [30]. The values of $\kappa$ that lead to a continuous derivative of $u$ at $r=R$ define allowed bound states. In the limit $R\to 0$, $\kappa$ is a root of the equation

[TABLE]

where $\gamma$ is Euler’s constant, and $\psi$ is the digamma function [30]. Note that for neutral particles ($\eta=0$) Eq. (2) depends only on $\kappa R$, hence, the result of taking the limit $R\to 0$ with fixed $\kappa$ is identical to that with $\kappa\to 0$ and fixed $R$. In other words, for weakly-bound states of neutral particles one may always rely on the zero-range limit. For charged particles with $\kappa\to 0$, by contrast, the zero-range limit is not necessarily accurate.

We rewrite Eq. (3) by following Ref. [29]. Collecting the energy-independent terms, we require their sum to approach the Coulomb modified scattering length,

[TABLE]

where

[TABLE]

is the length scale associated with the Coulomb potential (i.e., the generalized Bohr radius). This leads to a model-independent relation for $R\to 0$,

[TABLE]

which connects one observable (the binding energy) to another (the scattering length). Previously, Eq. (6) was derived for proton-proton states [31, 32, 33, 34]. It is related to the poles of the scattering amplitude defined as in Refs. [35, 36].

In this Letter, we consider Eq. (6) in the context of universal bound states which require $D\gg R$. We end up with the three different regions in Fig. 1 for the location of the scattering length $a_{C}$ with respect to the range, $R$, of $V_{S}$ and the Coulomb length, $D$. We consider the regions I ($a_{C}\gg D$) and II ($R\ll a_{C}\ll D$) where universal physics can be expected. (In region III this is not likely because the system probes the shape of $V_{S}$.) The left-hand-side of Eq. (6) is a monotonic negative function. Therefore, if $a_{C}>0$ a zero-range potential can support at most one bound state111 This statement is correct only for zero-range potentials. A finite range potential can support any finite number of bound states. Depending on the prescrition for taking the limit $R\to 0$, those states either vanish or become infinitely deep., and if $a_{C}<0$ there can be no bound states. In region II (“weak Coulomb”), we obtain

[TABLE]

This equation features a logarithmic correction to the standard expression $\kappa=1/a_{C}$ for neutral particles [1]. Logarithmic dependence on $D$ is typical for “weak Coulomb’ [32, 37]. In region I (“strong Coulomb”) we have

[TABLE]

which describes shallow bound states. The fact that $\kappa^{2}\sim 1/a_{C}$ in the limit $a_{C}\to\infty$ will be of utmost importance for finite-range corrections. We now focus on this new class of shallow states.

To calculate other observables, we note that for $R\to 0$ the particles move almost exclusively in the classically forbidden region. Indeed, the probability to find particles with $r>R$ is approaching unity: $P(r>R)=1-\int_{0}^{R}u^{2}\mathrm{d}r\xrightarrow[]{R\to 0}1\,$. To derive this limiting value, we notice that ${\int_{0}^{R}u^{2}\mathrm{d}r}<{u^{2}(R)R}$ for $R\to 0$. Therefore, observables for zero-range interactions are described with the wave function $W_{-\eta,1/2}(2\kappa r)$, defined by $\eta$ and $\kappa$. As an example, we use the root-mean-square (rms) radius, $\langle r^{2}\rangle\equiv\int u^{2}r^{2}\mathrm{d}r$, – a standard observable in few-body physics – given by

[TABLE]

where the subscript [math] refers to the zero-range limit. The right-hand-side of Eq (9) is a monotonically increasing function of $\eta$. The maximum value is attained at $1/\eta=0$ where $\sqrt{\langle r^{2}\rangle_{0}}/D=0.507\dots$; in this limit the size of the bound state is fully determined by $D$. The boundedness of the rms radius is relevant for charged halo nuclei [23]; it also supports the predicted discontinuous behavior of the mean distance between two polarons across the unbound-polarons to bipolaron transition [18, 38]. For $\eta\to 0$ the rms radius is determined from the standard relation $2\langle r^{2}\rangle_{0}\kappa^{2}\to 1$ [5].

3 Finite-range corrections

The rms radius in Eq. (9) does not diverge for $\kappa\to 0$, indicating that the inclusion of finite-range corrections is unavoidable for charged systems. These corrections must be small if $\sqrt{\langle r^{2}\rangle}\gg R$, which, according to Eq. (9), is satisfied for weakly-bound systems if $D\gg R$. For comparison, $D\simeq 57.6$ fm for two free protons and $D\simeq 0.1$ nm for free electrons. Therefore, if two protons (or a proton with a light nucleus) formed a shallow bound state it would be universal,222 One might speculate that proton-proton correlations in a nuclear medium (e.g., in the outer core of a neutron star) can potentially lead to states relevant for our results. Since these correlations are not fully understood [39], we omit this discussion here. since natural values of $R$ in this case are around $1$ fm. Two dressed electrons in solids could represent another universal system where the effective mass, the strength of the Coulomb potential (hence $D$), as well as $R$ depend on the material. In contrast, a shallow bound state of two free atomic ions ($D\simeq 57.6$ fm and $R\sim 0.1$ nm) cannot be universal.

To estimate finite-range corrections to the energy, we notice that the right-hand-side of Eq. (6) is the first term of the Coulomb-modified effective-range expansion. To account for the next term, one must use $-1/(a_{C}\kappa)-r_{\mathrm{eff}}\kappa/2$ instead of $-1/(a_{C}\kappa)$, where $r_{\mathrm{eff}}$ is the effective range [40, 41]

[TABLE]

The effective-range correction enters at leading order in the expansion of the energy for weakly bound states

[TABLE]

because the leading contribution to $\kappa^{2}$ in Eq. (8) is proportional to $1/a_{C}$. We note that Eq. (11) was previously derived in the context of connecting asymptotic normalization constants of charged bound states to scattering parameters [42]. Higher order finite-range corrections (e.g., due to the shape parameter) are not as important in the limit $1/a_{C}\to 0$ because they are convoluted with $\kappa^{n}$ where $n>2$. Note that the factor $(1-3r_{\mathrm{eff}}/D)$ in the denominator of $\kappa^{2}$ for $1/a_{C}\to 0$ implies that for potentials with weakly-bound states $r_{{\mathrm{eff}}}<D/3$ must hold, in agreement with the causality constraints of [43].

Our result can be used to define the leading order of an effective field theory for shallow bound states of charged particles where $r_{\mathrm{eff}}$ contributes at leading order, while higher effective range parameters can be included perturbatively. The fact that range corrections are enhanced in systems with strong Coulomb interactions was already observed for 17F [16] and 7Be [44], and attributed to an additional fine tuning. The importance of finite-range effects in systems with strong Coulomb interactions was also observed in [45, 46]. Here, we show that this enhancement is generic for “strong Coulomb”. Effective field theories designed for bound systems in region II that use only $a_{C}$ as two-body input in the leading order (see, e.g., [47, 48]) must be extended to describe shallow-bound nuclei close to the proton dripline where $r_{\mathrm{eff}}/D$ is not small.

To illustrate finite-range effects, we use the Gaussian potential

[TABLE]

where $R_{G}$ defines the range of the potential, and $g_{G}$ is used to fix the Coulomb-modified scattering length for a fixed value of $R_{G}$. For the sake of discussion, we use parameters of two $\alpha$-particles $\hbar^{2}/\mu=20.73$ MeV$\times$fm2 and $kQ_{1}Q_{2}=5.76$ MeV$\times$fm [$D\simeq 3.6$ fm]. We employ the Gaussian Expansion Method [24] to calculate $\kappa$. The result is plotted in Fig. 2a) as a function of $[a_{C}(1-3r_{\mathrm{eff}}/D)]^{-1}$ (for consistency, we take $r_{\mathrm{eff}}$ for $V_{S}^{G}$ with $\kappa=0$), which is the only relevant parameter for $1/a_{C}\to\infty$; see Eq. (10). Figure 2 shows that even though the universal prediction does not describe finite values of $r_{\mathrm{eff}}/D$ (see Fig. 2b)), it is still useful: Finite-range corrections for shallow bound states are captured by rescaling $a_{C}$ with $1-3r_{\mathrm{eff}}/D$, see Fig. 2a).

Expectation values of other observables $(\langle O\rangle\equiv\int u^{2}O\mathrm{d}r)$ also acquire finite-range corrections when $r_{\mathrm{eff}}\neq 0$. It is particularly easy to calculate these corrections for an observable $O$, which in the limit $R\to 0$ satisfies $\int_{0}^{R}Ou^{2}\mathrm{d}r\sim R^{1+\delta}$, where $\delta>0$, e.g., the rms radius. After straightforward but tedious calculations we derive (in the leading order in $r_{\mathrm{eff}}/D$)

[TABLE]

where $\langle O\rangle_{0}$ is the universal prediction for the same $\kappa$ and $\eta$. As anticipated, the universal value $\langle O\rangle_{0}$ is accurate if $R\ll D$ (note that $r_{\rm eff}\sim R$ for $D/R\gg 1$). The left-hand-side depends weakly on the energy ($\eta$) and for $\eta>2$ it can be accurately written as $1+\frac{r_{\mathrm{eff}}}{2D}(6+\frac{1.1}{\eta^{2}})$. Note also that the correction in Eq. (13) is independent of $O$. Therefore, it can, in principle, be used to relate different measurements.

Finally, we discuss what happens to a shallow bound state if there is a small perturbation to $V_{S}$ that changes the sign of the scattering length, i.e., $1/a_{C}=0^{-}$ in a new potential. Such a perturbation modifies the nature of a discrete state. The bound state turns into a resonance whose width is determined by the probability for a particle to overcome the Coulomb barrier. This probability is given by the Gamow-Sommerfeld factor, $P_{g}=\exp[-2\pi/(Dk)]$, where $k>0$ is the scattering wave vector, ($\kappa=ik$ in Eq. (1)). To make the discussion more quantitative, we expand the Coulomb-modified phase shift [29] for $k\to 0$

[TABLE]

If $a_{C}<0$ and large, then low-energy scattering shows a resonance feature that extends Eq. (11) to $a_{C}<0$. If $a_{C}>0$ or $a_{C}<0$ and small, then $\delta_{0}\sim P_{g}$, which is exponentially small for $k\to 0$. In this case, the scattering states are decoupled from the short-range potential; all scattering observables are accurately determined by $V_{C}$.

4 Three-Body System, Efimov effect, and Thomas collapse

Now we consider a three-body system of charged particles with mass $m_{i}$, charge $Q_{i}$, and coordinates $\mathbf{r}_{i}$, $i=1,2,3$, interacting via Coulomb and short-range pair interactions as in Eq. (1). For $Q_{i}\equiv 0$, this system features two hallmarks of few-body physics: the Efimov effect [49] and the Thomas collapse [50]. Both are connected in the hyperspherical formalism [51] to a super-attractive $\rho^{-2}$-potential in the hyperradius, $\rho^{2}=(r_{12}^{2}+r_{13}^{2}+r_{23}^{2})/3$, with $r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$. The Thomas collapse occurs due to the divergence of $1/\rho^{2}$ at the origin [29], whereas the infinite tower of Efimov states is supported by the scale invariance of this potential. The $1/\rho^{2}$ form of the potential strongly suggests that the Thomas effect is weakly modified by the Coulomb potential, but the shallow Efimov states must disappear. All in all, this indicates that universal systems of charged particles with zero-range interactions do obey universal scaling laws [5, 1] different from neutral particles.

In order to elucidate the nature of the universal scaling laws in the presence of Coulomb interactions, we investigate the fate of the Efimov effect and Thomas collapse which are hallmarks of universality for neutral systems. We first show the absence of the Efimov effect in a system of two identical heavy charged particles interacting with a light neutral particle, i.e., $m_{1}=m_{2},m_{3}/m_{2}\ll 1,Q_{1}=Q_{2}$ and $Q_{3}=0$. This system is conveniently studied within the Born-Oppenheimer approximation: First, the Schrödinger equation with heavy particles fixed at $\mathbf{r_{1}}$ and $\mathbf{r_{2}}$ is solved,

[TABLE]

which gives the energy $\epsilon(|\mathbf{r}_{1}-\mathbf{r}_{2}|)$. The superscript $lh$ is used to emphasize that the equation describes the light-heavy subsystem. Then, the energy spectrum is found from the two-body equation

[TABLE]

where $\mathbfcal{R}\equiv\mathbf{r}_{1}-\mathbf{r}_{2}$, and $\Phi^{hh}$ is the wave function that describes the relative motion of the two heavy particles. We are interested in the behavior of the potential, $\epsilon+V_{C}+V^{hh}_{S}$, at $\mathbfcal{R}\to\infty$. Therefore, without loss of generality, we assume that (i) $V^{hh}_{S}$ is an infinite barrier for $|\mathbfcal{R}|<R$ and zero otherwise; (ii) $V^{lh}_{S}$ is a separable s-wave interaction. The assumption (ii) allows us to write the function $\epsilon(\mathbfcal{R})$ analytically [52].

For an infinite heavy-light scattering length this function reads: $\epsilon(\mathbfcal{R}\to\infty)\simeq-\frac{\hbar^{2}A}{m_{1}\mathbfcal{R}^{2}}$ where $A>1/4$ generates infinitely many bound states if $V_{C}=0$. For charged particles the Schrödinger equation,

[TABLE]

cannot support infinitely many-bound states. It can support at most $N$ bound states. $N$ can be estimated using the Bargmann inequality [53]:

[TABLE]

where $D$ is the Coulomb length for two heavy particles. It is clear that only if $D/R\gg 1$ there can be many bound states. For example, for the parameters as in 9Be described as an $\alpha+\alpha+n$ system [54, 55] we have $N\lesssim 1.5$, where, for simplicity, we used natural values: $A=1.25$ and $R=1$ fm.

We now move on to charged particles with equal masses. It can be shown that the Coulomb potential also dominates the long-range behavior of the lowest adiabatic potential in the hyperspherical formalism (cf. [56]), which leaves no room for the Efimov effect with identical charged bosons. However, the low-lying Efimov states survive if the Coulomb interaction is sufficiently weak [57]. Our interest is in these states.

We study the Thomas collapse numerically (see B for more details) and observe that the ground state energy behaves similarly to that for neutral particles, i.e., $E\sim-1/R^{2}$, in the vicinity of the zero-range limit ($R\to 0$). This is consistent with the findings of Ref. [47] in their non-perturbative treatment of 3He. Therefore, to study three-body states in the universal limit, we introduce a three-body potential

[TABLE]

where $g_{{\mathrm{3b}}}$ is chosen to fix the three-body energy. As we show below, this three-body force also allows us to study a four-body problem without any additional parameters. This is similar to neutral systems [58, 59, 60],

5 Universal three- and four-body states

We use the Gaussian Expansion Method [24] to study few-body states in the zero-range limit. We obtain the energies by performing a sequence of calculations with small values of $R_{G}$ and extrapolate to $R_{G}\to 0$, which gives the value at $R_{G}=0$ and the error bars. Details of the extrapolation procedure are given in C. Figure 3 reports on the energies of two, three and four charged bosons. As before, the mass and charge of a boson are those of an alpha particle. The energies are fully determined by $D$, $a_{C}$ and an additional three-body parameter. The latter can be characterized either by $a_{C}^{-}$ which determines the three-body binding threshold or by the three-body energy at $1/a_{C}=0$, $\kappa_{3}$ (see Fig. 3). The energy of the three-body state at $1/a_{C}=0$, and, hence, $\kappa_{3}$, is fixed by the three-body force. For neutral particles another value of $\kappa_{3}$ would simply rescale the $y$-axis and $x$-axis due to the discrete scale invariance. For charged particles the discrete scale invariance is broken (cf. Eq. (17)). Thus, we also should investigate the effect of the three-body force; see the inset of Fig. 3. We use two different three-body forces whose $a_{C}^{-}$ differ by a factor of 10, and then rescale the $x$-axis and $y$-axis using $a_{C}^{-}$ and $\kappa_{3}$, correspondingly. We see that the effect of the three-body parameter leads to merely a rescaling of the axes for the considered cases. Therefore, we refrain from showing energies for other values of $a_{C}^{-}$.

We study finite-range corrections numerically; see Fig. 4. Even values of $R_{G}\ll D$ immediately lead to significant corrections to the energy; see the inset of Fig. 4. Similarly to two particles, these corrections can be accounted for by rescaling $a_{C}$. To demonstrate this, we define $a_{C}^{-}(R_{G})$ that determines the three-body binding threshold for a given value of $R_{G}$. We use $a_{C}^{-}(R_{G})$ to rescale the $x$-axis. The rescaled curves coincide for $a_{C}<0$, which suggests that the universal limit is a good starting point for studying Borromean three-body charged systems.

6 Outlook

We have demonstrated that the universal physics of shallow bound states of charged particles is different from neutral particles. The characteristic scale of the system is set by the Coulomb length scale $D$. Unless the Coulomb interaction is extremely weak ($D\to\infty$), the effective range is needed to determine the universal properties of shallow bound states. More detailed studies are required to see in which systems the features discussed here can be observed. To this end, charged quasi-particles and nuclei close to the proton drip-line [21, 61, 18, 62, 63, 15, 16, 23, 64, 65, 66] must be investigated. The concept of universality for neutral particles has also been explored in low and mixed spatial dimensions [67, 68, 69, 70, 71, 72], and with higher angular momenta [73, 74, 75]. It will be interesting to study the effect of the Coulomb potential on those universal states, especially in connection with low-dimensional bipolarons [76, 77, 78], p-wave halo nuclei such as 8B [79, 80], and exotic states of $\alpha$ particles [81]. Finally, it would be interesting to formulate an effective field theory for shallow bound states of charged particles based on our findings and calculate corrections from higher effective range parameters perturbatively.

We thank Wael Elkamhawy for comments on the manuscript, and Daniel Phillips for pointing us to [45]. We also thank Michael Birse, Wael Elkamhawy, Dmitri Fedorov, Aksel Jensen, Daniel Phillips, and Karsten Riisager for useful discussions about proton halos, and Emiko Hiyama for many inspiring conversations about the Gaussian Expansion Method. This work has been supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under project numbers 413495248 – VO 2437/1-1 and 279384907 – SFB 1245 and by the Bundesministerium für Bildung und Forschung (BMBF) through contract 05P18RDFN1.

Appendix A Numerical Methods

We employed the Gaussian expansion method [24] for numerical calculations in the main text. It is a variational method in which the wave function is expanded as a sum of Gaussians. In this section we briefly illustrate the method for a three-body system; see [82, 57] for a more detailed presentation. A variational wave function is written as

[TABLE]

where ${\bf x}={\bf r}_{1}-{\bf r}_{2}$ and ${\bf y}={\bf r}_{3}-({\bf r}_{1}+{\bf r}_{2})/2$ are the Jacobi coordinates, $\mathcal{S}$ is a symmetrization operator [the main text considers only spinless bosons], and $M$ is the basis size. The parameters $\alpha_{i}$ and $\beta_{i}$ are chosen in a form of a geometric progression, i.e., $\alpha_{i}=\alpha_{1}A^{i-1}$ and $\beta_{i}=\beta_{1}B^{i-1}$, where $A$ and $B$ are input parameters. Once the parameters $A,B$ and $M$ are given, the coefficients $a_{i}$ are found by minimizing the expectation value of the Hamiltonian, $E_{\Psi_{GEM}}=\langle\Psi_{GEM}|H|\Psi_{GEM}\rangle/\langle\Psi_{GEM}|\Psi_{GEM}\rangle$. We vary the parameters $A,B$ and $M$ to find the minimal value of $E_{\Psi_{GEM}}$, which is used in the main text as $E$. To benchmark our numerical calculations, we used known results for charged and neutral systems [25, 27, 83]. In addition, we cross-checked some of the energies presented in the main text using the stochastic variational method with correlated Gaussians (SVM) [26, 28]. In our implementation, the SVM assumes the variational wave function in the form

[TABLE]

where the parameters $\alpha^{S}_{i}$, $\beta^{S}_{i}$ and $\gamma^{S}_{i}$ [$\alpha^{S}_{i}\beta^{S}_{i}-\left(\gamma^{S}_{i}\right)^{2}/4>0$] are found by a stochastic search (cf. [26, 28]); $a^{S}_{i}$ are chosen to minimize the expectation value of the Hamiltonian. For weakly-bound neutral states one might restrain these parameters to better reproduce the tails of the wave function (cf. [84]). In our explorations, we saw that this option does not drastically improve convergence for charged systems, probably, because the corresponding tails decay faster with distances than those of neutral systems.

Finally, we briefly explain why our results do not gain any systematic errors due to the choice of the basis. To this end, we show that an eigenstate with $L=0$ ($L$ for total angular momentum) of the Hamiltonian can be accurately approximated by the variational wave function (20). This is not a trivial observation, since $\Psi_{SVM}$ depends only on the three variables: $x,y$, and $\theta_{xy}$ – the angle between ${\bf x}$ and ${\bf y}$. In general, an eigenstate might also depend on other combinations of angles that determine ${\bf x}$ and ${\bf y}$ (e.g., $\theta_{x}$,$\theta_{y}$, $\phi_{x}$, $\phi_{y}$ in spherical coordinates). A suitable angular basis for our discussion is $Y_{l_{1}m_{1}}(\theta_{y},\phi_{y})Y_{l_{2}m_{2}}(\theta_{x},\phi_{x})$, where $Y$ is a real spherical harmonic. Any suitable variational function can be written as

[TABLE]

where $f$ is the function that determines the expansion. We note that the Hamiltonian, $H_{3}$, does not mix the subset $\{P_{l}({\bf x}\cdot{\bf y})\}$ with the rest of the basis, $P_{l}$ is the Legendre polynomial [$P_{l}({\bf x}\cdot{\bf y})=\frac{4\pi}{2l+1}\sum_{m=-l}^{l}Y_{lm}(\theta_{y},\phi_{y})Y_{lm}(\theta_{x},\phi_{x})$]. Indeed, since two-body potentials $V$ depend only on $x,y$ and $\theta_{xy}$, we derive

[TABLE]

where $F_{l^{\prime}}$ is an irrelevant for our discussion function that depends on $x$ and $y$. Therefore, there are eigenstates of $H_{3}$ that can be written as

[TABLE]

It is expected that a square integrable function $\phi$ can be well represented by a suitable $\Psi_{SVM}$ [26]. This becomes intuitively clear after writing Gaussian functions in the form of the corresponding Maclaurin series. One can confirm that $\Psi_{SVM}$ describes the ground state by checking numerically that it does not change sign. For the considered cases, variational results using $\Psi_{SVM}$ and $\Psi_{GEM}$ [the Gaussian expansion method] agree well, which means that the form of $\Psi_{GEM}$ can approximate accurately the function $\phi$. The dependence on $\theta_{xy}$ in $\Psi_{GEM}$ appears due to the presence of the symmetrization operator $\mathcal{S}$.

Appendix B Thomas Collapse

If a zero-range two-body potential is set to reproduce a binding energy of two neutral particles, then the ground state of the corresponding three-body system is infinitely deep. This phenomenon is called the Thomas collapse [50]. A way to deal with this peculiarity in models with zero-range potentials is to introduce a three-body parameter. Here we show numerically that calculations with three charged particles, even in spite of a repulsive Coulomb potential, also require a three-body force. To this end, we compute three-body energies for two-body Gaussian potentials that have different values $R_{G}$ but lead to the same value of the two-body binding energy. For the sake of discussion, we use masses and charges of $\alpha$-particles, and assume that the binding energy is $1$ MeV. The results are presented in Fig. 5. The three-body bound state becomes unphysically deep for $R_{G}\to 0$ (the right panel of the figure suggests that $E_{3}\sim 1/R_{G}^{2}$ for $R_{G}\to 0$), which shows the necessity of a three-body force.

Appendix C Extrapolation Procedure

To obtain results in the limit $R_{G}\to 0$, we first calculate energies for a sequence of small values of $R_{G}$ and then use the extrapolation procedure described below. For trimers, the finite-range corrections to the energy for the smallest numerically accessible values of $R_{G}$ are negligible. Therefore, we present the procedure only for tetramers. The tetramer energies for the four smallest $R_{G}$ ($R_{G}=0.0075\,\text{fm},0.0150\,\text{fm},0.0212\,\text{fm}$ and $0.0374\,\text{fm}$) are presented in Fig. 6c). The three-body forces are taken from the corresponding trimer calculations. The figure shows that even a marginal change from $R_{G}=0.0075\,\text{fm}$ to $R_{G}=0.015\,\text{fm}$ changes noticeably the energies. To obtain energies for $R_{G}\rightarrow 0$ and error bars in the $x$ and $y$ direction, we extract data from Fig. 6c) at fixed values of $1/a_{C}$ and at fixed values of $E$. These data are shown in Figs. 6a) and b) for the four smallest values of $R_{G}$. Then we employ the fit function

[TABLE]

which assumes an analytical functional dependence on $r_{\text{eff}}$ near $r_{\text{eff}}=0$. In addition, we use the fit function

[TABLE]

which includes an additional $r_{\text{eff}}^{3}$ term in comparison to $f$. The fits with $f$ and $\tilde{f}$ yield $b_{1}$ and $b_{2}$, and the error bars on $b_{1}$ and $b_{2}$, which we refer to as $e_{1}$ and $e_{2}$. Note that $e_{2}=0$ because we use only four points for the fit. The average between $b_{1}$ and $b_{2}$ is used as a central value for the dots in Fig. 6c) (see also Fig. 2 of the main text). The error bars correspond to the smallest/largest value of $b_{i}-e_{i}$ / $b_{i}+e_{i}$, respectively. The result of this procedure is shown in Fig. 6c), where each point corresponds either to a vertical or horizontal extrapolation. The figure shows that the two extrapolation directions agree well.