Engineering momentum profiles of cold-atom beams

TL;DR

This paper presents a method to engineer cold-atom beams with specific momentum profiles by using a filter potential, enabling precise control over atomic beam properties for quantum probing applications.

Contribution

It introduces an algorithm to design filter potentials for tailored momentum profiles in cold-atom beams, including a narrow band-pass filter example.

Findings

Successfully designed a narrow band-pass quantum filter

Demonstrated potential for probing Bose polarons' properties

Provided a framework for custom momentum profile engineering

Abstract

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.

Tables

Table 1. Table 1: The explicit constraints on the potential parameters and the rationale for each constraint. The values of σ min subscript 𝜎 min \sigma_{\mathrm{min}} , σ max subscript 𝜎 max \sigma_{\mathrm{max}} , A min subscript 𝐴 min A_{\mathrm{min}} , and A max subscript 𝐴 max A_{\mathrm{max}} must be determined from the experimental context.

Constraints	Experimental Rationale
${\sum }_{i=1}^N{\mu }_i=0$	The cost function has a continuous degeneracy associated with overall translations of the link potential.
${\sigma }_{\min }\le {\sigma }_j\le {\sigma }_{\max }$	Laser beam widths fall between a minimum and maximum value.
$A_{\min }\le \left|A_j\right|\le A_{\max }$	Laser amplitudes fall between a minimum and maximum value.

Equations

$$-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}\psi +V_0(x)\psi =\frac{{\hslash }^2{k}^2}{2m}\psi ,$$

$$J_{\mathbf{\theta }}=\sum _{0<k<k_F}{w}_k{\left|T_0(k)-T_{\mathbf{\theta }}(k)\right|}^2,$$

$$w(k;r)=w_{\mathrm{bg}}(k)+r{T}_0(k)$$

$$w_{\mathrm{bg}}(k)={\left[\frac{{\sinh }^2(\sqrt{2}\pi k/k_0)}{{\sinh }^2(\sqrt{2}\pi k/k_0)+{\cosh }^2(\sqrt{7}\pi /2)}\right]}^{1/2},$$

$$V_i(x;A_i,{\mu }_i,{\sigma }_i)=\frac{A_i}{\sqrt{2\pi {\sigma }_i^2}}\exp \left[-\frac{(x-{\mu }_i{)}^2}{2{\sigma }_i^2}\right];$$

$$J_{\mathbf{\theta }}^{\mathrm{aug}}=J_{\mathbf{\theta }}+\alpha \sum _{i=1}^N\underset{|x|>x_0}{\int }dx|V_i(x;A_i,{\mu }_i,{\sigma }_i){|}^2,$$

$$T(k;k_0,b)={\left[1+\frac{(k-k_0{)}^2}{b^2}\right]}^{-1}$$

$$J_{\mathrm{boundary}}=\alpha \sum _i^N{J}_{\mathrm{boundary}}^i,$$

$$J_{\mathrm{boundary}}^i=\underset{|x|>x_0}{\int }dx|V_i(x;A_i,{\mu }_i,{\sigma }_i){|}^2.$$

$$J_{\mathrm{boundary}}^i=\frac{\sqrt{\pi }}{2}{A}_i^2{\sigma }_i\left[\mathrm{erfc}\left(\frac{x_0+{\mu }_i}{{\sigma }_i}\right)+\mathrm{erfc}\left(\frac{x_0-{\mu }_i}{{\sigma }_i}\right)\right],$$

$$H=-\frac{{\hslash }^2}{2m}\sum _{i=1}^N\frac{{\partial }^2}{\partial x_i^2}-\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial y^2}+\lambda \sum _{i>j=1}^N\delta (x_i-x_j)+g\sum _{i=1}^N\delta (x_i-y),$$

$$-\frac{{\hslash }^2}{2m}\sum _i\frac{{\partial }^2\psi }{\partial z_i^2}-\frac{{\hslash }^2}{2M}{\left(\sum _i\frac{\partial }{\partial z_i}\right)}^2\psi +i\frac{\hslash \mathcal{P}}{M}\sum _i\frac{\partial \psi }{\partial z_i}+\lambda \sum _{i>j}\delta (z_i-z_j)\psi =\left(E-\frac{{\mathcal{P}}^2}{2M}\right)\psi ,$$

