Engineering momentum profiles of cold-atom beams
D Hudson Smith Artem G Volosniev





Clemson University, Clemson, South Carolina 29634, USA


Institut für Kernphysik, Technische Universität Darmstadt, 64289 Darmstadt, Germany





We describe a procedure for engineering beams of cold atoms by selectively draining particles from a trapped gas based on momenta. Atoms escape through a filter potential that only transmits atoms with the desired momenta. We outline an algorithm that outputs a filter potential that produces a pre-specified beam momentum profile. We illustrate this procedure for the case of a narrow band-pass (NBP) quantum filter. Lastly, we discuss the application of the NBP filter for probing the self-energy and effective mass of Bose polarons, as well as the corresponding Landau criterion.



Engineering momentum profiles of cold-atom beams
    Artem G Volosniev
    December 30, 2023
================================================





Introduction. –   
In this article we discuss a procedure for creating cold-atom beams with momentum transport profiles that can be selected for the matter at hand.
Such beams would enable novel scattering experiments with quantum gases. In particular, 
they could be used to measure parameters that define
few- and many-body physics of cold-atoms systems, e.g., the
scattering length, three-body parameter <cit.> or
the self-energy and the effective mass of a polaron <cit.>.
The beams can also be used as the initial state for atom interferometry <cit.>.



Figure  <ref> summarizes our proposal. Analogous to a quantum switch device (`transistor') <cit.>, the flux of particles from the `source' (reservoir) is determined by the `gate' (link potential). However, rather than simply controlling the overall transmission rate, our proposal allows one to design the momentum profile of the outgoing flux.
 For simplicity, we illustrate our idea using a one-dimensional geometry (1D) (though our formalism also applies directly to a cylindrically symmetric three-dimensional geometry). We assume that the particles in the reservoir are non-interacting [Strong interactions can lead to the transmission behavior that is not captured by the one-body Schrödinger equation, e.g.,to a collective resonant transport <cit.>.], so that their scattering properties may be calculated using the one-body Schrödinger equation:

    -ħ^2/2m∂^2/∂ x^2ψ+V_0(x)ψ=ħ^2k^2/2mψ,

where m is the mass of a particle from the reservoir, ħ^2k^2/2m is its energy. The link potential V_0(x) produces the transmission coefficient, T_0(k). By carefully tuning V_0(x) one produces a T_0(k) that allows only particles with desired momenta to tunnel through the barrier into the flux region as required by the experimental application. The tunability of V_0(x) in cold-atom set-ups has been recently demonstrated <cit.>, 
suggesting that our proposal relies only on the toolbox available in current cold-atom laboratories.
Note that the momentum profile of the outgoing beam can be measured using single-atom momentum resolution techniques (e.g., <cit.>), allowing one to confirm that the beam has the desired flux profile.

In this article, we discuss how to determine an appropriate V_0(x) for a given desired flux profile T_0(k). Furthermore, we briefly discuss the 1D Bose polaron problem as a possible application of cold-atom beams. 

Other studies related to the engineering of atomic beams have focused on  the  outcoupling  of  atoms  from a Bose-Einstein condensate <cit.> or on controlling the transmission of atoms with careful manipulation of an optical lattice <cit.>. In contrast, the current article proposes a method for engineering atomic beams by coupling a trapped gas of atoms to a combination of individually tuned lasers. Moreover, the proposed scheme is generic in the sense that one can approximately produce arbitrary desired momentum profiles of the transmitted atom beam.





Procedure for Finding a Link Potential. – 
We find an appropriate link potential V_0(x) ≡ V(x,θ^*) by performing a global search over a family of possible potentials V(x,θ) for the parameters θ^* that reduce the k-integrated squared error between the desired transmission-momentum profile T_0(k) and the actual profile T_θ(k) produced by a sample potential V(x, θ). Concretely, we minimize the cost

    J_θ = ∑_0<k<k_F w_k|T_0(k) - T_θ(k)|^2,

where the k-integral has been approximated (up to a constant factor) by a sum over a discrete set of momentum values, the interval of integration is [0,k_F], and w_k is the weight given to momentum k. The weights are chosen to emphasize or de-emphasize special regions of k during the minimization. For instance, for a narrow band-pass filter that forbids transmission for all k except in the neighborhood of a chosen value k_0 (see Figs. <ref> and <ref>), it is appropriate to increase the weights in the region of k_0. The convergence of our approach is somewhat sensitive to the choice of w_k.

We developed a strategy for choosing the weights that takes into account the following considerations: i) the transport coefficient is zero for k=0 
(there can be no transmission for k=0 and V≠ 0 in 1D), so the weight can be smaller for small k, ii) the transport coefficient approaches 1 for large k so higher weights are required in the large-k region if one wishes to suppress flux at large k, and iii) to reproduce narrow features in the target transport profile, it may be helpful to increase the weight in the k-region of these features. Following these principles, we arrive at the following formula for the weight function:

    w(k;r) = w_bg(k) + rT_0(k)

