Systematic interpolatory ansatz for one-dimensional polaron systems

TL;DR

This paper introduces a new variational approach using a truncated basis to accurately study one-dimensional Fermi polaron systems at any interaction strength, improving computational efficiency and accuracy.

Contribution

The authors develop a systematic interpolatory ansatz that combines zero and infinite repulsion states for better modeling of 1D polaron systems.

Findings

Accurately describes polaron energies, densities, and correlations at finite interactions.

Validates the method against matrix product states and analytical solutions.

Applicable to various trapping potentials and particle numbers.

Abstract

We explore a new variational principle for studying one-dimensional quantum systems in a trapping potential. We focus on the Fermi polaron problem, where a single distinguishable impurity interacts through a contact potential with a background of identical fermions. We can accurately describe this system at arbitrary finite repulsion by constructing a truncated basis containing states at both the limits of zero and infinite repulsion. We show how to construct this basis and how to obtain energies, density matrices and correlation functions, and provide results both for a harmonic well and a double well for various particle numbers. The results are compared both with matrix product states methods and with the analytical result for two particles in a harmonic well.

Equations

$$H=\sum _{i=0}^N\left(\frac{{\mathbf{p}}_i^2}{2}+V({\mathbf{x}}_i)\right)+g\sum _{i=1}^N{V}_{0,i}({\mathbf{x}}_0-{\mathbf{x}}_i),$$

$${\vec{k}}_i=[k_0^{(i)};k_1^{(i)},\dots ,k_N^{(i)}].$$

$$|{\Phi }_{\vec{v}}⟩=\frac{1}{\sqrt{M!}}\sum _{\sigma }\text{sign}(\sigma )|f_{v_{\sigma (1)}}⟩\cdots |f_{v_{\sigma (M)}}⟩.$$

$$|{\varphi }_i⟩=|f_{k_0^{(i)}}⟩|{\Phi }_{{\vec{k}}_i[0]}⟩.$$

$${\vec{q}}_{\mu }=[q_0^{(\mu )},\dots ,q_N^{(\mu )}],$$

$${\vec{a}}_{\mu }=[a_0^{(\mu )},\dots ,a_N^{(\mu )}].$$

$$\sum _i(a_i^{(\mu )}{)}^2=N+1,\sum _i{a}_i^{(\mu )}{a}_i^{(\nu )}=0,\mu \ne \nu .$$

$${\psi }_{\mu }(x_0,\dots ,x_N)=a_l^{(\mu )}{\Phi }_{{\vec{q}}_{\mu }}(x_0,\dots ,x_N),\text{when }(x_0,\dots ,x_N)\in {\mathcal{M}}_l,$$

$${\varphi }_i(x_0,\dots ,x_N)=f_{k_0}(x_0){\Phi }_{{\vec{k}}_i[0]}(x_1,\dots ,x_N)$$

$${\Phi }_{{\vec{k}}_i[0]}=\frac{1}{\sqrt{N!}}\sum _{\pi }sign(\pi )f_{k_{\pi (1)}^{(i)}}\cdots f_{k_{\pi (N)}^{(i)}}.$$

$${\psi }_{\mu }={\alpha }_l{\Phi }_{{\vec{q}}_{\mu }}(x_0,x_1,\dots ,x_N),x_1,\dots ,x_N\in {\mathcal{M}}_l(x_0),$$

$$C_{i\mu }=\sum _{l=0}^N{\alpha }_l{I}_l,$$

$$I_l$$

$$=\frac{1}{\sqrt{N+1}}\sum _{j=0}^N(-1{)}^j{\int }_{-\infty }^{\infty }\mathrm{d}{x}_0{f}_{k_0^{(i)}}(x_0)f_{q_j^{(\mu )}}(x_0){\int }_{{\mathcal{M}}_l(x_0)}{\Phi }_{{\vec{k}}_i[0]}{\Phi }_{{\vec{q}}_{\mu }[j]}$$

$$=\frac{1}{l!\sqrt{N+1}}\sum _{j=0}^N(-1{)}^j{\int }_{-\infty }^{\infty }\mathrm{d}{x}_0{f}_{k_0^{(i)}}(x_0)f_{q_j^{(\mu )}}(x_0){\partial }_{ϵ}^l\det (A^j+ϵB^j{)}_{ϵ=0},$$

$$|{\chi }_{\mu }⟩=N_{\mu }\left(|{\psi }_{\mu }⟩-\sum _{i=0}^n{C}_{i\mu }|{\varphi }_i⟩-\sum _{\rho =0}^{\mu -1}{W}_{\rho \mu }|{\chi }_{\rho }⟩\right),$$

$$N_{\mu }={\left(1-\sum _{i=0}^n{C}_{i\mu }^2-\sum _{\rho =0}^{\mu -1}{W}_{\rho \mu }^2\right)}^{-1/2}.$$

$$W_{0\nu }=N_0\left({\delta }_{0\nu }-\sum _{i=0}^n{C}_{i0}{C}_{i\nu }\right)$$

$$W_{\mu \nu }=N_{\mu }\left({\delta }_{\mu \nu }-\sum _{i=0}^n{C}_{i\mu }{C}_{i\nu }-\sum _{\rho =0}^{\mu -1}{W}_{\rho \mu }{W}_{\rho \nu }\right)$$

$$H=H_0+gV,$$

$$L_{\mu \nu }=⟨{\psi }_{\mu }|H|{\chi }_{\nu }⟩=N_{\nu }\left(E_{\mu }{\delta }_{\mu \nu }-\sum _{i=0}^n{C}_{i\mu }{C}_{i\nu }{E}_i-\sum _{\rho =0}^{\nu -1}{W}_{\rho \nu }{L}_{\mu \rho }\right)$$