$$\psi(z_{i}=0)=\psi(z_{i}=L);\qquad\frac{\partial\psi}{\partial z_{i}}\bigg{|}^{z_{i}=0^{+}}_{z_{i}=L^{-}}=\frac{2g\kappa}{\hbar^{2}}\psi(z_{i}=0),$$

$$-\frac{{\hslash }^2}{2\kappa }\frac{{\partial }^2\Phi }{\partial z^2}+i\frac{\hslash \mathcal{P}}{M}\frac{\partial \Phi }{\partial z}-i\frac{{\hslash }^2(N-1)A}{M}\frac{\partial \Phi }{\partial z}+\lambda (N-1)|\Phi {|}^2\Phi =\mu \Phi ,$$

$$-\frac{{\partial }^2\Phi }{\partial z^2}+iv\frac{\partial \Phi }{\partial z}+\tilde{\lambda }(N-1)|\Phi {|}^2\Phi =\tilde{\mu }\Phi ,$$

$$\Phi(z=0)=\Phi(z=L);\qquad\frac{\partial\Phi}{\partial z}\bigg{|}^{z=0^{+}}_{z=L^{-}}=\tilde{g}\Phi(0),$$

$$\Phi =\sqrt{\frac{\tilde{\mu }}{\tilde{\lambda }(N-1)}}{\left(1-\beta {\mathrm{sech}}^2\left[\sqrt{\frac{\tilde{\mu }\beta }{2}}(z+z_0)\right]\right)}^{\frac{1}{2}}{e}^{i\varphi (z)},$$

$$\varphi (z)=-\pi \theta (z+z_d)+\arctan \left(\frac{\sqrt{\frac{2v^2}{\tilde{\mu }}\beta }}{\exp \left[\sqrt{2\tilde{\mu }\beta }(z+z_0)\right]-2\beta +1}\right),$$

$$\Phi =\sqrt{\frac{\tilde{\mu }}{\tilde{\lambda }(N-1)}}{\left(1-\beta {\mathrm{sech}}^2\left[\sqrt{\frac{\tilde{\mu }\beta }{2}}(z\pm z_0)\right]\right)}^{\frac{1}{2}}{e}^{i\varphi (z)},$$

$$\varphi (z)=\delta \varphi \theta (-z)+\arctan \left(\frac{\sqrt{\frac{2v^2}{\tilde{\mu }}\beta }}{\exp \left[\sqrt{2\tilde{\mu }\beta }(z\pm z_0)\right]-2\beta +1}\right),$$

$$\tilde{\mu }=\gamma {\rho }^2\frac{N-1}{N}\left(1-2\sqrt{2{\beta }_0}\frac{(\tanh (d)-1)}{\sqrt{\gamma }N}\right),$$

$$\frac{\tilde{g}}{\rho \sqrt{2\gamma }}=\frac{{\beta }_0^{\frac{3}{2}}\tanh (d)}{-{\beta }_0+{\cosh }^2(d)}.$$

$${\tanh }^2\left(\sqrt{\frac{\gamma {\beta }_0}{2}}\rho z_m\right)=\frac{\sqrt{1+\frac{4v^2}{\gamma {\rho }^2}}-(1+\frac{v^2}{\gamma {\rho }^2})}{2{\beta }_0},$$

$$\frac{\tilde{g}}{\rho \sqrt{\gamma }}=\frac{3-\sqrt{1+\frac{4v_c^2}{\gamma {\rho }^2}}}{-1+\sqrt{1+\frac{4v_c^2}{\gamma {\rho }^2}}}\sqrt{\sqrt{1+\frac{4v_c^2}{\gamma \rho }}-1-\frac{v_c^2}{\gamma {\rho }^2}}.$$

$$\mathcal{E}\equiv \underset{N\to \infty ,\frac{N}{L}\to \rho }{\lim }\left[E(c,P)-E(c=0,P=0)\right],$$

$$E(c,P)=\frac{{\mathcal{P}}^2}{2M}+\mu N-\frac{{\hslash }^2{A}^{2}N(N-1)}{2M}-\lambda N(N-1){\int }_0^{L/2}|\Phi {|}^4\mathrm{d}z.$$