where w_bg(k) is chosen to account for considerations i) and ii) above, and the term proportional to the positive constant r accounts for consideration iii). The form  of w_bg(k) can be inferred from typical transmission coefficients. For convenience, we use the analytic form for the transmission coefficient produced by
the Morse potential ħ^2k_0^2/[2m cosh^2(k_0x/√(2))] (cf. <cit.>):

    w_bg(k) = [sinh^2(√(2)π k/k_0)/sinh^2(√(2)π k/k_0) + cosh^2(√(7)π/2)]^1/2,

where the parameter k_0 determines a typical energy scale (for an example see our illustration below).
The parameter r was chosen by trial and error for each target profile T_0(k).

With the goal of discovering experimentally viable solutions, we parameterize the family of link potentials V(x, θ), as a sum of N Gaussians, each of the form

    V_i(x; A_i, μ_i, σ_i) = A_i/√(2πσ_i^2)exp[-(x-μ_i)^2/2σ_i^2];

the parameter θ denotes the parameter space {A_1, μ_1, σ_1,...,A_N, μ_N, σ_N}. While minimizing Eq. (<ref>), we enforce the parameter constraints listed in Table <ref>. In addition to these constraints on the parameters, we enforce a constraint requiring that the link potential should not extend beyond the region of potential support x∈[-x_0,x_0]. 
To accomplish this, we minimize the boundary-augmented cost function

    J_θ^aug = J_θ + α∑_i=1^N∫_|x|>x_0dx |V_i(x; A_i,μ_i,σ_i)|^2,

where α is a tuning parameter chosen empirically to aid in convergence. 
 If the potential V_i is outside [-x_0,x_0] then the boundary-augmented 
cost function evaluates to J_θ≃ J_θ +α A_i, which means that α
should be chosen of the order J_θ/A_min to dictate the constraint.
 The integral in Eq. (<ref>) evaluates to the complementary error function. The full form of J_θ^aug is given in Appendix <ref>.





For a particular choice of θ (and hence V(x;θ)), we solve for T_θ(k) by integrating the Schrödinger equation (<ref>) across the region of the potential and calculating the ratio of the transmitted to the incident flux. In order to do this efficiently, we discretize the second derivative in Eq. (<ref>), which transforms Eq. (<ref>) into
a banded linear system of equations solvable in O(M) time where M is the number of x-steps. Using these techniques we are able to evaluate the transmission coefficient 3.7 thousand times per second on a 7th generation Intel Core i7 processor for a test involving 1,000 randomly generated two-Gaussian potentials and 100 different scattering momenta. 






We minimize J_θ^aug for θ^* using the global optimization routine called Differential Evolution (DE) <cit.>. This evolutionary-based search algorithm is suitable given the non-convex (multiple local minima) nature of the optimization problem and the continuity of the parameter space. Despite its simplicity, DE does a good job of balancing exploration of the space of link potentials against the need to efficiently learn from each sample with little tuning of the model settings. Empirically, we found DE to converge to good solutions much more quickly than random search perhaps due to DE's ability to incorporate information from the previous iteration of the optimization algorithm. 

Narrow band-pass (NBP) filter transport profile. –
To illustrate the method described above we optimize for a NBP filter transport profile sharply-peaked near k=k_0. 
For our target transport profile, we use a Lorentz profile (see Figure <ref>)

    T(k; k_0,b) = [1 + (k-k_0)^2/b^2]^-1

where k_0 determines the peak position and b determines the width. 
For simplicity, in this subsection we adopt the units k_0=ħ=2m=1, which scales k_0
out of the problem. The value of b must be much smaller than k_0 
to have a well-defined peak, but not too small to have realistic time scales
for a one-body tunneling. We set b=0.03 k_0, which for a reasonable assumption 
ħ^2k_0^2/(2m)=k_B×μK, where k_B is the Boltzmann constant leads to 
the time scale associated with the resonance width 2m/ħ b^2∼ 10ms.

 For the constraints shown in Table <ref>, we use σ_min=0.2, σ_max=3, A_min=5, and A_max=30. We set k_F=2. We further simplify the optimization by searching over two-Gaussian link-potentials with equal amplitudes and widths. We anticipate that this potential might be the easiest to realize in the laboratory. 