$$⟨{\chi }_{\mu }|H_0|{\chi }_{\nu }⟩=N_{\mu }\left(L_{\mu \nu }-\sum _{\rho =0}^{\mu -1}{W}_{\rho \mu }⟨{\chi }_{\rho }|H_0|{\chi }_{\nu }⟩\right)$$

$$V_{i\mu }=N_{\mu }\left(-\sum _{j=0}^n{V}_{ij}{C}_{j\mu }-\sum _{\sigma =0}^{\mu -1}{W}_{\sigma \mu }{V}_{i\sigma }\right)$$

$$V_{\mu \nu }=N_{\mu }\left(-\sum _{j=0}^n{C}_{j\mu }{V}_{j\nu }-\sum _{\sigma =0}^{\mu -1}{W}_{\sigma \mu }{V}_{\sigma \nu }\right).$$

$$⟨x_1,x_2|V_2|y_1,y_2⟩=\delta (x_1-y_1)\delta (y_1-y_2)\delta (x_2-y_2)$$

$$⟨n_1,n_2|V_2|m_1,m_2⟩={\int }_{-\infty }^{\infty }{f}_{n_1}(x)f_{n_2}(x)f_{m_1}(x)f_{m_2}(x)$$

$$⟨n|V|m⟩=0,$$

$$⟨n|V|m⟩=⟨n_0,n_i|V_2|m_0,m_j⟩(-1{)}^{i-j},$$

$$⟨n|V|m⟩=\sum _{i=1}^N⟨n_0,n_i|V_2|m_0,m_i⟩.$$

$$|\Psi ⟩=\sum _{i=0}^n{C}_i|{\varphi }_i⟩+\sum _{\mu =0}^m{D}_{\mu }|{\psi }_{\mu }⟩,$$

$$\rho (x_0,y_0)=$$

Full text

1. 1 Introduction

2. 2 Hamiltonian

3. 3 Variational method

4. 3.1 States at zero interaction

5. 3.2 States at infinite interaction

6. 3.3 Overlaps between the zero and infinite interaction states

7. 3.4 Constructing the basis

8. 3.5 The Hamiltonian expressed in the basis

9. 4 Observables

10. 4.1 Single particle minority density matrix

11. 4.2 Majority particle density matrix

12. 4.3 Momentum distributions

13. 4.4 Minority-Majority correlation function

14. 5 Examples

15. 5.1 Two particles in a harmonic potential

16. 5.1.1 Basis

17. 5.1.2 Energies

18. 5.1.3 Position space densities

19. 5.1.4 Momentum space densities

20. 5.2 2+1

21. 5.2.1 Energies

22. 5.2.2 Position space densities

23. 5.2.3 Momentum space densities

24. 5.3 6+1

25. 6 Conclusions

26. 7 Acknowledgements

27. A Some integral formulas

28. A.1 Integral 1

29. A.2 Integral 2

30. B Two particle system

31. B.1 Odd states

32. B.2 Absolute coordinates

33. C Wavefunctions and energies for a smooth double well potential

34. D Matrix Product States

35. E Polynomial interpolation for computing determinants

1 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] and mathematical physics [5, 6, 7] 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 [8, 9, 10, 11] and trapping geometries opens up the possibility of using such experiments as quantum simulators for a multitude of interesting models [12], even beyond usual condensed matter models [13, 14, 15].

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 context of condensed matter, such systems can present interesting phenomena such as the Kondo effect [16] and the orthogonality catastrophe [17]. Theoretical and experimental studies with ultracold atomic setups have extensively explored both the bosonic [18, 19, 20, 21, 22, 23, 24, 25, 26, 27] and fermionic [28, 29, 30, 31] 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 [32, 33], which is exactly solvable through the Bethe ansatz approach [34]. Other approaches have later generalized the study of static properties to mixed fermionic systems in harmonic potentials [35, 36, 37, 38, 39, 40, 41]. On the dynamical side, impurity models have been shown to present exotic effects such as Bloch oscillations [42] and quantum flutter [43, 44].

In this work we present an original way to obtain the static properties of a Fermi polaron system, 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-Schmidth 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. [45], 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).

2 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

[TABLE]

where the potential $V$ is the background potential (which in this paper is either a harmonic potential or a double well) and $V_{0,i}$ is the interaction between the minority and majority, namely we have $\langle x_{0},x_{1},\ldots,x_{N}|V_{0,i}|x_{0}^{\prime},x_{1}^{\prime},\ldots,x_{N}^{\prime}\rangle=v(x_{0}-x_{i})\delta(x_{0}-x_{0}^{\prime})\cdots\delta(x_{N}-x_{N}^{\prime})$. In our case of a contact interaction, we have $v(x)=\delta(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].

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

3.1 States at zero interaction

The states at zero interaction are denoted by $|\phi_{i}\rangle$, for $0\leq i\leq n$. Each state $|\phi_{i}\rangle$ is defined by a collective index $\vec{k}_{i}$ of $N+1$ single particle states, namely

[TABLE]

Note that for the $k_{1}^{(i)},\ldots,k_{N}^{(i)}$, different orders correspond to the same state up to a sign since they correspond to the majority particles, while the single $k_{0}^{(i)}$ corresponds to the quantum number of the minority particle. We assume that $k_{1}^{(i)}<\ldots<k_{N}^{(i)}$. We define the totally antisymmetric state of a number of $M$ (ordered) quantum states $\vec{v}$ by

[TABLE]

Let us further denote $\vec{v}[i,j,\ldots]$ as the (ordered) set $\vec{v}$ with $v_{i}$, $v_{j}$, $\ldots$, removed. This notation will be used throughout this article. The zero interaction state is then given by

