Experimentally accessible invariants encoded in interparticle correlations of harmonically trapped ultra-cold few-fermion mixtures

TL;DR

This paper introduces a numerical method to identify a specific invariant in many-body correlations of ultra-cold fermion mixtures, aiding experimental analysis of interparticle interactions without affecting center of mass excitations.

Contribution

A novel numerical approach to extract an invariant from correlations in ultra-cold fermion mixtures, facilitating experimental identification of eigenstates.

Findings

Identification of a correlation-encoded invariant.

Invariant enables labeling of center of mass excitations.

Method applicable for arbitrary mass ratios.

Abstract

System of a two-flavor mixture of ultra-cold fermions confined in a one-dimensional harmonic trap is studied in the frame of the center of mass. We present a numerical method of obtaining energetic spectra in this frame for an arbitrary mass ratio of fermionic species. We identify a specific invariant encoded in many-body correlations which enable one to determine an eigenstate of the Hamiltonian and to label excitations of the center of mass. The tool presented may be particularly useful in experimental analysis of the interparticle interactions which do not affect the center of mass excitations in a harmonic potential.

Figures

Equations

$$\hat{H}=\sum _{i=1}^{N_{\downarrow }}\left[-\frac{{\hslash }^2}{2m_{\downarrow }}\frac{{\partial }^2}{\partial x_i^2}+\frac{m_{\downarrow }{\omega }^2}{2}{x}_i^2\right]+\sum _{j=1}^{N_{\uparrow }}\left[-\frac{{\hslash }^2}{2m_{\uparrow }}\frac{{\partial }^2}{\partial y_j^2}+\frac{m_{\uparrow }{\omega }^2}{2}{y}_j^2\right]+g_{1D}\sum _{i=1}^{N_{\downarrow }}\sum _{j=1}^{N_{\uparrow }}\delta (x_i-y_j),$$

$$\hat{H}=\sum _{i=1}^{N_{\downarrow }}\left[-\frac{1}{2}\frac{{\partial }^2}{\partial x_i^2}+\frac{1}{2}{x}_i^2\right]+\sum _{i=1}^{N_{\uparrow }}\left[-\frac{1}{\mu }\frac{1}{2}\frac{{\partial }^2}{\partial y_i^2}+\mu \frac{1}{2}{y}_i^2\right]+g\sum _{i=1,j=1}^{N_{\downarrow },N_{\uparrow }}\delta (x_i-y_j),$$

$$\hat{H}=\sum _{\sigma }\int \mathrm{d}x{\hat{\Psi }}_{\sigma }^{†}(x){\hat{H}}_{\sigma }(x){\hat{\Psi }}_{\sigma }(x)+g\int \mathrm{d}x{\hat{\Psi }}_{\uparrow }^{†}(x){\hat{\Psi }}_{\downarrow }^{†}(x){\hat{\Psi }}_{\downarrow }(x){\hat{\Psi }}_{\uparrow }(x),$$

$$H_{\downarrow }(x)$$

$$H_{\uparrow }(x)$$

$$R_{\mathrm{CM}}=\frac{{\sum }_{i=1}^{N_{\downarrow }}{x}_i\cdot m_{\downarrow }+{\sum }_{j=1}^{N_{\uparrow }}{y}_j\cdot m_{\uparrow }}{N_{\downarrow }\cdot m_{\downarrow }+N_{\uparrow }{m}_{\uparrow }}=(N_{\downarrow }+\mu N_{\uparrow }{)}^{-1}\left(\sum _{i=1}^{N_{\downarrow }}{x}_i+\sum _{j=1}^{N_{\uparrow }}\mu \cdot y_j\right).$$

$${\hat{R}}_{\mathrm{CM}}=(N_{\downarrow }+\mu N_{\uparrow }{)}^{-1}\int \mathrm{d}x\left({\hat{\Psi }}_{\downarrow }^{†}(x)x{\hat{\Psi }}_{\downarrow }(x)+\mu {\hat{\Psi }}_{\uparrow }^{†}(x)x{\hat{\Psi }}_{\uparrow }(x)\right).$$

$${\hat{\Psi }}_{\sigma }(x)=\sum _i{\varphi }_{i\sigma }(x){\hat{a}}_{i\sigma },$$

$$X_{ij,\sigma }={\mu }_{\sigma }\int \mathrm{d}x{\varphi }_{i\sigma }^{\ast }(x)x{\varphi }_{j\sigma }(x),$$

$${\hat{R}}_{\mathrm{CM}}=(N_{\downarrow }+\mu N_{\uparrow }{)}^{-1}\sum _{ij}\sum _{\sigma }{X}_{ij,\sigma }{\hat{a}}_{i\sigma }^{†}{\hat{a}}_{j\sigma }.$$

