A systematic interpolatory method for an impurity in a one-dimensional fermionic background

We explore a new numerical method for studying one-dimensional quantum systems in a trapping potential. We focus on the setup of an impurity in a fermionic background, where a single distinguishable particle interacts through a contact potential with a number of identical fermions. We can accurately describe this system, for various particle numbers, different trapping potentials and arbitrary finite repulsion, by constructing a truncated basis containing states at both zero and infinite repulsion. The results are compared with matrix product states methods and with the analytical result for two particles in a harmonic well.


Introduction
The investigation of one-dimensional quantum systems of interacting particles has, in the last decades, attracted renewed interest due to striking advances in experiments with cold atoms in optical traps [1]. Paradigmatic models extensively explored in the fields of condensed matter [2][3][4][5] and mathematical physics [6][7][8] are now within reach of experiments, and their exotic properties can be measure with great precision. Moreover, the degree of control over several experimental parameters, including interactions between the atoms [9][10][11][12] and trapping geometries opens up the possibility of using such experiments as quantum simulators for a multitude of interesting models [13], even beyond usual condensed matter systems [14]. One particular problem which has attracted interest in this context is that of a single distinct atom (or impurity) embedded in a background of identical particles. In the field of condensed matter, such systems can present interesting phenomena, such as the Kondo effect [15] and the orthogonality catastrophe [16]. Theoretical and experimental studies with ultracold atomic setups have extensively explored both the bosonic [17][18][19][20][21][22][23][24][25][26][27] and fermionic [28][29][30][31][32][33][34][35][36][37] manifestations of these models -the so-called Bose and Fermi polarons, respectively. The one-dimensional fermionic case, in particular, dates back to McGuire's impurity model in a homogeneous geometry [38,39], which is exactly solvable through the Bethe ansatz approach [40]. Other approaches have later generalized the study of static properties to mixed fermionic systems in harmonic potentials [41][42][43][44][45][46][47]. On the dynamical side, impurity models have been shown to present exotic effects such as Bloch oscillations [48] and quantum flutter [49,50].
In this work we present an original way to obtain the static properties of an impurity in a fermionic background, where the number of background fermions can be arbitrarily modified. We employ a variational principle where our ansatz for the wavefunction is a combination of states at zero and infinite interaction, relying on the fact that the analytical expressions for these limits are known. In practice, we construct a truncated basis by choosing a certain number of states at each limit and then employ the Gram-Schmidt ortonormalization process to construct an orthonormal basis. By diagonalizing the Hamiltonian in this basis, we obtain an approximation for the wavefunctions and eigenvalues. While it may be difficult to reach a regime of strong interactions with usual methods, our approach is exact in the zero and infinite interaction limits. This method is an extension of Ref. [51], where only two basis states were used. It can be applied to systems in different trapping geometries, and the repulsive interactions can be tuned from weak to strong. To validate our method, we compare our results for spatial densities and momentum distributions to simulations of the continuum obtained with Matrix Product States (MPS) [52,53], as well as the known analytical results for two atoms in a harmonic trap [54].

Hamiltonian
We focus on a one-dimensional system of N identical fermions (majority) which interacts with a single distinct particle (minority) with the same mass in the presence of a trapping potential. The Hamiltonian can be written as where the potential V is the background potential (which in this paper is either a harmonic potential, V (x) = mω 2 x 2 /2, or the double well constructed in Appendix C with parameters as shown in Figure 1) and V 0,i is the interaction between the minority and majority, namely we have x 0 , x 1 , . . . , x N |V 0,i |x 0 , x 1 , . . . , In our case of a contact interaction, we have v(x) = δ(x). Since all particles have the same mass, we can interpret the single impurity as a fermionic atom in a different internal state than the remaining majority atoms. Such systems can be realized in the lab with ultracold Li atoms in different hyperfine states [4].
For simplicity we will work in units whereh = 1 and m = 1, but we will reinstate these units in all figures.

Variational method
Our variational method consists of constructing a suitable truncated basis of states. The basis states are constructed by using both the analytically known eigenstates at zero interactions as well as the analytically known solutions at infinite interaction.