[TABLE]

3.2 States at infinite interaction

At infinite interaction, the states are denoted as $|\psi_{\mu}\rangle$, for $0\leq\mu\leq m$. 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 $|\psi_{\mu}\rangle$ corresponds to a collective index $\vec{q}_{\mu}$ of $N+1$ single particle states corresponding to a completely antisymmetric state $\Phi_{\vec{q}_{\mu}}$ built from the quantum states

[TABLE]

as well as a set of $N+1$ coefficients $\vec{a}_{\mu}$,

[TABLE]

Note that different orders of the $q_{i}^{(\mu)}$ correspond to the same state up to a sign, and we will assume that $q_{0}^{(i)}<\ldots<q_{N}^{(i)}$, and we will assume that the $\vec{a}_{\mu}$ satisfy

[TABLE]

The state at infinite interaction is then defined in the coordinate representation as

[TABLE]

where we denote $\mathcal{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},\ldots,x_{N}$.

These exact solutions of the Hamiltonian (1) at $g=+\infty$ [37, 38], 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.

3.3 Overlaps between the zero and infinite interaction states

In this section we will compute the overlaps between the infinite interaction states $\psi_{\mu}$ and the zero interaction states $\phi_{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 $\chi_{\mu}$. We will denote the overlaps between the states at zero interaction and at infinite interaction by $C_{i\mu}$, where by convention $i\in(0,\ldots,n-1)$ corresponds to the index for the zero interaction states and $\mu\in(0,\ldots,m-1)$ corresponds to the index for the states at infinite interaction. The zero interaction state is on the form

[TABLE]

where $\Phi_{\vec{k}_{i}[0]}$ is again the totally antisymmetric wave function

[TABLE]

Again, we use the notation where $\vec{k}[j]$ is the set $(k_{0},\ldots,k_{N})$ with $k_{j}$ removed. Recall that at $g=+\infty$, an eigenstate can be specified by a sequence of $N+1$ numbers $\alpha_{j}$, $0\leq j\leq N$, as well as a set $\vec{q}_{\mu}=(q_{0},\ldots,q_{N})$ of single particle quantum numbers, and is constructed by

[TABLE]

where $\mathcal{M}_{l}(x_{0})$ is the set where $x_{0}$ is larger than exactly $l$ of the $x_{j}$ with $j\geq 1$. The overlaps is thus given by

[TABLE]

where

[TABLE]

where in the last step we have used the formula in Appendix A.1. The matrix $A^{j}$ is defined by $A^{j}_{kl}=\int_{-\infty}^{x_{0}}f_{k^{(i)}_{k+1}}f_{q^{(\mu)}_{l}}$ for $l<j$ and $A^{j}_{kl}=\int_{-\infty}^{x_{0}}f_{k^{(i)}_{k+1}}f_{q^{(\mu)}_{l+1}}$ for $l\geq j$ while $B^{j}$ is defined by $B^{j}_{kl}=\int_{x_{0}}^{\infty}f_{k^{(i)}_{k+1}}f_{q^{(\mu)}_{l}}$ for $l<j$ and $B^{j}_{kl}=\int_{x_{0}}^{\infty}f_{k^{(i)}_{k+1}}f_{q^{(\mu)}_{l+1}}$ for $l\geq j$. To compute the derivatives efficiently we evaluate the determinant for several values of $\epsilon$, linearly spaced in $(-1,1)$, and fit a polynomial. We will encounter similar, but more involved, calculations when we compute the densities.

3.4 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 $|\psi_{\mu}\rangle$ corresponds to a state $|\chi_{\mu}\rangle$, which is a linear combination of all the $|\phi_{i}\rangle$ and the $|\chi_{\nu}\rangle$ with $\nu<\mu$ such that it is orthogonal to all of these states. This procedure will be explained below.

We define $C_{i\mu}=\langle\phi_{i}|\psi_{\mu}\rangle$ and $W_{\mu\nu}=\langle\chi_{\mu}|\psi_{\nu}\rangle$. Note that neither of these matrices are symmetric. The states $|\chi_{\mu}\rangle$ are then given by

[TABLE]

where the normalization constant is given by

[TABLE]

The $W_{\mu\nu}$ can be computed inductively. We first have that

[TABLE]

Then, for any $\mu$, assuming knowledge of $W_{\rho\sigma}$ where $\rho\leq\mu-1$, we can compute $W_{\mu\nu}$ as

[TABLE]

where $N_{\mu}$ is also given in terms of known $W_{\rho\sigma}$. Given the $W_{\mu\nu}$, $C_{i\mu}$ and $N_{\mu}$ we now know our truncated basis $(|\phi_{0}\rangle,\ldots,|\phi_{n-1}\rangle,|\chi_{0}\rangle,\ldots,|\chi_{m-1}\rangle)\equiv(\alpha_{0},\ldots,\alpha_{m+n-1})$. We will then express our Hamiltonian in this basis and numerically diagonalize it to find approximations to the eigenstates and energies.

3.5 The Hamiltonian expressed in the basis

We will now express the Hamiltonian in the $|\alpha_{i}\rangle$ basis by computing $\langle\alpha_{i}|H|\alpha_{j}\rangle$. We will write the Hamiltonian as

[TABLE]

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 $\langle\phi_{i}|H_{0}|\phi_{j}\rangle=\delta_{ij}E_{i}$, $\langle\phi_{i}|H_{0}|\chi_{\mu}\rangle=0$ due to the orthogonality propery of the basis and also $\langle\psi_{\mu}|H_{0}|\psi_{\nu}\rangle=\delta_{\mu\nu}E_{\mu}$. Let us define the quantity

[TABLE]

These can be computed recursively, starting with the known $L_{\mu 0}$. The matrix elements $\langle\chi_{\mu}|H_{0}|\chi_{\nu}\rangle$ are then given by

[TABLE]

Which can also be calculated recursively starting with the known $\langle\chi_{0}|H_{0}|\chi_{\nu}\rangle$.

Now let us look at the interaction operator $V$. Note that $\langle\alpha_{i}|V|\psi_{\mu}\rangle=0$ since $V$ is a contact interaction and $\langle x|\psi_{\mu}\rangle$ vanishes when $x_{0}=x_{j}$ for $1\leq j\leq N$. Let us define $V_{ij}=\langle\phi_{i}|V|\phi_{j}\rangle$. We then have

[TABLE]

which can be computed recursively starting with $V_{i,\mu=0}$. Given these quantities, the remaining matrix elements can be calculated as

[TABLE]

To compute the matrix elements $V_{ij}$, note that the interaction operator between two particles is defined as

[TABLE]

Thus the matrix elements between some discrete set of eigenstates given by $\langle x|n\rangle=f_{n}(x)$ are

[TABLE]

Now we would like to know the matrix elements of the total interaction operator between two many-body states $n=(n_{0};n_{1},\ldots,n_{N})$ and $m=(m_{0};m_{1},\ldots,m_{N})$ (and we again denote $n(0)=(n_{1},\ldots,n_{N})$, and we assume $n_{1}>\ldots>n_{N}$). The total interaction operator is given as $V=\sum_{j=1}^{N}V_{0j}$, where $V_{0j}$ is the interaction operator between particle 0 (the impurity) and particle with index $j$. 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$). Since $V_{0j}$ is diagonal in all other particles with index $i\neq 0,j$, we obtain that