$${\hat{P}}_{\mathrm{CM}}=-i\int \mathrm{d}x\left({\hat{\Psi }}_{\downarrow }^{†}(x){\partial }_x{\hat{\Psi }}_{\downarrow }(x)+{\hat{\Psi }}_{\uparrow }^{†}(x){\partial }_x{\hat{\Psi }}_{\uparrow }(x)\right),$$

$${\hat{P}}_{\mathrm{CM}}=\sum _{ij}\sum _{\sigma }{P}_{ij,\sigma }{\hat{a}}_{i\sigma }^{†}{\hat{a}}_{j\sigma },$$

$$P_{lm,\sigma }=-i\int \mathrm{d}x{\varphi }_{l\sigma }^{\ast }(x){\partial }_x{\varphi }_{m\sigma }(x).$$

$${\hat{H}}_{\mathrm{CM}}=\frac{{\hat{P}}_{\mathrm{CM}}^2}{2(N_{\downarrow }+\mu N_{\uparrow })}+\frac{(N_{\downarrow }+\mu N_{\uparrow })}{2}{\hat{R}}_{\mathrm{CM}}^2.$$

$$E({\nu }_1,{\nu }_2,\dots ,{\nu }_N)=\underset{({\nu }_1+\frac{1}{2})}{\underbrace{E({\nu }_1)}}+E({\nu }_2,{\nu }_3,\dots ,{\nu }_N),$$

$${\hat{H}}_{\mathrm{REL}}^{N=2}=-\frac{1}{2}\frac{{\partial }^2}{\partial r^2}+\frac{1}{2}{r}^2+\sqrt{\frac{\mu }{1+\mu }}g\delta (r).$$

$$C(x,y)=⟨i|{\hat{\Psi }}_{\downarrow }^{†}(x){\hat{\Psi }}_{\uparrow }^{†}(y){\hat{\Psi }}_{\uparrow }(y){\hat{\Psi }}_{\downarrow }(x)|i⟩.$$

$$\Psi (x_1,\dots ,x_{N_{\downarrow }},y_1,\dots ,y_{N_{\uparrow }})=\Psi (R)\Psi (r_1,r_2,\dots ,r_{N-1}),$$

$$\mathcal{I}=\int C(x,y)(x-y{)}^2\mathrm{d}x\mathrm{d}y.$$

$$\mathcal{J}=⟨i|{\hat{R}}_{\mathrm{CM}}^2|i⟩.$$

$$\mathcal{J}=\sum _{\sigma }{N}_{\sigma }{\mu }_{\sigma }^2\int {\rho }_{\sigma }(x)x^2\mathrm{d}x+\sum _{\sigma }{N}_{\sigma }(N_{\sigma }-1)\int C_{\sigma \sigma }(x,x^{\prime })x{x}^{\prime }\mathrm{d}x\mathrm{d}{x}^{\prime }+\mu N_{\downarrow }{N}_{\uparrow }\int C(x,y)\mathrm{d}x\mathrm{d}y,$$

$$C_{\sigma \sigma }(x,x^{\prime })=⟨i|{\hat{\Psi }}_{\sigma }^{†}(x){\hat{\Psi }}_{\sigma }^{†}(x^{\prime }){\hat{\Psi }}_{\sigma }(x^{\prime }){\hat{\Psi }}_{\sigma }(x)|i⟩.$$

Full text

††journal: Few-Body Systems

1 Introduction

Few-body problems has been recently vastly explored theoretically Blume2010TwoThreeToMany ; Blume2012FewBodyTraps ; VolosnievNC2014 . Partially it is due to the experiments on ultra-cold atoms making possible to verify theoretical concepts that had been infeasible before lewenstein2012ultracold . Recent developments in the field have shown that it is possible to control many parameters such as the number of atoms or the strength of mutual interactions with an exceptional precision wenz2013fewToMany ; zurn2012fermionization ; serwane2011tunableFermionSystem ; MurmannPRL2015a ; Kaufman2015Entangling ; Onofrio2016Mixtures . Apart from adjusting these parameters, also the geometry of an external trap can be tuned. In consequence, it is even possible to study effectively one-dimensional systems having very interesting and unique properties lewenstein2012ultracold .

Because of their complexity, few-body systems are usually very challenging to treat analytically. Also numerical studies may be very time-consuming due to a quite large number of degrees of freedom, particularly for fermionic species. Also the straightforward mean-field methods do not work well for a small number of particles. Therefore, in general, systems with mesoscopic number of particles are very difficult to treat theoretically.

