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

We introduce and describe the multiconfigurational time-depenent Hartree for indistinguishable particles (MCTDH-X) software, which is hosted, documented, and distributed at http://ultracold.org. This powerful tool allows the investigation of ground state properties and dynamics of interacting quantum many-body systems in different spatial dimensions. The MCTDH-X software is a set of programs and scripts to compute, analyze, and visualize solutions for the time-dependent and time-independent many-body Schrödinger equation for indistinguishable quantum particles. As the MCTDH-X software represents a general solver for the Schrödinger equation, it is applicable to a wide range of problems in the fields of atomic, optical, molecular physics, and condensed matter systems. In particular, it can be used to study light–matter interactions, correlated dynamics of electrons in the solid state as well as some aspects related to quantum information and computing. The MCTDH-X software solves a set of nonlinear coupled working equations based on the application of the time-dependent variational principle to the Schrödinger equation. These equations are obtained by using an ansatz for the many-body wavefunction that is a expansion in a set of time-dependent, fully symmetrized bosonic (X = B) or fully anti-symmetrized fermionic (X = F) many-body basis states. It is the time-dependence of the basis set that enables MCTDH-X to deal with quantum dynamics at a superior accuracy as compared to, for instance, exact diagonalization approaches with a static basis, where the number of basis states necessary to capture the dynamics of the wavefunction typically grows rapidly with time. Herein, we give an introduction to the MCTDH-X software via an easy-to-follow tutorial with a focus on accessibility. The illustrated exemplary problems are hosted at http://ultracold.org/tutorial and consider the physics of a few interacting bosons or fermions in a double-well potential. We explore computationally the position-space and momentum-space density, the one-body reduced density matrix, Glauber correlation functions, phases, (dynamical) phase transitions, and the imaging of the quantum systems in single-shot images. Although a few particles in a double well potential represent a minimal model system, we are able to demonstrate a rich variety of phenomena with it. We use the double well to illustrate the fermionization of bosonic particles, the crystallization of fermionic particles, characteristics of the superfluid and Mott-insulator quantum phases in Hubbard models, and even dynamical phase transitions. We provide a complete set of input files and scripts to redo all computations in this paper at http://ultracold.org/data/tutorial_input_files.zip, accompanied by tutorial videos at https://tinyurl.com/tjx35sq. Our tutorial should guide the potential users to apply the MCTDH-X software also to more complex systems.


I. INTRODUCTION
To date, the time-dependent many-body Schrödinger equation (TDSE) is a fundamental equation at the heart of many different fields of science: quantum chemistry, condensed matter, and atomic and molecular physics. Exact solutions to the TDSE exist only for model systems, like the time-dependent harmonic interaction model [1][2][3]. Even for the time-independent Schrödinger equation (TISE), exact solutions are scarce [4][5][6][7][8][9][10][11]. These exact solutions, however, are not generalizable to current experiments or theoretical studies. To obtain solutions to the TISE and TDSE, numerical methods as well as their implementation in software are therefore a must. Due to the fundamental nature of the problem, many methods have been put forward -each with their advantages and shortcomings: Matrix-product-state and density matrix renormalization group theory approaches employ a hierarchical partitioning of the many-body Hilbert space which is particularly well-suited for one-dimensional lattice systems [12,13]. Mean-field approaches like the time-dependent Hartree-Fock [14] and the time-dependent Gross-Pitaevskii methods [15][16][17] employ a drastically simplified approximation to the wavefunction of the state that ignores correlation effects.
Since the 1990s, the multiconfigurational time-dependent Hartree approach [18][19][20] (MCTDH) was put forward and applied very successfully in the field of theoretical chemistry, where systems involve coupled, distinguishable degrees of freedom. MCTDH enables the description of correlated wavefunctions; its ansatz for the wavefunction is a sum of all possible configurations of distinguishable degrees of freedom or particles in a set of time-dependent variationallyoptimized single-particle functions. In 2003 MCTDH was formulated for indistinguishable fermions [21][22][23] (MCTDH-F) and in 2007 for indistinguishable bosons [24,25] (MCTDH-B). MCTDH-F can be used to describe the correlated dynamics of electrons in atoms and molecules [26][27][28][29][30], or to describe ultracold atomic fermions [3]. MCTDH-B is chiefly used to describe the many-body properties of ultracold atomic bosons with a focus on fragmentation [31,32], where the reduced density matrix of the many-boson state attains several significant eigenvalues [33][34][35][36][37][38] as well as quantum fluctuations [39][40][41][42][43]. MCTDH-F and MCTDH-B can be formulated in a unified manner [44], where they share the same equations of motion; henceforth we will use the acronym MCTDH-X to refer to this unified prescription.
In this article, we provide an introduction to the software implementation of MCTDH-X hosted and distributed at http://ultracold.org [3,45,46]. Below, we merely provide a tutorial-type introduction to the MCTDH-X software with a focus on simplicity and instructiveness. The MCTDH-X software, however, can do way more than the examples we introduce below. The MCTDH-X software can deal with indistinguishable particles with internal degrees of freedom like spin [46], indistinguishable particles immersed in a high-finesse optical cavity [41,[47][48][49][50], indistinguishable particles with long-range dipolar interactions [42,43,51,52], and Hubbard models [53]. Moreover, the MCTDH-X software provides the possibility for an in-depth analysis of the computed solutions of the TISE and TDSE via full distribution functions [43], variances and quantum fluctuations [38,39,43] of observables, and correlation functions [2,38,43,54]; the MCTDH-X has been benchmarked against exact results [1,2], verified against experimental predictions [39], and was recently reviewed [55].
The objective, workflow and usage of the MCTDH-X software will be introduced in Sec. II and exemplified by a detailed tutorial in Sec. III, where ground states and dynamics of both bosons and fermions are inspected. Our focus is on introducing the usage of the software; details about the MCTDH-X theory are only briefly discussed where necessary.

