Strongly interacting impurities in a dilute Bose condensate

Let us first examine the case of weak attractive potential with a small scattering length $a<0$. We solve Eq. (21) in the interval $0\leq y\leq 1$ neglecting its right hand side as it contains a small parameter $\epsilon$. Then Eq. (21) reduces to a Schrödinger equation in the potential $U$ at zero energy,

$$-\frac{d^{2}\chi}{dy^{2}}+2mr_{c}^{2}U\chi=0,$$

(23)

whose solution $\chi_{0}$ must satisfy $\chi_{0}(0)=0$. At $y>1$ the potential $U$ is zero, so $\chi$ must be a linear function. Bethe-Peierls boundary conditions fix the form of this function to be

$$\chi_{0}(y)=\frac{\alpha}{1-\frac{r_{c}}{a}}\left(1-\frac{r_{c}}{a}y\right),\ y>1,$$

(24)

where $a$ is the scattering length in the potential $U$ and $\alpha$ is the normalization of the solution chosen in such a way that

$$\chi_{0}(1)=\alpha.$$

(25)

We can rewrite Bethe-Peierls boundary conditions in a convenient way by observing that they are equivalent to

$$\chi^{\prime}_{0}(1)=\frac{\alpha}{1-\frac{a}{r_{c}}}.$$

(26)

We could continue to solve Eq. (21) by means of successive approximations. The next correction $\chi_{1}$ to the solution satisfies

$$-\frac{d^{2}\chi_{1}}{dy^{2}}+2mr_{c}^{2}U\chi_{1}=\epsilon^{2}\left(\chi_{0}-\frac{\chi^{3}_{0}}{y^{2}}\right).$$

(27)

$\chi_{1}$ is a small correction to $\chi_{0}$ since $\epsilon\ll 1$. We will see that $\chi_{1}$ becomes important close to the unitary point where $a\rightarrow\infty$. Then $\chi^{\prime}_{0}(1)=0$ and it is $\chi_{1}$ which produces a nonzero contribution to the derivative. But right now we are interested in small $a$. We will see a little later that we do not need to construct $\chi_{1}$ in the lowest approximation in $\epsilon$ in that case.

Now we solve Eq. (22) neglecting its right hand side to find

$$u=Ae^{-\sqrt{2}z}.$$

(28)

We can construct corrections to it by means of successive approximations. Those will be small as long as $A$ is small, which as we will see in a second is the case here. So for now we neglect those corrections.

Matching amplitudes and derivatives of $\chi(y)$ and $u(z)$ at $z=\epsilon$, or correspondingly $y=1$, produces

$$A=\epsilon(\alpha-1),\ -\sqrt{2}A=\frac{\alpha}{1-\frac{a}{r_{c}}}-1.$$

(29)

Taking into account that $\epsilon\ll 1$, we find

$$A=-\frac{a}{\xi},\ \alpha=1-\frac{a}{r_{c}}.$$

(30)

Let us now examine if the terms neglected to arrive at this solution are indeed small. We solve (21) by successive approximations, plugging $\chi_{0}$ into the right hand side of Eq. (21) and producing a correction $\chi_{1}$. If $\left|a\right|<r_{0}$ then both $\chi_{0}(1)$ and $\chi^{\prime}_{0}(1)$ are of the order of $1$ while $\chi_{1}$ will be of the order of $\epsilon^{2}\ll 1$ and can be neglected. It gets more interesting if $\left|a\right|>r_{c}$. Then $\chi_{0}(1)=\alpha\sim a/r_{0}$, while $\chi^{\prime}_{0}(1)\sim 1$. At the same time, $\chi_{1}\sim\epsilon^{2}(a/r_{c})^{3}$. The magnitude of this had better be smaller than $1$, so that the contribution of $\chi^{\prime}_{1}(1)$ to the derivative of $\chi$ could be neglected. This condition gives $\epsilon^{2}{\left|a\right|^{3}}/{r_{c}^{3}}\ll 1$, or equivalently

$$\left|a\right|^{3}\ll\xi^{2}r_{c}.$$

(31)

This is the condition to neglect the right hand side of Eq. (21). We will see later that this condition signifies the transition from weak impurity potential, where Eq. (31) holds true, to the strong potential including the unitary point $a\rightarrow\infty$, where it is violated.

