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

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 to 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.


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].
Many of the mentioned works focus on ground state properties and effective mass of the one-dimensional Bose polaron with either zero or very small total momentum. Although the dynamics of impurities has 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 lowest energy at given momentum. Yrast states are stable as long as momentum is conserved, while adiabatic passage through the yrast states of different momentum 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 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 sprit 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 like 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 lowestenergy 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-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 analysing 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 Sec. 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 Sec. 3. Numerical results on properties of the yrast states and their interpretation in terms of the different physical regimes are presented in Sec. 4 starting with the yrast dispersion in Sec. 4.1. This is followed by the impurity momentum and the impurity-boson correlation function in Secs. 4.2 and 4.3, respectively. The effective mass in the soliton/depleton region of the dispersion is reported in Sec. 4.4 and the spin-flip energies in Sec. 4.5. Finally, we draw conclusions in Sec. 5 and we outline 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.

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 where x i (i = 1, · · · , N) and x 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 BB while the interaction of the impurity particle with the bosons is described by a contact potential of strength g IB . We consider repulsive interactions with g BB > 0 and g IB ≥ 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 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.

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 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 ∈ (−L/2, L/2] and the lattice constant is defined as α = L/M for M lattice points. In this representation, the Hamiltonian reads whereâ † k (â k ) are the boson creation (annihilation) operators; the corresponding operators for the impurity areb † k (b k ). The plane-wave eigenstates x|â † k |vac = e −ikx/α of momentumhkα are indexed with the dimensionless quantum numbers where j ∈ {1, 2, . . . , M} is an integer. The kinetic energy dispersion is the same for bosons and impurity (as they have equal mass) where we have introduced the unit of energy that will be used throughout this work 1 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]: where g 0 = π 2h2 /mα.

Connection to mean-field theory and choice of parameters
In the weakly interacting regime where γ 1, nonlinear phenomena in the Bose gas like 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 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 l 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 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 Eq. (3). As a compromise we 1 Note that Refs. [26,40] choose the Fermi energy ε F = π 2h2 n 2 2m = π 2 N 2 2 ε 0 as the energy unit. choose to work with N 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 h ≡ √ 2γN ≈ 12 for γ = 0.2 and L/l h ≈ 4 for γ = 0.02 for the two values of the Bose-gas interaction strength that we are using in this work.

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.

Bosonic Full Configuration Quantum Monte Carlo
In FCIQMC a basis of Fock states (occupation number basis) for N = ∑ M i=1 n i particles in M lattice sites is used Within this basis the Hamiltonian is represented as a matrix H and the quantum state (manybody wave function) as a vector 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 where the parameter δτ controls the size of the time step, c (n) is the approximation of the eigenvector at the n-th time step, and 1 is the identity matrix. The shift S (n) is a real number used to keep the norm of c (n) under control. It is adjusted by the following scheme: where N (n) w ≡ c (n) 1 is the 1-norm of c (n) , N t the parameter for the target norm, and ξ and ζ parameters that control the dynamics of the shift. In the steady state, the instantaneous norm N (n) w fluctuates around the value of N 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 H and c quickly become prohibitively large. To get around this problem, we replace the matrixvector multiplication in Eq. (10) with a stochastic sampling process. The sampling process is designed to reproduce the right hand side of Eq. (10) by 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 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 H.
The spawning from the configuration q to the configuration r = q can be described as where 1 p spawn is the inverse probability of picking r, i.e. the number of nonzero off-diagonal entries in the q-th column of H. If the occupation number 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 Eq. (10) is performed exactly. After a step is complete, we stochastically project the entries v i of the vector c to a threshold t; values |v i | < t are removed from the vector with probability p = 1 − |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 c (n) can be much smaller than the dimension of the Hilbert space while the expectation value of c (n) still approaches the exact eigenvector of the ground state of H. At the same time, the shift S (n) 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 Eq. (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].

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, 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 setup 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 highperformance 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 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 Eq. (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 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]. E.g., the number of particles N and modes M, and the length of a bitstring 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 compile 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 compilation is experienced in interactive use, but for the computationally intensive Monte Carlo calculations, the benefits from optimized code compilation are appreciable.