II. OBJECTIVE, STRUCTURE, AND WORKFLOW OF THE MCTDH-X SOFTWARE
A. Objective and main functionality of MCTDH-X The objective of the MCTDH-X software is to numerically solve the TISE or TDSE with a many-body Hamiltonian, which describes N interacting, indistinguishable bosons or fermions subject to a confining potential, and to analyze the computed solutions. The Hamiltonian has the form Here and in the following, all the quantities are given in dimensionless units ( = m = 1): x i is the coordinate of the i-th particle, − 1 2 ∂ 2 x is the kinetic energy, V (x) is a general, possibly time-dependent, one-body potential, and W (x, x ) is a general, possibly time-dependent interparticle interaction.
The HamiltonianĤ can in general be dependent or independent on time. The corresponding time-dependent or time-independent many-particle Schrödinger equation (TDSE/TISE) read respectively, andĤ whereĤ 0 is of the same type asĤ, but explicitly time-independent; |Ψ E is an eigenstate ofĤ 0 with eigenvalue, i.e. energy, E, and |Ψ(t) stands for the solution of the TDSE at time t. In the TISE, we used the symbolĤ 0 for the Hamiltonian to indicate that it is generally distinct from the Hamiltonian used for the TDSE and, moreover,Ĥ 0 has to be time-independent in order for the TISE to be a time-independent equation. Technically, MCTDH-X is currently only capable of computing a few lowest-in-energy eigenstates using the Davidson or short iterative Lanczos routine from the Heidelberg MCTDH package.
The MCTDH-X theory [25,44] uses an ansatz for the wavefunction that is a time-dependent superposition of time-dependent many-body basis states: |Ψ(t) = n C n (t)| n; t ; n = (n 1 , ..., n M ) T ; Here, the C n (t) are referred to as coefficients and the | n; t as configurations. Each configuration is a fully symmetric or fully anti-symmetric many-body basis state built from M orthonormal time-dependent single-particle states, {φ k (x, t); k = 1, ..., M }, or orbitals. To fully specify the solution of the TISE or TDSE, the MCTDH-X software computes and stores the coefficients C n (t) and the orbitals {φ k (x, t); k = 1, ..., M } at times t that are specified by the user.

B. Observables in MCTDH-X
For the sake of simplicity and brevity, we will omit the dependence of quantities on time where it is not imperative to explicitly mention it. Once the coefficients C n (t) and the orbitals φ k (x, t) are computed, the MCTDH-X software can analyze the solution and calculate a series of observables. These include, respectively, the real-space and momentum space density distributions the Glauber one-body and two-body correlation functions, and the orbital occupations ρ i , which are the eigenvalues of the reduced one-body density matrix, Moreover, the MCTDH-X software can also perform simulations of single-shot images [39][40][41][42][43]. In each single-shot image s, all N particles are drawn at the positions (x 1 ,x 2 ,...,x N ) distributed according to the probability A collection of these single-shot images is used to provide information on the phases. These single-shot images are used in two different ways in this tutorial.
In the first way, it is used to measure the particle correlations between the two wells of a double-well potential. For each shot, the number of particles in one well is calculated n i = x B i (x), where Σ indicates that the summation is within a certain well, B i (x) is the value of single-shot image i at position x, and N shot is the total number of single-shot images considered. We then calculate the probability of n i among all single shot images P (n) = N (n i = n)/N shot , n ∈ {0, 1, 2, . . . , N } (10) where N (n i = n) denotes the total number of shots where n i = n. The distribution of this probability depends on the particle correlations.
In the second way, it is used to extract the quantum fluctuations of the considered system [39][40][41][42][43]. To quantify the position-dependent quantum fluctuations of the particle number, we will calculate the variance V(x) from single-shot simulations: Many of these observables are accessible in experiments. For example, the momentum space density distribution is accessible by time-of-flight measurements, the one-body particle correlations are accessible by thermodynamic quantities like kinetic energy [56], and the single-shot images are accessible by single-shot experiments [39][40][41][42][43].