Under the condition (31) $A\ll 1$, so obviously we can indeed neglect the right hand side of (22) as we did above.

We can now plug our solution into Eq, (18). The integral in this formula needs to be split into two parts, $0<r<r_{c}$ and $r_{c}<r$. It is easy to see that the contribution of the interval $0<r<r_{c}$, under the condition (31), is negligible, so we only need to integrate over $r>r_{c}$ with $u=A\exp(-\sqrt{2}z)$. Performing the integration and again taking into account Eq. (31) to get rid of some terms suppressed relative to the main contribution we find

$$E=\frac{2\pi n_{0}a}{m}.$$

(32)

This is the well known result for the energy of the polaron at weak scattering length $|a|\ll\xi$. We now see that the validity condition for this to hold true is Eq. (31).

Suppose now that the potential $U$ is made more attractive so that its scattering length increases, violating the condition (31) and eventually reaching infinity at the unitary point. We can follow the same strategy to obtain the solution in this case. The new element is that Bethe-Peierls boundary conditions now imply $\chi^{\prime}_{0}(1)=0$, so we need to solve Eq. (21) perturbatively, using its right hand side as a perturbation, to find nonzero $\chi^{\prime}_{1}(1)$ contributing to $\chi^{\prime}(1)$. Same goes for Eq. (22).

In case when $a$ was small we found earlier that the amplitude $\chi_{0}(1)=\alpha$, $\alpha=1-a/r_{c}$. This is of the order of $1$ when $a$ is very small, but starts growing as $|a|$ is increased. When $|a|^{3}\sim\xi^{2}r_{c}$, $\alpha\sim\epsilon^{-2/3}$. We will see that $\chi(1)$ remains of the order of $\epsilon^{-2/3}$ even as $a$ is taken all the way to infinity, so it is convenient to introduce the notation

$$\beta=\alpha\epsilon^{2/3},$$

so that

$$\chi_{0}(y)=\frac{\beta}{\epsilon^{2/3}}v(y).$$

(33)

Here $v$ is the solution of the Schrödinger equation

$$-\frac{d^{2}v}{dy^{2}}+2mr_{c}^{2}Uv=0$$

(34)

normalized so that $v(1)=1$. $v^{\prime}(1)=0$ since $U$ is tuned to unitarity. $\beta$ is thus a yet unknown $\epsilon$-independent normalization coefficient.

We will need a correction to this which satisfies

$$-\frac{d^{2}\chi_{1}}{dy^{2}}+2mr_{c}^{2}U\chi_{1}=-\epsilon^{2}\frac{\chi_{0}^{3}}{y^{2}}.$$

(35)

The term $\epsilon^{2}\chi_{0}$ from the right hand side of Eq. (27) goes as $\epsilon^{4/3}$ and can be neglected. At the same time, we see that $\chi_{1}$ is of the order of $\epsilon^{0}$. Solving Eq. (35) gives

$$\chi_{1}=\beta^{3}v(y)\int_{0}^{y}\frac{ds}{v^{2}(s)}\int_{0}^{s}\frac{dt\,v^{4}(t)}{{t}^{2}}.$$

Putting it together produces

$$\chi=\frac{\beta}{\epsilon^{2/3}}v(y)+\beta^{3}v(y)\int_{0}^{y}\frac{ds}{v^{2}(s)}\int_{0}^{s}\frac{dt\,v^{4}(t)}{{t}^{2}}+\dots.$$

(36)

The next term $\chi_{2}$ which can be obtained by continuing successive approximations goes as $\epsilon^{2/3}$. We will not need it here, but note that it will have an even more complicated dependence on $v$ and by extension on features of the potential $U(r)$ than the already obtained term $\chi_{1}$.

From this solution we find that

$$\chi(1)=\frac{\beta}{\epsilon^{2/3}}+\mathcal{O}(1),\quad\chi^{\prime}(1)=\beta^{3}c+\mathcal{O}(\epsilon^{2/3}),$$

(37)

where the dimensionless coefficient $c$ is defined via

$$\quad c=\int_{0}^{1}\frac{dy\,v^{4}(y)}{y^{2}}.$$

(38)

Now we turn our attention to Eq. (22). Its solution $u(z)$ needs to be matched with the boundary conditions (37). Easy to verify that these boundary conditions imply