One way to overcome these problems is to simplify the description by introducing variables directly related to the natural degrees of freedom of the system. It is known that these generalized variables, in which description of the system is very simple, do exist. However, even for small number of particles one needs to perform not obvious transformations Busch1998 ; LiuHuDrummond2010ThreeAttractive ; KestnerDuan2007ThreeBody ; WernerCastin2006UnitaryThreeBody ; loft2015variational ; Koscik2012 ; Koscik2015vonNeuman ; Koscik2016Bipartite ; Olshanii20154Body ; Harshman2016I ; Harshman2016II . For example, for four particles it is convenient to introduce Jacobi coordinates Dehkharghani2016Impenetrable ; Garcia-MarchPRA2014 . The situation is even worse for larger number of particles and in practice these transformations are very hard to implement. Apart from these problems, in the case of quantum description, additional complexity appears because of the indistinguishability of particles and quantum uncertainty of the center of mass position. If this uncertainty is relatively large, what can happen in attractive systems, then experimental separtion of relative dynamics from the center of mass motion becomes not trivial and has to be related to higher order correlations Gajda2006 . On top of this, the symmetrization or anti-symmetrization of the wave function, when expressed in the language of generalized variables leads to the highly non-trivial constrains in the many-body picture.

Moreover, the generalized coordinates are not well suited for experiments. Since atomic microscopes detect directly positions of particles. Therefore, generalized variables are not helpful in the interpretation of experimental outputs. It would be most convenient if a theoretical description of the system could be done in the most natural variables, i.e., positions of particles measured in experiments. Unfortunately, in this case a lot of redundant information comes from trivial dynamics of the centre of mass of the system which is completely insensitive to interparticle interactions. This general problem can be overcome in the case of the harmonic confinement, where the centre-of-mass degree of freedom can be separated out from the relative motion Busch1998 ; LiuHuDrummond2010ThreeAttractive ; KestnerDuan2007ThreeBody ; WernerCastin2006UnitaryThreeBody ; Garcia-MarchPRA2014 ; loft2015variational ; lawson1974 ; FedorovMikkelsen2015Recombination ; Amin2016CorrelationsHO ; Gajda2000COM ; Bienias2014Rydberg ; Brandt2016COMconvergence ; Blakie2016PCS ; Bradley2016Brownian . Here, we study one-dimensional fermionic mixtures and we show how to distillate the relevant information about relative excitations encoded in the full many-body wave functions. Particularly, we identify experimentally accessible quantity which are invariant under the change of the centre of mass excitations.

Additionally, we present a simple numerical method of extracting the eigenstates of the centre of mass and the relative motion from particular eigenstate of the full many-body Hamiltonian. In consequence, it is possible to simplify the description of the many-body spectrum. The method is based on a trivial properties of the centre-of-mass motion, however the procedure may shed some light on the problem of few-body systems. We also show a technique of obtaining information about the relative excitations from the simple two-body correlation, namely the two-body density profile. We believe that it may be a useful tool in the field of experimental few-body physics where the two-body correlations are experimentally available to observe. We show that the procedure is very general and can be used for fermionic mixtures of different masses.

The article is composed as following. In the next Sec. 2 we describe the system under consideration. Then we introduce the frame of the centre of mass in Sec. 3 which is followed by the description of our numerical methods in Sec. 4. The subsequent Sec. 5 contains the analysis of the spectra in the centre-of-mass frame. Then in Sec. 6 a method of decoding an information from the interparticle correlations is presented. Finally, the last Sec. 7 concludes our paper.

2 The model

Let us consider a one-dimensional system of two distinguishable types of fermions confined in an external harmonic potential of a frequency $\omega$. We will label them with the quantum number $\sigma\in\{\downarrow,\uparrow\}$. The quantum number $\sigma$ does not necessary mean that the particles differ with the spin projection. In fact they may not have the spin degree of freedom at all. We assume that $\sigma$ is conserved and it is not a dynamical variable. Experimentally, this can be realized for example by choosing two proper hyperfine states of a nucleus wenz2013fewToMany or two different elements Wille6Li40K ; tiecke2010Feshbach6Li40K . We model interactions between particles of different flavours with a $\delta$-like contact potential PethickSmith which is well justified in the ultra-cold regime and causes no mathematical problems in one dimension (there is no need of regularization). In this approximation, the interaction within a given flavour vanishes due to the Pauli exclusion principle. The Hamiltonian of the system reads:

[TABLE]

where $N_{\downarrow},N_{\uparrow}$ are the numbers of fermions with masses $m_{\downarrow}$ and $m_{\uparrow}$, respectively. The Hamiltonian of the system  (1) commutes with the operators of particle numbers, $[\hat{H},\hat{N}_{\uparrow}]=[\hat{H},\hat{N}_{\downarrow}]=0$, what implies that the number of particles of each kind is conserved and fermions cannot change their flavour. In real experiments this condition is met. In the realizations using two different hyperfine states the interactions channels that change the quantum number $\sigma$ are practically disabled. For a mass imbalanced system, there is an additional superselection rule resulting from the conservation of the mass that forbids changing atoms of one species to another. An effective one-dimensional coupling parameter $g_{\mathrm{1D}}$ can be obtained from the full three-dimensional theory of scattering by integrating out two dimensions as shown in Olshanii1998 . What is worth noticing, the coupling $g_{\mathrm{1D}}$ may be tuned experimentally, for example, by varying an external magnetic field Wille6Li40K ; tiecke2010Feshbach6Li40K .