Moreover, it allows us to give a physical interpretation in terms of quasi-discrete energy levels supported by the link.  Even though we work here with a very simple example with only three unknown parameters, a method for globally searching the space of possible link potentials is still required, because the cost function for the optimization has many local minima corresponding to the many ways to produce a resonance state near the scattering energy k_0^2. Moreover this global search technique extends to more complicated transport profiles which necessitate more complicated families of link potentials.

Figure <ref> shows the link potential and transport profile resulting for the NBP filter optimization. The solution suppresses transport except near k=1 as set by our Lorentz target profile and near k=2 resulting from a second resonance in the scattering potential. 
In general, the transmission coefficient will always be non-zero to the right of the target profile. 
This, however, need not be a problem if the atoms are sourced from a thermal reservoir with sufficiently low population to the right of the target profile [For a reservoir at a finite temperature it is beneficial to include the energy distribution, n(k),
in the cost function. To this end, one should simply convolute T_θ(k) in Eq. (<ref>) with n(k). 
The function n(k) is determined by the Bose-Einstein or Fermi-Dirac distributions, temperature and the shape of the trap. For simplicity, 
we assume in current Eq. (<ref>) a one-dimensional Fermi gas at zero temperature in a harmonic trap, i.e., n(k<k_F)=1 and zero otherwise.].
Moreover, the transmission coefficient can be further shaped by using a second NBP filter.
 
It may be experimentally problematic if the transport profile shown in Fig. <ref> were highly sensitive to the potential parameters. Such sensitivity would require extremely fine control over the laser amplitudes, positions, and widths in order to produce the desired transport profile. To test this sensitivity, we generate 2000 perturbed potentials by varying the six parameters of the two-Gaussian solution shown in Fig. <ref> by a random-normal multiplicative factor with mean 1 and standard deviation 0.05. The distributions of the peak positions and heights for the 2000 perturbed potentials are shown in Fig. <ref>. Both the peak positions, k_max, and the peak heights, T_max, undergo perturbations on the scale of the 5%  potential perturbations suggesting that transport properties are relatively insensitive to slight errors in the potential parameters. 

Though we have demonstrated our optimization method in a very simple scenario, it is possible to apply this technique to more complicated scenarios such as a double band-pass filter or step transport profiles. These more complicated transport profiles require more than two Gaussian potentials,
because they rely on multipath interference that must suppress tunneling for certain values of momenta. In our explorations, we found that link potentials made of 3- or 4-Gaussian potentials tended to be more sensitive to random variations of potential parameters. If such potentials are needed to produce the desired transport profile, it may be possible to further augment the cost function in order to preference solutions that are less sensitive to potential perturbations. We leave the thorough exploration of these ideas to future work. In the next section, we discuss a possible experimental application of the NBP filter.








Application. – A Bose or Fermi gas can be placed in the flux region  (see Fig. <ref>) to study quantum environments with neutral, mobile impurities – an important research venue promoted by cold-atom simulators <cit.>. With our proposal, it is possible to investigate the dynamics of impurities that initially have a known momentum profile, thus allowing for a direct measurement of the effective mass and the critical momentum. A detailed discussion of these concepts is beyond the scope of this article. Still, we find it useful to briefly explain them in connection to our proposal. To this end, we consider a degenerate one-dimensional Bose gas with an impurity of momentum P. To model this system, we employ a non-linear Schrödinger equation for a Bose gas with an impurity atom (See Appendix <ref>). This equation was solved analytically in the context of a nonlinear flow past an obstacle <cit.>, which allows us to work out all properties of the dressed particle in a simple manner; note that this (or a similar) non-linear equation was discussed in Refs. <cit.>; see, also, <cit.> for other relevant studies. 

The lowest energy state of the non-linear Schrödinger equation with a given P is a combination of two solitons. They make a dissipationless defect in the Bose gas, which accompanies the impurity. The corresponding lowest energy is given by E≃ E_B+ϵ+P^2/(2m_eff), where E_B is the energy of the gas without an impurity, ϵ is the self-energy of the dressed particle, and m_eff is its effective mass. The solution is stable only for P<P_c; impurities with P>P_c generate grey solitons (cf. <cit.>). Note that quantum fluctuations lead to a finite dissipation (cf. <cit.>) even for P<P_c. We do not consider this effect as it does not change our qualitative presentation.