$$\mathcal{E}=\frac{P^2}{2M}+\frac{{\hslash }^2{\rho }^2}{2\kappa }\frac{\sqrt{2\gamma \beta }}{3}\left[4b+(-4b+\beta {\mathrm{sech}}^2(d))\tanh (d)\right]+\frac{\hslash P}{M}\underset{N\to \infty }{\lim }AN,$$

$$\mathcal{E}\simeq ϵ+\frac{P^2}{2m_{\mathrm{eff}}},$$

Full text

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 Braaten and Hammer (2006); Bloch et al. (2008) or the self-energy and the effective mass of a polaron Massignan et al. (2014); Schmidt et al. (2018). The beams can also be used as the initial state for atom interferometry Impens et al. (2006); Cronin et al. (2009).

Figure  1 summarizes our proposal. Analogous to a quantum switch device (‘transistor’) Micheli et al. (2004); Marchukov et al. (2016); Thuberg et al. (2017), 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 111Strong 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 Paul et al. (2005)., so that their scattering properties may be calculated using the one-body Schrödinger equation:

[TABLE]

where $m$ is the mass of a particle from the reservoir, $\frac{\hbar^{2}k^{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 Krinner et al. (2015); Lebrat et al. (2018, 2019), 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., Bergschneider et al. (2018)), 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 Guerin et al. (2006); Gattobigio et al. (2009); Büning et al. (2010) or on controlling the transmission of atoms with careful manipulation of an optical lattice Carusotto and La Rocca (2000); Fabre et al. (2011); Damon et al. (2015). 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)\equiv V(x,\bm{\theta}^{*})$ by performing a global search over a family of possible potentials $V(x,\bm{\theta})$ for the parameters $\bm{\theta}^{*}$ that reduce the $k$-integrated squared error between the desired transmission-momentum profile $T_{0}(k)$ and the actual profile $T_{\bm{\theta}}(k)$ produced by a sample potential $V(x,\bm{\theta})$. Concretely, we minimize the cost

[TABLE]

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. 1 and 2), 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\neq 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:

[TABLE]

