MCTDH-X: The multiconfigurational time-dependent Hartree method for indistinguishable particles software

The objective of the MCTDH-X software is to numerically solve the TISE or TDSE for a given many-body Hamiltonian, which describes N interacting, indistinguishable bosons or fermions subject to a confining potential, and to analyze the computed solutions. A general Hamiltonian has the form

$$\begin{eqnarray}&&\hat{H}=\sum _{i=1}^{N}\left[-\displaystyle \frac{1}{2}{\partial }_{{{\bf{x}}}_{i}}^{2}+\hat{V}({{\bf{x}}}_{i};t)\right]+\sum _{i\lt j}^{N}\hat{W}({{\bf{x}}}_{i},{{\bf{x}}}_{j};t),\end{eqnarray} \tag{ 1 }$$

that can be either explicitly dependent on time or not. Here, x_i is the coordinate of the ith particle, $-\tfrac{1}{2}{\partial }_{{\bf{x}}}^{2}$ is the kinetic energy operator, $\hat{V}$(x; t) is a general, possibly time-dependent, one-body potential, and $\hat{W}$(x, x'; t) is a general, possibly time-dependent interparticle interaction operator. The symbol ˆ denotes operators. All the quantities are given in dimensionless units. The length scale L can be chosen to appropriately represent the physical problem; the corresponding time and energy scales are then determined as ${{mL}}^{2}/{\hslash }$ and ${{\hslash }}^{2}/({{mL}}^{2})$, respectively. In particular, in the presence of a harmonic confinement potential, it is natural to choose the time scale as the inverse of the harmonic trapping frequency ω, i.e. $L=\sqrt{{\hslash }/(m\omega )}$.

The time-independent many-particle Schrödinger equation (TISE) corresponding to the Hamiltonian of equation (1) is

$$\begin{eqnarray}&&\hat{H}| {{\rm{\Psi }}}_{E}\rangle =E| {{\rm{\Psi }}}_{E}\rangle ,\end{eqnarray} \tag{ 2 }$$

while the time-dependent Schrödinger equation (TDSE) is

$$\begin{eqnarray}&&\hat{H}| {\rm{\Psi }}(t)\rangle ={\rm{i}}{\partial }_{t}| {\rm{\Psi }}(t)\rangle .\end{eqnarray} \tag{ 3 }$$

Note that $\hat{H}$ in the TISE needs to be a time-independent Hamiltonian. In equation (2), $| {{\rm{\Psi }}}_{E}\rangle$ is an eigenstate of $\hat{H}$ with eigenvalue (energy) E. $| {\rm{\Psi }}(t)\rangle$ stands for the solution of the TDSE at time t. Technically, the MCTDH-X software is currently capable of accurately computing few lowest-in-energy eigenstates using Davidson or short iterative Lanczos routine from the Heidelberg MCTDH package [67].

The MCTDH-X theory [26, 27] uses an ansatz for the wavefunction that is a time-dependent superposition of time-dependent many-body basis states:

$$\begin{eqnarray}\begin{array}{rcl}| {\rm{\Psi }}(t)\rangle & = & \sum _{\vec{n}}{C}_{\vec{n}}(t)| \vec{n};t\rangle ;\ \ \vec{n}={\left({n}_{1},\ldots ,{n}_{M}\right)}^{T};\\ | \vec{n};t\rangle & = & { \mathcal N }\prod _{i=1}^{M}{\left[{\hat{b}}_{i}^{\dagger }(t)\right]}^{{n}_{i}}| \mathrm{vac}\rangle ;\qquad {\phi }_{j}({\bf{x}};t)=\langle {\bf{x}}| {\hat{b}}_{j}(t)| 0\rangle .\end{array}\end{eqnarray} \tag{ 4 }$$

Here, the ${C}_{\vec{n}}(t)$ are referred to as coefficients, the $| \vec{n};t\rangle$ as configurations, and the normalization factor is ${ \mathcal N }=\tfrac{1}{\sqrt{{\prod }_{i=1}^{M}{n}_{i}!}}$ for bosons and ${ \mathcal N }=\frac{1}{\sqrt{N!}}$ for fermions. Each configuration is a fully symmetric or fully anti-symmetric many-body basis state built from M orthonormal time-dependent single-particle states, or orbitals, $\{{\phi }_{k}({\bf{x}},t);k=1,\ \ldots ,\ M\}$. The orbitals form an incomplete basis of the Hilbert space which is optimized to represent the many-body wavefunction with a minimal error. The orbitals may be transformed to the eigenfunctions of the reduced one-body density matrix, i.e., the so-called natural orbitals (see below in section 2.2), via a unitary transformation. To fully specify the solution of the TISE or TDSE, the MCTDH-X software computes and stores the coefficients ${C}_{\vec{n}}(t)$ and the orbitals $\{{\phi }_{k}({\bf{x}},t);k=1,\ \ldots ,\ M\}$ at times t that are specified by the user.

The set of equations of motion for the parameters in equation (4) comprises a coupled set of first-order differential equations for time-dependent coefficients ${C}_{\vec{n}}(t)$ and nonlinear integro-differential equations for the orbitals ϕ_j(x; t). The details about these equations of motion and their derivation can be found, for instance, in [3, 4, 26, 27]. The MCTDH-X software solves these equations of motion using the so-called constant mean-field integration scheme (see, for instance, [21, 26]). The constant mean-field scheme features an adaptive time step for which the equations of motion for the coefficients and the orbitals are decoupled.

We note that the equation for the orbitals contains the inverse of the matrix elements of the one-body density matrix. In cases where the reduced one-body density matrix has zero eigenvalues and is not invertible, the orbital equations are therefore undefined and problematic. In almost all practical cases and, particularly, for the computations we present below, the regularization strategy documented in [21]—a posteriori adding negligibly small eigenvalues to make the inversion possible—is sufficient. More elaborate schemes to improve or avoid this regularization have been developed [68, 69].

Once the coefficients ${C}_{\vec{n}}(t)$ and the orbitals ϕ_k(x, t) are computed, the MCTDH-X software can analyze the solution $| {\rm{\Psi }}(t)\rangle$ and calculate several quantities of interest. These include, respectively, the real-space and momentum-space density distributions

$$\begin{eqnarray}&&\rho ({\bf{x}})=\langle {\rm{\Psi }}| {\hat{{\rm{\Psi }}}}^{\dagger }({\bf{x}})\hat{{\rm{\Psi }}}({\bf{x}})| {\rm{\Psi }}\rangle /N,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\tilde{\rho }({\bf{k}})=\langle {\rm{\Psi }}| {\hat{{\rm{\Psi }}}}^{\dagger }({\bf{k}})\hat{{\rm{\Psi }}}({\bf{k}})| {\rm{\Psi }}\rangle /N,\end{eqnarray} \tag{ 6 }$$

the Glauber one-body and two-body correlation functions

$$\begin{eqnarray}&&{g}^{(1)}({\bf{x}},{\bf{x}}^{\prime} )=\displaystyle \frac{\langle {\rm{\Psi }}| {\hat{{\rm{\Psi }}}}^{\dagger }({\bf{x}})\hat{{\rm{\Psi }}}({\bf{x}}^{\prime} )| {\rm{\Psi }}\rangle }{N\sqrt{\rho ({\bf{x}})\rho ({\bf{x}}^{\prime} )}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{g}^{(2)}({\bf{x}},{\bf{x}}^{\prime} )=\displaystyle \frac{\langle {\rm{\Psi }}| {\hat{{\rm{\Psi }}}}^{\dagger }({\bf{x}}){\hat{{\rm{\Psi }}}}^{\dagger }({\bf{x}}^{\prime} )\hat{{\rm{\Psi }}}({\bf{x}}^{\prime} )\hat{{\rm{\Psi }}}({\bf{x}})| {\rm{\Psi }}\rangle }{{N}^{2}\rho ({\bf{x}})\rho ({\bf{x}}^{\prime} )},\end{eqnarray} \tag{ 8 }$$

and the natural orbitals ${\phi }_{i}^{(\mathrm{NO})}$ and the orbital occupations ρ_i, which, respectively, are the eigenfunctions and eigenvalues of the reduced one-body density matrix

$$\begin{eqnarray}&&{\rho }^{(1)}({\bf{x}},{\bf{x}}^{\prime} )=\displaystyle \frac{1}{N}\langle {\rm{\Psi }}| {\hat{{\rm{\Psi }}}}^{\dagger }({\bf{x}}^{\prime} )\hat{{\rm{\Psi }}}({\bf{x}})| {\rm{\Psi }}\rangle =\sum _{i}{\rho }_{i}{\phi }_{i}^{(\mathrm{NO}),* }({\bf{x}}^{\prime} ){\phi }_{i}^{(\mathrm{NO})}({\bf{x}}).\end{eqnarray} \tag{ 9 }$$

Here ${\hat{{\rm{\Psi }}}}^{\dagger }({\bf{x}})(\hat{{\rm{\Psi }}}({\bf{x}}))$ is a field operator that creates (annihilates) a particle at position x. The natural orbitals are ranked in the order that the first orbital has the highest occupation while the last orbital has the lowest ${\rho }_{1}\,\geqslant ...\geqslant \,{\rho }_{M}$. Another important quantity, the correlation order parameter (COP), is defined as the sum of the squares of the orbital occupations ρ_i [58, 59]

$$\begin{eqnarray}&&{\rm{\Delta }}=\sum _{i=1}^{M}{\rho }_{i}^{2}.\end{eqnarray} \tag{ 10 }$$

The Glauber correlation functions and the orbital occupations are related to each other. When only the first orbital is macroscopically occupied and thus Δ ≈ 1, the system is in a coherent state. In this case, both the one-body and two-body correlation functions for large N are unity, ${g}^{(1)}({\bf{x}},{\bf{x}})={g}^{(2)}({\bf{x}},{\bf{x}})=1$. The many-body wavefunction Ψ of the system is a product of the same single-particle wavefunction ψ, ${\rm{\Psi }}({{\bf{x}}}_{1},\ldots ,{{\bf{x}}}_{N})={\prod }_{i=1}^{N}\psi ({{\bf{x}}}_{i})$. When more than one orbital is macroscopically occupied and thus Δ < 1, the system is fragmented. In this case, the Glauber correlation functions become non-trivial and differ from unity. For the sake of brevity, above and in the following we omit the dependence of quantities on time unless when imperative or instructive.

Even more importantly, the MCTDH-X software can simulate single-shot images [55–59]. These single-shot images are the standard way to measure quantum many-body systems of ultracold atoms [70–73]. We note here that single-shot images were also recently simulated, chiefly for mixtures of different species of ultracold atoms, with ML-MCTDH-X, see [74] and references by Mistakidis et al. therein. By simulating such single-shot images, the MCTDH-X software can reproduce the quantum measurements in the laboratory. In real space, a single-shot image can be obtained by drawing random positions $({\tilde{{\bf{x}}}}_{1}$, ${\tilde{{\bf{x}}}}_{2}\,...{\tilde{{\bf{x}}}}_{N})$ distributed according to the probability

$$\begin{eqnarray}&&P({{\bf{x}}}_{1},\ \ldots ,\ {{\bf{x}}}_{N})={\left|{\rm{\Psi }}({{\bf{x}}}_{1},{{\bf{x}}}_{2},...,{{\bf{x}}}_{N})\right|}^{2}.\end{eqnarray} \tag{ 11 }$$

A similar definition holds in momentum-space too. Due to the presence of quantum fluctuations and correlations, a single-shot image, which is distributed according to $| {\rm{\Psi }}{| }^{2}$, can be drastically different than the density distributions ρ(x), $\tilde{\rho }({\bf{k}})$ in equations (5) and (6). The deviation of single-shot images from the density distributions is especially significant when the particle number is small and/or the quantum correlations are large.

A large collection of single-shot images can be used to provide information on the system and particularly on its correlations and phases. Below, we consider a total number of N_shot images, with the value of the ith image at position x given by ${{ \mathcal B }}_{i}({\bf{x}})$, and provide two types of single-shot analyses:

• 1.  
  Single shots for particle correlations between the two wells of a double-well potential: for each shot, the number of particles in one well, ${n}_{i}={\sum }_{{\bf{x}}}^{{\prime} }{{ \mathcal B }}_{i}({\bf{x}})$, is calculated. Here, Σ' indicates that the summation is within a certain well (left or right). We then calculate the probability of finding n = n_i particles in the considered well among all single-shot images

  $$\begin{eqnarray}&&P(n)=\displaystyle \frac{{N}_{n}({n}_{i}=n)}{{N}_{\mathrm{shot}}},\quad n\in \{0,1,2,\ \ldots ,\ N\},\end{eqnarray} \tag{ 12 }$$

  where N_n(n_i = n) denotes the total number of shots with n_i = n. We will show that the distribution of this probability depends on the correlations between particles.
• 2.  
  Quantum fluctuations can also be extracted from single-shot images [55–59]. To quantify the position-dependent quantum fluctuations of the particle number, we calculate the variance ${ \mathcal V }(x)$ from single-shot simulations as:

  $$\begin{eqnarray}&&{ \mathcal V }({\bf{x}})=\displaystyle \frac{1}{{N}_{\mathrm{shot}}}\sum _{i=1}^{{N}_{\mathrm{shot}}}{{ \mathcal B }}_{i}^{2}({\bf{x}})-{\left[\displaystyle \frac{1}{{N}_{\mathrm{shot}}}\sum _{i=1}^{{N}_{\mathrm{shot}}}{{ \mathcal B }}_{i}({\bf{x}})\right]}^{2}.\end{eqnarray} \tag{ 13 }$$

All these quantities above are chosen because of their accessibility in experiments and direct comparability to experimental results. For example, the momentum-space density distribution is accessible by time-of-flight measurements, the one-body particle correlations are accessible from thermodynamic quantities like kinetic energy [75], and single-shot images are the standard way of measuring cold-atom systems [70–73]. Other quantities of interest currently available in the MCTDH-X software include, for instance, the Glauber one-body and two-body correlation functions in momentum space [76] and many-body entropies of the system [77]. We note that the MCTDH-X software provides the full many-body wavefunction and therefore, in principle, any desired and computationally feasible quantity of interest can be computed with it.

The MCTDH-X software mainly consists of two programs:

• 1.  
  The main program, MCTDHX, computes the numerical solution $| {\rm{\Psi }}\rangle$ of the TISE or TDSE.
• 2.  
  The analysis program, MCTDHX_analysis, analyzes the found solution $| {\rm{\Psi }}\rangle$.

Here and in the following, we use the verbatim font to refer to code, including executable commands, files, statements, and variables.

To set up a numerical task, the user modifies and chooses the parameters via the text input file MCTDHX.inp. A detailed description of the available options is given in the manual [78] of the MCTDH-X software. The file Get_1bodyPotential.F is used to specify custom one-body potentials [$\hat{V}$(x_i; t) in equation (1)], while the file Get_Interparticle.F allows for custom two-body potentials [$\hat{W}$(x_i, x_j; t) in equation (1)]. For custom (initial) states, the Get_Initial_Coefficients.F and Get_Initial_Orbitals.F files can be used to specify the (initial) coefficients ${C}_{\vec{n}}(t={t}_{0})$ [see ${C}_{\vec{n}}(t)$ in equation (4)] and orbitals $\{{\phi }_{k}({\bf{x}};t={t}_{0}),k=1,\ \ldots ,\ M\}$ [see ϕ_j(x; t) in equation (4)], respectively.

The workflow and structure of the MCTDH-X software follows naturally from its main objectives to determine a numerical solution to the TISE or the TDSE and then to extract desired quantities of interest from the solution. This workflow can be summarized in the following steps, which are also visualized in figure 1:

Zoom In Zoom Out Reset image size

• 1.  
  Determine the initial state.

  • (a)  
    Default: compute the ground state of some Hamiltonian $\hat{H}$ by running MCTDHX and configuring the numerical task ('relaxation mode', Hamiltonian, integration procedure etc) in the input file MCTDHX.inp.
  • (b)  
    Advanced: manually set the coefficients and orbitals that determine the initial wavefunction via the files Get_Initial_Coefficients.F and Get_Initial_Orbitals.F
• 2.  
  Analyze the initial state by choosing the desired quantities of interest in analysis.inp and running the analysis program MCTDHX_analysis

  • (a)  
    Default: call supplied visualization scripts or gnuplot to obtain plots or videos of the results
  • (b)  
    Advanced: create custom visualizations or customize supplied visualization scripts
• 3.  
  Compute the time-evolution of the initial state with given Hamiltonian $\hat{H}$. The dynamics of the system is obtained by choosing the numerical task ('propagation mode', Hamiltonian, integration procedure etc.) in the input file MCTDHX.inp and running MCTDHX.
• 4.  
  Analyze the computed time-evolving state by choosing the desired quantities of interest in analysis.inp and running the analysis program MCTDHX_analysis. Visualization of the results as in step 2.(a) and 2.(b).

The scripts in the MCTDH-X software generally fall into two categories: either scripts to automate computations like parameter scans or scripts to visualize data in plots and videos. A series of tutorial videos illustrating the workflow step-by-step is available at https://tinyurl.com/tjx35sq. Convergence tests on different variables, including orbital number and grid size, are needed to ensure the accuracy of the simulations. A discussion on the convergence, particularly in orbital number, can be found in the supplementary material available online at stacks.iop.org/QST/5/024004/mmedia ⁹ (including [79–84]). For support and documentation, the website of the MCTDH-X software, http://ultracold.org [37], the MCTDH-X manual [78], and the email address mctdhx@ultracold.org are available. Feature requests should be directed towards the support email address and/or be discussed on the web forum.

We first solve the relevant TISE to find the ground state of N = 6 bosons or fermions in the double-well potential of equation (14) by propagating the MCTDH-X coefficients ${C}_{\vec{n}}(t)$ and basis states ϕ_i(x; t) in imaginary time [see equation (4)]. We vary the barrier heights E_dw and interaction strengths g to investigate where the superfluid-Mott insulator transition happens. This is done by adapting the input MCTDHX.inp and running the program MCTDHX. The number of orbitals is chosen to be M = 10. The density distributions in position space and momentum space and the one-body correlation function of the ground states in different potentials are shown in figure 2.

A bosonic superfluid state and a bosonic Mott-insulator state are shown in figures 2(a)–(f), respectively. Compared to the Mott-insulator state, the superfluid state has two extra peaks in momentum space and a non-zero one-body correlation between the two wells g⁽¹⁾(x, −x) > 0. The comparison between these two kinds of states has been investigated thoroughly with the MCTDH-X software [57, 62, 63, 95] and the results are consistent with other theoretical and experimental results [88, 96, 97]. In this tutorial, we use the terms 'Mott insulating' or 'Mott insulation' to refer to situations where the one-body correlation function [equation (7)] has vanishing values in off-diagonal blocks, i.e. situations where $| {g}^{(1)}(x,x^{\prime} )| \approx 0$ is true for positions where x and x' are in distinct wells or at the position of distinct peaks of the one-body density, see figure 2(f), for instance.

As the interaction between the bosons further increases, it induces a self-organized lattice-like structure of the density and correlations even within each of the double-well sites [figures 2(g)–(i) for g = 20]. This emergence of structure heralds the onset of so-called fermionization; the real-space densities of bosons with large contact interactions approach those of non-interacting fermions [98]. The reason for the emergence of fermionization is that two bosons with infinitely large repulsive contact interactions, like fermions, cannot pass through each other in a one-dimensional system [98]. Such similarities lay the foundation of useful tools like bosonization [99, 100].

However, fermionization and fermionic nature are driven by different mechanisms. This can be intuitively understood through their many-body wavefunctions. The wavefunction of a fermionized state of bosons, $| {{\rm{\Psi }}}_{B}\rangle$, is still symmetric, while the wavefunction of a non-interacting fermionic state, $| {{\rm{\Psi }}}_{F}\rangle$, is given by an anti-symmetric Slater determinant built from the first N single-particle eigenstates. In position-space representation $| {{\rm{\Psi }}}_{B}{| }^{2}\equiv | {{\rm{\Psi }}}_{F}{| }^{2}$ holds true in the fermionization limit; the wavefunctions of bosons and fermions, however, are completely different ${{\rm{\Psi }}}_{B}({x}_{1},{x}_{2},\ \ldots ,\ {x}_{N})={\prod }_{i\lt j}\mathrm{sgn}({x}_{i}-{x}_{j}){{\rm{\Psi }}}_{F}({x}_{1},{x}_{2},\ \ldots ,\ {x}_{N})$ [98]. This difference is reflected in distinct features of the momentum distributions of bosons in the fermionization limit and non-interacting fermions [see figure 2(b), (e), (h), (k), (n) and discussion below].

The density distribution at g = 20 in figure 2(g) is similar to the one of non-interacting fermions in figure 2(j). The one-body correlation function of the fermions has a complicated pattern of significant and non-trivial correlations [figure 2(l)]. The bosons have a slightly different correlation function than the fermions, but the distinct fermionic pattern is already visible. We thus refer to such states as 'en route to fermionization'. As discussed in [95, 101] and the supplementary material, if an extremely large repulsive interaction and an adequate number of orbitals are used, the fermionization limit can be accurately captured by the MCTDH-X approach [see stacks.iop.org/QST/5/024004/mmedia and footnote ⁹].

The difference between the bosons en route to fermionization and the fermions, however, appears most explicitly in momentum space. The momentum distribution for the bosons en route to fermionization [figure 2(h)] has the same single-peak structure as the normal Mott insulating bosons [figure 2(e)], where both widths and heights of the peaks are similar. This confirms the similarity between fermionized bosons and Bose–Einstein condensates in momentum space [98]. In contrast, the fermionic state has three peaks in its momentum distribution due to the Pauli principle [figure 2(k)]. These peaks correspond to the three fermions in each well since the two wells are Mott insulating [figure 2(l)].

For long-range dipolar interactions crystallization emerges; for this case, bosons and fermions have been compared, for instance, in [102]. A crystallized bosonic state and a crystallized fermionic state have similar real-space density distributions and they both have many contributing eigenvalues in the reduced one-body density matrix. This indicates that the long-range interaction dominates over the particle statistics. We note that crystallized bosons and fermionized bosons can be distinguished via the spread of their density matrices as a function of the interparticle interaction strength [101].

A fermionic state crystallizes in the presence of sufficiently strong long-range interactions g = 5 [see equation 15(b)]. As expected, the repulsive interaction increases the distance between the fermions in real space, see figure 2(m). The interactions also have a pronounced impact in the momentum space distribution [figure 2(n)] and the particle correlations, see figure 2(o). The peaks in the real-space density distribution within each of the wells, which are induced by the Pauli exclusion principle in the absence of interaction [see figures 2(j), (l)], become Mott insulating for crystallized fermions [see figures 2(m), (o)].

The correlations and fluctuations of particles can also be revealed by (simulations of) single-shot images [55–59]. For an illustration of what can be extracted from simulated single-shots, we fix the barrier height E_dw = 20 and compute bosonic and fermionic ground states for various interactions g. For every computed state, we generate 10 000 simulated single-shots (adapting analysis.inp and running MCTDHX_analysis). For each set of single-shot images, we calculate the frequency of n = 0, ..., 6 particles being in the left well, P(n), and show it in figure 3. For the fermions, each well contains exactly half of the particles regardless of the presence of interactions. This can be seen as a consequence of the Pauli principle. For the bosons, in the superfluid limit g = 0, the system is in a coherent state [g⁽¹⁾ ∼ 1 for all x, x' in figure 2(c)]. In such a coherent state, the distributions of all bosons are independent of each other. As a result of this and the symmetry of the potential, P(n) should be given by the binomial distribution $B(N,1/2)=\left(\displaystyle \genfrac{}{}{0em}{}{N}{n}\right){/2}^{N}$ and this is confirmed by the single-shot results. As superfluidity gives way to Mott insulation, the distributions of different bosons become interdependent. Due to the repulsive interaction, the Hubbard model predicts that the particles tend to be distributed evenly in each well and P(3) gradually increases while $P(n\ne 3)$ gradually decreases. For our MCTDH-X results in the Mott-insulator limit, g = 1, E_dw = 20, the repulsion between particles becomes so strong that there are exactly three particles in each well P(3) = 1 as predicted by the Hubbard model.

Zoom In Zoom Out Reset image size

The behavior of the system in the superfluid, Mott-insulating crystallized and phases can be summarized and represented by the COP Δ [equation (10)] [58, 59]. The dependence of the COP on the barrier height E_dw and interaction strength g in bosonic and fermionic systems is shown in figure 4. In the bosonic system, the COP decreases from almost unity in the superfluid phase to 0.5 in the Mott insulating phase, as expected for a fragmented state with ρ₁ ≈ ρ₂ ≈ 0.5 [figures 4(a) and (b)]. As the interaction strength increases and the system is en route to fermionization, the COP drops further to Δ ≈ 0.2. In the fermionic system, Δ is much less sensitive to the interaction strength. It only drops very slightly when the fermions crystallize [figure 4(c)]. The value Δ ≈ 0.167 indicates each fermion occupies a single orbital, agreeing with the Pauli principle.

Zoom In Zoom Out Reset image size

MCTDH-X is capable of solving the TDSE to capture the dynamics of a system as a reaction to a time-dependent Hamiltonian or to a quench of a parameter. As an example, we prepare the system in the ground state of a harmonic trap, ${V}_{\mathrm{har}}(x)=\tfrac{1}{2}{x}^{2}$ and subsequently, we ramp up a barrier at x = 0. We thus smoothly transform the harmonic trap into the double-well potential given in equation (14). We use a linear ramp with a time scale τ

$$\begin{eqnarray}{E}_{\mathrm{dw}}(t)=\left\{\begin{array}{ll}{E}_{\max }t/\tau , & t\leqslant \tau \\ {E}_{\max }, & t\gt \tau .\end{array}\right.\end{eqnarray} \tag{ 16 }$$

In the following, we compare the states that result for different τ. There are two situations where the resulting state can be different from the ground state of the instantaneous Hamiltonian. (i) If the system undergoes a first-order phase transition—while the parameters of the Hamiltonian change smoothly, the system can get stuck in an excited state of the final Hamiltonian in the absence of a strong perturbation. (ii) With a fast and non-adiabatic change of parameters in the Hamiltonian, excitations are inevitably triggered. Such a quench thus may result in a finite occupation of a large number of eigenstates of the final Hamiltonian [103]. To investigate these two scenarios, we now discuss two examples where the system contains N = 5 or N = 6 weakly-interacting bosons or fermions.

A first-order transition between a superfluid and a Mott insulator has already been observed in the presence of a three-body interaction [104] in a bosonic system. Such a first-order transition is also observed in our simulations even without three-body interactions, when the total number of bosons or fermions is odd (N = 5 in our simulations). There is a large residual correlation $| {g}^{(1)}(x,-x)| \approx 0.6$ in the time-evolving state—even in the slow evolution limit with ramping time τ = 100 [figures 5(b), (d)] for both the bosonic and the fermionic cases, compared to the ground state where the one-body correlations between the two wells vanish completely [x, x' where $| {g}^{(1)}| \approx 0$ in figures 5(a), (c)].

Zoom In Zoom Out Reset image size

The interaction between particles is an essential ingredient in this first-order transition. It lifts the degeneracy between the Mott-insulating state and the state with large residual correlations. However, as the barrier of the double well increases, one of the particles has no preference for entering either of the two wells and thus straddles across both wells. This straddling particle builds up the correlations between the two wells that we observe.

To justify the claim of a first-order transition, we investigate the hysteretical behavior of one-body correlations $| {g}^{(1)}(1,-1)|$ for N = 5 bosons, as shown in figure 5(e). In the ground state, there is clearly a jump at roughly E_dw ≈ 11. To study the dynamical effects, we choose a ground state in the superfluid limit, E_dw = 0 and a ground state in the Mott-insulator limit, E_dw = 20. We propagate them in 'opposite directions' across the phase boundary by slowly ramping up (down) the barrier from E_dw = 0 to E_dw = 20 (E_dw = 20 to E_dw = 0) in a time of τ = 100. It is clear that there is a hysteresis in the particle correlations. Strikingly, the hysteresis covers an infinite area [see the yellow marked region in figure 5(e)], as the residual correlations will never vanish in the Mott limit. The jump in the correlation functions in the ground states and the hysteresis in the time-evolving states are not seen with an even number of particles [figure 5(f)]. In this even-number case, there is thus a second-order phase transition.

Since the system follows the ground state in the slow limit (large τ) with an even number of particles, we now investigate the effect of a quench (smaller τ) using N = 6 bosons or fermions. We ramp up the barrier from zero to E_dw = 20 in three different ramping times, τ = 1, τ = 10 and τ = 50, to which we refer as a fast, an intermediate, and a slow ramp, respectively, in the following. Their one-body correlation functions at a given time, t = 80, are shown in figures 6(a)–(c) for bosons and in figures 6(g)–(i) for fermions.

Zoom In Zoom Out Reset image size

With a slow ramp (τ = 50), the ground state of the final Hamiltonian is recovered—the correlation between the two wells vanishes [see areas where x is to the left (right) of the barrier and x' to its right (left) and $| {g}^{(1)}(x,x^{\prime} )| \approx 0$ in figures 6(c), (i)]. As the ramping time becomes shorter τ = 10, fluctuations appear on the correlation function, indicating excitations [see areas where x is to the left (right) of the barrier and x' to its right (left) and $| {g}^{(1)}(x,x^{\prime} )| \gt 0$ in figures 6(b), (h)]. These excitations accumulate rapidly as the ramping time shortens and, eventually, with a fast ramping τ = 1, we obtain a strongly fluctuating one-body correlation function [figures 6(a), (g)].

To analyze how the quantum fluctuations in the time-evolving states are affected by the speed of the ramp, we quantify the density fluctuations using the position-dependent variance ${ \mathcal V }(x)$ [equation (13)] extracted from 10 000 single-shot simulations for bosons in figures 6(d)–(f) and for fermions in figures 6(j)–(l).

Generally, we find the fluctuations are large where the density is large. For bosons, a slower ramp results in an increase of the magnitude of fluctuations [compare figures 6(d)–(f)]. For fermions, the magnitude of fluctuations is practically unaffected by the pace of the ramp [compare figures 6(j)–(l)]. We infer that the behavior of the magnitude of the fluctuations is analogous to the COP Δ [equation (10)]: for bosons, the COP increases a lot as the ramping rate decreases in long time t [see figure 5(g)], whereas the COP varies only very slightly at different ramping rates for the case of fermions [see figure 5(h)]. This analogous behavior of the quantum fluctuations and the COP has previously been used to extract the phases of dipolar ultracold atoms in lattices [58, 59].

The excitations generated by the quench thus also manifest themselves in the eigenvalues of the reduced one-body density matrix or orbital occupations as represented by the COP and its time-dependence [figures 5(g), (h)]. For bosons, in the case of a slow ramp, only two orbitals are macroscopically occupied as in the ground state, while in the case of a fast ramp, the contribution of the first 10 orbitals are of the same order of magnitude and the system is thus highly correlated. The extra fragmentation comes from the dynamics of the system.

For bosons, in the slow case, the COP decreases as the barrier increases and converges to Δ = 0.5 as in the ground state [see figure 4(a)]. As the ramping becomes faster, the COP starts to oscillate and the amplitude becomes larger and larger. However, the dynamical behavior of the COP becomes qualitatively different as the ramping time becomes shorter than τ = 10. It no longer oscillates in a regular manner but rather decreases rapidly and fluctuates. This implies that many high-energy eigenstates of the final Hamiltonian are populated in the dynamics of the system.

For fermions, the absolute values of the orbital occupations are not drastically influenced by the system dynamics and neither is the COP. However, the dynamical behavior of the COP is completely different in the fast and slow ramping limits—in contrast to the bosonic case, the COP increases slightly for a slow ramp, moving closer to the value for non-interacting fermions Δ = 0.167. This indicates that the effects of the interparticle interactions become weaker since the fermions are now divided into two well-separated groups.

The results for the one-body correlation, the COP, and the single-shot simulations point towards a dynamical phase transition at roughly τ_C ∼ 10 for fermions and bosons alike. The system behaves like the ground state for ramps slower than τ_C, but becomes highly excited for ramps faster than τ_C.