$$u(\epsilon)=\beta\epsilon^{1/3}+\mathcal{O}(\epsilon),\quad u^{\prime}(\epsilon)=-1+\beta^{3}c+\mathcal{O}(\epsilon^{2/3}).$$

(39)

Eq. (22) differs from Eq. (21) in that its nonlinear terms do not have an explicit factor of $\epsilon$ in front of them. We will nonetheless solve Eq. (22) by means of successive approximations, and verify later that this is a legitimate approach. Without its right hand side, the solution to Eq. (22) reads as before,

$$u_{0}(z)=A\,e^{-\sqrt{2}z}.$$

(40)

We can plug it into the right hand side of Eq. (22), however we will follow a slightly different procedure. We use Eq. (40) to rewrite Eq. (22) as an integral equation via a standard procedure. This involves solving the auxiliary equation

$$\frac{d^{2}u}{dz^{2}}-2u=g(z)$$

with arbitrary given $g(z)$, then substituting the actual right hand side of Eq. (22). We find

$$u(z)=A\,e^{-\sqrt{2}z}+\frac{e^{-\sqrt{2}z}}{2\sqrt{2}}\int_{z}^{\infty}\frac{ds\,e^{\sqrt{2}s}}{s}\left(3u^{2}(1+u^{\prime})+\sqrt{2}u^{3}\right)-\frac{e^{\sqrt{2}z}}{2\sqrt{2}}\int_{z}^{\infty}\frac{ds\,e^{-\sqrt{2}s}}{s}\left(3u^{2}(1+u^{\prime})-\sqrt{2}u^{3}\right).$$

(41)

We now use this equation to calculate $u(\epsilon)$ and $u^{\prime}(\epsilon)$ in perturbative expansion in powers of $\epsilon$. Anticipating that the leading behavior $A\sim\epsilon^{1/3}$, as should be clear from comparing Eq. (40) and Eq. (39), we iterate Eq. (41) by plugging $u_{0}(z)$ into the right hand side of Eq. (41). The resulting integrals can be computed in terms of Gamma functions and expanded in powers of $\epsilon$. We omit rather lengthy algebra involved, and just state that this allows us to evaluate $u(\epsilon)$ to find

$$u(\epsilon)=A+\frac{3\ln 3}{2\sqrt{2}}A^{2}-A^{3}\ln\epsilon+\mathcal{O}(\epsilon).$$

(42)

We also evaluate $u^{\prime}(\epsilon)$. Differentiating Eq. (41) gives

$$u^{\prime}(z)=-\sqrt{2}Ae^{-\sqrt{2}z}-\frac{u^{3}}{z}-\frac{e^{-\sqrt{2}z}}{2}\int_{z}^{\infty}\frac{ds\,e^{\sqrt{2}s}}{s}\left(3u^{2}(1+u^{\prime})+\sqrt{2}u^{3}\right)-\frac{e^{\sqrt{2}z}}{2}\int_{z}^{\infty}\frac{ds\,e^{-\sqrt{2}s}}{s}\left(3u^{2}(1+u^{\prime})-\sqrt{2}u^{3}\right).$$

(43)

We can again substitute $u_{0}(z)$ into the integrals on the right hand side of Eq. (43), to find

$$u^{\prime}(\epsilon)=-\sqrt{2}A-\beta^{3}+3A^{2}\ln\epsilon+\mathcal{O}(\epsilon^{2/3}).$$

(44)

Here we took advantage of the boundary conditions (39) which tell us that $u(\epsilon)=\beta\epsilon^{1/3}$ within the accuracy that we work with.

Combining Eq. (42) and Eq. (44) with Eq. (39) gives

$$\begin{matrix}A+\frac{3\ln 3}{2\sqrt{2}}A^{2}-A^{3}\ln\epsilon&=&\beta\epsilon^{1/3},\cr-\sqrt{2}A-\beta^{3}+3A^{2}\ln\epsilon&=&-1+\beta^{3}c.\end{matrix}$$

(45)

We now need to solve these equations for $A$ and $\beta$ perturbatively, in powers of $\epsilon$. Introduce

$$\delta=\frac{\epsilon}{1+c}.$$