[TABLE]

if $|n(0)-m(0)|>1$. If $|n(0)-m(0)|=1$, we obtain

[TABLE]

where $i$ and $j$ are the unique indices such that $n_{i}\neq n_{j}$ and $n(0)(i)=n(0)(j)$. If $n(0)=m(0)$, we obtain

[TABLE]

This concludes our construction of the Hamiltonian $H=H_{0}+gV$, and all that remains is diagonalizing the matrix $\langle\alpha_{i}|H|\alpha_{j}\rangle$ to find the energies and wavefunctions.

4 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 $|\Psi\rangle$, or $\Psi(x_{0},\ldots,x_{N})=\langle x_{0},\ldots,x_{N}|\Psi\rangle$ in the coordinate basis. This state is expressed as a linear combination of the zero interaction states and infinite interaction states

[TABLE]

which can be obtained easily given the expansion of $\Psi$ in the basis $\{\phi_{i},\chi_{\mu}\}$. Note that the $\{\phi_{i},\psi_{\mu}\}$ 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.

4.1 Single particle minority density matrix

The single particle density matrix is defined by integrating out the coordinates of the majority particles as

[TABLE]

where the integral is short for $\int=\int_{-\infty}^{\infty}dx_{1}\cdots\int_{-\infty}^{\infty}dx_{N}$. 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

[TABLE]

Here we are again using the notation that $\vec{k}[0]$ is equal to $\vec{k}$ with $k_{0}$ removed, namely the set $\{k_{1},\ldots,k_{N}\}$, and the Kronecker delta is thus equal to one if and only if the sets $\{k_{1}^{(i)},\ldots,k_{N}^{(i)}\}$ and $\{k_{1}^{(j)},\ldots,k_{N}^{(j)}\}$ are the same.

For the cross terms $\beta_{i,\mu}$, it will be useful to split up the integral into several regions, and we write

[TABLE]

where $\mathcal{M}_{l}$ as the set of points where $y_{0}$ is smaller than exactly $l$ of the $x_{1},\ldots,x_{N}$. We then split up the term $\beta_{i,\mu}(x_{0},y_{0})$ as

[TABLE]

The cross term is then given by

[TABLE]

where $\Phi$ represents a totally antisymmetric state. The matrix $A^{J}$ is defined by $A^{J}_{ab}=\int_{-\infty}^{x_{0}}f_{k^{(i)}_{a+1}}f_{q^{(\mu)}_{b}}$ for $b<J$ and $A^{J}_{ab}=\int_{-\infty}^{x_{0}}f_{k^{(i)}_{a+1}}f_{q^{(\mu)}_{b+1}}$ for $b\geq J$ while $B^{J}$ is defined by $B^{J}_{ab}=\int_{x_{0}}^{\infty}f_{k^{(i)}_{a+1}}f_{q^{(\mu)}_{b}}$ for $b<J$ and $B^{J}_{ab}=\int_{x_{0}}^{\infty}f_{k^{(i)}_{a+1}}f_{q^{(\mu)}_{b+1}}$ for $b\geq J$. Here $a$ and $b$ take the values $0,\ldots,N-1$. For a derivation of this equation see A.1.

For the density where $x_{0}=y_{0}$, this works also for the infinite interaction terms, namely we can write

[TABLE]

where now the matrices $A^{J,J^{\prime}}$ and $B^{J,J^{\prime}}$ are defined by $A^{J,J^{\prime}}_{ab}=\int_{-\infty}^{x_{0}}f_{q^{(\mu)}_{a+\sigma}}f_{q^{(\nu)}_{b+\delta}}$ and $B^{J,J^{\prime}}_{ab}=\int_{x_{0}}^{\infty}f_{q^{(\mu)}_{a+\sigma}}f_{q^{(\nu)}_{b+\delta}}$ where $\sigma=0$ for $a<J$, $\sigma=1$ for $a\geq J$, $\delta=0$ for $b<J^{\prime}$ and $\delta=1$ for $b\geq J^{\prime}$. Here $a$ and $b$ take the values $0,\ldots,N-1$ and we refer again to A.1 for a derivation of the determinant formulas.