To measure ϵ, m_eff and P_c, one can use a narrow band-pass filter as shown in Fig. <ref> to create a flux of particles with momenta close to P. 
The width (value of b) of the target profile must be chosen such that the current is weak, i.e., there is a negligible probability to find two flux particles at distances smaller than the healing 
length of the Bose gas. Then the impurity-in-a-gas picture is applicable by construction. 
For simplicity, we assume that initially the impurity is in a hyperfine state 
that does not interact with the Bose gas. To transfer to a strongly interacting hyperfine state 
one has to deposit enough energy to compensate for the interaction effects; see Fig. <ref>a). 
Therefore, the radio-frequency responce (e.g., the transferred fraction) at different momenta directly measures 
the self-energy and effective mass of the polaron.

A measurement such as this would be similar to the measurement in a recent experiment with a three-dimensional Fermi polaron <cit.> but 
with a superior control over the impurity momentum. 
Moreover, the overlap between the non-interacting and interacting states, i.e., the residue, can
be measured, allowing one to test different theoretical methods <cit.>
that, while qualitatively agreeing on m_eff and ϵ, contradict each other on the residue.
Since the momentum of the impurity is known, not only the effective parameters but also the limits of applicability 
of the polaron model will be seen in the radio-frequency responce, in particular, P_c.
In Fig. <ref>b) we present P_c for impurities whose mass M is much larger than m <cit.>.
For weak interactions (g→ 0) the critical momentum is determined by the speed of sound, c, in accordance with the Landau criterion.
In the opposite limit, g→∞, the critical momentum goes to zero as 1/g (cf. <cit.>):   
P_c is limited by the timescale for a two-body exchange. 
 



Summary. –  We outline a procedure for engineering beams of particles with desired momentum profiles using a filter potential connected to a reservoir (see Fig. <ref>). Such a beam can be used to probe cold-atom systems. It can also be used for quantum simulations, as we illustrated with a narrow band-pass filter and a one-dimensional Bose gas in the flux region. Polarons in two-, three- and mixed-dimensional geometries can similarly be created. 






We thank Peter Schlagheck for referring to <cit.>,
and Joachim Brand and Volodymyr Pastukhov for useful discussions.
A. G. V. gratefully acknowledges the support of the Humboldt Foundation and the Deutsche Forschungsgemeinschaft
(VO 2437/1-1).




The authors contributed equally to this work.




apsrev4-1






§ BOUNDARY-AUGMENTED COST FUNCTION

Equation (6) shows the boundary-augmented cost function which includes a term that increases the cost for link-potential solutions that extend beyond the support region x∈[-x_0,x_0]. This added terms is

    J_boundary = α∑_i^N J_boundary^i,

where 

    J_boundary^i = ∫_|x|>x_0dx |V_i(x; A_i,μ_i,σ_i)|^2.

Assuming the Gaussian potential form as in Eq. (5), this evaluates to

    J_boundary^i = √(π)/2A_i^2σ_i[
        (x_0+μ_i/σ_i) +
        (x_0-μ_i/σ_i)
      ],

where  is the complementary error function.




§ IMPURITY IN A BOSE GAS



To model one impurity atom that moves through a one-dimensional environment made of N cold bosonic atoms, we employ the following Hamiltonian

    H=-ħ^2/2m∑_i=1^N∂^2/∂ x_i^2-ħ^2/2M∂^2/∂ y^2+λ∑_i>j=1^Nδ(x_i-x_j)+g ∑_i=1^N δ(x_i-y),

where M is the mass of the impurity atom, and m is the mass of a bosonic particle. The position of the impurity is y, bosons are at the coordinates {x_i}. 
We assume that the realistic boson-boson and boson-impurity interactions are well-described by the zero-range potentials of strengths λ and g, respectively. 
The environment is large by assumption. To describe it, the periodic boundary conditions are used: The particles move in a ring of the circumference L, such that 0<x_i<L and 0<y<L.
We are interested in the thermodynamic limit: N, L→∞ with a fixed value of the density ρ=N/L.