The solution to Eq. (45) reads

	$\displaystyle A$	$\displaystyle=$	$\displaystyle\delta^{1/3}-\left(\frac{3\ln 3}{2\sqrt{2}}+\frac{\sqrt{2}}{3}\right)\delta^{2/3}+2\delta\ln\delta+\mathcal{O}(\delta),$		(46)
	$\displaystyle\beta\epsilon^{1/3}$	$\displaystyle=$	$\displaystyle\delta^{1/3}-\frac{\sqrt{2}}{3}\delta^{2/3}+\delta\ln\delta+\mathcal{O}(\delta).$		(47)

We can now use the parameters we obtained in this way to calculate the energy and the particle number of the polaron. It turns out to be technically easier to calculate the particle number first and then use Eq. (3) to find the energy, which is the strategy we will follow here.

The excess number of particles due to the polaron is given by

$$N=\int d^{3}x\left[\left|\psi\right|^{2}-n_{0}\right]=4\pi n_{0}\xi^{3}\int_{0}^{\infty}z^{2}dz\left[\phi^{2}-1\right].$$

It is natural to split the integral over $z$ into two intervals, from $0$ to $\epsilon$ and from $\epsilon$ to infinity. Now the contribution of the first interval can be safely neglected. Indeed, it gives

$$\int_{0}^{\epsilon}z^{2}dz\left[\phi^{2}-1\right]=\epsilon^{3}\int_{0}^{1}dy\left(\frac{\chi^{2}}{y^{2}}-1\right)\sim\epsilon^{5/3}.$$

Here we used that $\chi(y)\sim 1/\epsilon^{2/3}$. This contribution is very small and exceeds the accuracy in $\epsilon$ with which we were doing our calculations. This also indicates that the bulk of the particles bound by the impurity are located farther than distance $r_{c}$ away from the impurity.

The contribution of the second interval gives

$$\int_{\epsilon}^{\infty}z^{2}dz\left(\left(1+\frac{u}{z}\right)^{2}-1\right).$$

(48)

To evaluate this integral we again iterate Eq. (41) once, to find $u$ up to the terms of the order of $A^{2}$, and substitute that into Eq. (48). The result of this evaluation is

$$N=4\pi n_{0}\xi^{3}\left(\delta^{1/3}-\frac{5}{3\sqrt{2}}\delta^{2/3}+2\delta\ln\delta+\dots\right).$$

(49)

Thus we evaluated the number of particles trapped in the polaron up to terms of the order of $\delta\ln\delta$. To go beyond this order, starting from terms of the order of $\delta$ and beyond represented by the dots above, we would need to go beyond the terms presented in Eq. (36). We expect that this will produce terms which depend on the features of the potential other than the coefficient $c$.

To construct the energy of the polaron, it is easiest at this stage to take advantage of Eq. (3). The subtlety in evaluating the derivative there is that the particle number $n_{0}$ as well as $\xi$ have to be traded for $\mu$ before differentiating. Doing the algebra we arrive at

$$E=-\frac{\pi n_{0}\xi}{m}\left(3\delta^{\frac{1}{3}}-2\sqrt{2}\delta^{\frac{2}{3}}+4\delta\ln\delta+\dots\right).$$

(50)

This constitutes the answer that we seek, the energy of the polaron at unitarity as an expansion in powers of $\delta\sim r_{c}/\xi$.

Finally, let us examine $\delta=\epsilon/(1+c)$ in a little more detail. From the definition of $c$ given in Eq. (37) we can write

$$1+c=1+r_{c}\int_{0}^{r_{0}}\frac{dr\,v^{4}}{r^{2}}=r_{c}\left(\int_{r_{c}}^{\infty}\frac{dr}{r^{2}}+\int_{0}^{r_{c}}\frac{dr\,v^{4}}{r^{2}}\right).$$

$v$ is the solution of the Schrödinger equation with the potential tuned to unitarity, so that $v^{\prime}(r_{c})=0$. Since it is normalized such that $v(r_{c})=1$, it will naturally satisfy $v(r)=1$ for all $r\geq r_{c}$. Therefore we can rewrite this as

$$1+c=r_{c}\int_{0}^{\infty}\frac{dr\,v^{4}}{r^{2}}.$$

Now

$$\delta=\frac{\epsilon}{1+c}=\frac{1}{\xi\int_{0}^{\infty}\frac{dr\,v^{4}}{r^{2}}}.$$