Data structures and distributed Computation
For representing the coefficient vector 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 c, which is realized as a hash table [78] and thus provides access times that are nearly independent of the number on nonzero 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 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.

The Initiator Approximation
With some Hamiltonians, FCIQMC exhibits the sign problem. The problem manifests itself when the number of walkers, which is equal to c 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 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.

Simulation Details
All simulations were performed with the Rimu.jl [77] code (version v0.6.0) 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 comparing to the shift estimator, provided 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 expectation 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 γ and η. 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 App. A.
For calculations with small η, the equilibration can take a very long time. To overcome this problem, we used equilibrated wave functions from a system with much larger η as the starting vector, and re-equilibrated the wave function with the desired small η. This procedure speeds up the equilibration process significantly.

Results
Yrast states are the lowest energy states at a given non-zero momentum. We denote the energy of the yrast state |Ψ P as E N,N imp (P) [or E(P) for short], where N is the particle number of the Bose gas, N 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 → ∞ while the density n = N/L is finite, the momentum becomes a continuous variable.
Yrast dispersions for a finite system with N tot = N + N imp = 20 particles are shown in Fig. 1. Special points on the dispersion occur at integer multiples of the "umklapp" momentum P = 2πhN tot /L = N tot P 0 , where P 0 = 2πh/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(P) = E(0) + P 2 /(2N tot m), as indicated by the dash-dotted (green) parabola in Fig. 1. For a pure one-dimensional Bose gas (where N imp = 0), the yrast states and their energy can be generated exactly via the Bethe ansatz [84] 2 . The yrast dispersion for N = 20 bosons from FCIQMC is shown in Fig. 1 by the blue data (solid line). The umklapp points at P = N tot P 0 (and integer multiples) have the meaning of a superfluid ring current [12]. The rest of the yrast dispersion are 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(P)/dP, given by the slope of the yrast dispersion [51]. At the half umklapp momentum P = N tot P 0 /2, a dark soliton forms with a π phase step. It is associated with a negative effective mass m * = (d 2 E/dP 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].
An yrast dispersion in the presence of a spin impurity (N imp = 1) is shown with orange data (dashed line) in Fig. 1. The excitation energy of yrast states E(P) − E(0) is generally smaller in the presence of a spin impurity compared to a Bose gas with the same number of particles N 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] 3 . The yrast dispersion in the presence of the impurity is approximately quadratic at 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 Fig. 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 half umklapp momentum at P = N 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 Fig. 1 depicts data for N tot = 20 particles, the situation is generic (upon adjusting the scales) for finite particle number where the umklapp momentum is found at P = N tot P 0 . Figure 1 displays strong finite size effects in terms of the center-of-mass dispersion P 2 /2N 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 like 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, Ω(P), as

Yrast dispersion with weak and strong boson-impurity coupling strength
where E(P) 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 tot m that a classical particle of mass N 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πhN tot /L = N 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 tot P 0 /2. In Fig. 2, we present two sets of finite-size corrected yrast dispersions with boson-boson coupling strengths of γ = 0.02 and γ = 0.2, which are both considered to be weak interactions. The boson-impurity coupling is chosen in the range from η = 0.01 to 1.0, which covers both η > γ and η < γ 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 η, as previously predicted [19,43]. The quadratic polaron part of the dispersion near P = 0 and the umklapp points reduces its curvature with increasing η, 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 Fig. 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 η γ. A cusp for weak coupling in an infinite system was predicted in Ref. [43] based on 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 (η = 0) immersed in a weakly interacting Bose gas at γ = 0.2 in Fig. 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 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 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.

