Polaron-Depleton Transition in the Yrast Excitations of a One-Dimensional Bose Gas with a Mobile Impurity

Abstract

We present exact numerical data for the lowest-energy momentum eigenstates (yrast states) of a repulsive spin impurity in a one-dimensional Bose gas using full configuration interaction quantum Monte Carlo (FCIQMC). As a stochastic extension of exact diagonalization, it is well suited for the study of yrast states of a lattice-renormalized model for a quantum gas. Yrast states carry valuable information about the dynamic properties of slow-moving mobile impurities immersed in a many-body system. Based on the energies and the first and second-order correlation functions of yrast states, we identify different dynamical regimes and the transitions between them: The polaron regime, where the impurity’s motion is affected by the Bose gas through a renormalized effective mass; a regime of a gray soliton that is weakly correlated with a stationary impurity, and the depleton regime, where the impurity occupies a dark or gray soliton. Extracting the depleton effective mass reveals a super heavy regime where the magnitude of the (negative) depleton mass exceeds the mass of the finite Bose gas.

1. Introduction

The study of a single quantum impurity in a surrounding many-body medium has fascinated scientists for many decades [1,2]. Beyond the historical interest around the influence of the crystal lattice on the motion of an electron—the original “polaron” [3], or impurity atoms in superfluid helium [4]—there has recently been a surge of interest in the field of ultracold atoms, where interactions can be readily tuned with the help of Feshbach resonances [5] and excitation spectra probed with spectroscopic methods [6]. A particular focus of experimental scrutiny has been the Bose polaron, where an impurity atom is coupled with a bosonic bath [7,8,9,10].

Restricting the dimensionality to one spatial dimension provides access to the special physics of one-dimensional quantum liquids [11,12], where impurities have been predicted to undergo Bloch oscillations [13,14]: Due to the periodicity of the dispersion relation, an impurity experiencing a weak force periodically alters its excitation state without contributing to transport in real space, as originally predicted [15] and later observed [16,17] for a particle in an external lattice potential. The prediction of Bloch oscillations in a one-dimensional quantum liquid, even in the absence of a periodic potential [13,14], was debated [18,19] but eventually confirmed in an experiment with spin impurities in a one-dimensional gas of cesium atoms [20]. Other experiments probing impurity physics in one-dimensional quantum gases also employed spin impurities (where the impurity atoms have the same mass and only differ in a spin quantum number) [21,22], or different types of atoms [23,24].

While experimental studies have started, the quantitative understanding of impurity physics in a one-dimensional Bose gas is far from complete. The different theoretical approaches to the problem range from mean-field theory [25,26] and related variational theory [27,28,29] via the path-integral approach [30,31], the renormalization group [32,33,34] and flow equation [35] method to multiconfiguration time-dependent Hartree [36] and quantum Monte Carlo methods [37,38,39,40].

Many of the mentioned works focus on ground-state properties and effective mass of the one-dimensional Bose polaron with either zero or a very small total momentum. Although the dynamics of impurities have also been actively studied [28,29,36,41], limited understanding has been achieved on the full dispersion relation of a Bose gas coupled with a mobile impurity. There are analytical results on the dispersion relations restricted to specific models, such as the Yang–Gaudin model [42], and the Luttinger liquid [43].

In a homogeneous one-dimensional gas, e.g., in a ring geometry with periodic boundaries, translational invariance makes the total momentum a good quantum number. This allows for the study of yrast states, which are eigenstates with the lowest energy at given momentum. Yrast states are stable as long as momentum is conserved, while adiabatic passage through the yrast states of different momentums is responsible for the Bloch oscillation phenomena of Refs. [13,20]. The yrast states of a bosonic superfluid in the absence of impurities are intimately connected [44,45,46,47,48,49,50,51] to localized nonlinear waves, known as dark solitons [52]. Dark solitons are ubiquitous features of superfluids, which can be characterized by a localized density depression and a phase jump [51,53,54].

When a repulsive impurity is introduced into the Bose gas, two different low-energy configurations can exist depending on the momentum: At a lower momentum, the impurity moves relative to the quantum gas forming a polaron. At higher momentum, the bulk of the momentum is taken by the Bose gas forming a gray or dark soliton modified by the presence of the impurity. This situation was named the “depleton” in Refs. [14,19].

For strongly correlated impurity problems that are outside the reach of analytically solvable models, quantum Monte Carlo (QMC) methods have proven to be invaluable tools [37,38,39,40,55]. In this work, we employ the full configuration interaction quantum Monte Carlo (FCIQMC) [56,57] method. It can be seen as a natural stochastic extension to the exact diagonalization method, which allows one to treat a larger Hilbert space that could otherwise not fit into the computer memory. While different in detail, it is similar in spirit to earlier versions of projector Monte Carlo methods [58] in sampling the ground state wave function. When applied to a translationally invariant Hamiltonian in momentum space, FCIQMC has the advantage over other QMC methods, such as diffusion Monte Carlo (which is formulated in real space) or auxiliary field QMC, that momentum is strictly conserved in each elementary stochastic operation. Thus yrast states can be obtained easily by projection onto the lowest-energy state within a total-momentum sector starting from an initial state with the same momentum. Moreover, FCIQMC mitigates the sign problem by walker annihilation in many systems when a sufficient number of walkers is present [59]. Additionally, the initiator approximation [57] can be applied to suppress the sign problem with the trade-off that a small initiator bias is introduced.

The FCIQMC method was originally developed for fermionic many-body problems. It has been applied to the electronic structure of molecules and solids [60,61,62] and the Hubbard model [63,64,65]. Recently, it was used to study the yrast states in a superfluid of spin-$\frac{1}{2}$ fermions [66]. Here, we use FCIQMC for the first time to quantitatively study the physics of a bosonic many-body problem, while previously, bosonic Hamiltonians were employed when developing and analyzing the FCIQMC procedures [67,68].

In this work, we use the FCIQMC method to obtain numerical results for the yrast states of a one-dimensional Bose gas containing a repulsive spin impurity. We characterize the polaron and depleton regimes of the yrast dispersion, as well as the transitions between them, by examining the energies and the first and second-order correlation functions of yrast states. The extracted depleton effective mass reveals a super-heavy regime where the magnitude of the (negative) depleton mass exceeds the mass of the finite Bose gas. The results also show that the depleton picture becomes inadequate for smaller impurity–boson interactions where the impurity and Bose-gas motion decouples.