It is now convenient to define

$$R^{-1}=\int_{0}^{\infty}\frac{dr\,v^{4}}{r^{2}}=\int\frac{d^{3}x}{4\pi}\,\left|\psi_{0}\right|^{4},$$

(51)

where $\psi_{0}=v/r$ is the solution of the Schrödinger equation

$$-\frac{\Delta}{2m}\psi_{0}+U\psi_{0}=0.$$

$R$ constitutes a properly defined range of the potential, finite even for potentials which extend all the way to infinity, which correctly captures its extent for the purposes of solving the polaron problem.

This gives us a definition

$$\delta=\frac{R}{\xi}.$$

(52)

At this stage $r_{c}$ drops from the equations and no longer needs to be finite. It can be taken to infinity if desired, with the answer for the energy of the polaron (50) as well as the number of particles trapped in the polaron (49) expressed entirely in terms of $\delta=R/\xi$.

We can go further and generalize the above analysis to account for the small deviations away from unitarity using a $1/a$ expansion. To accomplish this, it is convenient to parametrize $a$ by a new variable $\eta$ according to

$$\eta=\frac{r_{c}^{1/3}\xi^{2/3}}{a}.$$

(53)

Declaring $\eta\ll 1$ is equivalent to saying that we are in the strongly interacting regime (compare with the condition (31) which is violated once we cross over from weak to strong potentials as explained earlier). We then proceed to work out the expansion in powers of $\eta$. The advantage of this reparametrization is due to form that Bethe-Peierls boundary conditions take. Indeed, Eq. (26) at large $a$ where close to unitarity we should write $\alpha=\beta/\epsilon^{2/3}$ imply

$$\frac{\beta}{\epsilon^{2/3}(1-\frac{a}{r_{c}})}\approx-\frac{\beta r_{c}}{\epsilon^{2/3}a}=-\beta\eta.$$

(54)

Because $r_{c}/a\sim\epsilon^{2/3}$ the correction to the nonlinear term and to the RHS of the top equation in (45) that arise from correction to $v(y)$ in Eq. (36) are of the higher order in power of $\epsilon$ and thus can be neglected. Thus slighly away from unitarity the only modification to the formalism presented in the previous section is the addition of the $-\beta\eta$ term to RHS of the second equation in (45):

$$\begin{matrix}A+\frac{3\ln 3}{2\sqrt{2}}A^{2}-A^{3}\ln\epsilon&=&\beta\epsilon^{1/3},\cr-\sqrt{2}A-\beta^{3}+3A^{2}\ln\epsilon&=&-1-\beta\eta+\beta^{3}c.\end{matrix}$$

(55)

Notice that now in order to find the leading term in the expansion of $\beta$ in powers of $\epsilon$, we have to solve the cubic equation:

$$\beta_{0}^{3}(1+c)-\beta_{0}\eta-1=0.$$

The number of real roots of this equation depends on the value of the discriminant $\Delta=4(\frac{\eta}{1+c})^{3}-\frac{27}{(1+c)^{2}}$. If $\Delta>0$, there are three real roots, and if $\Delta<0$, only a single real root. The critical value of $\eta$ is obtained by setting $\Delta=0$ and gives $\eta_{c}=(\frac{27(1+c)}{4})^{1/3}>0$. A positive value of $\eta$ corresponds to a positive value of the scattering length, which means that we are in the regime where the potential has a single bound state. It is known that when a potential admits $\nu$ number of bounds states, then the Gross-Pitaevskii equation can have up to $2\nu+1$ solutions [22]. This means that the approach outlined above can in principle be used to construct other solutions for the Gross-Pitaevskii equation, but since we are interested only in the ground state physics, we do not consider them here. Instead, just as already elaborated above we are going to assume that $\left|\eta\right|\ll 1$, so that we are in the regime where there is only a single solution exists and we can do expansion in $\eta$ on top of the expansion in $\epsilon$. The first order correction in $1/a$ corresponds to the first order correction in $\eta$. Note that the case of large and positive scattering length corresponds to the presence of the bound state with energy $E_{\text{binding}}=-\frac{1}{2ma^{2}}$. Since $r_{c}/a=\eta\epsilon^{2/3}$, in order to capture $1/a^{2}$ effects we would need to construct the further corrections to the matching equations that would include higher powers of $\epsilon$ and therefore depend on other details of the potential. Therefore, in the lowest order approximation in $\eta$ that we are going to employ, we are not sensitive to the presence of the bound state even if $a$ is positive.