We assume the same trapping frequency $\omega$ for both species which implies the same energy scale $\hbar\omega$ independently on their masses. In the following, we choose to measure all the masses in the unit of $m_{\downarrow}$ and consequently the unit of length is $\sqrt{\hbar/(m_{\downarrow}\omega)}$. The Hamiltonian (1) in a dimensionless form reads:

[TABLE]

where we introduced the mass ratio of the particles, $\mu=m_{\uparrow}/m_{\downarrow}$, and the dimensionless interaction strength coupling, $g=g_{\mathrm{1D}}\sqrt{m_{\downarrow}/(\hbar^{3}\omega)}$.

Before we apply our numerical techniques, it is extremely convenient to rewrite the Hamiltonian (3) in the second quantization form:

[TABLE]

where appropriate single-particle Hamiltonians $\hat{H}_{\sigma}(x)$ are defined as:

[TABLE]

A field operator $\hat{\Psi}_{\sigma}(x)$ annihilates a particle of a type $\sigma$ in a position $x$. Within the same flavour, for the indistinguishable fermions, the anticommutation relations hold $\{\hat{\Psi}_{\sigma}(x),\hat{\Psi}_{\sigma}^{\dagger}(x^{\prime})\}=\delta(x-x^{\prime})$ and $\{\hat{\Psi}_{\sigma}(x),\hat{\Psi}_{\sigma}(x^{\prime})\}=0$. For the particles of different flavours, the choice of commutation relation is not unique, since the physical results do not depend on this choice (as long as the relations are used consistently).

3 The centre-of-mass frame

In the classical mechanics as well as in the quantum mechanics IBB1985COM ; IBB2002COM ; Sowinski2007Linear there is a particular class of problems which are described by a Hamiltonians having biquadratic form of positions and momenta. In this case, if the interaction between particles depends only on their relative positions, it is possible to canonically transform dynamical variables in such a way that the centre-of-mass motion separates from the relative motion of the system. The motion of the centre of mass does not depend on the interaction which influences only the motion of the remaining degrees of freedom. The two-body problem Busch1998 is an example but the separation can also be done for larger number of particles LiuHuDrummond2010ThreeAttractive ; KestnerDuan2007ThreeBody ; WernerCastin2006UnitaryThreeBody ; FedorovMikkelsen2015Recombination ; Amin2016CorrelationsHO at the cost of using more sophisticated methods of describing the relative motion loft2015variational ; Dehkharghani2016Impenetrable . In practice, the mentioned transformation of variables is difficult to perform straightforwardly, however it is a standard tool in a theoretical analysis of few-body problems Brandt2016COMconvergence ; Blakie2016PCS ; Bradley2016Brownian . In the following we will separate the the centre-of-mass motion within the scope of the second quantisation formalizm.

Position of the centre of mass for classical system of $N_{\uparrow}+N_{\downarrow}$ particles has a form:

[TABLE]

Since the quantity (5) is a linear function of the positions of particular particles, the corresponding quantum-mechanical counterpart, written in the second quantization picture, is defined as:

[TABLE]

The field operators can be decomposed into the basis of the eigenfunctions of an appropriate harmonic oscillator $\phi_{i\sigma}(x)$:

[TABLE]

where the operator $\hat{a}_{i\sigma}$ annihilates a fermion of type $\sigma$ from the orbital $i$ and the following condition is met $(H_{\sigma}(x)-\hbar\omega(i+1/2))\phi_{i\sigma}(x)=0$. The sum (7) runs from zero to infinity. For convenience we introduce the coefficients:

[TABLE]

where $\mu_{\downarrow}=1$ and $\mu_{\uparrow}=\mu$. Then, the centre-of-mass operator gets simplified form:

[TABLE]

Analogously, one introduces the operator of the total momentum:

[TABLE]

which can be rewritten in a decomposed form as:

[TABLE]

where appropriate coefficients are defined as:

[TABLE]

It can be shown that the operator of the total momentum $\hat{P}_{\mathrm{CM}}$ and the operator of the position of the centre of mass $\hat{R}_{\mathrm{CM}}$ are canonically conjugate variables and natural commutation relations between them hold, $[\hat{R}_{\mathrm{CM}},\hat{P}_{\mathrm{CM}}]=i\hbar$ and $[\hat{R}_{\mathrm{CM}},\hat{R}_{\mathrm{CM}}]=[\hat{P}_{\mathrm{CM}},\hat{P}_{\mathrm{CM}}]=0$. The Hamiltonian of the centre of mass written in a dimensionless units reads:

[TABLE]

As suspected, the Hamiltonian (13) describes a particle of the mass equal to the total mass of the system $(N_{\downarrow}+\mu N_{\uparrow})$ trapped in the harmonic potential. As a consequence, the total Hamiltonian can be separated to $\hat{H}=\hat{H}_{\mathrm{CM}}+\hat{H}_{\mathrm{REL}}$. $\hat{H}_{\mathrm{REL}}$ represents the relative Hamiltonian and only this part contains all information about internal interactions between ultra-cold fermions. It is also easy to check that the total Hamiltonian of the system $\hat{H}$ and the Hamiltonian of the centre of mass $\hat{H}_{\mathrm{CM}}$ as well as the interaction Hamiltonian $\hat{H}_{\mathrm{REL}}$ do commute, i.e., $[\hat{H},\hat{H}_{\mathrm{CM}}]=[\hat{H},\hat{H}_{\mathrm{REL}}]=0$. It means that there exists a basis in the Hilbert space in which all three Hamiltonians have simultaneously diagonal form.

Every eigenstate of the Hamiltonian (3) can be denoted by the set of $N=N_{\downarrow}+N_{\uparrow}$ quantum numbers. Since the centre of mass is separated out, the remaining $N-1$ quantum numbers describe the relative degrees of freedom of the system. Therefore, the total energy $E$ of a given eigenstate of the Hamiltonian (3) can be written in a form:

[TABLE]

where ${\bm{\nu}}=(\nu_{1};\nu_{2},\ldots,\nu_{N})$ is a set of quantum numbers describing a given eigenstate of the Hamiltonian $\hat{H}$ and $E(\nu_{1})$ corresponds to the centre-of-mass energy present in this state. Naturally, for a given set of quantum numbers $(\nu_{2},\nu_{3},\ldots,\nu_{N})$ describing the relative motion excitations, there are infinitely many states with different values of $\nu_{1}$ with energies separated by the excitations of the centre of mass.

The method of getting the energy spectrum of the Hamiltonian (3) filtered out from the centre-of-mass excitations is based on a quite trivial observation that the eigenstates of the Hamiltonian $\hat{H}$ are composed as products of eigenstates of $\hat{H}_{\mathrm{CM}}$ and $\hat{H}_{\mathrm{REL}}$. Therefore, we can find the eigenvalue $E(\nu_{1})$ of any eigenstate of $\hat{H}$ by calculating an expectation value of the Hamiltonian of the centre of mass $\hat{H}_{\mathrm{CM}}$ in this state. Then, by collecting all eigenstates of the total Hamiltonian $\hat{H}$ with the same eigenvalue $E(\nu_{1})$, we get the spectrum that contains only one specific excitation of the centre of mass. The choice of the quantum number $\nu_{1}$ is arbitrary, but for technical reasons (our method has the best accuracy in the low energy regime) the best choice is to gather all eigenvalues for $\nu_{1}=0$ corresponding to the ground state of the centre of mass with the eigenenergy $\hbar\omega/2$.

4 Method

The Hamiltonian (3) does not couple the Hilbert subspaces of different particle numbers, thus it can be diagonalized in each of those subspaces independently. The numerical method of the exact diagonalization of the Hamiltonian $\hat{H}$ returns the eigenenergies $E_{i}$ as well as the many-body eigenstates $|i\rangle$. Of course to apply the method, one has to introduce in the sum (7) a cutoff $K$ at sufficiently large index. It is justifiable because for any finite interaction $g$, only the lowest excitations play a role. In principle, we assume that the cutoff is appropriate if its further increase does not affect the results. In practice, it is done by calculating consecutive fidelities $\mathbf{F}_{K}(|i\rangle)$ defined for chosen many-body eigenstate $|i\rangle$ as the overlap between states $|i\rangle_{K}$ and $|i\rangle_{K+1}$, i.e., the eigenstate calculated $|i\rangle$ for cutoffs $K$ and $K+1$. We accept the cutoff as a sufficiently large if the fidelity deviates from one by less than $10^{-3}$. In Fig. 1 we show example adaptation of this procedure to the many-body ground-state of the system for different number of particles, different mass ratios $\mu$ and interaction $g=3$.

Since the two Hamiltonians $\hat{H}$ and $\hat{H}_{\mathrm{CM}}$ do commute, one can calculate the energy in the centre-of-mass frame using the Hamiltonian of the centre of mass $\hat{H}_{\mathrm{CM}}$ for any eigenstate of the Hamiltonian $\hat{H}$. Then by collecting all the states in the ground state of the centre-of-mass motion, i. e. $\hat{H}_{\mathrm{CM}}|i\rangle=\hbar\omega/2|i\rangle$, one finally gets all states with excited internal motion only.