In Section 2, we introduce the lattice discretized model Hamiltonian with renormalized coupling constants to be studied in this work. The FCIQMC algorithm and modifications made for treating bosonic Fock states, along with implementation details and parameter choices, are discussed in Section 3. The numerical results on the properties of the yrast states and their interpretation in terms of the different physical regimes are presented in Section 4, starting with the yrast dispersion in Section 4.1. This is followed by the impurity momentum and the impurity–boson correlation function in Section 4.2 and Section 4.3, respectively. The effective mass in the soliton/depleton regions of the dispersion is reported in Section 4.4, and the spin-flip energies in Section 4.5. Finally, we draw conclusions in Section 5, and we outline the possible future prospects of this work. Appendix A presents data on the elimination of systematic biases in FCIQMC, which is relevant for validating the computational method.

2. The Model

2.1. The Hamiltonian in One-Dimensional Real Space

We consider a single impurity particle immersed in a one-dimensional interacting Bose gas of N identical particles. The Hamiltonian reads

$$H=-\frac{{\hslash }^2}{2m}{\displaystyle \sum }_{i=1}^N\frac{{\partial }^2}{\partial x_i^2}+g_{\mathrm{BB}}{\displaystyle \sum }_{i<j}\delta (x_i-x_j)-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x_{\mathrm{I}}^2}+g_{\mathrm{IB}}{\displaystyle \sum }_{i=1}^N\delta (x_i-x_{\mathrm{I}}),$$

(1)

where $x_i(i=1,\cdots ,N)$ and $x_{\mathrm{I}}$ are the coordinates of the bosons and the impurity, respectively. We have already assumed that the impurity has the same mass m as the bosons do and will continue to do so throughout this work. This is adequate for a spin impurity where the impurity atom becomes distinguishable from the remaining bosons by changing a spin quantum number, e.g., changing a hyperfine quantum number of an ultracold atom. The N bosons are interacting with a contact potential of strength $g_{\mathrm{BB}}$ while the interaction of the impurity particle with the bosons is described by a contact potential of strength $g_{\mathrm{IB}}$. We consider repulsive interactions with $g_{\mathrm{BB}}>0$ and $g_{\mathrm{IB}}\ge 0$ in order to access the physics of the polaron–depleton transition. We leave the detailed study of attractive impurities, which bind to bosons rather than to the hole-like dark soliton excitations, to future work.

Following the definitions from previous works on polaron problems [26,37,40], we introduce the dimensionless coupling parameters

$$\gamma =\frac{m{g}_{\mathrm{BB}}}{{\hslash }^{2}n},\eta =\frac{m{g}_{\mathrm{IB}}}{{\hslash }^{2}n},$$

(2)

to represent the boson–boson and boson–impurity interaction strengths, respectively. The density of the Bose gas is $n=N/L$. Note that with the impurity and the bosons having equal mass, the reduced mass, $m_r$, used in other works, becomes $m/2$.

2.2. Lattice Discretized Continuum Model

For our numerical simulations, we consider a finite system in a one-dimensional box of length L with periodic boundary conditions. We discretize the model using a lattice with M lattice sites with renormalized contact interactions [69,70]. In order to access the yrast states numerically with FCIQMC, we use a momentum-space representation of the Hamiltonian. In FCIQMC, individual stochastic sampling steps then conserve momentum, which allows us to access the yrast states with this projector QMC method. The spatial domain is $x\in \left(-L/2,L/2\right]$ and the lattice constant is defined as $\alpha =L/M$ for M lattice points. In this representation, the Hamiltonian reads

$$H^{\mathrm{mom}}={\displaystyle \sum }_k{ϵ}_k{\widehat{a}}_k^{†}{\widehat{a}}_k+{\displaystyle \sum }_k{ϵ}_k{\widehat{b}}_k^{†}{\widehat{b}}_k+\frac{U}{2M}{\displaystyle \sum }_{spqr}{\widehat{a}}_s^{†}{\widehat{a}}_p^{†}{\widehat{a}}_q{\widehat{a}}_r{\delta }_{s+p,q+r}+\frac{V}{M}{\displaystyle \sum }_{spqr}{\widehat{a}}_s^{†}{\widehat{b}}_p^{†}{\widehat{b}}_q{\widehat{a}}_r{\delta }_{s+p,q+r},$$

(3)

where ${\widehat{a}}_k^{†}\left({\widehat{a}}_k\right)$ are the boson creation (annihilation) operators; the corresponding operators for the impurity are ${\widehat{b}}_k^{†}\left({\widehat{b}}_k\right)$. The plane-wave eigenstates $\langle x|{\widehat{a}}_k^{†}|\mathrm{vac}\rangle =e^{-ikx/\alpha }$ of momentum $\hslash k\alpha$ are indexed with the dimensionless quantum numbers