Section VII further elaborates on what we expect to happen in the regime where $a>0$ and a bound state is present once we move away from the unitary point. In particular, we expect that in this regime the behavior of the polaron strongly depends on the details of the impurity-boson potential $U$.

Let us go back to staying in the vicinity of the unitary point by demanding $|\eta|\ll 1$. Since $\frac{\eta}{(1+c)^{1/3}}=\frac{\delta^{1/3}\xi}{a}$, solution to the system (55) reads:

	$\displaystyle A$	$\displaystyle=$	$\displaystyle\delta^{1/3}\left(1+\frac{\delta^{1/3}\xi}{3a}\right)-\left(\frac{3\ln 3}{2\sqrt{2}}+\frac{\sqrt{2}}{3}\right)\delta^{2/3}-\frac{\ln 3}{\sqrt{2}}\cdot\frac{\xi\delta}{a}+2\left(1+\frac{2\delta^{1/3}\xi}{3a}\right)\delta\ln\delta+\mathcal{O}(\delta),$		(56)
	$\displaystyle\beta\epsilon^{1/3}$	$\displaystyle=$	$\displaystyle\delta^{1/3}\left(1+\frac{\delta^{1/3}\xi}{3a}\right)-\frac{\sqrt{2}}{3}\delta^{2/3}+\left(1+\frac{\delta^{1/3}\xi}{3a}\right)\delta\ln\delta+\mathcal{O}(\delta).$		(57)

Note that at this step $r_{c}$ dropped out just like it did in the previous subsection of this paper, getting replaced by $R$ via the definition (52). $R$ is defined even with for potentials which extend all the way to infinity, so the results obtained here, just as elsewhere in this paper, are valid as long as $R$ can be defined via Eq. (51) and is finite.

Repeating the same steps that lead to (49) and (50) we finally obtain the expression for the number of trapped bosons and the energy of the polaron:

	$\displaystyle N$	$\displaystyle=$	$\displaystyle 4\pi n_{0}\xi^{3}\left[\delta^{\frac{1}{3}}-\frac{5}{3\sqrt{2}}\delta^{\frac{2}{3}}+2\delta\ln\delta+\dots+\frac{\xi\delta^{\frac{1}{3}}}{3a}\left(\delta^{\frac{1}{3}}-\sqrt{2}\delta^{\frac{2}{3}}+4\delta\ln\delta+\dots\right)\right],$		(58)
	$\displaystyle E$	$\displaystyle=$	$\displaystyle-\frac{\pi n_{0}\xi}{m}\left[3\delta^{\frac{1}{3}}-2\sqrt{2}\delta^{\frac{2}{3}}+4\delta\ln\delta+\dots+\frac{\xi\delta^{\frac{1}{3}}}{a}\left(2\delta^{\frac{1}{3}}-\frac{4\sqrt{2}}{3}\delta^{\frac{2}{3}}+4\delta\ln\delta+\dots\right)\right].$		(59)

We verify the above expressions by numerically solving the Gross-Pitaevskii equation (17), where we choose $U$ to be the square well potential. The analytic expression for the scattering length in this potential is well known, so it is easy to tune the strength of the potential to get the desired value of $a$. Far away from impurity one can use the asymptotic solution $\phi=1+\frac{u}{z}$, where $u$ is given by Eq. (28). With two values of parameter $A$, $A_{\text{min}}$ and $A_{\text{max}}$, our algorithm performs the Newton bisection untill we find the value of $A$ such that the derivative of the solution at the origin is zero. The graphs of the polaron energy and the number of trapped bosons slighly away from unitarity are presented in Fig. 1 and Fig. 2.