However, when $x_{0}\neq y_{0}$, it is necessary to split the integral in more regions. We then write

[TABLE]

where $\mathcal{M}_{l,s}$ is the region where $x_{0}$ and $y_{0}$ are smaller than exactly $l$ respectively $s$ of the $x_{1},\ldots,x_{N}$. We then split up the terms $\gamma_{\mu,\nu}(x_{0},y_{0})$ as

[TABLE]

The term only involving infinite interaction states is then given by

[TABLE]

Now the matrices are defined as $A_{ab}=\int_{-\infty}^{\text{min}(x,x^{\prime})}f_{q^{(\mu)}_{a+\sigma}}(x^{\prime\prime})f_{q^{(\nu)}_{b+\delta}}(x^{\prime\prime})\mathrm{d}x^{\prime\prime}$, $B_{ab}=\int_{\text{min}(x,x^{\prime})}^{\text{max}(x,x^{\prime})}f_{q^{(\mu)}_{a+\sigma}}(x^{\prime\prime})f_{q^{(\nu)}_{b+\delta}}(x^{\prime\prime})\mathrm{d}x^{\prime\prime}$ and $C_{ab}=\int_{\text{max}(x,x^{\prime})}^{\infty}f_{q^{(\mu)}_{a+\sigma}}(x^{\prime\prime})f_{q^{(\nu)}_{b+\delta}}(x^{\prime\prime})\mathrm{d}x^{\prime\prime}$ where $\sigma=0$ for $a<J$, $\sigma=1$ for $a\geq J$, $\delta=0$ for $b<J^{\prime}$ and $\delta=1$ for $b\geq J^{\prime}$. The indices $a$ and $b$ take the values $0,\ldots,N-1$. See A.2 for a derivation of this formula.

4.2 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

[TABLE]

where in all but the last line the integral is short for $\int=\int_{-\infty}^{\infty}dx_{2}\cdots\int_{-\infty}^{\infty}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

[TABLE]

The latter delta function means that this expression is zero unless the set $\vec{k}^{(i)}$ with $k_{0}^{(i)}$ and $k_{I}^{(i)}$ removed and the set $\vec{k}^{(j)}$ with $k_{0}^{(j)}$ and $k_{J}^{(j)}$ removed, are equal, in which case it is equal to one. For the cross term $\beta^{\text{maj}}_{i,\mu}(x_{1},y_{1})$, we split it up into $N$ terms $\beta^{\text{maj},l}$, corresponding to $x_{0}$ being smaller than exactly $l$ of the $x_{2},\ldots,x_{N}$. We have

[TABLE]

where the integral is again short for $\int=\int_{-\infty}^{\infty}dx_{0}dx_{2}\cdots\int_{-\infty}^{\infty}dx_{N}$. The sign $sgn(J,J^{\prime})$ is defined as $(-1)^{J+J^{\prime}}$ if $J^{\prime}<J$ and $-(-1)^{J+J^{\prime}}$ otherwise. We have defined $a^{(\mu)}(l,x_{0},y_{1})$ as being equal to $a^{(\mu)}_{l}$ if $x_{0}<y_{1}$ and equal to $a^{(\mu)}_{l+1}$ otherwise. This formula does not simplify much for the density where $x_{1}=y_{1}$. We have defined the matrices $A^{0,I;J,J^{\prime}}_{ab}=\int_{-\infty}^{x_{0}}f_{[\vec{k}^{(i)}(0)(I)]_{a}}(x^{\prime})f_{[\vec{q}^{(\mu)}(J)(J^{\prime})]_{b}}(x^{\prime})dx^{\prime}$ and $B^{0,I;J,J^{\prime}}_{ab}=\int_{x_{0}}^{\infty}f_{[\vec{k}^{(i)}(0)(I)]_{a}}(x^{\prime})f_{[\vec{q}^{(\mu)}(J)(J^{\prime})]_{b}}(x^{\prime})dx^{\prime}$, and simplified the notation by assuming that $\vec{S}(I)$ is the (ordered) set $S$ with the element with index $I$ removed.

Now let’s consider the term $\gamma$. We now have the expression

[TABLE]

The matrices are now analogously defined, namely

$A^{I,I^{\prime};J,J^{\prime}}_{ab}=\int_{-\infty}^{x_{0}}f_{[\vec{q}^{(\mu)}(I)(I^{\prime})]_{a}}(x^{\prime})f_{[\vec{q}^{(\nu)}(J)(J^{\prime})]_{b}}(x^{\prime})dx^{\prime}$ and

$B^{I,I^{\prime};J,J^{\prime}}_{ab}=\int_{x_{0}}^{\infty}f_{[\vec{q}^{(\mu)}(I)(I^{\prime})]_{a}}(x^{\prime})f_{[\vec{q}^{(\nu)}(J)(J^{\prime})]_{b}}(x^{\prime})dx^{\prime}$.

4.3 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 $\rho_{\text{min}}(x,y)$ and $\rho_{\text{maj}}(x,y)$ for the minority respectively majority species. The momentum distributions are then defined as

[TABLE]

and have the same normalization as the coordinate space densities.

4.4 Minority-Majority correlation function

The last observable we will consider is the coordinate space minority-majority correlation function, which is defined as

[TABLE]