The situation is slightly more complicated if the eigenstates are $j$-fold degenerated. Then, the whole set of the eigenstates $\{|m_{1}\rangle,\ldots,|m_{j}\rangle\}$ of the full Hamiltonian $\hat{H}$ has the same energy. In these cases we calculate all matrix elements of the centre-of-mass Hamiltonian $\hat{H}_{\mathrm{CM}}$ in the basis of the degenerated eigenstates $\langle i|\hat{H}_{\mathrm{CM}}|i^{\prime}\rangle$, where $|i\rangle\in\{|m_{1}\rangle,\ldots,|m_{j}\rangle\}$, and diagonalize it. The lowest eigenstate we got after diagonalization is obviously an eigenstate of $\hat{H}_{\mathrm{CM}}$ and it is a superposition of some eigenstates of $\hat{H}$. This eigenstate has the lowest energy of the centre-of-mass excitation. Therefore, this combination is the state we are interesting in.

Note, that in the procedure described above any assumptions about interactions, mass ratios or any additional symmetries are not made. Thus the technique is very general. In the following, the described method of filtering out the spectrum of the higher excitations of the centre of mass, is applied to the systems of equal masses $\mu=1$ as well as to the systems with mass imbalance $\mu\neq 1$, bearing in mind that nowadays experiments on mixtures of different atomic species are also possible Wille6Li40K ; tiecke2010Feshbach6Li40K .

5 Spectrum of the Hamiltonian

First, we apply the procedure described in the previous section to the two-flavoured systems of several particles of the same mass $\mu=1$. The results are shown in Figure 2. The top panels show the spectrum in the centre-of-mass frame (the top left panel reproduces the result obtained analytically in Busch1998 ). The eigenstates in the top panels are chosen to be in the ground state of the centre of mass, thus by shifting this plot by $\hbar\omega,2\hbar\omega,\ldots$ one can obtain the full spectrum of the system (shown on the bottom panels of Figure 2) reproducing results well-known for the equal mass mixtures of fermions SowinskiGrass2013FewInteracting ; Gharashi2013UpperBranchCorrelations . By increasing the number of particles, the spectrum becomes more complicated. Even after simplification to the centre-of-mass frame, the spectrum is complex due to increasing number of degrees of freedom in a relative motion (right top panel of Figure 2).

Similar analysis can be done also for the mass-imbalanced mixtures as shown in Figure 3. The top panels show the spectrum in the centre-of-mass frame. Shifting this spectrum by $\hbar\omega,2\hbar\omega,\ldots$ leads directly to the full many-body spectrum (bottom panels). Apart from the trivial two-body case, the spectra in the centre-of-mass frame differ qualitatively from spectra of the equal mass systems. For example, totally antisymmetric states do not exist in the mass-imbalanced case. This is a direct consequence of different shapes of the single-particle orbitals (for more details see Pecak2016Separation ; Pecak2016Transition ).

It is worth noticing that the two-particle problem is basically the same for equal mass and mass imbalanced cases. In this particular case $N_{\uparrow}=N_{\downarrow}=1$, changing the mass of one component leads only to an effective change in the interaction strength. It is quite obvious when one introduces the centre-of-mass position, $R_{\mathrm{CM}}=(x_{1}+\mu x_{2})/(1+\mu)$, and the relative position, $r=(x_{1}-x_{2})\sqrt{\mu/(1+\mu)}$. Then, the centre-of-mass Hamiltonian $\hat{H}_{\mathrm{CM}}$ gets a form of a Hamiltonian of the standard harmonic oscillator. Whereas, the relative Hamiltonian $\hat{H}_{\mathrm{REL}}$, due to a specific normalization of the relative position, has the form:

[TABLE]

It means that indeed, for any $\mu$ one can rescale the interaction strength, $g\rightarrow g\sqrt{(1+\mu)/\mu}$, to map the different mass spectrum to the equal mass spectrum (compare top panels of Figure 2 and Figure 3). Note however, that although the energies are rescaled properly, other quantities, like for example single-particle densities, are different Pecak2016Separation .

6 Two-body correlations

All the information about the state of the system is encoded in its wave function. However, this mathematical object is very complex, especially for many-body states, due to a large number of degrees of freedom. To get more insight into properties of a given quantum state, other simpler quantities, such as energy, single-particle density profile or correlations Manz2010G2 , may be considered to describe the properties of the system. However, recent experiments allow one not only to take a ’photograph’ of all ultra-cold atoms but also to measure all their positions simultaneously bakr2010probing ; sherson2010single ; OmranBloch2015PauliBlocking ; Cheuk2015FermionicMicroscope ; Cheuk2016MottMicroscope ; Cheuk2016Correlations ; ParsonsPRL2015 ; EdgePRA2015ImagingFermions ; haller2015single ; ParsonsPRL2015 . It means that, in principle, it is possible to measure many-body correlations Gajda2016Imaging . This new tool gives another opportunity to probe theoretically predicted correlations in quantum systems and motivates us to describe closely the properties of higher correlations in few-body systems.

