{'fraction_non_alphanumeric': 0.031275198709826056, 'fraction_numerical': 0.008178781246400185, 'mean_word_length': 4.7550546900894926, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Several fermions strongly interacting with a heavy mobile impurity in a one-dimensional harmonic trap

We propose a numerically exact method for a mixture with a single impurity immersed in several majority fermions, confined in a harmonic potential. We separate one of the degrees of freedom through an appropriately tailored canonical transformation and perform exact diagonalization on the simplified Hamiltonian. This method is especially effective for a heavy impurity, where it outmatches the typical exact diagonalization approach. We used our method to calculate energy and density profiles of the first few eigenstates for the mixture with up to ten majority fermions.


I. INTRODUCTION
Recent progress in experiments with ultracold atoms has allowed for physical realizations of seemingly purely theoretical models. Notable examples are few-body fermionic mixtures produced in Heidelberg [1][2][3]. In these experiments, 6 Li atoms, in two different hyperfine states, have been trapped in a quasi-one-dimensional harmonic potential. Atomic species interact by contact interactions, which can be tuned using Feshbach resonance. The number of atoms in the mixture has also been precisely controlled. These experiments are the motivation for the further theoretical studies of few-body ultracold mixtures in one dimension (see [4] for a review).
In the present work, we analyze a special case of a fermionic mixture with a single particle immersed in N majority fermions. In the many-body limit (an impurity in a Fermi sea), this is known as a Fermi polaron problem [5]. In the few-body case, this problem has been extensively studied for equal masses of components. For a mixture in uniform potential, the eigenproblem has been solved analytically, both for repulsive [6] and attractive interactions [7]. For the trapped systems, especially important are mixtures in harmonic confinement since they were obtained experimentally [3]. In the simplest twobody case (N = 1), there exist an analytical solution [8]. For larger mixtures, the problem has been analyzed numerically [9][10][11][12][13][14]. Other studied confinements include double-well potential [15], optical lattice [16], and a combined harmonic and periodic potential [17]. Moreover, the mixture has been analyzed for an arbitrary confining potential in the regime of strong interactions [18][19][20]. In the case of mass imbalance,the problem has been addressed for the uniform potential [21], mainly for the light impurity [22,23].
In our analysis, we focus on the strongly interacting mixture in the harmonic confinement with the heavy impurity [24,25]. We use a numerical method specifically designed to handle this problem. We perform dedicated canonical transformation on the Hamiltonian. The resulting Hamiltonian is split into the center-of-mass part, which has an analytical solution, and the relative part described as a system of weakly interacting fermions. The latter is solved using the exact diagonalization method. This paper is organized as follows. The analyzed model is described in Section II. In Section III, we introduce the coordinate transformation used to simplify the problem. In Section IV we present the numerical method used to solve the problem and compare it to the typical approach. The results for the first few stationary states of the system are presented in Section V. Finally, Section VI contains our conclusions.

II. THE SYSTEM STUDIED
In this paper, we consider a system of N spinless fermions with mass m, interacting with a single impurity through contact interactions. The mixture is trapped in a one-dimensional harmonic trap with the same frequency Ω for both components. The analysis is performed in natural units of the harmonic oscillator for the majority fermions. Therefore units of energy, length, and interaction strength are Ω, mΩ , and m 3 Ω , respectively. The Hamiltonian of the system has a form: where M is the mass of the impurity expressed as a ratio to the mass of the other component. The Hamiltonian is partially written in the second quantization. The majority fermions are described with fermionic field operatorŝ ψ † (x) andψ(x), which obey an anti-commutation rela- . Therefore their fermionic statistic is incorporated into the Hamiltonian. For the impurity, we use the first quantization because it simplifies the coordinates transformation presented in the next section. The Hamiltonian in this mixed formalism acts on a tensor product of a single-particle Hilbert space in the position representation and a Fock space in the occupation number representation. In the analysis, we focus on the case of strong but finite repulsive interactions (g ≥ 1) and a heavy impurity (M > 1).
Calculating properties of this system is a numerical challenge due to the strong interactions. One can arXiv:2206.01086v2 [cond-mat.quant-gas] 17 Nov 2022 approach this problem using the exact diagonalization method [26]. It allows for high accuracy but requires large amounts of computational resources. Therefore it is limited to few-body systems. On the other hand, this approach has a wide range of applications -it can be used for any mass ratio M , any trapping potential, as well as for the systems with more than one impurity. In this paper, we propose an alternative solution, designed specifically for systems with a single heavy impurity in uniform or quadratic potential.

III. THE TRANSFORMATION
In the proposed method, we start by changing the coordinate system. Instead of the particle positions, we use the center-of-mass coordinate and the relative distances between the fermions and the impurity: New coordinates are obtained through the canonical transformation governed by the operator: are single-particle position and momentum operators for the majority fermions, respectively. By using this transformation, we separate the Hamiltonian (1) into two parts:Ĥ where:Ĥ describes a harmonic oscillator with eigenstates: (H n () is n-th Hermite polynomial) and eigenenergies E n = n + 1 2 . The second term describes a system of N interacting fermions: with single-particle part: and an effective interaction: This transformation can be viewed as an extension of the so-called Lee-Low-Pines (LLP) transformation [27,28], which also can be applied to the Fermi polaron problem [29]. The LLP transformation separates relative motion, but the applicability is limited to the uniform potential. Our transformation generalizes this idea to any potential with a separable center of mass.

IV. THE EXACT DIAGONALIZATION
After the transformation is performed, the problem is reduced to the system of N fermions described by the HamiltonianĤ x . The effective interaction between these fermions depends both on their positions and momenta. Moreover, the interaction term is inversely proportional to the mass ratio M . Therefore, in the case of heavy impurity, it is a weakly interacting system. In order to solve this problem, we use the exact diagonalization method. In this approach, we want to represent the Hamiltonian (8) as a matrix in the Fock basis of the noninteracting system. We begin by decomposing the field operator in a single-particle basis: whereâ m are annihilation operators, and φ m (x) are eigenfunctions of the single-particle part of the hamil-tonianĤ x : Fortunately, this eigenproblem can be solved analytically [8]. Using an effective mass µ = , the eigenfunctions can be written as: if m is even (13) where U () is the confluent hypergeometric function of the second kind. The eigenvalues are equal to m = ν m , where ν m = m+ 1 2 for odd m, while for even m, they are solutions of the equation: After applying decomposition (11), we rewrite the Hamiltonian (8) to the form In this matrix representation of the hamiltonianĤ x , every element corresponds to the Slater determinant of N single-particle states. The matrix has the infinite number of elements, but high-energy states are numerically irrelevant when we calculate the low-energy many-body eigenstates. Therefore we cut the sum in (11) on some index K. It is clear that, for a big enough cutoff K, the problem can be solved with any desired accuracy. As a result, the hamiltonianĤ x is reduced to a finite matrix that can be diagonalized. The overall spectrum of hamil-tonianĤ is obtained by combining the eigenvalues of this matrix with the center-of-mass energies E n . To evaluate how effective the new method is, we compare its predictions with the exact diagonalization used directly on the original Hamiltonian (1). In this standard approach, both components of the mixture are described in the second quantization with the field operators. The operators are decomposed in a single-particle basis (with the same cutoff K) and used to rewrite the Hamiltonian into a matrix which is then diagonalized. We check how accurately these two methods obtain the ground state of the system depending on the value of the cutoff. The results are presented in Fig 1. The top figures show the ground state energy as a function of the cutoff K. The comparison is made for the system with N = 3 fermions and four different sets of parameters M and g. As expected, the results tend to decrease with increasing cutoff, as they approach the actual groundstate energy. It is also clear that, for all considered cases, the proposed method is much more accurate. That is the case both for heavy impurity and for the mass-balanced system. The reason is that accuracy of this method is less affected by the strength of the contact interactions since they are included directly in single-particle functions and not in the interaction term. For the system with a very light impurity, the effective interaction term (10) would be large and, therefore, the proposed method would be less effective.
For further comparison, we also calculate the fidelity: where |Ψ K is the eigenvector of the system with the cutoff K. This quantity describes how much the resulting state change when the cutoff increase. It is clear that as the eigenvector approaches the actual ground state, the fidelity should approach unity. For convenience, the results in Fig 1b are presented as 1 − F K such that they can be displayed in the logarithmic scale. The plot shows that the fidelity (as a function of the cutoff) approaches unity much faster in the proposed method. This again shows the advantage of this method. Moreover, the results become even better for larger impurity mass, which is not the case for the conventional exact diagonalization. On the other hand, stronger contact interactions decrease the accuracy of both methods. It should also be recalled that, in the proposed approach, one degree  of freedom is calculated analytically. Therefore even for the same cutoff, the new method requires less computational resources than the typical exact diagonalization method. Fidelity can be used to determine whether the cutoff is large enough for accurate numerical calculations. Here we assume that cutoff K is sufficiently large when for every k ≥ K, fidelity F k is larger than 0.999. Such sufficient cutoffs are plotted in Fig 2 for a system with heavy impurity and strong interactions (M = 5 and g = 5). Results were calculated for the ground state and the first few excitations of the relative motion. As one can see from the graph, both the ground state and the lowest excited state require a relatively low cutoff for accurate calculations (thanks to significant mass imbalance). It should be noted that choosing higher center-of-mass excitations does not influence the accuracy of the method, since this part of the Hamiltonian is always solved analytically.

V. STATIONARY PROPERTIES
Using our method, we calculate some stationary properties of the few-body fermi gas interacting with the impurity. We start with the ground state energy for the mass-balanced systems (M = 1). In this case, there are available experimental data [3], which can be compared to the numerical results. In the experiment, the researchers measured the interaction energy instead of the ground state energy. Interaction energy is defined as a difference between the energy of the system with and without the inter-component interactions. The ground state energy of a non-interacting system is known analytically to be equal to N 2 +1 2 . Therefore the interaction energy is directly related to the ground state energy of the interacting system. A comparison of numerical results with experimental data is presented on the left panel of Fig 3. Three lines correspond to weak (g = 0.36), medium (g = 1.14), and strong (g = 2.80) interactions. As one can see from the graph, our results are in good agreement with the experimental data. It may seem that there is a significant relative error for N = 5 and strong interactions. That is because we are comparing the interaction energy, which, in this case, is only a small part of the ground state energy that we calculate. For the ground state energy, the relative error is smaller than 1% for all compared data.
On the right panel of Fig 3, there is interaction energy calculated for the mass-imbalanced systems in the strong interaction regime (g = 2.80). As can be seen from the graph, the interaction energy, as a function of the particle's number, changes qualitatively for a strong mass imbalance, from a smooth curve to a step function. Increasing the number of fermions by one to the even values affects the interaction energy significantly less than with the odd values. This behavior can be understood using our transformation. In the limit of infinitely heavy impurity, the effective interaction term (10) disappears. Therefore, the relative part of the hamiltonian reduces to:Ĥ which describes non-interacting fermions in the harmonic trap with the additional delta potential in the center. In the ground state, the fermions occupy the first N single-particle eigenstates. Every second such state is an odd function with respect to the spatial inverse x → −x.
Thus it is not affected by the Dirac potential in the middle. These single-particle states do not depend on the interaction strength g. Adding a new particle to the system with an even number of fermions increases groundstate energy in the same way regardless of whether there are interactions. Therefore interaction energy does not increase. In this limit, steps are ideally horizontal. Besides energies, one can also calculate density profiles. For the coordinates after the transformation, the density profiles of the eigenstate |Φ n ⊗ |Ω m (center-ofmass and relative part respectively) are equal: where Φ n (y) are functions from eq (7). Finding this quantity for the original Hamiltonian is not as simple as for the energy. The eigenstate has to be transformed back using the operatorŜ: whereΦ n () is an extension of the scalar function Φ n () to the function acting on matrices. Since Φ n () is an analytical function, this extension can be obtained rigorously by replacing scalars with matrices in the Taylor series. For the eigenstate in the original coordinates (20), density profiles are equal: whereψ(x − y) = e iypψ (x)e −iyp . This can be generalized further to many-body density profiles: From the numerical standpoint, one can calculate these quantities through diagonalization of the operatorx in the Fock basis. The single-particle density profiles for the system with N = 10, interaction strength g = 5 and for gradually increasing mass ratio M are presented in Fig 4. In the ground state (first column), the density profile for fermions doesn't change significantly as mass imbalance increases. Even for M = 5, it is qualitatively the same as for M = 100. Quantitatively, the most noticeable difference is that the local minimum at the center decreases as the mass imbalance increases. This is expected since the impurity becomes more and more localized in the middle of the trap and repulses the fermions from there. Surprisingly, the results for the excited states (second and third column) change qualitatively as mass imbalance increases. For M = 10 there is still a local maximum in the middle of the trap, which is not the case in the heavy impurity limit. Even for M = 20 shape of the density profile is significantly different than for M = 100. Therefore finite-mass corrections are important even for the system with relatively heavy impurity, at least for the excited states.

VI. CONCLUSION
In this paper, we've presented a new method to calculate the properties of the system of a few fermions interacting with a heavy impurity and confined in a harmonic trap. By comparing it to the standard approach of the exact diagonalization, we've shown that this new method allows for high accuracy of the results with smaller restrictions on the size of the system. Next, the method has been used to calculate basic stationary properties for the systems with up to 10 fermions, namely the energies and the density profiles of the eigenstates.
It is worth noticing that this method can also be used to calculate the dynamic properties of the analyzed system. Once the initial state is transformed to the new coordinates, it can be decomposite into the Fock basis and evolve using any standard numerical method (for example Runge-Kutta algorithm). For the evolved state, the energy and density profiles can be calculated analogously to the eigenstates.