where the integral is short for $\int=\int_{-\infty}^{\infty}dx_{2}\cdots\int_{-\infty}^{\infty}dx_{N}$. The computation of the minority-majority 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 $\beta$ and $\gamma$ terms is trivial, and the $\alpha$ term is given by

[TABLE]

5 Examples

5.1 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. As explained in Section 3, when describing the basis we will write $[k_{0};k_{1},\ldots,k_{N}]$ to denote a zero interaction state with the single particle in state $k_{0}$ and the majority particles in the antisymmetric state with quantum numbers $k_{1}>\ldots>k_{N}$. The states at infinite interaction will be denoted by two sets of numbers, $[q_{0},\ldots,q_{N}]_{\infty}$ and $[a_{0},\ldots,a_{N}]$, such that the wavefunction is a totally antisymmetric wavefunction built from $q_{0}>\ldots>q_{N}$ and which is multiplied with the coefficient $a_{i}$ if $x_{0}$ is greater than exactly $i$ of $x_{1},\ldots,x_{N}$.

5.1.1 Basis

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 $\vec{q}=[q_{0},q_{1}]_{\infty}$ 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}\Phi_{\vec{q}}(x_{0},x_{1})$ when $x_{0}>x_{1}$ and $a_{1}\Phi_{\vec{q}}(x_{0},x_{1})$ when $x_{0}<x_{1}$ where $\Phi$ is the totally antisymmetric state. A basis for such states is given by all antisymmetric states and the coefficients $\vec{a}^{(1)}=[1,1]$ and $\vec{a}^{(2)}=[1,-1]$. However, note that $\vec{a}=[1,1]$ just corresponds to the totally antisymmetric state and is thus included among (a linear combination of) the zero interaction states. 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 $\vec{a}^{(1)}=[1,1]$ and thus for two particles it is enough to specify a state at infinite interaction only by the quantum numbers $\vec{q}=[q_{0},q_{1}]$ and we leave the coefficients $[1,-1]$ implicit. When building our basis, we will typically increase the size by increasing the maximum energy of our states (above the lowest state). For example, if we say that we include all zero interaction states with an energy not greater than 2 (above the lowest energy state), we have the basis $[0,0]$, $[0,1]$, $[1,0]$, $[2,0]$, $[0,2]$ and $[1,1]$, and if we include all states at infinite interaction with energy not greater than 2 (above the lowest energy state), we have the infinite interaction states $[1,0]_{\infty}$, $[2,0]_{\infty}$, $[3,0]_{\infty}$ and $[2,1]_{\infty}$ (with the implicit coefficients $[1,-1]$). For simplicity we will restrict to having the same energy cutoff on both the infinite interaction states and the zero interaction states, but it is possible that an optimal scheme with different energy cutoffs for zero and infinite interaction exists.

5.1.2 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]_{\infty}$ at $g=+\infty$ (with the implicit coefficients $[1,-1]$).The first and third excited states are totally antisymmetric states (unaffected by the interaction) with the quantum numbers $[1,0]$ and $[2,0]$.

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 all basis states at zero and infinite interaction with a total energy (above the lowest energy state) not greater than some $E$, where $E$ is increased in steps of two. Thus the data points are for $E=0,2,4,\ldots$. 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.

5.1.3 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 well with the analytical result. To compute the density from the analytical result, we need to perform an integral transforming from Jacobi coordiantes to absolute coordinates (see Appendix B), and it turned out that the most accurate approach was to perform this integral numerically for various grid sizes $N$ and then fit a function of the form $f(N)=a+b/N+c/N^{2}$ to extrapolate to a final value of the density. For basis sizes larger than 40 we don’t see much improvement, but we do not claim to have that high numerical precision in neither our method nor in the numerical integral used for the analytical formula.

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

5.2 2+1

5.2.1 Energies

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 an energy of [math], $2$, $4$ and $6$ above the ground states, and the solid lines are computed using an energy 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.

5.2.2 Position space densities

Figure 8 and 9 shows the position space minority and majority density at $g=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.

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

5.2.3 Momentum space densities

In Figure 11 we compare the momentum space densities at $g=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.

5.3 6+1

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=1$ and $g=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.

6 Conclusions

In this paper we explored a new method for studying strongly coupled one-dimensional polaron systems, a method that generalizes that of [45]. Our results compare well both with analytical methods for two particles and with numerical methods based on matrix product states. The method converges well (exceptionally well for two particles), but does get worse when the number of particles increase.

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 (and in particular density matrices or momentum distributions) is however a significantly time consuming step, but again this part must only be carried out once for each chosen basis and we can then easily obtain the densities for any interaction strength $g$.

7 Acknowledgements

We would like to thank Artem Volosniev and Molte Andersen for useful discussions. REB acknowledges funding from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES). NTZ would like to acknowledge funding from Aarhus University Research Foundation under the JCS Fellowship program.

Appendix A Some integral formulas

A.1 Integral 1

In this section we will derive an expression for

[TABLE]

where $\mathcal{M}_{k}(x)$ is the set where $x$ is smaller than exactly $k$ of the coordinates $x_{1},\ldots,x_{n}$ and $\Phi_{\vec{v}}$ is the (normalized) totally antisymmetric wave function of the states corresponding to the quantum numbers in $\vec{v}=(v_{1},\ldots,v_{n})$. We will use induction to show that

[TABLE]

where $A$ is the matrix defined by $A_{ab}(x)=\int_{-\infty}^{x}f_{q_{a}}(x^{\prime})f_{p_{b}}(x^{\prime})\mathrm{d}x^{\prime}$ and $B_{ab}=\int_{x}^{\infty}f_{q_{a}}(x^{\prime})f_{p_{b}}(x^{\prime})\mathrm{d}x^{\prime}=\delta_{ab}-A_{ab}$. For $k=0$ we easily obtain