Impurity momentum
In order to understand the physical nature of an 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, it's momentum is no longer a good quantum number. Thus, we calculate the expectation value of the impurity momentum with respect to the yrast state |Ψ P 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 P imp P as a function of the total momentum of the yrast state for different values of the impurity coupling strength η 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 P Bg P = P − P imp P , wherê P Bg = ∑ k M k=k 1 kP 0â † kâ k is the operator for the momentum of the Bose gas alone. In the case of weak impurity coupling (η = 0.01) shown in Fig. 4(a) the situation is very close to the non-interacting limit of Fig. 3 Figure 4. The expectation value of the impurity momentum P imp P against the total momentum P in the system. The dots (•) are directly calculated data with FCIQMC. The crosses (×) show the m dE dP computed numerically from the yrast spectum in Fig. 2. The dotted line shows the value P imp P = 0.5P 0 as a guide to the eye. The boson-boson coupling is γ = 0.2 for all cases, and N = 19 and N imp = 1, which means that P = 10P 0 corresponds to half umklapp and P = 20P 0 is the full umklapp point. total momentum P = P 0 the impurity carries (almost) the full momentum of the system as this is energetically favourable. 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 quantised 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 Fig. 2(b). An interesting situation occurs directly at the half umklapp point where P imp ≈ 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 P imp ≈ 0 where the Bose gas carries the full (half-umklapp) momentum 10P 0 .
At larger interaction strengths η 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 P imp ≈ 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 P imp ≈ 1.25 at η = γ = 0.2. Further increase of the impurity coupling then leads to a reduced expectation value for the impurity momentum going against P imp ≈ 0.5P 0 over the whole P range.
The panels of Fig. 4 also show m dE dP ≡ 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 weaklyinteracting impurity, and over the whole dispersion relation when η γ. When η γ 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 η γ.

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 (2) P (d) for the yrast state |Ψ P in real space as g (2) 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 into momentum space usingâ † k = e ikx/αψ † (x) dx andb † k = e ikx/αψ † imp (x) dx, to obtain the equivalent representation The chosen normalization ensures g P (d) = 1 in a non-interacting system for any yrast state. In an interacting system g (2) P still obeys a reflection symmetry g  P (−d) and is a periodic function with period L. Furthermore, as a function of the yrast momentum P, the correlation function g (2) P of a finite system is periodic in P with reflection symmetry around P = 0 and around the half-umklapp point, as does the yrast dispersion relation in the thermodynamic limit. Due to these symmetries we show the correlation functions only in the nontrivial intervals 0 ≤ d ≤ L/2 and 0 ≤ P ≤ N tot P 0 /2. (d) P/P 0 = 0 P/P 0 = 1 P/P 0 = 2 P/P 0 = 3 P/P 0 = 4 P/P 0 = 5 P/P 0 = 6 P/P 0 = 7 P/P 0 = 8 P/P 0 = 9 P/P 0 = 10 P Figure 5. The impurity-boson correlation function g  For the smallest interaction strength η = 0.01 in Fig. 5(a) we can identify clear evidence of the two transitions discussed in Sec. 4.1: For P ≤ P 0 the very small deviations of g There is evidence for a weak correlation hole, and the shape of g (2) P (d) (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 (2) P (d) changes with displaying positive curvature at small P to negative at 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 Sec. 4.2.
Increasing the impurity-boson coupling strength η in panels 5(b)-(d) the changes in the correlation function for different P values become smaller. For η = 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, comparing with Fig. 4(b) we see that this is not yet completely the case and still requires larger η to become fully accurate.
The physics of the polaron regime is more resilient and survives to larger values of η up to η ≈ 0.2 as seen in panels 5(b) and (c). We note that the impurity carries almost the full momentum of the system at P = P 0 for η ≤ 0.2 as seen in Fig. 4, which is consistent with the polaron picture. Seeing the (anti-)correlation with the Bose gas strengthened with increasing η in panels 5(b) and (c) is consistent with the decrease in the curvature of the dispersion observed in Fig. 2 and associated increase in polaron mass.
At the largest value of η = 1 shown in Fig. 5(d) 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 h ) 2 , consistent with a healing length l h ≡ L/ √ 2γN ≈ 4.2α.
At the half umklapp point P = N 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/dP 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 Ω(P) of Eq. (13) near the half umklapp point. Figure 6 shows the extracted effective mass as a function of the impurity coupling strength η for two different values of the Bose gas interaction constant γ. In the regime η > γ our data shows a linear trend with η. A linear dependence of the effective mass on η is expected from exact results for an equal-mass impurity in a Tonks-Girardeau gas (γ = ∞) [19]. The effective mass becomes particularly heavy for small γ and large η, up to several times the mass of the dark soliton at the same value of γ, as seen in Fig. 6(b). 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 Eq. (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 Fig. 6(b) is that the impurity effective mass is approximately equal to the soliton mass for η = γ. In the regime where η < γ (seen for γ = 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 Sec. 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 η = γ, thus shedding some light on the question of the critical coupling which remains unsolved from Ref [43].

Spin-flip Energy
So 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 SF (P) as E SF (P) = E N tot −1,1 (P) − E N tot ,0 (P).
(17) Figure 7 shows the spin-flip energy as a function of the impurity coupling strength η for yrast states at different total momentum P. The two panels refer to different values of the boson-boson interaction strength γ. The brown data 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 η and crosses zero, P/P 0 = 0 P/P 0 = 1 P/P 0 = 2 P/P 0 = 3 P/P 0 = 4 P/P 0 = 5 P/P 0 = 6 P/P 0 = 7 P/P 0 = 8 P/P 0 = 9 P/P 0 = 10 meaning that for large η flipping the spin becomes energetically unfavorable. The vertical lines indicate where η = γ. 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 Fig. 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 η < γ as seen in Fig.  7(b). This behavior can be rationalized from the quantitative changes in the yrast dispersions shown in Fig. 2.

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 η = γ in all calculated quantities. For the η > γ 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, η < γ, with the impurity not longer hybridizing with the soliton. This behavior is consistent with an observed break-down of the effective mass concept below η = γ.
In this work, the FCIQMC method is 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 η or γ 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 question. 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 impurityboson 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 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,[90][91][92][93], at non-zero momentum is also interesting due to its connection to high-temperature superconductivity [91,92].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Eliminating Biases
While bosonic systems can often be described by stoquastic Hamiltonians characterized by having only non-positive, real off-diagonal elements, the momentum-space Hamiltonian of Eq. 3 considered here, is non-stoquastic. As a consequence, one has to deal with the QMC sign problem, that originates from the fact that different configurations can spawn into the same configuration with incoherent signs.
In FCIQMC, the initiator approach can be used to mitigate this sign problem by restricting the walker spawning process to the dominant configurations. This enforces a better coherence in the sign structure of the wave function. Albeit typically small [57,78], an initiator bias can be observed as a consequence of the initiator approximation when an insufficient number of walkers is used to sample a much larger Hilbert space. Furthermore, the population control bias [68,94,95] is a stochastic bias that appears as a result of sampling noise. It is typically much smaller than the initiator bias (where the latter is applicable) and scales like a power law with the number of walkers [68]. Both biases can be reduced below the size of statistical error bars by increasing the walker population. To make sure our calculated energies are bias free and a sufficiently large walker population is used, one can check the ground-state energy as a function of the equilibrated walker number, N w , as shown in Fig. A1. It can be seen that when N w > 10 4 the biases in the shift and projected energies are smaller than the statistical errors, and converged to the same energy throughout. This convergence check is carried out with a larger interaction strength (η = 2, γ = 0.2) than used for any of the data presented in the main article, and thus should overestimate the bias for the presented data. For all energies presented in Sec. 4, a walker population of N w = 10 6 is used. Hence we are confident that the presented data is free of both, the initiator and the population control bias.