If the system is non-interacting (λ=g=0), the eigenstates are written as e^2π i n_1x_1+…+n_N x_N+my/L, where n_1,…,n_N and m are arbitrary integers. 
For non-vanishing interactions we use these functions to write an eigenfunction of the Hamiltonian as Ψ=∑_{n_j},m a_{n_j},me^2π i ∑ n_jx_j+my/L. 
Because all interactions are pairwise, 
the total (angular) momentum of the system must be conserved, and we write it as P=2πħ/L(∑_j n_j+m).
A conserved quantity (P) allows us to exclude one variable from the consideration. 
We write the function Ψ as Ψ=e^i  P y/ħ∑_{n_j},m a_{n_j},me^2π i ∑ n_jz_j/L≡ e^i  P y/ħψ(z_1,...,z_N) 
with z_i=Lθ(y-x_i)+x_i-y, where θ(x) is the Heaviside step function, i.e., 
θ(x>0) = 1 and zero otherwise. The variables z_i are defined such that 0≤ z_i ≤ L and the impurity is placed at z=0 (z=L). Now if we insert this 
function into the Schrödinger equation, HΨ=EΨ, we obtain the following equation for ψ(0<z_i<L)

    -ħ^2/2m∑_i∂^2 ψ/∂ z_i^2-ħ^2/2M(∑_i∂/∂ z_i)^2ψ
    + i ħ P/M∑_i∂ψ/∂ z_i+λ∑_i>jδ(z_i-z_j)ψ=(E- P^2/2M)ψ,

which must be supplemented with the boundary conditions:

    ψ(z_i=0)=ψ(z_i=L);     ∂ψ/∂ z_i|^z_i=0^+_z_i=L^-= 2 g κ/ħ^2ψ(z_i=0),

where κ=mM/(m+M) is the reduced mass.

By assumption the bosons interact weakly, such that the ansatz ψ=∏_i Φ(z_i) can be used to approximate the system. 
To minimize the expectation value of the Hamiltonian the function Φ(z) must satisfy the following non-linear Schrödinger equation

    -ħ^2/2κ∂^2Φ/∂ z^2+iħ P/M∂Φ/∂ z
    -i ħ^2 (N-1) A/M∂Φ/∂ z + λ (N-1)|Φ|^2Φ=μΦ,

where A=-i∫Φ(x)^*∂/∂ xΦ(x)dx defines the momentum of a boson, and μ is the Lagrange multiplier. 
We rewrite this equation as

    -∂^2Φ/∂ z^2+i v ∂Φ/∂ z + λ̃(N-1)|Φ|^2Φ=μ̃Φ,

where μ̃=2 κμ/ħ^2, λ̃=2 κλ/ħ^2, and 
v= 2 κ P/M ħ with P= P-ħ A(N-1). P defines the momentum of the impurity in the thermodynamic limit;
note that A is determined by P, and there is a unique value of P for a given P.
The boundary conditions for Eq. (<ref>) read

    Φ(z=0)=Φ(z=L);     ∂Φ/∂ z|^z=0^+_z=L^-= g̃Φ(0),

where g̃= 2κ g/ħ^2. The non-linear equation (<ref>) 
has an analytic steady solution <cit.>, which determines the properties 
of the dressed impurity in our problem. Let us first consider the non-interacting impurity g=0. 
In this case the solution for v>0 is <cit.>

    Φ=√(μ̃/λ̃(N-1))(1-βsech^2[√(μ̃β/2)(z+z_0)])^1/2e^iϕ(z),


    ϕ(z)=-πθ(z+z_d)+arctan(√(2 v^2/μ̃β)/exp[√(2μ̃β)(z+z_0)]-2β+1),

where  β=1- v^2/(2μ̃), z_0 is some parameter that determines the origin and z_d is the point where arctan reaches π/2.
 It is worthwhile noting that the solution for v<0 is Φ^*. 
The solution from Eqs. (<ref>) and (<ref>) is plotted in Fig. <ref>; 
for simplicity it is plotted in the interval -L/2<z<L/2, the region 0<z<L easily follows. 



 
 