[TABLE]

which proves the base case. Now assume that $I_{\vec{q},\vec{p}}^{j}(x)=\frac{1}{j!}\partial_{\epsilon}^{j}\det(A+\epsilon B)_{\epsilon=0}$ for $j<k$. We then have

[TABLE]

where we have use the notation that $N(i)(j)$ is the matrix $N$ with row $i$ and column $j$ removed and similarly $\vec{q}(i)$ is the ordered set with the element indexed $i$ removed. We also used the formula $\mathrm{tr}[M\mathrm{adj}N]=\sum(-1)^{i+j}M_{ij}\det N(i)(j)=\partial_{\epsilon}\det(M+\epsilon N)_{\epsilon=0}$ and he factor $1/k={N\choose k}/\left(N{N-1\choose k-1}\right)$ can be inferred from combinatorics and the normalization of the wavefunctions. Thus our formula is proven by induction.

A.2 Integral 2

Let us now consider the integral

[TABLE]

where $\mathcal{M}_{k,l}(x,x^{\prime})$ is the set where $x$ is smaller than exactly $k$ of the coordinates $x_{1},\ldots,x_{n}$ and $x^{\prime}$ is smaller than exactly $l$ of the coordinates $x_{1},\ldots,x_{n}$. $\Phi_{\vec{v}}$ is the (normalized) totally antisymmetric wave function of the states corresponding to the quantum numbers in $\vec{v}=(v_{1},\ldots,v_{n})$. The result is

[TABLE]

where $A_{ij}=\int_{-\infty}^{min(x,x^{\prime})}f_{n_{i}}(x^{\prime\prime})f_{m_{j}}(x^{\prime\prime})\mathrm{d}x^{\prime\prime}$, $C_{ij}=\int_{min(x,x^{\prime})}^{max(x,x^{\prime})}f_{n_{i}}(x^{\prime\prime})f_{m_{j}}(x^{\prime\prime})\mathrm{d}x^{\prime\prime}$ and $B_{ij}=\int_{max(x,x^{\prime})}^{\infty}f_{n_{i}}(x^{\prime\prime})f_{m_{j}}(x^{\prime\prime})\mathrm{d}x^{\prime\prime}$. We can also prove this by induction. Note that if we assume $x>x^{\prime}$ and $k=0$, the formula is the same as (46) if the upper integral limit is changed from $\infty$ 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^{\prime}$, 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

[TABLE]

Since we assumed that $k<l$ and $x>x^{\prime}$, and the exact same proof can be done for $k>l$ and $x<x^{\prime}$, formula (50) follows.

Appendix B Two particle system

In this section we review the analytical solution of two particles in a harmonic trap, with a delta function interaction [46]. The full Hamiltonian is

[TABLE]

where

[TABLE]

By introducing Jacobi coordinates $x=(x_{1}-x_{2})/\sqrt{2}$, $p=(p_{1}-p_{2})/\sqrt{2}$, $X=(x_{1}+x_{2})/\sqrt{2}$ and $P=(p_{1}+p_{2})/\sqrt{2}$ we can split this Hamiltonian into two parts, namely

[TABLE]

where $H_{\mathrm{CM}}=X^{2}/2+P^{2}/2$ is just a harmonic oscillator corresponding to the center-of-mass motion, and

[TABLE]

The hard part, which will occupy most of this appendix, is solving for the eigenstates of $H_{\mathrm{rel}}$. The full set of eigenstates and eigenenergies are then obtained by tensor product with the eigenstates of $H_{\mathrm{CM}}$.

We will solve for the wavefuctions by first expanding in a harmonic oscillator basis. The Harmonic oscillator eigenfunctions are given by

[TABLE]

where $H_{n}$ are the Hermite polynomials. The energy is given by $E_{n}=n+1/2$. Let $|\Phi\rangle$ be an eigenstate for $H_{\mathrm{rel}}$. We have

[TABLE]

Solving for $c_{n}\equiv\langle n|\Phi\rangle$ and defining the quantity $A=\sum f_{n}(0)c_{n}$ we obtain

[TABLE]

Now multiplying both sides by $f_{n}(0)$ and summing over $n$, we can cancel $A$ from both sides to obtain

[TABLE]

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

[TABLE]

It thus makes sense to treat these simultaneously, so let us define

[TABLE]

To compute this function, we use the following relation between Hermite polynomials and Laguerre polynomials

[TABLE]

We thus obtain

[TABLE]

Now we use the integral representation

[TABLE]

to obtain

[TABLE]

Now we can recognize the generating function $e^{-tx/(1-t)}(1-t)^{-\alpha-1}=\sum t^{n}L_{n}^{\alpha}(x)$ to obtain

[TABLE]

where we have used a standard representation for the confluent hypergeometric function $U$. At $x=0$, we can use the relation $U(-\nu,1/2,0)=\Gamma(1/2)/\Gamma(1/2-\nu)=\sqrt{\pi}/\Gamma(1/2-\nu)$, to obtain

[TABLE]

Thus for the energy, we must solve the equation

[TABLE]

For the wavefunction, we instead have

[TABLE]

To find the normalization constant $A$ we can consider the normalization constraint

[TABLE]

Defining $\psi(x)=\Gamma^{\prime}(x)/\Gamma(x)$ and using the energy formula (68), we can simplify this to

[TABLE]

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_{\mathrm{rel}}+H_{\mathrm{CM}}$ is given by

[TABLE]