States at zero interaction
The states at zero interaction are denoted by |φ i , for 0 ≤ i ≤ n − 1. Each state |φ i is defined by a collective index k (i) of N + 1 single particle states, namely Note that for the k N , different orders correspond to the same state up to a sign since they correspond to the majority particles, while the single k (i) 0 corresponds to the quantum number of the minority particle. We assume that k We define the totally antisymmetric state of a number of M (ordered) quantum states v by Let us further denote v[i, j, . . .] as the (ordered) set v with v i , v j , . . ., removed. This notation will be used throughout this article. A zero interaction state is then of the form

States at infinite interaction
At infinite interaction, the states are denoted as |ψ µ , for 0 ≤ µ ≤ m − 1. Note that the number of states at zero interaction, n, is not necessarily the same as the number of states at infinite interaction. Each state |ψ µ corresponds to a collective index q (µ) of N + 1 single particle states corresponding to a completely antisymmetric state Φ qµ built from the quantum states as well as a set of N + 1 coefficients a (µ) , Note that different orders of the q (µ) i correspond to the same state up to a sign, and we will assume that q (i) 0 < . . . < q (i) N , and we will assume that the a (µ) satisfy N l=0 (a (µ) The state at infinite interaction is then defined in the coordinate representation as where we denote M l as the set of points where x 0 , the coordinate for the minority particle, is smaller than exactly l of the x 1 , . . . , x N .
These exact solutions of the Hamiltonian (1) at g = +∞ [43,44], are orthogonal and properly normalized to unity provided that (7) holds. However, they are not orthogonal to the zero interaction eigenstates, and in Section 3.4 we will apply the Gram-Schmidt process to construct an orthonormal basis.

Overlaps between the zero and infinite interaction states
In this section we will compute the overlaps between the infinite interaction states ψ µ and the zero interaction states φ i , which is a necessary input for the construction in Section 3.4 and for computing the matrix elements and overlaps that include the states χ µ . We will denote the overlaps between the states at zero interaction and at infinite interaction by C iµ , where by convention i ∈ (0, . . . , n−1) corresponds to the index for the zero interaction states and µ ∈ (0, . . . , m − 1) corresponds to the index for the states at infinite interaction. The zero interaction state is on the form where Φ k (i) [0] is again the totally antisymmetric wave function Again, we use the notation where J removed. Recall that at g = +∞, an eigenstate can be specified by a sequence of N + 1 numbers a l , 0 ≤ l ≤ N , as well as a set q (µ) = (q (µ) 0 , . . . , q (µ) N ) of single particle quantum numbers, and is constructed by where M l (x 0 ) is the set where x 0 is smaller than exactly l of the x j with j ≥ 1. The overlaps are thus given by where where in the last step we have used the formula in Appendix A.
for l ≥ J. To compute the derivatives efficiently we evaluate the determinant for several values of , linearly spaced in (−1, 1), and fit a polynomial. We will encounter similar, but more involved, calculations when we compute the densities.

Constructing the basis
We will construct our basis by starting with the n states at zero interaction. We then add the infinite interaction states one by one, and orthonormalize after each added state. In other words, each state at infinite interaction |ψ µ corresponds to a state |χ µ , which is a linear combination of all the |φ i and the |χ ν with ν < µ such that it is orthogonal to all of these states. This procedure will be explained below.
We define C iµ = φ i |ψ µ and W µν = χ µ |ψ ν . Note that neither of these matrices are symmetric. The states |χ µ are then given by where the normalization constant is given by The W µν can be computed inductively. We first have that Then, for any µ, assuming knowledge of W ρσ where ρ ≤ µ − 1, we can compute W µν as where N µ is also given in terms of known W ρσ . Given the W µν , C iµ and N µ we now know our truncated basis (|φ 0 , . . . , |φ n−1 , |χ 0 , . . . , |χ m−1 ) ≡ (ξ 0 , . . . , ξ m+n−1 ). We will then express our Hamiltonian in this basis and numerically diagonalize it to find approximations to the eigenstates and energies.
Note that the construction of the basis ξ i in this subsection assumes that none of the states are linearly dependent. The only states at infinite interaction that can be linearly dependent with the states at zero interaction are the totally antisymmetric states, namely the states with coefficients a (µ) 0 = . . . = a (µ) N = 1, and thus one may have to exclude some of these states (or exclude some states at zero interaction) such that the final basis only contains states that are linearly independent. However, note that totally antisymmetric states are anyway already eigenstates to the Hamiltonian at finite g, and thus can be safely removed since they will be orthogonal to any non-trivial eigenstates. We will thus always exclude totally antisymmetric states at infinite interaction.

The Hamiltonian expressed in the basis
We will now express the Hamiltonian in the |ξ i basis by computing ξ i |H|ξ j . We will write the Hamiltonian as where V is the contact interaction between the majority and minority particles and H 0 is the Hamiltonian at zero interaction. We will treat these two terms individually.
For the zero interaction states, we have φ i |H 0 |φ j = δ ij E i , φ i |H 0 |χ µ = 0 due to the orthogonality propery of the basis and also ψ µ |H 0 |ψ ν = δ µν E µ . Let us define the quantity These can be computed recursively, starting with the known L µ0 . The matrix elements χ µ |H 0 |χ ν are then given by Which can also be calculated recursively starting with the known χ 0 |H 0 |χ ν . Now let us look at the interaction operator V . Note that ξ i |V |ψ µ = 0 since V is a contact interaction and x|ψ µ vanishes when x 0 = x j for 1 ≤ j ≤ N . Let us define V ij = φ i |V |φ j . We then have which can be computed recursively starting with V i,µ=0 . Given these quantities, the remaining matrix elements can be calculated as To compute the matrix elements V ij , note that the interaction operator between two particles, which we denote by V 2 , is defined as where x 1 , x 2 | is a position space eigenstate for two particles. Now consider some discrete set of single particle states |i , i = 0, 1, . . ., defined by the wavefunction by x|i = f i (x). The matrix element for the interaction operator between two particles in this basis is then Now we would like to know the matrix elements of the total interaction operator between two many-body states k (i) = [k where V 0l is the interaction operator between particle 0 (the impurity) and particle with index l. For two sets A and B with equal size, let us define |A − B| be the number of elements that only appear in A (or equivalently in only B). For a particle with index p, with p = 0, l, V 0l is diagonal and thus we have where l and l are the unique indices such that k This concludes our construction of the Hamiltonian H = H 0 +gV , and all that remains is diagonalizing the matrix ξ i |H|ξ j to find the energies and wavefunctions.

Observables
In this section we explain how to compute several important observables. They will all be computed starting with a specific eigenstate, which we denote by |Ψ , or Ψ(x 0 , . . . , x N ) = x 0 , . . . , x N |Ψ in the coordinate basis. This state is expressed as a linear combination of the zero interaction states and infinite interaction states which can be obtained easily given the expansion of Ψ in the basis {φ i , χ µ }. Note that the {φ i , ψ µ } is not an orthonormal basis.
We will start by computing the single particle density, which is the easiest observable presented in this section. The equations for the other observables are similar in nature but with varying extra degrees of complexity and subtleties, and thus it is recommended to understand the single particle density computation in detail first.

Single particle minority density matrix
The single particle density matrix is defined by integrating out the coordinates of the majority particles as where the integral is short for = The density matrix is useful since it is related to the momentum distribution by a simple Fourier transform. For just the particle density in coordinate space, we set x 0 = y 0 . We will comment on how the computations simplify for this special case.
The simplest term, namely between the zero interaction states is given by Here we are again using the notation that k[0] is equal to k with k 0 removed, namely the set {k 1 , . . . , k N }, and the Kronecker delta is thus equal to one if and only if the sets {k For the cross terms β i,µ , it will be useful to split up the integral into several regions, and we write where M l as the set of points where y 0 is smaller than exactly l of the x 1 , . . . , x N . We then split up the term β i,µ (x 0 , y 0 ) as The cross term is then given by where Φ represents a totally antisymmetric state. The matrix A J is defined by For the density where x 0 = y 0 , this works also for the infinite interaction terms, namely we can write where now the matrices A J,J and B J,J are defined by A J,J ab = However, when x 0 = y 0 , it is necessary to split the integral in more regions. We then write where M l,s is the region where x 0 and y 0 are smaller than exactly l respectively s of the x 1 , . . . , x N . We then split up the terms γ µ,ν (x 0 , y 0 ) as The term only involving infinite interaction states is then given by Now the matrices are defined as

Majority particle density matrix
The single particle majority density matrix is defined by integrating out the coordinate of the single minority particle and the coordinates of N − 1 of the majority particles. We thus write (38) where in all but the last line the integral is short for = ∞ −∞ dx 2 · · · ∞ −∞ dx N and in the last line we have separated out the dx 0 integral in all but the first term. In this case there are not many simplifications when x 1 = y 1 . The zero interaction term is given by The latter delta function means that this expression is zero unless the set k (i) with k J removed, are equal, in which case it is equal to one. For the cross term β maj i,µ (x 1 , y 1 ), we split it up into N terms β maj,l , corresponding to x 0 being smaller than exactly l of the x 2 , . . . , x N . We have where the integral is again short for = as being equal to a dx , and simplified the notation by assuming that S(I) is the (ordered) set S with the element with index I removed. Now let's consider the term γ. We now have the expression The matrices are now analogously defined, namely

Momentum distributions
The momentum distributions are obtained as a Fourier transform of the single particle density matrices. Let us denote the single particle density matrices by ρ min (x, y) and ρ maj (x, y) for the minority respectively majority species. The momentum distributions are then defined as (42) and have the same normalization as the coordinate space densities.

Minority-Majority correlation function
The last observable we will consider is the coordinate space minority-majority correlation function, which is defined as where the integral is short for = The computation of the minoritymajority correlation functions is very similar to the computation of the majority density; take the formulas for the majority density matrix and set x 1 = y 1 , and drop the integrals over x 0 . Dropping the x 0 integral in the β and γ terms is trivial, and the α term is given by

Examples
In this section we will consider several examples with different number of particles and in different background potentials. For two particles, we will compare with the analytic result for a harmonic trap. For more particles, we will compare with results computed using matrix product states as explained in Appendix D.
Recall from Section 3, that when describing the basis we will write [k When truncating the infinite basis (of both zero and infinite interaction states), some choice must be made on how this truncation should be done. In this paper, we will simply truncate when the number of particle excitations above the lowest state is above a certain threshold, namely we will construct the basis from all zero interaction states that satisfy a ≤ and all infinite interaction states that satisfy a q (µ) a − a q (0) a ≤ . We will refer to this cutoff as simply the basis threshold and denote it by in this section. For the harmonic oscillator potential, this is equivalent to an energy cutoff on the states, since the energy is linear in the quantum numbers, but for other potentials this is not the Double well potential Harmonic potential Figure 1: The two potentials used in this paper, the harmonic well with angular frequency ω and the double well described in Appendix C with parameters x 2 = −x 0 = 2 h ωm , ω 0 = ω 1 ≡ ω, ∆ 0 = 0, ∆ 1 = 1.5 ×hω and ∆ 2 = 0.8 ×hω.
case. We will leave a deeper investigation on the optimal selection of states for future work.
Note that in this paper we work in absolute coordinates and thus our energy spectra, such as figures 2 and 7, show all energy levels. This must be taken into account when comparing with for instance [51], where only the internal energy levels are shown and thus does not include excitations of the center of mass of the particles.

Two particles in a harmonic potential
In this section we will make detailed comparisons between the methods in this paper and the analytically known formula for two particles. A full derivation of the two particle system can be found in Appendix B. The basis is built from states at zero interaction and from states at infinite interaction. The states at zero interaction are specified by two quantum numbers, denoted [k 0 ; k 1 ]. There are no constraints on these two quantum numbers as we are dealing with two distinguishable particles. At infinite interaction, the states are built by taking a totally antisymmetric state, denoted by q = [q 0 , q 1 ] ∞ with q 0 > q 1 , but by multiplying with different coefficients a 0 and a 1 depending on the position space coordinates. In other words, the wave function is given by a 0 Φ q (x 0 , x 1 ) when x 0 > x 1 and a 1 Φ q (x 0 , x 1 ) when x 0 < x 1 where Φ is the totally antisymmetric state. A basis for such states is given by all antisymmetric states and the coefficients a (1) = [1, 1] and a (2) = [1, −1]. However, note that a = [1, 1] just corresponds to the totally antisymmetric state and is thus included among (a linear combination of) the zero interaction states. As was mentioned in the end of Section 3.4, it is important to exclude such linearly dependent states to avoid singular behaviour in the Gram-Schmidt orthogonalization process when constructing the basis. We will thus exclude the states with coefficients a (1) = [1, 1] and thus for two particles it is enough to specify a state at infinite interaction only by the quantum numbers q = [q 0 , q 1 ] and we leave the coefficients [1,

Energies
In Figure 2 we show the energy of the lowest six states computed using our variational approach and the analytic formula, for various values of the coupling g. The ground state interpolates between the state [0, 0] at g = 0 to the state [1, 0] ∞ at g = +∞ (with the implicit coefficients [1, −1]).The first and third excited states are totally antisymmetric states (unaffected by the interaction) with the quantum numbers [0, 1] and [0, 2].
As we can see from this plot, there is an agreement between the results, but it is difficult to appreciate exactly how well they agree. In Figure 3 we therefore plot the energy difference of the analytic result and the variational result for g = 1/2, 1, 2 for the ground state and one of the excited states for various basis sizes. As we can see, they agree to an extraordinary accuracy (note the logarithmic scale). Each data point corresponds to a basis constructed with basis threshold = 0, 2, 4, . . . The x-axis then shows the total size of the basis.
The reason why we look at the fourth excited state is that this is the first "nontrivial" excited state when computed using the analytical formula. As explained in Appendix B, the non-trivial part of the analytical derivation is computing the eigenstates of the relative motion Hamiltonian, and to get the full spectrum we also need to add the energy for the  center of mass Hamiltonian which is just a free harmonic oscillator. In the variational method, where we work directly in absolute coordinates, we automatically get all states. It turns out that the first excited state is just a totally antisymmetric state, the second excited state is just the first state plus a center of mass excitation, and the third excited state is then also just a totally antisymmetric state (actually the first excited state plus center of mass motion). The fourth excited state is then the first excited state which corresponds to a non-trivial eigenstate to the relative motion Hamiltonian and where the center of mass energy is zero.

Position space densities
In Figure 4 we show the position space density for the ground state at g = 1 compared to the analytical result. We see that when we only use two states in the basis there is a small discrepancy between the two methods, but when we use larger basis sizes the methods agree very well. A more detailed comparison can be seen in Figure 5, where we plot the density at three arbitrary values of x as a function of the basis size. We see that they agree exceptionally well with the analytical result. To compute the density from the analytical result, we need to perform an integral transforming from Jacobi coordinates to absolute coordinates (see Appendix B).

Momentum space densities
Finally we will compare the momentum densities, which is computed from the density matrix by equation 42. The comparison between the analytical and the variational methods is shown in Figure 6. Again, there is a discrepancy with the analytical result when only using 2 basis states, but when using 10 basis states the results agree very well.

Three particles
In this section we will consider the case of three particles, namely one single particle and two identical fermions. The states at zero interaction are now given by three quantum numbers [k    : Energies for the lowest seven states for the 2+1 system computed using the variational method. The result using the matrix product states method for increasing accuracy is shown in red dots, with the black star being the extrapolated value. The energies are computed using basis states with basis threshold 8. The ground state is computed also using basis threshold 0, 2, 4 and 6 to show the convergence, which are shown with dashed lines. Figure 7 shows the lowest seven energies for the 2+1 system with harmonic potential. We compare with the matrix product states (MPS) result at g = 1.0 for the ground state. Note that to obtain good agreement, we need to compute the energy for several numerical accuracies and then extrapolate the result. The MPS computations for the different accuracies are given by the red dots, and the extrapolated value is the black cross. The dashed lines are the ground state computed with the variational method using basis states with a basis threshold of = 0, 2, 4 and 6, and the solid lines are computed using a basis cutoff of = 8. These correspond to basis sizes of 1+2,7+8, 22+22, 50+46 and 95+82 respectively, where the first (second) number is the number of zero (infinite) interaction states the basis is constructed from. Our vatiational method easily gives us the energies of several states at many different values of g, which is one of the main advantages of the method compared to for example the MPS method where each computation only yields the energy and wavefunction at one particular interaction. : Position space density for the ground state at g m h 3 ω = 1 for the 2+1 system in the harmonic potential, for different basis sizes compared to the matrix product states method. The density profile localized in the center is the minority density and the other one is the majority density. Figure 8 and 9 shows the position space minority and majority density at g m h 3 ω = 1 for the 2+1 system for different basis sizes compared with MPS method. In Figure 8 we assume a harmonic potential, while in Figure 9 we consider the double-well geometry shown in Figure 1. In both cases we have good agreement with the MPS result. The computations are for energy cutoffs of 0, 2 and 4.

Position space densities
In Figure 10 we plot the integral of the squared difference of the densities for different basis sizes to better compare the convergence.

Momentum space densities
In Figure 11 we compare the momentum space densities at g m h 3 ω = 1 in the harmonic well with the MPS result. We see that they agree quite well already for the lowest possible number of basis states, and we again see that the discrepeancy goes to zero as we increase the basis size.      Figure 1 for the 6+1 system. The density profile localized in the left well is the minority density and the other one is the majority density.

Seven particles
In this section we study the 6+1 system. Figure 12 and Figure 13 shows the position space density profiles for the ground state in the double well potential for g m h 3 ω = 1 and g m h 3 ω = 10. We see that for large number of majority particles the system starts to look like a single impurity in a homogeneous bath. Moreover, when the interaction increases the minority particle density clearly gets deformed, which is reproduced with both methods, and we see that our method does work well both for intermediate and strong interactions. However, the discrepancy with the MPS result is clearly larger compared to the 2+1 system. Figure 14 shows the position space density profiles for the ground state in the harmonic trap for g m h 3 ω = 1 for the 10+1 system, and we see that the method still works well even for higher particle numbers.  Figure 1 for the 6+1 system. The density profile localized in the left well is the minority density and the other one is the majority density.

Conclusions
In this paper we explored a new method for studying strongly coupled one-dimensional systems where an impurity interacts with a background of identical fermions, a method that generalizes that of [51]. Our results compare well both with analytical methods for two particles and with numerical methods based on matrix product states. However, our method has the fundamental advantage of allowing calculations for arbitrary values of the interaction strength by only constructing the basis once. Generally, numerical approaches would require a full calculation for every value of the interaction strength. To compute the eigenstates and energies, we just need to change the interaction parameter g in the Hamiltonian before diagonalizing. Moreover, most numerical methods would perform worse the stronger the interaction strength is, but our method is exact at infinite interaction and thus works well both for small and strong interactions, with a peak of slower convergence at some intermediate interaction strength. Since our states are chosen such as to well approximate a state at finite interaction, the basis size is also relatively small and the computational power needed for the diagonalization is negligible. In particular, the method does not require sophisticated diagonalization algorithms or high performance computing tools, which is often the case for exact diagonalization methods. Note, moreover, that the matrix we diagonalize is not a particularly sparse matrix. Computing densities is more computationally heavy since they need to be evaluated for each position space coordinate separately, and thus scales linearly with the number of grid points. For majority densities, a double numerical integral must be performed which makes them heavier than the minority density. Furthermore, each position space coordinate for the minority density can be computed in parallel which can further decrease the computation time. While it is possible that our numerical techniques can be significantly improved, such as finding a way to efficiently parallelize the majority density calculations or optimize the grids on which we evaluate the numerical integrals, we will leave such investigations for future work. We also want to again stress that the time consuming calculations only have to be performed once for each chosen basis and we can then easily obtain the densities for any interaction strength g.
In this paper we only considered repulsive interactions. The method will not generalize trivially to attractive interaction since on the attractive sides, there are bound states that do not converge to any state when the interaction goes to infinite interaction strength, but instead diverges with infinitely negative energy. Such states would not be well approximated by a linear combination of exact solutions at infinite interaction and states at zero interaction. It should be possible to extend the basis considered in this paper to states that would also capture these bound states, similarly to what was done in [51], but we will leave the study of attractive interactions for future work.
We can also prove this by induction. Note that if we assume x > x and k = 0, the formula is the same as (46) if the upper integral limit is changed from ∞ to x and the same proof goes through. We will thus use this as a base case for our induction proof and thus assuming without loss of generality that x > x , we can prove the formula for l, k with k < l by assuming that it holds for k − 1, l − 1. Following the exact same reasoning as in the proof in A.1, we have Since we assumed that k < l and x > x , and the exact same proof can be done for k > l and x < x , formula (50) follows.

B Two particle system
In this section we review the analytical solution of two particles in a harmonic trap, with a delta function interaction [54]. The full Hamiltonian is where By introducing Jacobi coordinates x = (x 1 − x 2 )/ √ 2, p = (p 1 − p 2 )/ √ 2, X = (x 1 + x 2 )/ √ 2 and P = (p 1 + p 2 )/ √ 2 we can split this Hamiltonian into two parts, namely where H CM = X 2 /2 + P 2 /2 is just a harmonic oscillator corresponding to the center-ofmass motion, and The hard part, which will occupy most of this appendix, is solving for the eigenstates of H rel . The full set of eigenstates and eigenenergies are then obtained by tensor product with the eigenstates of H CM .
We will solve for the wavefuctions by first expanding in a harmonic oscillator basis. The Harmonic oscillator eigenfunctions are given by where H n are the Hermite polynomials. The energy is given by E n = n + 1/2. Let |Φ be an eigenstate for H rel . We have Solving for c n ≡ n|Φ and defining the quantity A = f n (0)c n we obtain Now multiplying both sides by f n (0) and summing over n, we can cancel A from both sides to obtain For the case where A = 0, for which we can not cancel it from both sides to obtain equation (59), see Appendix B.1. For the Hermite polynomials, we have H n (0) = 0 if n is odd, and H 2n (0) = (−1) n (2n!)/n!, which is the reason why we have omitted the odd terms. The wavefunction is given by a similar formula, namely It thus makes sense to treat these simultaneously, so let us define To compute this function, we use the following relation between Hermite polynomials and Laguerre polynomials H 2n (x) = (−1) n 2 2n n!L −1/2 n (x 2 ).
We thus obtain Now we use the integral representation to obtain Now we can recognize the generating function e −tx/(1−t) (1 − t) −α−1 = t n L α n (x) to obtain where we have used a standard representation for the confluent hypergeometric function U . At x = 0, we can use the relation U (−ν, 1/2, 0) = Γ(1/2)/Γ(1/2 − ν) = √ π/Γ(1/2 − ν), to obtain Thus for the energy, we must solve the equation For the wavefunction, we instead have (69) To find the normalization constant A we can consider the normalization constraint Defining ψ(x) = Γ (x)/Γ(x) and using the energy formula (68), we can simplify this to

B.1 Odd states
What we have obtained so far are all even parity states where the wavefunction in position space is an even function. The odd parity states are just odd harmonic oscillator states and they are unaffected by the interaction since they vanish at x = 0. These states would have A = 0 and thus the step to obtain equation (59) would be illegitimate.

B.2 Absolute coordinates
The full eigenstates are then obtained by also multiplying by the center of mass states. The complete wave function for H = H rel + H CM is given by where we have labeled all eigenstates of H rel (both even and odd) by Φ k for k = 0, 1, . . . and f n are just the standard harmonic oscillator wavefunctions. The energy is likewise E k,n = E Φ k + E n where E n = n + 1/2 is the nth harmonic oscillator energy.
To compare with the variational method in this paper, we would also like to compute the coordinate and momentum densities. Recall that x = (x 1 − x 2 )/ √ 2 and X = (x 1 + x 2 )/ √ 2. The single particle density matrix is just the square of the wavefunction in absolute coordinates, namely ρ(x 1 , x 2 ) = Φ 2 k,n ( and the density is thus given by The momentum density can then be obtained by C Wavefunctions and energies for a smooth double well potential In this appendix we give details on energies and wavefunctions of the double well potential. The double well potential is defined as where we require x 0 < x 1 < x 2 and x L < x R . Continuity of the potential as well as its derivatives at two points x L and x R implies the equations This system is uniquely solved for ω 1 , x 1 , x L and x R given the physically relevant quantities ω 0 , ω 2 , x 0 , x 2 , ∆ 0 , ∆ 1 and ∆ 2 . The solution is given by and then x L and x R are given by Extra care for these formulas must be taken when evaluating these expressions for ∆ 0 = ∆ 2 . In this case we have For the symmetric case (symmetric around x 1 = (x 0 + x 2 )/2) where we also have ω 0 = ω 2 , we have We can compute an upper limit on the parameter ∆ 1 . The highest value is the value such that ω 1 = ∞, namely we have the more well known double well potential which has a discontinuous derivative between the wells. For such a potential the discontinuity is at the intersection of the left and right wells, namely we solve which results in the solution Then the upper limit of ∆ 1 is given by ∆ 1,max = 1 2 (x M − x 0 ) 2 ω 2 0 + ∆ 0 . We will now work out the wavefunctions and energies. We will work in units wherē h = 1 and m = 1 (the mass of the particle) and for simplicity and we will define ν 0 and ν 2 by E = ω 0 (ν 0 + 1 2 ) + ∆ 0 = ω 2 (ν 2 + 1 2 ) + ∆ 2 . The eigenfunctions are now uniquely given by for x > x R and for some constants C 0 , C 2 (this follows since these are the only solutions with the correct falloffs at x → ±∞). The function D is the parabolic cylinder function given by where 1 F 1 is the confluent hypergeometric function. Note that this function is a linear combination of the two linearly independent solution of the Schrödinger equation in a harmonic well, and the relative coefficient has been fixed by requiring falloff at infinity. In the intermediate region we need to solve the Schrödinger equation for an inverted harmonic well. It can be showed that the solution then is and where we have parametrized the energy as E = ω 1 (ν 1 + 1 2 ) + ∆ 1 (which we recall is also equal to ω 0 (ν 0 + 1 2 ) + ∆ 0 = ω 2 (ν 2 + 1 2 ) + ∆ 2 ). Despite the complex arguments, these are real functions. These solutions should now be glued smoothly across the points x L and x R such that ψ and ψ are continuous. To simplify the equations, we will define r = ω 2 /ω 1 , R = ω/ω 1 , ∆ =hω 1 δ, C =hω 1 c and we work in units where µω 1 /h = 1. This gives the equations √ If we are given ν 0 , ν 1 , ν 2 (which are all determined by the energy E), this is a linear system of equations for C 0 , C 2 , C 1 , C 1 . For this system to have a non-trivial solution, the determinant of the corresponding matrix must vanish and this condition is what determines the energy (or equivalently the parameters ν 0 , ν 1 , ν 2 ). This system of equations, supplemented with normalization of the wave function, then fixes all the constants C 0 , C 2 , C 1 , C 1 . In general, if we piece together N different quadratic (or other analytically solvable) potentials, the energy will be obtained by solving the equation resulting from enforcing zero determinant of a 2(N − 1) × 2(N − 1) matrix.

D Matrix Product States
Throughout this work we compare our analytical method with simulations performed with Matrix Product States (MPS), using the Open Source MPS (OSMPS) libraries [52]. In these calculations, we employ the Hubbard model as an approximation to the continuum in order to obtain static properties of a fermionic polaron system. Thus the spinful lattice Hamiltonian is written as H = −t j,σ (c † j+1,σ c j,σ + H.c.) + U j n j,↑ n j,↓ + j,σ j n j,σ , where c † and c are the creation and annihilation operators, respectively, t is the hopping parameter and U denotes the strength of the on-site interactions between fermions with different spin projections. We denote the internal states as |↑ for the background fermions and |↓ for the impurity. Since we consider only a single |↓ fermion, we have naturally j n j,↓ = 1, with j n j,↑ also being normalized to the number of background fermions. We include additionally the trapping potential as the position-dependent j parameter.
We simulate the continuum by taking a total of L = 256 sites. We thus obtain a lattice spacing a = l/L where l is the total length assumed for the trapping potential. The hopping parameter is related to the kinetic term in the continuum as t = 1/(2ma 2 ), where m is the atomic mass, which we take to be 1. The continuum and discrete interaction parameters are related as U = g/a. To obtain matching energies, we must include an additional term in the Hamiltonian given by j 1/a 2 . In some cases, to improve the accuracy we compute the results for several increasing values of L and then extrapolate to a final value using a function of the form f (L) = A + B/L + C/L 2 .

E Polynomial interpolation for computing determinants
At several stages in the technique used in this paper we have to compute derivatives of determinants of the form ∂ i D( )| =0 = ∂ i 1 i! det(M ( ))| =0 , where M ( ) is some n × n matrix and i = 0, . . . , n. We evaluate these derivatives by computing the function D( ) on n + 1 values with i = −1 + 2i/n, i = 0, . . . , n, and then fitting a polynomial to these values and extracting the coefficients. These coefficients can be obtained by multiplying the vector D( i ) with the inverse of the matrix K ij ≡ j i .