$$k_j=\left\{\begin{array}{cc}-\pi +j\frac{2\pi }{M}\hfill & \mathrm{if}M\mathrm{is}\mathrm{even}\hfill \\ -\pi \frac{M+1}{M}+j\frac{2\pi }{M}\hfill & \mathrm{if}M\mathrm{is}\mathrm{odd}\hfill \end{array}\right.$$

(4)

where $j\in \{1,2,\dots ,M\}$ is an integer. The kinetic energy dispersion is the same for bosons and impurity (as they have equal mass)

$${ϵ}_k=\frac{{\hslash }^2{k}^2}{2m{\alpha }^2}=\frac{1}{2}{M}^2{k}^2{\epsilon }_0,$$

(5)

where we have introduced the unit of energy that will be used throughout this work

$${\epsilon }_0=\frac{{\hslash }^2}{m{L}^2}.$$

(6)

Note that Refs. [26,40] choose the Fermi energy ${\epsilon }_F=\frac{{\pi }^2{\hslash }^2{n}^2}{2m}=\frac{{\pi }^2{N}^2}{2}{\epsilon }_0$ as the energy unit.

The parameters U and V are the lattice on-site interaction strengths for boson–boson and boson–impurity, respectively. They are renormalized to generate the correct scattering length for a two-particle scattering problem at zero energy [69,70]:

$$U\alpha =\frac{g_{\mathrm{BB}}}{1+\frac{g_{\mathrm{BB}}}{g_0}},V\alpha =\frac{g_{\mathrm{IB}}}{1+\frac{g_{\mathrm{IB}}}{g_0}},$$

(7)

where $g_0={\pi }^2{\hslash }^2/m\alpha$.

2.3. Connection to Mean-Field Theory and Choice of Parameters

In the weakly interacting regime where $\gamma \ll 1$, nonlinear phenomena in the Bose gas, such as dark and gray solitons, are accurately described by the Gross–Pitaevskii equation [51,52]. A similar mean-field treatment is also available for a Bose gas with an impurity [19,32]. The relevant length scale in this theory is the healing length, which is the shortest length scale on which the superfluid order parameter can change

$$l_{\mathrm{h}}=\frac{\hslash }{\sqrt{2m{g}_{\mathrm{BB}}n}}=\frac{L}{\sqrt{2\gamma }N}.$$

(8)

In order to obtain insights into the physics of solitons and their interaction with impurities, we need to choose the parameters of our model system such that $L\gg l_{\mathrm{h}}$. In order to obtain results relevant for the thermodynamic limit, it would be desirable to choose both the particle number N and the box size L to be large. However, we are constrained by the fact that FCIQMC has a sign problem, which limits the maximum number of particles and modes that can be accurately computed. The sign problem also grows more severe with stronger interaction strength, which affects the off-diagonal matrix elements in the Hamiltonian of Equation (3). As a compromise, we choose to work with $N_{\mathrm{tot}}=20$ particles and $M=50$ modes and explore the weakly interacting regime. This yields a ratio of box size to healing length of $L/l_{\mathrm{h}}\equiv \sqrt{2\gamma }N\approx 12$ for $\gamma =0.2$ and $L/l_{\mathrm{h}}\approx 4$ for $\gamma =0.02$ for the two values of the Bose-gas interaction strength that we are using in this work.

3. Computational Method and Simulation Details

Full configuration interaction quantum Monte Carlo is a projector quantum Monte Carlo method that can be used to determine the ground-state energies of quantum many-body systems. It was originally formulated to solve problems in quantum chemistry [56]. In this section, we describe the algorithm and some of the modifications we implemented to extend the algorithm to describe bosonic systems.

3.1. Bosonic Full Configuration Quantum Monte Carlo

In FCIQMC, a basis of Fock states (occupation number basis) for $N={{\displaystyle \sum }}_{i=1}^M{n}_i$ particles in M lattice sites is used

$$|n_1,n_2,\dots n_M\rangle =\prod _{i=1}^M\frac{1}{\sqrt{n_i!}}{\left({\widehat{a}}_i^{†}\right)}^{n_i}|\mathrm{vac}\rangle .$$

(9)

Within this basis, the Hamiltonian is represented as a matrix $\mathbf{H}$ and the quantum state (many-body wave function) as a vector $\mathbf{c}$ containing the signed weights of the individual Fock states as coefficients.

The ground-state coefficient vector is then found in an iterative manner by repeatedly applying the equation

$${\mathbf{c}}^{(n+1)}={\mathbf{c}}^{\left(n\right)}+\delta \tau \left(\mathbf{1}{S}^{\left(n\right)}-\mathbf{H}\right){\mathbf{c}}^{\left(n\right)},$$

(10)

where the parameter $\delta \tau$ controls the size of the time step, ${\mathbf{c}}^{\left(n\right)}$ is the approximation of the eigenvector at the n-th time step, and $\mathbf{1}$ is the identity matrix. The shift $S^{\left(n\right)}$ is a real number used to keep the norm of ${\mathbf{c}}^{\left(n\right)}$ under control. It is adjusted by the following scheme

$$S^{(n+1)}=S^{\left(n\right)}-\frac{\zeta }{\delta \tau }ln\left(\frac{N_{\mathrm{w}}^{(n+1)}}{N_{\mathrm{w}}^{\left(n\right)}}\right)-\frac{\xi }{\delta \tau }ln\left(\frac{N_{\mathrm{w}}^{(n+1)}}{N_{\mathrm{t}}}\right),$$

(11)

where $N_{\mathrm{w}}^{\left(n\right)}\equiv {\parallel {\mathbf{c}}^{\left(n\right)}\parallel }_1$ is the 1-norm of ${\mathbf{c}}^{\left(n\right)}$, $N_{\mathrm{t}}$ the parameter for the target norm, and $\xi$ and $\zeta$ are parameters that control the dynamics of the shift. In the steady-state, the instantaneous norm $N_{\mathrm{w}}^{\left(n\right)}$ fluctuates around the value of $N_{\mathrm{t}}$ [67]. It is important to control the vector norm in FCIQMC, as it is a proxy for the number of (stored) non-zero elements of the coefficient vector and thus for both memory and runtime requirements of the simulation.

Because the size of the Hilbert space grows exponentially with system size, both $\mathbf{H}$ and $\mathbf{c}$ quickly become prohibitively large. To get around this problem, we replace the matrix-vector multiplication in Equation (10) with a stochastic sampling process. The sampling process is designed to reproduce the right-hand side of Equation (10) by the expected value while at the same time replacing most coefficients in the vector with zero, such that the values do not have to be stored. Concretely, we divide the values of the entries in $\mathbf{c}$ into integer units called “walkers”. At each time step, each walker attempts to “spawn” to a configuration connected by a non-zero entry in the corresponding column of the matrix $\mathbf{H}$.

The spawning from configuration q to configuration $r\ne q$ can be described as

$$c_r\leftarrow c_r-\frac{\delta \tau }{p_{\mathrm{spawn}}}{H}_{r,q}{c}_q,$$

(12)

where $\frac{1}{p_{\mathrm{spawn}}}$ is the inverse probability of picking r, i.e., the number of non-zero off-diagonal entries in the q-th column of $\mathbf{H}$. If the occupation number ${\mathbf{c}}_q$ is greater than the number of non-zero entries in this column, the spawns can be performed exactly. In addition to the off-diagonal spawns, the diagonal part of the matrix-vector multiplication in Equation (10) is performed exactly. After a step is complete, we stochastically project the entries $v_i$ of the vector $\mathbf{c}$ to a threshold t; values $|v_i|<t$ are removed from the vector with probability $p=1-\frac{|v_i|}{t}$. Otherwise, their value becomes $v_i=t$. In practice, we usually set $t=1$.

By using this approach, the length of the vectors ${\mathbf{c}}^{\left(n\right)}$ can be much smaller than the dimension of the Hilbert space while the expectation value of ${\mathbf{c}}^{\left(n\right)}$ still approaches the exact eigenvector of the ground state of $\mathbf{H}$. At the same time, the shift $S^{\left(n\right)}$ equilibrates to fluctuating around the ground state eigenvalue with a small stochastic bias [68,71]. The spawning process described above differs from the original one of Ref. [56] and is similar in spirit but more efficient than the modifications discussed in Refs. [72,73,74]. It will be described in greater detail elsewhere [75].

While Equation (3) defines the Hamiltonian used in this study, the FCIQMC method is completely agnostic to the nature of the Hamiltonian, as long as it results in a sparse matrix where elements can be computed efficiently on the fly. Thus, it is possible to study multi-dimensional models, long-range interactions, or even complex-valued problems [62].

3.2. Implementation Details

We have implemented the FCIQMC algorithm in the high-level and high-performance programming language Julia [76]. Both the FCIQMC algorithm and all analysis tools are implemented as a library. This way, the calculations and all parameters can be defined in a concise script, written in the same language as the library, without the need for input files in a different format. The set-up is very flexible and makes it easy to experiment interactively with immediate visualization of data, e.g., in a notebook interface, or deploy code to a high-performance computer. The library code Rimu.jl used for all calculations in this work is available as an open-source software project [77].

While in practice, the matrix $\mathbf{H}$ is extremely large, it is also extremely sparse, and it is easy to compute its matrix elements on the fly. To facilitate this, we index the matrix and vectors with the Fock states of Equation (9) directly. To encode the occupation number representation of a bosonic Fock state, we use a bit string where a sequence of n ones encodes n particles in a mode (lattice site), and zeros are used as separators between the modes. As an example, the state $|0,0,3,0,1,2\rangle$ would be encoded as the bit string “00111001011”. Using this scheme, storing N particles in M modes requires a bit string of length $N+M-1$. This representation is both extremely compact and allows for efficient on-the-fly calculations through bit manipulations.

The Rimu.jl code makes extensive use of Julia’s type system and code optimization capabilities through the multiple-dispatch paradigm and just-in-time compilation [76]. For example, the number of particles N and modes M, and the length of a bit string are all encoded in the type of a Fock-state address as type parameters. This allows us to easily write generic, well-tested, and reusable library code for manipulating bit strings and matrix-element calculation for bit strings of arbitrary length and type. As the type information is available at compiling time, part of the computational workload related to specializing the code to a specific physical problem is off-loaded to the compiler. Julia’s just-in-time compiler can thus produce optimized code for the particular parameters of the physics problem, which is easily defined in the script that is used to initiate the computation. As a consequence of this approach, some lag from the compilation is experienced in interactive use, but for the computationally intensive Monte Carlo calculations, the benefits from optimized code compilation are appreciable.

3.3. Data Structures and Distributed Computation

For representing the coefficient vector $\mathbf{c}$ it is important to access the data quickly based on the Fock space address. This is important as spawns hitting the same configuration must be allowed to annihilate [56], but also to save memory by encoding all walkers on a single configuration in a single number. We thus use a dictionary data structure to store the non-zero elements of $\mathbf{c}$, which is realized as a hash table [78] and thus provides access times that are nearly independent of the number on non-zero vector elements.

Another benefit of this approach is that it is relatively easy to distribute the data and computations to be processed in parallel. Our approach to parallelization follows Ref. [78] and divides the vector $\mathbf{c}$ into approximately equally-sized chunks, which are assigned to different workers. The workers perform the spawning step independently. After each step, but before the vector compression, a communication step is performed, where the newly spawned entries are transferred to the correct workers. In our implementation, we use the Message Passing Interface (MPI) [79] through its Julia bindings [80] to handle the data distribution and communication between workers.

3.4. The Initiator Approximation

With some Hamiltonians, FCIQMC exhibits the sign problem. The problem manifests itself when the number of walkers, which is equal to ${\parallel \mathbf{c}\parallel }_1$, is too small. In such a regime, the energy estimates given by FCIQMC become completely unusable [59].

A well-known solution to the sign problem in FCIQMC is the initiator approximation [57], which trades the sign problem for a small bias. It works by suppressing spawns from configurations with low walker occupation. To be precise, it divides the entries of $\mathbf{c}$ into two classes: initiators and non-initiators. For the initiators, the algorithm is unchanged, while the non-initiators are only allowed to perform spawns to configurations that are themselves, initiators. A configuration is an initiator if its occupation number is strictly greater than a chosen initiator threshold. In our computations, the initiator threshold was always set to 1.

3.5. Simulation Details

All simulations were performed with the Rimu.jl [77] code written by the authors. Energy estimators are computed as averages from a time series collected from the simulation, discarding data from an initial equilibration phase. The projected energy is used throughout this work as it has a much smaller fluctuation compared to the shift estimator, provided a sufficient number of walkers occupies the reference configurations. Error bars were determined using the blocking analysis of Ref. [81] supplemented by hypothesis testing of Ref. [82].

When calculating an expectated value of an observable, such as the two-body correlation and the momentum of the impurity, the replica trick [83] is used. It uses two independent FCIQMC wave functions to avoid a bias that would appear if correlated data was used.

For most of the calculations, one million floating-point walkers are used. As mentioned previously, the initiator approach is applied to all calculations with a threshold value of 1. This is necessary for controlling the sign problem in our simulations in the parameter regimes of larger values of $\gamma$ and $\eta$. We have performed extensive tests to control the biases introduced by population control and the initiator approximation and present some exemplary data from these efforts in Appendix A.

For calculations with small $\eta$, the equilibration can take a very long time. To overcome this problem, we used equilibrated wave functions from a system with much larger $\eta$ as the starting vector and re-equilibrated the wave function with the desired small $\eta$. This procedure speeds up the equilibration process significantly.

4. Results

Yrast states are the lowest energy states at a given non-zero momentum. We denote the energy of the yrast state $|{\mathsf{\Psi }}_P\rangle$ as $E_{N,N_{\mathrm{imp}}}\left(P\right)$ (or $E\left(P\right)$ for short), where N is the particle number of the Bose gas, $N_{\mathrm{imp}}$ the number of impurities present, and P the total (conserved) momentum. We also refer to the energy as a function of momentum as the “dispersion”. In the thermodynamic limit where $N,L\to \infty$ while the density $n=N/L$ is finite, the momentum becomes a continuous variable.

Yrast dispersions for a finite system with $N_{\mathrm{tot}}=N+N_{\mathrm{imp}}=20$ particles are shown in Figure 1. Special points on the dispersion occur at integer multiples of the “umklapp” momentum $P=2\pi \hslash N_{\mathrm{tot}}/L=N_{\mathrm{tot}}{P}_0$, where $P_0=2\pi \hslash /L$. At these umklapp points, the system’s internal state is identical to the ground state with a Galilean boost applied, such that every particle gains a momentum of unit $P_0$. Thus, the umklapp points have the energy $E\left(P\right)=E\left(0\right)+P^2/\left(2N_{\mathrm{tot}}m\right)$, as indicated by the dash-dotted (green) parabola in Figure 1.

For a pure one-dimensional Bose gas (where $N_{\mathrm{imp}}=0$), the yrast states and their energy can be generated exactly via the Bethe ansatz [84]. Note that yrast states in the one-dimensional Bose gas were denoted as type II excitations by Lieb in Ref. [84]. The yrast dispersion for $N=20$ bosons from FCIQMC is shown in Figure 1 by the blue line (solid line). The umklapp points at $P=N_{\mathrm{tot}}{P}_0$ (and integer multiples) have the meaning of a superfluid ring current [12]. The rest of the yrast dispersion is associated with dark and gray soliton phenomena [52] characterized by a localized dip in the density and step in the superfluid phase. While the momentum eigenstates are translationally invariant and can be thought of as a superposition of the (localized) solitons at various positions [48], wave-packet-like superpositions of nearby momentum eigenstates reveal localized soliton solutions that move at the velocity of $v=dE\left(P\right)/dP$, given by the slope of the yrast dispersion [51]. At the half-umklapp momentum $P=N_{\mathrm{tot}}{P}_0/2$, a dark soliton forms with a $\pi$ phase step. It is associated with a negative effective mass $m^{*}={(d^{2}E/d{P}^2)}^{-1}$, and, as a consequence, will oscillate around localized density maxima created by trapping potentials [85,86]. At small momentum (and next to any umklapp point), the Bose gas dispersion is linear, and the slope becomes the Bogoliubov speed of sound in the thermodynamic limit [84].

A yrast dispersion in the presence of a spin impurity ($N_{\mathrm{imp}}=1$) is shown with the orange line (dashed line) in Figure 1. The excitation energy of yrast states $E\left(P\right)-E\left(0\right)$ is generally smaller in the presence of a spin impurity compared to a Bose gas with the same number of particles $N_{\mathrm{tot}}$ apart from the umklapp points. This can be attributed to the fact that the spin impurity is not a part of the superfluid and thus does not fully contribute to the energy cost of forming a soliton by creating a twist in the phase—a phenomenon that can be rationalized with the phase rigidity of a superfluid [87]. Note that phenomena associated to phase rigidity occur in a one-dimensional Bose gas even though Bose–Einstein condensation is absent in the thermodynamic limit [12,51]. The yrast dispersion in the presence of the impurity is approximately quadratic at a small momentum (and near the umklapp points), in contrast to the pure Bose gas, and thus can be assigned an effective mass. The fitted, idealized parabola is shown in Figure 1 as a dotted (red) line. The effective mass determined by the curvature may differ from the bare mass m of the impurity due to interactions with the bosonic superfluid. We refer to this quadratic part of the dispersion near the umklapp points as the polaron.

Near the umklapp momentum of $P=N_{\mathrm{tot}}{P}_0/2$ we may expect the physics of a depleton, i.e., a dark or gray soliton that is affected by the presence of the impurity [19]. While Figure 1 depicts data for $N_{\mathrm{tot}}=20$ particles, the situation is generic (upon adjusting the scales) for a finite particle number where the umklapp momentum is found at $P=N_{\mathrm{tot}}{P}_0$.

4.1. Yrast Dispersion with Weak and Strong Boson–Impurity Coupling Strength

Figure 1 displays strong finite-size effects in terms of the center-of-mass dispersion $P^2/2N_{\mathrm{tot}}m$, the classical kinetic energy associated with the translation of the whole system, which provides a lower limit for the yrast excitation energies (shown as a dash-dotted line). In the thermodynamic limit, this energy contribution vanishes due to the total system mass appearing in the denominator. The detailed relation between the yrast dispersion of a finite system and its thermodynamic limit has been worked out for the pure Bose gas in Ref. [51] in terms of quantities, such as the phase step, the associated superfluid backflow current, and the depleted particle number for a dark/gray soliton. Here we correct for the dominant finite-size effect in a simple way by subtracting the center-of-mass kinetic energy from the yrast excitation energy. We thus define the finite-size-corrected yrast dispersion, $\mathsf{\Omega }\left(P\right)$, as

$$\mathsf{\Omega }\left(P\right) =E\left(P\right)-E\left(0\right)-\frac{P^2}{2N_{\mathrm{tot}}m},$$

(13)

where $E\left(P\right)$ is the lowest energy at fixed momentum P. The finite-size energy correction is equivalent to a Galilean boost into a reference frame that moves with the velocity $P/N_{\mathrm{tot}}m$ that a classical particle of mass $N_{\mathrm{tot}}m$ would have at momentum P. By removing the center of mass kinetic energy from the total energy, the yrast dispersion becomes periodic in P with the umklapp momentum $2\pi \hslash N_{\mathrm{tot}}/L=N_{\mathrm{tot}}{P}_0$ as the period, as in the thermodynamic limit. Additionally, the finite-size-corrected yrast dispersion has reflection symmetry across the half-umklapp point $N_{\mathrm{tot}}{P}_0/2$.

In Figure 2, we present two sets of finite-size-corrected yrast dispersions with boson–boson coupling strengths of $\gamma =0.02$ and $\gamma =0.2$, which are both considered to be weak interactions. The boson–impurity coupling is chosen in the range from $\eta =0.01$ to 1.0, which covers both $\eta >\gamma$ and $\eta <\gamma$ scenarios. Our results show that the yrast excitation energy is consistently lower in the presence of the spin impurity compared to the pure Bose gas at any value of $\eta$, as previously predicted [19,43]. The quadratic polaron part of the dispersion near $P=0$ and the umklapp points reduces its curvature with increasing $\eta$, which is consistent with an increase of the polaron effective mass. Quantitative results for the polaron effective mass were previously reported from diffusion Monte Carlo calculations [40] and mean-field theory [26].

As shown in Figure 2, the shape of the yrast dispersion is smooth in general. However, although very subtle, a cusp at the half-umklapp point can be seen in panel (b), when $\eta \ll \gamma$. A cusp for weak coupling in an infinite system was predicted in Ref. [43] based on the Luttinger liquid theory. It was pointed out that in a Luttinger liquid, the cusp is expected to vanish discontinuously at some critical value of the coupling strength between the impurity and the quantum liquid, but the exact transition point is difficult to determine. While we have only finite system data available from our calculations, we examine this transition in the remainder of this work and present further data that provides insights into the physics at play.

It is easiest to understand the origin of the cusp from the situation where the impurity does not interact with the Bose gas at all. We thus show data for a non-interacting impurity ($\eta =0$) immersed in a weakly interacting Bose gas at $\gamma =0.2$ in Figure 3, which demonstrates two interesting transition points in the yrast dispersion. The transition points originate from a trade-off between the energy cost of depositing momentum into either the impurity or the Bose gas. While for small momentum it is favorable to deposit momentum into the impurity (dotted red line), at $P>P_0$ it becomes favorable to deposit additional momentum into the Bose gas instead (orange squares show data where $P_{\mathrm{imp}}=P_0$). This is the first transition point. The second transition happens at the half-umklapp point (red diamond), where the yrast state (indicated with a green line) switches again to one with zero impurity momentum $P_{\mathrm{imp}}=0$. This switch generates a cusp in the yrast dispersion, connecting segments with a gray soliton moving to the right (at $P<10P_0$) and a gray soliton moving to the left ($P>10P_0$).

A symmetrical scenario to the first transition happens near the umklapp point. Both transitions become sharp quantum phase transitions (cusp with discontinuous derivative) in the thermodynamic limit. These phase transitions are first-order (level-crossing type) transitions without diverging quantum fluctuations.

4.2. Impurity Momentum

In order to understand the physical nature of a yrast state, it is of great interest to understand how the momentum is distributed between the impurity and the Bose gas. In the case of an interacting impurity, its momentum is no longer a good quantum number. Thus, we calculate the expectation value of the impurity momentum with respect to the yrast state $|{\mathsf{\Psi }}_P\rangle$

$${\langle {\widehat{P}}_{\mathrm{imp}}\rangle }_P={\displaystyle \sum }_{k=k_1}^{k_M}k{P}_0\langle {\mathsf{\Psi }}_P|{\widehat{b}}_k^{†}{\widehat{b}}_k|{\mathsf{\Psi }}_P\rangle .$$

(14)

Unbiased estimators for such a symmetric expectation value can be obtained from FCIQMC using the replica trick: two propagating independent stochastic representations of the quantum state are used for the bra and the ket state, respectively [83].

Figure 4 shows the expectation value of the impurity momentum ${\langle {\widehat{P}}_{\mathrm{imp}}\rangle }_P$ as a function of the total momentum of the yrast state for different values of the impurity coupling strength $\eta$ as blue dots. Because the total momentum is fixed to the value P for each state, the expectation value of momentum in the Bose gas is given by the difference ${\langle {\widehat{P}}_{\mathrm{Bg}}\rangle }_P=P-{\langle {\widehat{P}}_{\mathrm{imp}}\rangle }_P$, where ${\widehat{P}}_{\mathrm{Bg}}={{\displaystyle \sum }}_{k=k_1}^{k_M}k{P}_0{\widehat{a}}_k^{†}{\widehat{a}}_k$ is the operator for the momentum of the Bose gas alone.

In the case of weak impurity coupling ($\eta =0.01$) shown in Figure 4a, the situation is very close to the non-interacting limit of Figure 3 discussed in the previous section: For small total momentum $P=P_0$ the impurity carries (almost) the full momentum of the system, as this is energetically favorable. At larger values of P, additional momentum is taken up by the Bose gas while the impurity momentum stays at about $P_0$, before switching abruptly to approximately zero at the half-umklapp point. At the full-umklapp point $P=20P_0$ the impurity momentum jumps back to $P_0$, consistent with the expectation that every particle, including the impurity, carries a single unit of quantized momentum at the umklapp point. The abrupt change near the half-umklapp point $P=10P_0$ is consistent with the cusp observed in the yrast dispersion in Figure 2b. An interesting situation occurs directly at the half-umklapp point, where $\langle {\widehat{P}}_{\mathrm{imp}}\rangle \approx 0.5P_0$, which indicates that this state is an entangled superposition of a state where the impurity has momentum $P_0$ and the Bose gas $9P_0$, and a state with $\langle {\widehat{P}}_{\mathrm{imp}}\rangle \approx 0$ where the Bose gas carries the full (half-umklapp) momentum $10P_0$.

At larger interaction strengths $\eta$ shown in panels (b) to (d), the curves keep the inversion symmetry around the half-umklapp point. At this point, we find the entangled superposition state as described, with $\langle {\widehat{P}}_{\mathrm{imp}}\rangle \approx 0.5P_0$. Strong changes are found in the polaron regions. At intermediate interactions, additional momentum is deposited in the impurity with a maximum of $\langle {\widehat{P}}_{\mathrm{imp}}\rangle \approx 1.25$ at $\eta =\gamma =0.2$. The further increase of the impurity coupling then leads to a reduced expectation value for the impurity momentum going against $\langle {\widehat{P}}_{\mathrm{imp}}\rangle \approx 0.5P_0$ over the whole P range.

The panels of Figure 4 also show $m\frac{dE}{dP}\equiv mv$, where v is the group velocity of the system with orange crosses. This data indicates that the impurity moves with the group velocity in the polaron part of the dispersion relation (close to $P=0$ or umklapp points) for a weakly interacting impurity and over the whole dispersion relation when $\eta ⪆\gamma$. When $\eta \ll \gamma$ and outside of the polaron section, the impurity is rather transparent to the Bose gas and does not follow the group velocity. This indicates that the depleton picture where the impurity hybridizes with a dark or gray soliton is only valid when $\eta ⪆\gamma$.

4.3. Two-Body Correlation Function

The two-body correlation function contains important information about how particles interact with one another. In particular, the impurity–boson correlation provides direct evidence for the transition from a polaron to a depleton. Here, we define the dimensionless impurity–boson correlation function $g_P^{\left(2\right)}\left(d\right)$ for the yrast state $|{\mathsf{\Psi }}_P\rangle$ in real space as

$$g_P^{\left(2\right)}\left(d\right)=\frac{L}{N}{\int }_0^L\langle {\mathsf{\Psi }}_P|{\widehat{\psi }}^{†}(x+d){\widehat{\psi }}_{\mathrm{imp}}^{†}\left(x\right){\widehat{\psi }}_{\mathrm{imp}}\left(x\right)\widehat{\psi }(x+d)|{\mathsf{\Psi }}_P\rangle \mathrm{d}x,$$

(15)

where d is the distance between the impurity and a boson. In order to evaluate this correlation function in the lattice discretized model, we transform it into momentum space using ${\widehat{a}}_k^{†}=\int e^{ikx/\alpha }{\widehat{\psi }}^{†}\left(x\right)\mathrm{d}x$ and ${\widehat{b}}_k^{†}=\int e^{ikx/\alpha }{\widehat{\psi }}_{\mathrm{imp}}^{†}\left(x\right)\mathrm{d}x$, to obtain the equivalent representation

$$g_P^{\left(2\right)}\left(d\right) =\frac{1}{M}{\displaystyle \sum }_{s,p,q,r=1}^{M}exp\left(-id(p-q)\frac{2\pi }{L}\right)\langle {\mathsf{\Psi }}_P|{\widehat{a}}_s^{†}{\widehat{b}}_p^{†}{\widehat{b}}_q{\widehat{a}}_r|{\mathsf{\Psi }}_P\rangle {\delta }_{s+p,q+r}.$$

(16)

The chosen normalization ensures $g_P^{\left(2\right)}\left(d\right) =1$ in a non-interacting system for any yrast state. In an interacting system, $g_P^{\left(2\right)}$ still obeys a reflection symmetry $g_P^{\left(2\right)}\left(d\right)$ = $g_P^{\left(2\right)}(-d)$ and is a periodic function with period L. Furthermore, as a function of the yrast momentum P, the correlation function $g_P^{\left(2\right)}$ of a finite system is periodic in P with a reflection symmetry around $P=0$ and around the half-umklapp point, as does the yrast dispersion in relation to the thermodynamic limit. Due to these symmetries, we show the correlation functions only in the nontrivial intervals $0\le d\le L/2$ and $0\le P\le N_{\mathrm{tot}}{P}_0/2$.

Figure 5 shows the correlation functions $g_P^{\left(2\right)}\left(d\right)$ for yrast states with different momentums over a range of impurity–boson coupling strengths for $\gamma =0.2$. Significant changes in the correlation functions with respect to different momentum values are clearly visible.

For the smallest interaction strength $\eta =0.01$ in Figure 5a, we can identify clear evidence of the two transitions discussed in Section 4.1: For $P\le P_0$ the very small deviations of $g_P^{\left(2\right)}\left(d\right)$ from the background value of 1 indicate very weak correlations consistent with the polaron regime. There is evidence for a weak correlation hole, and the shape of $g_P^{\left(2\right)}\left(d\right)$ (negative curvature) is consistent with an otherwise homogeneous Bose gas. In the intermediate momentum range $P_0<P<10P_0$, the correlations are stronger, and the shape of $g_P^{\left(2\right)}\left(d\right)$ changes, displaying a positive curvature at a small P to a negative on at a larger P. This is consistent with the impurity weakly correlating with a gray soliton forming in the Bose gas—explaining the shape of the correlation function. A much stronger correlation is observed at the half-umklapp point $P=10P_0$, consistent with the superposition state expected at the cusp of the dispersion, as discussed in Section 4.2.

Increasing the impurity–boson coupling strength $\eta$ in panels Figure 5b–d changes the correlation function for different P values, it becomes smaller. For $\eta =0.05$ in panel (b), the transition from the half-umklapp momentum $P=10P_0$ to smaller momentum values is less dramatic and smoother, which indicates that the depleton picture of the impurity being localized inside a (modified) gray soliton is becoming adequate. However, compared with Figure 4b, we see that this is not yet completely the case and still requires a larger $\eta$ to become fully accurate.

The physics of the polaron regime is more resilient and survives larger values of $\eta$, up to $\eta \approx 0.2$, as seen in Figure 5b,c. We note that the impurity almost carries the full momentum of the system at $P=P_0$ for $\eta \le 0.2$, as seen in Figure 4, which is consistent with the polaron picture. Seeing the (anti-)correlation with the Bose gas strengthened with increasing $\eta$ in Figure 5b,c is consistent with the decrease in the curvature of the dispersion observed in Figure 2 and the associated increase in polaron mass.

At the largest value of $\eta =1$ shown in Figure 5d, the correlation function drops close to zero at zero distance $d=0$ for any value of P consistent with the picture that the impurity now acts as a weak link in the Bose gas, almost severing the superfluid [19]. At the half-umklapp momentum $P=10P_0$, the shape of the correlation function now closely traces the shape of a dark soliton density ∼$tanh{(d/l_{\mathrm{h}})}^2$, consistent with a healing length $l_{\mathrm{h}}\equiv L/\sqrt{2\gamma }N\approx 4.2\alpha$.

4.4. Effective Mass at Half Umklapp $P=N_{\mathrm{tot}}{P}_0/2$

At the half-umklapp point $P=N_{\mathrm{tot}}{P}_0/2$, we may expect the physics of the yrast states to be dominated by a dark soliton in the interacting Bose gas and by a depleton if an interacting spin impurity is present. The effective mass $m^{*}= {(d^{2}E/d{P}^2)}^{-1}$ is negative due to the concave shape of the dispersion relation. We extract the effective mass by fitting a parabola to three points of the finite-size-corrected dispersion relation $\mathsf{\Omega }\left(P\right)$ of Equation (13) near the half-umklapp point.

Figure 6 shows the extracted effective mass as a function of the impurity–coupling strength $\eta$ for two different values of the Bose gas interaction constant $\gamma$. In the regime $\eta >\gamma$, our data shows a linear trend with $\eta$. A linear dependence of the effective mass on $\eta$ is expected from exact results for an equal-mass impurity in a Tonks–Girardeau gas ($\gamma =\infty$) [19]. The effective mass becomes particularly heavy for a small $\gamma$ and a large $\eta$, up to several times the mass of the dark soliton at the same value of $\gamma$, as seen in Figure 6b. As the magnitude of the extracted effective mass (from the finite-size-corrected dispersion relation) becomes larger than the total system mass of $20m$, we call this the super-heavy regime. Note that without the finite-size correction of Equation (13), the curvature of the yrast dispersion changes from concave to convex, which means that the uncorrected effective mass diverges and changes sign (not shown). The heavy effective mass regime has potential experimental relevance, as it is relevant for realizing physical phenomena, such as Bloch oscillations [14,19]. Furthermore, a recent study demonstrates that the dynamical phenomenon of temporal orthogonality catastrophe is exhibited, given the impurity–boson couplings are sufficiently stronger than the intra-species background ones [28].

Another interesting feature shown in Figure 6b is that the impurity effective mass is approximately equal to the soliton mass for $\eta =\gamma$. In the regime where $\eta <\gamma$ (seen for $\gamma =0.2$), the effective mass becomes very small in magnitude, trending towards zero. This is consistent with the establishment of a cusp in the dispersion relation at the half-umklapp point, a feature that was already discussed in Section 4.1. While the concept of an effective mass breaks down in the cusp regime, our data can be used to determine that the transition happens approximately where $\eta =\gamma$, thus shedding some light on the question of the critical coupling which remains unsolved from Ref [43].

4.5. Spin-Flip Energy

Thus far, we have considered the yrast excitation energies, which measure how much energy is required to deposit momentum into the system on top of the energy of the ground state at $P=0$. Now we want to examine the energy that is required to flip a spin in the Bose gas at fixed momentum. We define the spin-flip energy $E_{\mathrm{SF}}\left(P\right)$ as

$$\begin{array}{c}\hfill E_{\mathrm{SF}}\left(P\right) =E_{N_{\mathrm{tot}}-1,1}\left(P\right)-E_{N_{\mathrm{tot}},0}\left(P\right).\end{array}$$

(17)

Figure 7 shows the spin-flip energy as a function of the impurity coupling strength $\eta$ for yrast states at different total momentum P. The two panels refer to different values of the boson–boson interaction strength $\gamma$. The brown line shows the spin-flip energy for the $P=0$ ground state. It is separated by a gap from the spin-flip energies at other momentum values, which are all lower. The ground state spin-flip energy increases with $\eta$ and crosses zero, meaning that for a large $\eta$, flipping the spin becomes energetically unfavorable. The vertical lines indicate where $\eta =\gamma$. At this point, the impurity is distinguishable from the background Bose gas, but all physical properties, such as mass and interactions, are the same. The spin-flip energy is thus solely due to quantum statistics. From the data shown in Figure 7, we see that the spin-flip energies are all negative at this point, and thus the system with impurity has lower energy than the pure Bose gas, i.e., it is favorable to flip the spin for any value of the total momentum P. Furthermore, the energy gain is larger the higher the momentum (up to the half-umklapp value). This trend is remarkably not maintained when $\eta <\gamma$, as seen in Figure 7b. This behavior can be rationalized from the quantitative changes in the yrast dispersions, shown in Figure 2.

5. Conclusions

Using the FCIQMC method, we investigated the properties of the yrast states of Bose gases coupled with a mobile impurity in one spatial dimension. Based on the energies and the first and second-order correlation functions of yrast states, we identified the polaron and depleton regimes, as well as the transitions between them. The extracted depleton effective mass revealed a super-heavy regime where the magnitude of the (negative) depleton mass exceeds the mass of the finite Bose gas. We also observed a qualitative change in behavior crossing $\eta =\gamma$ in all calculated quantities. For the $\eta >\gamma$ regime, we can identify the formation of depletons around the half-umklapp point, where the impurity is more or less confined to the density hole of the gray/dark soliton of the Bose gas. The depleton picture becomes inadequate for smaller interactions between the impurity and the bosons, $\eta <\gamma$, with the impurity no longer hybridizing with the soliton. This behavior is consistent with an observed break-down of the effective mass concept below $\eta =\gamma$.

In this work, the FCIQMC method was applied to a bosonic many-body problem for the first time. Due to the non-stoquastic nature of the momentum-space Hamiltonian (3), the sign problem exists and becomes severe when either $\eta$ or $\gamma$ is large. Through this study, we demonstrated the effective suppression of the sign problem in FCIQMC by the application of the initiator approximation, showing the potential of FCIQMC for studying complex bosonic many-body systems.

Outlook: The demonstrated computational method is extremely versatile and can be applied to a wide range of physics questions. Possible future extensions of this study include extrapolating the results to the thermodynamic limit. A transcorrelated Hamiltonian [88,89] can be applied to accelerate the basis set convergence to the infinite limit. There are also many interesting set-ups that we wish to study further. In this work, we only focus on the cases where the impurity and bosons all have identical mass and repulsive interactions, which could be extended to unequal masses and attractively interacting impurities. Attractive impurity–boson coupling has been studied in the polaron regime [26,40] but not yet explored in the context of depleton physics. In addition, the case of an impurity in a strongly interacting Bose gas or with long-range interactions [90,91,92,93] is interesting to study, where perturbative and mean-field approaches are of limited use or invalid. A more complex system with two impurity atoms, known as the bipolaron problem [27,94,95,96,97], at non-zero momentum is also interesting due to its connection to high-temperature superconductivity [95,96].