C. Structure and Workflow of the MCTDH-X software
The two main programs in the MCTDH-X software are for (i) computing the numerical solution of the TISE or TDSE (MCTDHX) and for (ii) analyzing said solution (MCTDHX analysis). Here and in the following, we use the verbatim font to refer to code, commands/executables, or files and changes to them.
To set up the numerical task to be solved by running the command MCTDHX, the input file MCTDHX.inp can be used -a detailed description of the available options is given in the manual [57] of the MCTDH-X software. The file Get 1bodyPotential.F can be used to specify custom one-body potentials [V (x i ; t) in Eq. (1)], while the file Get Interparticle.F allows for custom two-body potentials [W (x i , x j ; t) in Eq. (1)]. For specifying custom (initial) states the Get Initial Coefficients.F and Get Initial Orbitals.F files can be used to specify the (initial) coefficients at time 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 observables from the solution. This workflow can be summarized in the following steps, which we also visualize in Fig. 1: 1. Determine the initial state.
(a) Default: computing the ground state of some HamiltonianĤ 0 by running MCTDHX and configuring the numerical task ("relaxation mode", Hamiltonian, integration procedure, etc.) in the input file MCTDHX.inp.
(b) Advanced: manually setting 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 observables in analysis.inp and running the analysis program MCTDHX analysis  The scripts in the MCTDH-X software generally fall into two categories: either scripts to automate computations, i.e., parameter scans, or scripts to visualize data in plots and videos. For support and documentation, the website of the MCTDH-X software, http://ultracold.org [45], the MCTDH-X manual [57], and the email address mct-dhx@ultracold.org are available. Feature requests should be directed towards the email address and discussed on the web forum. j e Q X j h P U i c i 1 5 y C k B I / n 2 w Z m f e c B G E I S Z z P M K V P E + 9 k w l A + 4 P c e X L G + X 7 G K q + X X Z q z h T 4 L 3 E L U k Y F L n z 7 z e v H N I 2 Y B C q I 1 l 3 X S a C X E Q W c C p a X v F S z h N A h u W Z d Q y W J m O 5 l 0 7 N y v G e U P g 5 j Z V o C n q r f N z I S a T 2 O A j M Z E R j o 3 9 5 E / M / r p h C e 9 D I u k x S Y p J 8 P h a n A E O N J R r j P F a M g x o Y Q q r j 5 K 6 Y D o g g F k 2 T J h O D + P v k v a d V r 7 m G t f n l U b p w W c S y h H b S L K s h F x 6 i B z t E F a i K K b t E 9 e k R P 1 p 3 1 Y D 1 b L 5 + j M 1 a x s 4 1 + w H r 9 A D r B n i Q = < / l a t e x i t >

Hamiltonian for eigenstatê
< l a t e x i t s h a 1 _ b a s e 6 4 = " J X + P X Q R + L g j q c t B n O 2 a l Y o G 5 o M 0 = " > A A A C G H i c b V D L S s N A F J 3 4 r P V V d e l m s A i u a l I F 3 Q h F K X R Z w T 6 g C W E y n b R D J 5 M w c y O U 2 M 9 w 4 6 + 4 c a G I 2 + 7 8 G 6 e P h b Y e u H A 4 5 1 7 u v S d I B N d g 2 9 / W y u r a + s Z m b i u / v b O 7 t 1 8 4 O G z q O F W U N W g s Y t U O i G a C S 9 Y A D o K 1 E 8 V I F A j W C g Z 3 E 7 / 1 y J T m s X y A Y c K 8 i P Q k D z k l Y C S / c O 7 2 C W S 1 k W / j J + z W N f e r 2 F V E 9 g T D N 7 i 6 J P q F o l 2 y p 8 D L x J m T I p q j 7 h f G b j e m a c Q k U E G 0 7 j h 2 A l 5 G F H A q 2 C j v p p o l h A 5 I j 3 U M l S R i 2 s u m j 4 3 w q V G 6 O I y V K Q l 4 q v 6 e y E i k 9 T A K T G d E o K 8 X v Y n 4 n 9 d J I b z 2 M i 6 T F J i k s 0 V h K j D E e J I S 7 n L F K I i h I Y Q q b m 7 F t E 8 U o W C y z J s Q n M W X l 0 m z X H I u S u X 7 y 2 L l d h 5 H D h 2 j E 3 S G H H S F K q i G 6 q i B K H p G r + g d f V g v 1 p v 1 a X 3 N W l e s + c w R + g N r / A P d D p 5 p < / l a t e x i t >