To describe an interacting impurity, we combine two moving solitons with ± z_0, which creates a singularity at z=0 <cit.>.
Therefore, a dressed impurity in our model is a topological defect with a dissipationless propagation. 
We write the corresponding `wave function' as

    Φ=√(μ̃/λ̃(N-1))(1-βsech^2[√(μ̃β/2)(z± z_0)])^1/2e^iϕ(z),

with 

    ϕ(z)=δϕθ(-z)+arctan(√(2 v^2/μ̃β)/exp[√(2μ̃β)(z± z_0)]-2β+1),

where z_0>0 is discussed below, the parameter δϕ is not important for the further derivations, it reassures that the phase is a continuous function;
the plus sign in ± corresponds to z>0 and the minus sign to z<0. This function is illustrated in Fig. <ref>. 
The density has a non-analytic derivative at z=0. The phase is a continuous function at z=0 (its derivative is also continuous).
Note that the wave function is not periodic (see Eq. (<ref>)). This non-periodicity is not important 
for our discussion, because we are interested in the behavior of the bosons close to the impurity. It suggests that a grey soliton
must be formed upon a change of interaction parameters to take care of the phase slip.

The parameter μ̃ is found from the normalization condition ∫Φ^2=1. For N→∞, we obtain

    μ̃=γρ^2 N-1/N(1-2√(2β_0) (tanh(d)-1)/√(γ)N),

where γ=λ̃/ρ, ρ=N/L, β_0=1-v^2/(2γρ^2), and d=√(γβ_0/2)ρ z_0.
The equation to determine z_0 is found by using the boundary conditions at z={0,L}

    g̃/ρ√(2γ)=β_0^3/2tanh(d)/-β_0+cosh^2(d).

This equation is cubic (in tanh(d)), hence, the solutions can be found in a closed form.
There are three solutions. However, only two will lead to the acceptable values of z_0. 
We will refer to these steady solutions as the `polaron' and the `polaron-soliton' pair, because in the limit g→ 0 the former 
corresponds to the ground state, and the latter to a gray soliton. 
The `polaron-soliton' pair is expected to be unstable (small perturbations lead to 
a decay of this steady solution <cit.>),
therefore, we do not consider it. The solutions merge for z_m 

    tanh^2(√(γβ_0/2)ρ z_m)=√(1+4v^2/γρ^2)-(1+v^2/γρ^2)/2β_0,

which is derived by taking a derivative of Eq. (<ref>) with respect to z_0 and equating the resulting expression to zero – this determines 
the maximum value of g for which (for a fixed β_0) there is a steady solution.
Equations (<ref>) and (<ref>) give the equation for the critical value of v_c:

    g̃/ρ√(γ)=3-√(1+4v_c^2/γρ^2)/-1+√(1+4v_c^2/γρ^2)√(√(1+4v_c^2/γρ)-1-v_c^2/γρ^2).

For v>v_c (see Fig. 4 of the main text) there are no steady solutions.


Now we can calculate the energy of the dressed impurity in the thermodynamic limit

    E≡lim_N→∞, N/L→ρ[E(c,P)-E(c=0,P=0)],

where

    E(c,P) = 𝒫^2/2M + μ N-ħ^2 A^2 N(N-1)/2M-λ N(N-1)∫_0^L/2|Φ|^4dz.

Using these expressions we derive 

    E=P^2/2M+ħ^2ρ^2/2κ√(2  γβ)/3[4 b + (-4b+βsech^2(d))tanh(d)]+ħ P/Mlim_N→∞ A N,

where b=1+v^2/4λ̃ρ=1+κ P^2/2M^2 λρ. This energy for v→ 0 can be written as 

    E≃ϵ+ P^2/2m_eff,

where ϵ is the effective energy of the dressed impurity, and m_eff is the effective mass.

The parameters m_eff and ϵ calculated using Eq. (<ref>)  
agree well with the results in the literature <cit.>, supporting 
the use of the non-linear Schrödinger equation for solving the problem. 
However, further work is required to understand other properties of a dressed impurity in this formalism.
First of all, it will be interesting to investigate the critical momentum,
which is supersonic in the employed model for g→ 0 and for not heavy impurities.
Indeed, the model we solve is equivalent to a heavy impurity moving in a gas of bosons with mass κ,  
which has a different speed of sound. 
Note that Fig. 4 of the main text reports on a heavy impurity (M/m≫ 1) for which this problem does not occur.
It will also be interesting to investigate the residue – the overlap between the wave function that describes 
a state with g=0 and the wave function that describes an interacting state. In the present model, the impurity changes the order parameter
only locally, which means a non-zero residue (see <cit.> for P=0), contradicting other 
studies on the topic <cit.>. To understand this disagreement, one could calculate the overlap using an exactly solvable model, e.g., 
a heavy impenetrable impurity in a Bose gas (solvable by the Bethe ansatz). 