where we have labeled all eigenstates of $H_{\mathrm{rel}}$ (both even and odd) by $\Phi_{k}$ for $k=0,1,\ldots$ and $f_{n}$ are just the standard harmonic oscillator wavefunctions. The energy is likewise $E_{k,n}=E_{\Phi_{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})/\sqrt{2}$ and $X=(x_{1}+x_{2})/\sqrt{2}$. The single particle density matrix is just the square of the wavefunction in absolute coordinates, namely

[TABLE]

and the density is thus given by

[TABLE]

The momentum density can then be obtained by

[TABLE]

Appendix 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

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

This system is uniquely solved for $\omega_{1}$, $x_{1}$, $x_{L}$ and $x_{R}$ given the physically relevant quantities $\omega_{0}$, $\omega_{2}$, $x_{0}$, $x_{2}$, $\Delta_{0}$, $\Delta_{1}$ and $\Delta_{2}$. The solution is given by

[TABLE]

[TABLE]

and then $x_{L}$ and $x_{R}$ are given by

[TABLE]

[TABLE]

Extra care for these formulas must be taken when evaluating these expressions for $\Delta_{0}=\Delta_{2}$. In this case we have

[TABLE]

[TABLE]

For the symmetric case (symmetric around $x_{1}=(x_{0}+x_{2})/2$) where we also have $\omega_{0}=\omega_{2}$, we have

[TABLE]

We can compute an upper limit on the parameter $\Delta_{1}$. The highest value is the value such that $\omega_{1}=\infty$, 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

[TABLE]

which results in the solution

[TABLE]

Then the upper limit of $\Delta_{1}$ is given by $\Delta_{1,\text{max}}=\frac{1}{2}(x_{M}-x_{0})^{2}\omega_{0}^{2}+\Delta_{0}$.

We will now work out the wavefunctions and energies. We will work in units where $\hbar=1$ for simplicity and we will define $\nu_{0}$ and $\nu_{2}$ by $E=\omega_{0}(\nu_{0}+\frac{1}{2})+\Delta_{0}=\omega_{2}(\nu_{2}+\frac{1}{2})+\Delta_{2}$. The eigenfunctions are now uniquely given by

[TABLE]

for $x<x_{L}$ and

[TABLE]

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\rightarrow\pm\infty$). The function $D$ is the parabolic cylinder function given by

[TABLE]

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

[TABLE]

where

[TABLE]

and

[TABLE]

and where we have parametrized the energy as $E=\omega_{1}(\nu_{1}+\frac{1}{2})+\Delta_{1}$ (which we recall is also equal to $\omega_{0}(\nu_{0}+\frac{1}{2})+\Delta_{0}=\omega_{2}(\nu_{2}+\frac{1}{2})+\Delta_{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 $\psi$ and $\psi^{\prime}$ are continuous. To simplify the equations, we will define $r=\omega_{2}/\omega_{1}$, $R=\omega/\omega_{1}$, $\Delta=\hbar\omega_{1}\delta$, $C=\hbar\omega_{1}c$ and we work in units where $\mu\omega_{1}/\hbar=1$. This gives the equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

If we are given $\nu_{0},\nu_{1},\nu_{2}$ (which are all determined by the energy $E$), this is a linear system of equations for $C_{0},C_{2},C_{1}^{(1)},C_{1}^{(2)}$. 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 $\nu_{0},\nu_{1},\nu_{2}$). This system of equations, supplemented with normalization of the wave function, then fixes all the constants $C_{0},C_{2},C_{1}^{(1)},C_{1}^{(2)}$. 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)\times 2(N-1)$ matrix.

Appendix 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 [47]. 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

[TABLE]

where $c^{\dagger}$ 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 $\lvert\uparrow\rangle$ for the background fermions and $\lvert\downarrow\rangle$ for the impurity. Since we consider only a single $\lvert\downarrow\rangle$ fermion, we have naturally $\sum_{j}n_{j,\downarrow}=1$, with $\sum_{j}n_{j,\uparrow}$ also being normalized to the number of background fermions. We include additionally the trapping potential as the position-dependent $\epsilon_{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 $\sum_{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}$.

Appendix 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 $\partial_{\epsilon}^{i}D(\epsilon)|_{\epsilon=0}=\partial_{\epsilon}^{i}\frac{1}{i!}\det(M(\epsilon))|_{\epsilon=0}$, where $M(\epsilon)$ is some $n\times n$ matrix and $i=0,\ldots,n$. We evaluate these derivatives by computing the function $D(\epsilon)$ on $n+1$ values with $\epsilon_{i}=-1+2i/n$, $i=0,\ldots,n$, and then fitting a polynomial to these values and extracting the coefficients. These coefficients can be obtained by multiplying the vector $D(\epsilon_{i})$ with the inverse of the matrix $K_{ij}\equiv\epsilon_{i}^{j}$.

For the single-particle density matrix, we also need to compute terms of the form $\partial_{\epsilon}^{i}\partial_{\delta}^{j}D(\epsilon,\delta)|_{\epsilon=0,\delta=0}$. This is done similarly be fitting a polynomial of two variables to the values $D(\epsilon_{i},\epsilon_{j})$ with $\epsilon_{i}=-1+2i/n$, $i=0,\ldots,n$. We carry out the polynomial fit by applying the $(n+1)^{2}\times(n+1)^{2}$ matrix $K_{IJ}\equiv\epsilon_{\lfloor I/(n+1)\rfloor}^{\lfloor J/(n+1)\rfloor}\epsilon_{I\textrm{ mod }(n+1)}^{J\textrm{ mod }(n+1)}$ on the $(n+1)^{2}$ vector $D_{I}=D(\epsilon_{\lfloor I/(n+1)\rfloor},\epsilon_{I\textrm{ mod }(n+1)})$.