where $w_{\mathrm{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_{\mathrm{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 $\hbar^{2}k_{0}^{2}/[2m\cosh^{2}(k_{0}x/\sqrt{2})]$ (cf. Landau and Lifshitz (1977)):

[TABLE]

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,\bm{\theta})$, as a sum of $N$ Gaussians, each of the form

[TABLE]

the parameter $\bm{\theta}$ denotes the parameter space $\{A_{1},\mu_{1},\sigma_{1},...,A_{N},\mu_{N},\sigma_{N}\}$. While minimizing Eq. (2), we enforce the parameter constraints listed in Table 1. 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\in[-x_{0},x_{0}]$. To accomplish this, we minimize the boundary-augmented cost function

[TABLE]

where $\alpha$ 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_{\bm{\theta}}\simeq J_{\bm{\theta}}+\alpha A_{i}$, which means that $\alpha$ should be chosen of the order $J_{\bm{\theta}}/A_{\mathrm{min}}$ to dictate the constraint. The integral in Eq. (6) evaluates to the complementary error function. The full form of $J_{\bm{\theta}}^{\mathrm{aug}}$ is given in Appendix A.

For a particular choice of $\bm{\theta}$ (and hence $V(x;\bm{\theta})$), we solve for $T_{\bm{\theta}}(k)$ by integrating the Schrödinger equation (1) 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. (1), which transforms Eq. (1) 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_{\bm{\theta}}^{\mathrm{aug}}$ for $\bm{\theta}^{*}$ using the global optimization routine called Differential Evolution (DE) Storn and Price (1997). 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 2)

[TABLE]

where $k_{0}$ determines the peak position and $b$ determines the width. For simplicity, in this subsection we adopt the units $k_{0}=\hbar=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.03k_{0}$, which for a reasonable assumption $\hbar^{2}k_{0}^{2}/(2m)=k_{B}\times\mu$K, where $k_{B}$ is the Boltzmann constant leads to the time scale associated with the resonance width $\frac{2m}{\hbar b^{2}}\sim 10$ms.

For the constraints shown in Table 1, we use $\sigma_{\mathrm{min}}=0.2$, $\sigma_{\mathrm{max}}=3$, $A_{\mathrm{min}}=5$, and $A_{\mathrm{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 2 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 222For 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_{\bm{\theta}}(k)$ in Eq. (2) 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. (2) 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. 2 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. 2 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. 3. Both the peak positions, $k_{\mathrm{max}}$, and the peak heights, $T_{\mathrm{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. 1) to study quantum environments with neutral, mobile impurities – an important research venue promoted by cold-atom simulators Schirotzek et al. (2009); Nascimbène et al. (2009); Kohstall et al. (2012); Spethmann et al. (2012); Catani et al. (2012); Fukuhara et al. (2013); Hu et al. (2016); Jørgensen et al. (2016); Scazza et al. (2017); Peña Ardila et al. (2019). 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 B). This equation was solved analytically in the context of a nonlinear flow past an obstacle Hakim (1997), 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. Schecter et al. (2016); Volosniev and Hammer (2017); Mistakidis et al. (2019); Dehkharghani et al. (2018); Pastukhov (2019); Panochko and Pastukhov (2019); see, also, Sacha and Timmermans (2006); Catani et al. (2012); Kain and Ling (2016); Parisi and Giorgini (2017); Grusdt et al. (2017); Pastukhov (2017); Kain and Ling (2018) 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\simeq E_{B}+\epsilon+P^{2}/(2m_{\mathrm{eff}})$, where $E_{B}$ is the energy of the gas without an impurity, $\epsilon$ is the self-energy of the dressed particle, and $m_{\mathrm{eff}}$ is its effective mass. The solution is stable only for $P<P_{c}$; impurities with $P>P_{c}$ generate grey solitons (cf. Hakim (1997)). Note that quantum fluctuations lead to a finite dissipation (cf. Astrakharchik and Pitaevskii (2004); Sykes et al. (2009); Cherny et al. (2012)) even for $P<P_{c}$. We do not consider this effect as it does not change our qualitative presentation.

To measure $\epsilon,m_{\mathrm{eff}}$ and $P_{c}$, one can use a narrow band-pass filter as shown in Fig. 2 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. 4a). 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 Scazza et al. (2017) 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 Volosniev and Hammer (2017); Pastukhov (2017); Grusdt et al. (2017); Pastukhov (2019) that, while qualitatively agreeing on $m_{\mathrm{eff}}$ and $\epsilon$, 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. 4b) we present $P_{c}$ for impurities whose mass $M$ is much larger than $m$ Hakim (1997). For weak interactions ($g\to 0$) the critical momentum is determined by the speed of sound, $c$, in accordance with the Landau criterion. In the opposite limit, $g\to\infty$, the critical momentum goes to zero as $1/g$ (cf. Schecter et al. (2016)): $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. 1). 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.

The authors contributed equally to this work.

Appendix A 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\in[-x_{0},x_{0}]$. This added terms is

[TABLE]

where

[TABLE]

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

[TABLE]

where $\operatorname{erfc}$ is the complementary error function.

Appendix B 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

[TABLE]

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 $\lambda$ 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\to\infty$ with a fixed value of the density $\rho=\frac{N}{L}$.

If the system is non-interacting ($\lambda=g=0$), the eigenstates are written as $e^{2\pi i\frac{n_{1}x_{1}+\ldots+n_{N}x_{N}+my}{L}}$, where $n_{1},\ldots,n_{N}$ and $m$ are arbitrary integers. For non-vanishing interactions we use these functions to write an eigenfunction of the Hamiltonian as $\Psi=\sum_{\{n_{j}\},m}a_{\{n_{j}\},m}e^{2\pi i\frac{\sum n_{j}x_{j}+my}{L}}$. Because all interactions are pairwise, the total (angular) momentum of the system must be conserved, and we write it as ${\cal P}=\frac{2\pi\hbar}{L}\left(\sum_{j}n_{j}+m\right)$. A conserved quantity (${\cal P}$) allows us to exclude one variable from the consideration. We write the function $\Psi$ as $\Psi=e^{i\frac{{\cal P}y}{\hbar}}\sum_{\{n_{j}\},m}a_{\{n_{j}\},m}e^{2\pi i\frac{\sum n_{j}z_{j}}{L}}\equiv e^{i\frac{{\cal P}y}{\hbar}}\psi(z_{1},...,z_{N})$ with $z_{i}=L\theta(y-x_{i})+x_{i}-y$, where $\theta(x)$ is the Heaviside step function, i.e., $\theta(x>0)=1$ and zero otherwise. The variables $z_{i}$ are defined such that $0\leq z_{i}\leq L$ and the impurity is placed at $z=0$ ($z=L$). Now if we insert this function into the Schrödinger equation, $H\Psi=E\Psi$, we obtain the following equation for $\psi(0<z_{i}<L)$

[TABLE]

which must be supplemented with the boundary conditions:

[TABLE]

where $\kappa=mM/(m+M)$ is the reduced mass.

By assumption the bosons interact weakly, such that the ansatz $\psi=\prod_{i}\Phi(z_{i})$ can be used to approximate the system. To minimize the expectation value of the Hamiltonian the function $\Phi(z)$ must satisfy the following non-linear Schrödinger equation

[TABLE]

where $A=-i\int\Phi(x)^{*}\frac{\partial}{\partial x}\Phi(x)\mathrm{d}x$ defines the momentum of a boson, and $\mu$ is the Lagrange multiplier. We rewrite this equation as

[TABLE]

where $\tilde{\mu}=\frac{2\kappa\mu}{\hbar^{2}}$, $\tilde{\lambda}=\frac{2\kappa\lambda}{\hbar^{2}}$, and $v=\frac{2\kappa P}{M\hbar}$ with $P={\cal P}-\hbar 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 ${\cal P}$ for a given $P$. The boundary conditions for Eq. (15) read

[TABLE]

where $\tilde{g}=\frac{2\kappa g}{\hbar^{2}}$. The non-linear equation (15) has an analytic steady solution Hakim (1997), 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 Tsuzuki (1971); Ishikawa and Takayama (1980)

[TABLE]

[TABLE]

where $\beta=1-v^{2}/(2\tilde{\mu})$, $z_{0}$ is some parameter that determines the origin and $z_{d}$ is the point where $\mathrm{arctan}$ reaches $\pi/2$. It is worthwhile noting that the solution for $v<0$ is $\Phi^{*}$. The solution from Eqs. (17) and (18) is plotted in Fig. 5; 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 $\pm z_{0}$, which creates a singularity at $z=0$ Hakim (1997); Schecter et al. (2016). Therefore, a dressed impurity in our model is a topological defect with a dissipationless propagation. We write the corresponding ‘wave function’ as

[TABLE]

with

[TABLE]

where $z_{0}>0$ is discussed below, the parameter $\delta\phi$ is not important for the further derivations, it reassures that the phase is a continuous function; the plus sign in $\pm$ corresponds to $z>0$ and the minus sign to $z<0$. This function is illustrated in Fig. 5. 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. (20)). 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 $\tilde{\mu}$ is found from the normalization condition $\int\Phi^{2}=1$. For $N\to\infty$, we obtain