First, let us focus mainly on pair correlations. Particularly, we analyze the properties of the system encoded in the pair-density profile, i.e. the diagonal part of the reduced two-body density matrix. For a state $|i\rangle$ it is defined as:

[TABLE]

This quantity is the simplest one that contains information about both components. Integrating out any of the two variables would reduce the information to just one component, i.e., to the single-particle density profile.

In Figure 4 we show the pair-density profile calculated for the ground state and one of the excited states of the system with $N_{\uparrow}=2,N_{\downarrow}=1$, and $g=3$ (upper left and bottom left panels). The excited state is excited only in the internal motion. In the right panels we show the same quantity calculated for the corresponding states with a single excitation of the centre of mass. It is seen that based on the shape of the pair-density profile shown in Figure 4, it is very hard to predict the set of quantum numbers ${\bm{\nu}}$. It means that if one obtains such correlations experimentally for some states, it may be difficult to distinguish and properly classify them, especially if they contain unknown centre of mass or the relative motion excitations.

The classification becomes even more confusing for the mass-imbalanced case. The pair-density profiles for $\mu=5$ are shown in Figure 5. Here, in contrast to the equal mass case, states being in the same state of the centre-of-mass motion and differing only by relative excitations have almost the same profile of the pair density (left panels correspond to the ground state and right panels to the first excited state of the centre-of-mass motion). The relative motion excitations change marginally the pair-density profile. As a consequence, these states are in practice indistinguishable. Additionally, the states (as seen in Figure 3) are almost degenerated and have similar single-particle density profiles. The only quantity that distinguishes these eigenstates is the centre-of-mass motion excitation. However, much more information is encoded in the pair-density profiles. As shown in the following, this hidden information on the centre of mass and the relative motion excitations can be extracted directly from the experimentally measured pair-density profiles.

The extraction is possible because one can separate some degrees of freedom in the system. The separation of the centre of mass from the internal degrees of freedom means that a wave function of any eigenstate of the system can be presented as a product:

[TABLE]

where $R$ denotes the centre of mass coordinate and $r_{1},r_{2},\dots,r_{N-1}$ are some general relative coordinates. From all possible combinations of coordinate differences of the form $x_{i}-y_{j},x_{i}-x_{j}$ and $y_{i}-y_{j}$, only $N_{\uparrow}+N_{\downarrow}-1$ are independent. This implies that one has a lot of freedom in choosing relative distances between particles. The relative distance $r_{1}=x_{1}-y_{1}$ between a pair of fermions of different species is extremely useful and obviously, due to indistinguishability, does not depend on the choice of particular coordinate $x_{i}$ or $y_{j}$.

To explore the information hidden in two-body correlations of a state $|i\rangle$, let us consider the second moment of the pair-density profiles:

[TABLE]

It can be shown straightforwardly that this quantity does not depend on the centre of mass position $R$, i.e., it is insensitive to the excitations of the centre of mass of the system. Therefore, as long as the eigenstates of the Hamiltonian (3) have the same value of relative degrees quantum numbers $(\nu_{2},\nu_{3}\ldots\nu_{N})$, they share the same value of ${\cal I}$. Since ${\cal I}$ is invariant under changes of the centre of mass excitation, it is possible to divide all eigenstates of the system with respect to excitations of the relative motion. It can be very useful when one is interested in preparing the state with a specific excitation of the internal motion but simultaneously is not able to control excitations of the centre of mass. By calculating the invariant ${\cal I}$ from the experimentally accessible pair-density profile an internal excitation can be determined. The method is quite useful since, as it is shown in Figure 4 and Figure 5, the pair-density profile may have completely different shapes for states with different centre of mass excitations. However, ${\cal I}$ extracted from these shapes is exactly the same.