Initial state
< l a t e x i t s h a 1 _ b a s e 6 4 = " Z Q S 3 q D H j f S e a 9 S l J z r k D x j d m 9 9 I = "

III. TUTORIAL APPLICATION OF THE MCTDH-X SOFTWARE
As a tutorial example, we use the MCTDH-X software to study one-dimensional bosons and fermions in singleand double-well one-body potentials. Importantly, seen as a lattice with two sites, the double-well potential can be used as a demonstration on how we can use the MCTDH-X software to reveal the static and dynamical properties of the building blocks of lattice systems, i.e., periodic structures of potential wells. This choice reduces computational difficulty and improves the clarity of our presentation. In the following, we use the notation x instead of x.
The TISE and TDSE for lattices are a cornerstone in many different physical systems in condensed matter physics and ultracold atoms. Typically, the Hamiltonian for lattices is approximated by the Hubbard Hamiltonian (see [58]); this approximation uses a basis of static, suboptimal site-localized Wannier functions. MCTDH-X goes beyond the Hubbard approximation and uses a basis set which is time-dependent, variationally-optimized, and not necessarily localized at sites. MCTDH-X thus enables a modelling of the many-body state with an improved accuracy as compared to the Hubbard-lattice-models which presuppose the single-particle states are static Wannier functions [see supplementary material [58]]. This improved accuracy is of particular importance in cases where the single-band Hubbard description fails, for instance, due to the action of a blue-detuned cavity [49] or long-range dipolar interactions [42,43,52]. Direct comparisons between the MCTDH-X and Hubbard models have been performed in Refs. [2,34,36,50,[59][60][61][62][63][64].
The double-well potential is a good demonstration of the MCTDH-X software because it has a simple Hamiltonian, is experimentally relevant [65][66][67], and captures well the superfluid-Mott insulator transition which is amply-studied in bosonic Hubbard models [68,69]. Moreover, the double-well is simple enough to showcase physical systems where a Hubbard description breaks down and the improved accuracy that MCTDH-X can deliver is of importance, see Refs. [34,36,64] and Fig. 3 below. However, the double-well model fails to capture many aspects of a lattice, especially those which are intrinsic to a lattice with a large number of sites, like the Peierls transition, the SSH model [70,71], etc. We defer the discussion of the Hubbard model to the supplemental material [58].
In a superfluid state, all the bosons occupy one single-particle state, i.e., |Ψ ∼ |N, 0, ...; t . In a Mottinsulating state, all bosons in one potential well occupy the same single-particle state, i.e., for a double well |Ψ ∼ |N/2, N/2, 0, ...; t (for even N ). In MCTDH-X, to correctly simulate a Mott insulator state each lattice site requires its own orbital and the required number M of orbitals is at least as large as the number of lattice sites. If the number of orbitals is inadequate and two lattice sites have to "share" one orbital, the emergence of Mott-type correlations between these two sites cannot be captured correctly.
In our examples, we follow and illustrate the workflow described in Sec. II. We investigate first the ground states and, subsequently, the time-evolution of a few bosonic or fermionic particles. All the input files and scripts in this section are hosted at http://ultracold.org/data/tutorial and a detailed description on how to do the computations is given in supplementary material [58].
For bosons, we consider contact interactions W = W B in Eq. (1). For fermions, we consider the regularized Coulomb interaction, W = W F in Eq. (1): where g > 0 is the repulsive interaction strength, δ(x) is the Dirac delta distribution, and the parameters α = 0.1 and β = 100 are used as in Refs. [3,72]. We model the double-well potential with an external harmonic confinement and a central Gaussian barrier: Here, E dw controls the barrier height. The two minima to the left and to the right of the Gaussian barrier correspond to two lattice sites. The hopping strength of an analogous Hubbard model, t, is mainly controlled by the barrier height, E dw , while the on-site interaction, U , is mainly determined by interaction strength, g (see supplementary material [58] for reference about the Hubbard model).
The detailed procedure of the modification of the input file for the simulations (MCTDHX.inp) and the extraction of observables (analysis.inp) is described in the supplementary material [58]. These observables contain sufficient information to extract useful properties of the quantum state and reveal phenomena like the superfluid to Mott-insulator transition, the fermionization of bosons, crystallization of fermions, the order and presence of phase transitions, and the excitations induced by a quench.

A. Ground state properties
We first compute the ground state of N = 6 bosons or fermions in the double-well potential of Eq. (13) via imaginary-time propagation of the MCTDH-X working equations. We vary the barrier heights E dw and interaction strengths g. The number of orbitals is chosen as M = 10 such that an accurate representation of the state is obtained. 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 Fig. 2.
A bosonic superfluid state and a bosonic Mott-insulator state are shown in Fig. 2(a-c) and Fig. 2(d-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 (1) (x, −x) > 0. The comparison between these two kinds of states has been investigated thoroughly with MCTDH-X [41,49,50,73], and the results are consistent with previous theoretical and experimental results [69,74,75]. In this tutorial, we use the terms "Mott insulating" or "Mott insulation" to refer to situations where the one-body correlation function [Eq. (6)] has vanishing values in off-diagonal blocks, i.e., situations where g (1) (x, x ) ≈ 0 is vanishing for positions where x and x are in distinct wells or at the position of distinct peaks of the one-body density, see Fig. 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 [ Fig. 2(g-i)]. 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 [76]. 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 [76]. Such similarities lay the foundation of useful tools like bosonization [77,78].
However, fermionization and the 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, |Ψ B , is still symmetric, while the wavefunction of a non-interacting fermionic state, |Ψ F , is given by an anti-symmetric Slater determinant built from the first N single-particle eigenstates. While in the position space representation |Ψ B | 2 ≡ |Ψ F | 2 is true in the fermionization limit [76], the wavefunctions of bosons and fermions are completely different |Ψ B = |Ψ F ; this difference is reflected in distinct features of the momentum distributions of bosons in the fermionization limit and non-interacting fermions [see Fig. 2(b),(e),(h),(k),(n) and discussion below].
The density distribution at g = 20 in Fig. 2(g) is similar to the one of non-interacting fermions in Fig. 2(j). The one-body correlation function of the fermions has a complicated pattern of significant and non-trivial correlations [ Fig. 6(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 Ref. [73,79] 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 [58].  Table S1 of the supplementary material [58].
The difference between the bosons en route to fermionization and the fermions, however, lies most essentially in the momentum space. The momentum distribution for the bosons en route to fermionization [ Fig. 2(h)] has the same single-peak structure as the normal Mott insulating bosons [ Fig. 6(e)], where both widths and heights of the peaks are similar. This confirms the similarity between fermionized bosons and Bose-Einstein condensate in the momentum space [76]. In contrast, the fermionic state has three peaks in its momentum distribution due to the Pauli principle [ Fig. 2(k)]. These peaks correspond to the three fermions in each well since the two wells are Mott insulating [ Fig. 2(l)].
For long-range dipolar interactions crystallization emerges; for this case, bosons and fermions have been compared, for instance, in Ref. [80]. 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. Nevertheless, crystallized bosons and fermions can be distinguished via the spread of their density matrices as a function of the interparticle interaction strength [79].
A fermionic state crystallizes in the presence of strong long-range interactions g = 5 [cf. Eq. (12b)]. As expected, the repulsive interaction increases the distance between the fermions in real space, Fig. 2(m). The interactions also have a pronounced impact in the momentum space distribution Fig. 2(n) and the particle correlations Fig. 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 [cf. Fig. 2(j),(l)], become Mott insulating for crystallized fermions [cf. Fig. 2(m),(o)].
The correlations and fluctuations of particles can also be revealed by (simulations of) single-shot images [39][40][41][42][43]. For an illustration of what can be extracted from single-shot simulations, 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 single-shot simulations (adapting analysis.inp and running MCTDHX analysis). For each set of single shot images, we calculate the frequency of n particles being in the left well P (n) and show it in Fig. 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) ∼ 1 for all x, x in Fig. 2(c)]. In such a coherent state, the distribution of all bosons are independent from each other. As a result, P (n) should be given by the binomial distribution B(N, 1/2) = N n /2 N , and this is confirmed by the single-shot results. As the superfluidity gives way to Mott insulation, the distribution 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 = 3) gradually decreases. In the Mott-insulator limit, g = 1, E dw = 20, the repulsion between particles become so strong that there are exactly three particles in each well P (3) = 1.
The behavior of the system in the superfluid, Mott-insulating, and crystal phases can be summarized and represented by the correlation order parameter (COP) [42,43].
The COP is the sum of the squares of the eigenvalues ρ i of the reduced one-body density matrix [cf. Eq. (8)] The dependence of ∆ on the barrier height E dw and interaction strength g in bosonic and fermionic systems is shown in Fig. 4. In the bosonic system, ∆ decreases from almost unity in the superfluid phase to 0.5 in the Mott insulating phase, as expected for a fragmented state with ρ 1 = ρ 2 ≈ 0.5 [ Fig. 4(a) and (b)]. As the interaction strength increases and the system is en route to fermionization, ∆ 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 [ Fig. 4(c)]. The value ∆ ≈ 0.17 indicates one fermion occupies a single orbital, agreeing with the Pauli principle.

B. Dynamical behavior
Apart from solving the TISE for the ground state, MCTDH-X is also capable of solving the TDSE to capture the dynamics of a system as a reaction to a time-dependent Hamiltonian or the a quench of a parameter. As an example, we prepare the system in the ground state of a harmonic trap, V har (x) = 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 Eq. (13). We use a linear ramp with a time scale τ , 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   Table S2 of the supplementary material [58]. (g,h) Correlation order parameter ∆ [Eq. (14)] of (g) N = 6 bosons and (h) N = 6 fermions as a function of time with different ramping times τ . The dynamical behavior is completely different between τ < 10 and τ > 10 for bosons and fermions alike. A dynamical phase transition is implied. The input files of the simulations are given as simulations #7 and #8 in Table S2 of the supplementary material [58].
without three-body interactions, when the total number of bosons or fermions is odd (N = 5 in our simulations). In 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) | ≈ 0 in Fig. 5(a,c)], there is a large residual correlation |g (1) (x, −x)| ≈ 0.6 in the time-evolved state -even in the slow evolution limit with ramping time τ = 100 [ Fig. 5(b,d)]. 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.
To justify the claim of a first-order transition, we investigate the hysteretical behaviors of one-body correlations |g (1) (1, −1)| for N = 5 bosons, as shown in Fig. 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 with 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 Fig. 5(e)], as the residual correlations will never vanish in the Mott limit. The jump in the orbital occupation and the hysteresis are not seen with an even number of particles [ Fig. 5(f)]. In this even-number case, it is a second-order phase transition.
Since the system follows the ground state in the slow limit with an even number of particles, we investigate the effect of a quench 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 investigate the effects of a fast, an intermediate, and a slow ramp. Their one-body correlation functions at a given time t = 80 are shown in Fig. 6(a-c) for bosons, and in Fig. 6(g-i) for fermions.
With a slow ramp (τ = 50), the ground state of the final Hamiltonian is recovered -the correlation between the two wells vanishes [cf. areas where x is to the left (right) of the barrier and x to its right (left) and |g (1) Fig. 6(c,i)]. As the ramping time becomes shorter τ = 10, fluctuations appear on the correlation function, indicating excitations [cf. areas where x is to the left (right) of the barrier and x to its right (left) and |g (1) (x, x )| > 0 in Fig. 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 [ Fig. 6(a,g)].
To analyze the quantum fluctuations in the time-evolved states are affected by the speed of the ramp, we quantify the density fluctuations using the position-dependent variance V(x) [Eq. (11)] extracted from 10, 000 single-shot simulations for bosons in Fig. 6(d),(e),(f) and for fermions in Fig. 6(j),(k),(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 Fig. 6(d),(e),(f)]. For fermions, the magnitude of fluctuations is practically unaffected by the pace of the ramp [compare Fig. 6(j),(k),(l)]. We infer that the behavior of the magnitude of the fluctuations is analogous to the correlation order parameter (COP): for bosons, the COP increases as the pace of the ramp decreases [cf. Fig. 5(g) for large times t], whereas the COP varies only very slightly for the case of fermions [cf. Fig. 5(h)]. This analogous behavior of the quantum fluctuations and the correlation order parameter has previously been used to extract the phases of dipolar ultracold atoms in lattices [42,43].
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 ∆ [Eq. (14)], and its time-dependence [ Fig. 5(g,h)]. For bosons, in the slow case, only two orbitals are macroscopically occupied as in the ground state, while in the fast ramping case, the contribution of the first 10 orbitals are in 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 [cf. Fig. 4(a)]. As the ramping becomes faster, the COP starts to oscillate and the amplitude becomes larger and larger. However, the dynamical behavior of ∆ 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 does the COP. However, the dynamical behaviors of the COP has completely different behavior in the fast and slow ramping limits. In contrast to the bosonic case, ∆ increases slightly for a slow ramp, moving closer to the value of non-interacting fermions ∆ = 0.167.
The results from the one-body correlation, COP, and single-shot simulations point towards a dynamical phase transition at roughly τ = 10 for fermions and bosons alike. The system behaves like the ground state for slower ramps, but becomes highly excited for faster ramps.

C. Conclusion and Discussions
In this work, we have described and demonstrated the workflow, usage, and structure of the MCTDH-X software hosted, documented, and distributed through http://ultracold.org. In a step-by-step tutorial, we show how to examine particle correlations and fluctuations using the MCTDH-X software. We exhibited applications of MCTDH-X to solutions of the time-dependent and time-independent many-body Schrödinger equation for systems of bosons and fermions in a double-well potential. We have computed the many-body wavefunction and calculated various observables, including correlations, real and momentum space density distributions, and single-shot experiments.
The double-well potential captures many salient features of the Hubbard model, like the superfluid-Mott insulator transition of bosons. Additionally, we find a crystal state emerging for fermions that bears some similarities with strongly interacting bosons that are en route to the fermionization limit.
In the ground state, we thus have already observed a wide range of states of fermions and bosons by tuning the interparticle interactions and the barrier height of the double well. Apart from the transition into the Mott phase, the bosonic system approaches the non-interacting fermionic one as the interparticle interactions increases towards the fermionization limit. Fermionized bosonic states have many similarities to non-interacting fermions in the real space representation, and similarities to a usual bosonic gas in the momentum space representation. We demonstrate that the superfluid, Mott insulating, fermionized, and crystallized phases of many-body states can be distinguished by their correlation functions and by a correlation order parameter that is a function of the eigenvalues of the reduced one-body density matrix. Dynamics of a system may drive it out of the ground state, and induce quantum fluctuations and correlations. We use MCTDH-X to solve the time-dependent Schrödinger equation for bosons and fermions in a one-body potential that is smoothly transformed from a single into a double well by ramping up a Gaussian barrier at different time scales. We find an unexpected hysteretic behavior that heralds a first-order phase transition for the case of an odd number of bosons or fermions. When the particle number is even, there is no hysteresis, heralding a second-order phase transition. We thus demonstrate that the order of the superfluid-Mott insulator transition depends on the number of particles in our double-well system. With an odd number of particles, strong residual correlations between  Table S3 of the  supplementary material the two wells prevail even in the limit of large barriers for any pace of the ramp of the barrier. This implies that the superfluidity of the system cannot be eliminated by increasing the barrier, and the system cannot enter the Mott insulator phase. For faster ramps, our tutorial example approaches a quench scenario: When the double well's barrier is quenched up, a considerable amount of eigenstates of the final Hamiltonian become relevantly occupied and participate in the emergent quantum dynamics. This suggests strong fluctuations, which can be confirmed by an analysis of the variance in single-shot images. The significant differences between the fast and slow ramps that we observe point to a dynamical phase transition.
Finally, we note that although we only show examples with a few particles in this tutorial, the MCTDH-X is able to provide many-body states of systems with a much larger number of particles [82,83].
In this supplemental material, we give details on how to obtain the results and visualizations presented in the main text in Sec. S1. A discussion on the convergence of orbital number is given in Sec. S2, where the physics of the Bose-Hubbard model is also discussed in detail.

S1. PREPARING THE INPUT FILE FOR THE MCTDH-X SIMULATIONS IN THE MAIN TEXT
In this section we demonstrate how to prepare the input file and mention some caveats while using the program. In Tables S1, S2, S3, S4 and S5, we show how to modify the default input files included in the MCTDH-X software in order to repeat the simulations that are presented in the main text.
For the relaxations, the program (command MCTDHX ) can be directly executed when the modified input file and the library file libmctdhx.so that contains the Hamiltonian are present in the current working directory. For relaxations, we use randomly generated initial states to circumvent a convergence to excited states. When "repetition" is given as a keyword in a table, the program should be executed multiple (∼10) times with the same input file, and the energies of the resulting states should be compared in order to choose the correct ground state.
In our simulations of time-evolved states, the initial states are all obtained from relaxations. The results of the relaxations are saved in the binary files PSI bin, CIc bin and Header. They should be copied into the folder where propagations are performed. With the input file, binary files, and library in one folder, we can execute the program with the command MCTDHX. We should always verify that the propagations and the corresponding relaxations have compatible input variables like orbital number, particle number, grid size, etc.. We collect all the input files and scripts necessary to reproduce the results in the main text. Its folder tree is shown in Fig. S1, and it is available for download at http://ultracold.org/data/tutorial_input_files.zip.
There are several input variables which need more detailed explanation.
• The potentials HO1D, HO1D+td gauss, and HO1D+td g r are defined respectively as where x is the position, t is the imaginary or real time, and p 1 to p 6 correspond to the the input variables parameter1 to parameter6, respectively.
• The interparticle interaction can be chosen in the following way. The contact interaction can be chosen by setting Interaction Type = 0. In this case, the variables which interaction, Interaction Parameter1, and Interaction Parameter2 are irrelevant. The regularized Coulomb interaction can be chosen by setting which interaction = 'regC', Interaction Type = 4, and Interaction Parameter1 and Interaction Parameter2 are respectively α and β in Eq.(10b) in the main text. Other type of interaction can be chosen by using different keywords for which interaction. Detailed descriptions of these interactions are given in the manual at http://ultracold.org/data/manual.
• The number of grids NDVR X is recommended to be an exponent of 2.

S2. DISCUSSION ON THE CONVERGENCE IN THE ORBITAL NUMBER
The convergence in the orbital number is an important issue in MCTDH-X. The choice of the orbital number is a trade-off between computation speed and simulation accuracy. Without enough orbitals, we will not only obtain incorrect particle correlations, but sometimes even incorrect density distribution. The simplest method to determine the convergence is by comparing the occupations. A higher order orbital is considered insignificant, when its occupation is significantly lower than the lower order orbitals. This criterion sometimes fails when the seemingly insignificant orbitals contain important information about particle correlations and fluctuations of the simulated many-body state.   A most intuitive example is the strongly-interacting bosons. Fig. S2(a) shows the real-space density distribution of N = 2 bosons with extremely strong interaction g = 500 simulated with different orbital number. When M = 1 orbital is used, the interaction can only be captured in the Gross-Pitaevskii mean-field manner, and the many-body state simply reduces to a Thomas-Fermi cloud. As soon as the second orbital is used, the double peak density distribution can be captured, though an obvious discrepancy with the fermionic distribution can be observed. Such discrepancy diminishes as more orbitals are used.
A more detailed discussion on the convergence of orbital number can be illustrated by the comparison between the standard single-band Bose-Hubbard model and the simulation results for bosons in a Mott insulator state.

A. Bose-Hubbard model and beyond
The Hubbard model usually serves as good approximation for a lattice system. It assumes only one Wannier function contributes in each lattice sites, and its Hamiltonian takes the following form in one dimension, where,ĉ i is the annihilation operator for a particle in the Wannier state w i (x) that is localized at site i of the lattice andn i =ĉ † iĉ i is the number operator. In the Hubbard model, the superfluid to Mott-insulator transition is a result of the competition between the hopping strength t ∝ dxw * i (x)∂ 2 x w i+1 (x + a) between lattice sites with lattice constant a, and the on-site interaction U ∝ dx|w i (x)| 4 between particles. We consider a system with a fixed total number of particles N and therefore neglect the chemical potential.
In the limit t/U → ∞, the system is in the superfluid phase, and the ground state is given by In this state, the U (1) symmetryĉ i → e iθĉ i of the system is broken spontaneously, such that all the particles have the same phase. Particles in different sites are coordinated; they are condensed in the same single-particle state. This condensation entails a Glauber one-body correlation between sites which is unity. In the limit U/t → ∞, the system in a Mott insulator phase where the ground state of the system is given by |Ψ ∝ i a † i ni |0 , where n i is the number of particles in the ith site and satisfies the constraint i n i = N . In this state, particles in different sites are disconnected and have no information on each other's behaviors. The phases of the particles in different sites lose coherence. Strong correlation effects build up and the one-body correlation between sites vanishes. The transition point given in units of t/U depends on the structure of the lattice. In particular, it is roughly proportional to the filling factor and the inverse of the dimensionality.
For illustration, we consider N = 6 weakly-interacting (g = 1) bosons in a double-well potential with high barrier E dw = 20, where the two wells are Mott insulating from each other. According to the single-band Bose-Hubbard model, the energy of this system is minimzed when the bosons are evenly distributed into the two wells, with three bosons in each well. The many-body state can thus be written as |Ψ ∝ ĉ † L 3 ĉ † R 3 |0 , whereĉ † L andĉ † R are the creation operators for the band in the left and right well, respectively. Although such picture is well supported by the momentum space density distribution and the one-body correlation function [g (1) ][cf. Fig. 2(e,f) in the main text], it contradicts with the two-body correlation function [g (2) ] obtained from the simulations. In a Mott insulator state where g (1) between lattice sites vanishes, the diagonal term of g (2) is directly related to the number of particles n i in the corresponding site through g (2) (x, x) = 1 − 1/n i . If the single-band Bose-Hubbard model holds true, the diagonal term of g (2) would be uniformly g (2) (x, x) = 2/3 since n i = 3. This is not the case in the simulations when sufficient orbitals are used. The landscape of the correlation function is uneven, with value ranging from 0.4 to 0.8 [ Fig. S2(a)]. On the contrary, if only M = 2 orbitals are used in the simulations, indeed g (2) (x, x) = 2/3 holds uniformly [ Fig. S2(b)]. These results show that the higher order (M > 3) orbitals are indispensable for correctly capturing the particle correlations, and indicate that more than one mode has contributions in each well. The singleband Bose-Hubbard model is an over-simplification for the system. We have seen that although the third orbital ρ 3 = 4.432 × 10 −3 is two order of magnitude lower than the second orbital ρ 2 = 0.493, it has significant impact on the simulated particle correlations.
In summary, to decide whether the convergence in orbital number has been achieved, we should perform simulations with more orbitals and compare the convergence in different observables.