Now that we have established the machinery for solving the Gross-Pitaevskii equation in both regions, let us revisit the case of the weak potential and find the next order correction to the energy (32) in powers of the scattering length $a$. Note that the contribution to the energy from the first interval goes in powers of $a/r_{c}$, while in the second interval in powers of $a/\xi=a\epsilon/{r_{c}}$, so in principle, the energy expansion should be in powers of $(\frac{a}{r_{c}})^{i}\epsilon^{j}$. Moreover, from the equation one can show that the contribution of the region $z\in[0,\epsilon]$ to the energy goes as $\sim\frac{n_{0}a}{m}(\epsilon+\frac{a}{\xi}\epsilon)$, while the region $z\in[\epsilon,\infty]$ contributes $\sim\frac{n_{0}a}{m}(1+\frac{a}{\xi})$, so in the limit $\epsilon\rightarrow 0$, only the second region is important. The $\sim a^{2}$ contribution to the energy comes from $\frac{a^{2}}{r_{c}^{2}}\epsilon^{2}$ term, so in order to find the energy up to order $a^{2}$, one needs to find the expression for the amplitude $A$ in the second region up-to order $\epsilon^{2}$. The structure of equations in (29) tells us that the matching of amplitudes determines $A$, while the matching of derivatives determines $\alpha$. Since we need to know only $A$ up to order $\epsilon^{2}$, we need to correct the first equation by the terms of order $\epsilon^{2}$, while the second only by the terms of order $\epsilon$. Notice that iteration of the RHS of Eq. (22) will produce terms of higher order than we need, so this tells us that the result should be independent on the further details of the potential. Having in mind that $A\sim\epsilon$, in the second interval we expand Eq. (41) up to $\epsilon^{2}$ terms to produce:

$$\begin{matrix}A+\frac{3\ln 3}{2\sqrt{2}}A^{2}-\sqrt{2}\epsilon A&=&\epsilon(\alpha-1),\cr-\sqrt{2}A&=&\frac{\alpha}{1-\frac{a}{r_{c}}}-1.\end{matrix}$$

(60)

Looking for the solution of this system in the form $\alpha=\alpha_{0}+\epsilon\alpha_{1}$ and $A=\epsilon A_{0}+\epsilon^{2}A_{1}$, we get

	$\displaystyle\alpha$	$\displaystyle=$	$\displaystyle 1-\frac{a}{r_{c}}+\sqrt{2}\frac{a}{r_{c}}\left(1-\frac{a}{r_{c}}\right)+\mathcal{O}(\epsilon^{3}\ln\epsilon),$		(61)
	$\displaystyle A$	$\displaystyle=$	$\displaystyle-\frac{a}{r_{c}}\epsilon-\frac{(4+3\ln 3)}{2\sqrt{2}}\left(\frac{a}{r_{c}}\epsilon\right)^{2}+\mathcal{O}(\epsilon^{3}\ln\epsilon).$		(62)

The energy of the polaron is given by

$$E=-\frac{\lambda}{2}\int d^{3}x\left(|\psi|^{4}-n_{0}^{2}\right)=-2\pi\lambda n_{0}^{2}\xi^{3}\int_{0}^{\infty}dzz^{2}\left(\phi^{4}-1\right).$$

(63)

The region from 0 to $\epsilon$ does not contribute, while the second interval produces:

$$\int_{\epsilon}^{\infty}dzz^{2}\left(\left(1+\frac{u}{z}\right)^{4}-1\right).$$

(64)

Once again, we use Eq. (41) with $A$ given by Eq. (61) to compute the integral in Eq. (64) and expand the result up-to $\epsilon^{2}$ terms. Plugging this into Eq. (63), we finally produce

$$E=\frac{2\pi n_{0}a}{m}\left(1+\frac{\sqrt{2}a}{\xi}+\dots\right).$$

(65)

The second term in the brackets was obtained before in Ref. [28]. This shows that the Gross-Pitaevskii approach discussed above is well suited for studying both weak potentials and potentials tuned to unitarity.

Having established the expressions for the solution of the Gross-Pitaevskii equation, now we can compute other key quasiparticle properties such as quasiparticle residue $Z$ and Tan’s contact $C$. This was already reported in Ref. [1]. Let us reproduce it here for completeness.

The residue $Z$ quantifies the overlap between the solutions in presence and absence of the impurity. Within the GP treatment, this is given by [29]

$$\ln Z=-\int d^{3}x\,|\psi(x)-\sqrt{n_{0}}|^{2}.$$

At unitarity, the above analysis shows that to leading order

$$\ln Z=-\sqrt{2}\pi n_{0}\xi^{3}\delta^{2/3}+\dots.$$