We calculated ${\cal I}$ for many different parameters (mass ratio $\mu$, number of particles $N_{\uparrow},N_{\downarrow}$ and the interaction strength $g$). The results are shown in Figure 6. The consecutive eigenstates of the centre of mass excitation $\nu_{1}=0,1$ and $2$ are marked by black crosses, red squares, and blue circles, respectively. We found, as expected, that for a given set of quantum numbers describing the relative motion ($\nu_{2},\nu_{3},\ldots,\nu_{N}$), the quantity ${\cal I}$ does not depend on the centre of mass excitations — for a given interaction $g$ cross, square and circle are centred at the same point. Some discrepancy may be seen for strong interactions $g$, especially for the equal mass case $\mu=1$ (top panels of Figure 6). This is an artefact of our numerical method caused by a finite size of the basis. This deviation may be reduced by introducing larger single-particle basis (larger cutoff parameter), increasing simultaneously a time-cost of calculations. Since the basis functions for mass-imbalanced systems are different for different species, therefore a numerical error for $\mu=5$ is smaller than in the equal mass case $\mu=1$ (bottom and top panels of Figure 6, respectively). We can also see that ${\cal I}$ grows monotonically with the interaction strength $g$. The growth is what one would expect, since the invariant ${\cal I}$, as seen in the definition (18), is a quantity that measures a width of the wave function in the relative coordinate (the width grows with increasing repulsions). There is one exception for the totally antisymmetric states (so-called Girardeau), which do not depend on the interaction $g$. They exist only in the equal mass ($\mu=1$) mixtures and are visible in the spectra of Figure 2 as horizontal lines. In these states, due to the antisymmetry of the wave function under exchange of any two particles, the interaction energy vanishes.

Finally, let us note that it is possible also to define other invariants calculated from the many-body correlations. For example, instead of ${\cal I}$, one can consider an existence of a different quantity ${\cal J}$ related to the width of the centre of mass distribution of a state $|i\rangle$:

[TABLE]

This quantity, in contrast to ${\cal I}$, is invariant under the change of excitations of the relative coordinates, i.e. it does not depend on the quantum numbers ($\nu_{2},\nu_{3},\ldots,\nu_{N}$). Therefore, it can be used as an indicator of excitation of the centre of mass regardless of the internal motion. It can be shown that, in contrast to ${\cal I}$, which depends only on the two-body pair-density profile, ${\cal J}$ depends also on single-particle densities and pair-density profiles of the fermions of the same type:

[TABLE]

where the pair-density profile for the same type fermions is defined analogously to (16):

[TABLE]

Since the invariant ${\cal J}$ depends on all three two-body density profiles, it is a less convenient quantity than ${\cal I}$ for a practical use.

7 Conclusions

We have studied the properties of the system of a mass imbalanced two-component fermionic mixture in the centre-of-mass frame. We show how to numerically construct the Hamiltonian of the centre of mass $\hat{H}_{\mathrm{CM}}$ and how to obtain eigenvectors of both Hamiltonians $\hat{H}$ and $\hat{H}_{\mathrm{CM}}$ expressed in a natural Fock basis. Then, by determining excitations of the centre-of-mass motion, we filtered out its excitations and obtained simplified spectrum of the system being in the ground state of the centre-of-mass motion. In this way we were able to perform analysis of many different properties of the system, regardless of the excitations of the centre of mass. The procedure has been applied to systems with different number of fermions and different mass ratio generalizing previous results Busch1998 ; loft2015variational .

To investigate the properties of a system confined in a harmonic trap we have shown that eigenstates are characterized by specific invariants that can be straightforwardly calculated from the observable available in the experiment. The fact that some important information about the eigenstates can be extracted from the experimentally accessible correlations comes straightforwardly from from the separation of the centre-of-mass motion.

By comparing the values of the invariant ${\cal I}$ calculated theoretically with these measured experimentally one can control and engineer output eigenstates in experiments. Moreover, both methods (filtering out the centre of mass excitations and using invariants encoded in pair correlations) can be very useful in distinguishing eigenstates labeled by different quantum numbers ${\bm{\nu}}$. Therefore, the method can be used for analysing experimental outputs as long as the eigenstates of the system are consdered. If the output state is a superposition or a mixed state, then the method is still valid, but might be not as helpful as in the case of eigenstates.

Let us note that the results presented may be important also in the case of deviations from a harmonic potential. Then the centre-of-mass motion is coupled to the relative motion of the particles and can lead to very interesting phenomena known as Inelastic Confinement-Induced Resonances that has been recently studied both theoretically and experimentally Sowinski2012Creation ; Saenz2012ICIR ; Saenz2013ICIR .

We believe that our results may shed some light on the strongly correlated few-body systems regardless the centre of mass excitations. In broader context, the above methods may be used in other fields, like nuclear physics Ozen2014mappingNuclear ; Zinner2013nucleiColdAtoms or chemistry. The giant numerical complexity of the many-body systems increases exponentially with the number of particles, therefore the limit for our scheme is set by the numerical power of computers. The method can be easily generalized to Bose-Bose and Fermi-Bose mixtures and can be used beyond the short-range contact potential as long as the form of the interatomic potential allows one to separate out the centre-of-mass motion.

Acknowledgements.

This work was supported by the (Polish) National Science Center Grants No. 2016/21/N/ST2/03315 (DP) and 2016/22/E/ST2/00555 (TS). MG acknowledges support from the EU Horizon 2020-FET QUIC641122.