[TABLE]

where $\gamma=\tilde{\lambda}/\rho$, $\rho=N/L$, $\beta_{0}=1-v^{2}/(2\gamma\rho^{2})$, and $d=\sqrt{\frac{\gamma\beta_{0}}{2}}\rho z_{0}$. The equation to determine $z_{0}$ is found by using the boundary conditions at $z=\{0,L\}$

[TABLE]

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\to 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 Hakim (1997)), therefore, we do not consider it. The solutions merge for $z_{m}$

[TABLE]

which is derived by taking a derivative of Eq. (22) with respect to $z_{0}$ and equating the resulting expression to zero – this determines the maximum value of $g$ for which (for a fixed $\beta_{0}$) there is a steady solution. Equations (22) and (23) give the equation for the critical value of $v_{c}$:

[TABLE]

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

[TABLE]

where

[TABLE]

Using these expressions we derive

[TABLE]

where $b=1+\frac{v^{2}}{4\tilde{\lambda}\rho}=1+\frac{\kappa P^{2}}{2M^{2}\lambda\rho}$. This energy for $v\to 0$ can be written as

[TABLE]

where $\epsilon$ is the effective energy of the dressed impurity, and $m_{\mathrm{eff}}$ is the effective mass.

The parameters $m_{\mathrm{eff}}$ and $\epsilon$ calculated using Eq. (28) agree well with the results in the literature Mistakidis et al. (2019), 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\to 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 $\kappa$, which has a different speed of sound. Note that Fig. 4 of the main text reports on a heavy impurity ($M/m\gg 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 Volosniev and Hammer (2017) for $P=0$), contradicting other studies on the topic Pastukhov (2017); Grusdt et al. (2017). 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).