1 Introduction

In recent years, significant experimental and theoretical advances have been made in the field of quantum manybody systems, driven by the quest to understand and manipulate strongly correlated states. These states play a crucial role in various physical phenomena and hold promise for applications such as quantum computingand quantum materials. However, the inherent complexity of strongly correlated systems often presents substantial challenges, particularly in their large Hilbert space which made exact calculation of the intractable. In addition to the the recent development of numerical algorithms, such as tensor network or DMRG, which posed assumptions to the many body system, we still wish to apply exact treatment to many body systems for the purpose of quantum simulation of strongly correlated systems in table top experiments, such as trapped ions, neutral atoms or superconducting qubits. The project is organized to: 1. study the basic techniques of exact diagonalization of quantum many body systems and 2. learn how this can be applied to models frequently used in quantum simulation, such as Bose Hubbard model.

The first part of this paper focuses on the essential prerequisites for exact diagonalization. We begin with a detailed examination of the Schrieffer-Wolff (SW) transformation, using a lambda system as an illustrative example. This transformation is pivotal in reducing the complexity of the Hamiltonian by decoupling high-energy and lowenergy states. Additionally, we explore the concepts of symmetries and symmetry operators, demonstrating their applications with the hydrogen atom. By exploiting the inherent symmetries of a system, we can simplify

The second part of the paper explores specific case studies of quantum many-body systems. We begin by examining a onedimensional spin model, with a focus on the degeneracy of the ground state energy and the process of block diagonalization. This model is crucial for understanding how symmetry considerations can expedite the exact diagonalization process.

Following this, we analyze the Bose-Hubbard model, which is key to describing the behavior of interacting bosons on a lattice and offers insights into phenomena such as the superfluid-Mott insulator transition in ultracold atomic systems. Specifically, we investigate the phase diagram of the 1D Bose-Hubbard model using both exact diagonalization and mean field theory.

2 Prerequisite to the study of quantum many body system

In the study of complex multilevel systems, such as atoms or molecules, leveraging symmetries to reduce the Hilbert space dimension is a fundamental step. Beyond this initial reduction, it is often possible to further isolate a subset of states that are only weakly and non-resonantly coupled to the rest of the system. The dynamics of this subset can then be approximated by a Hamiltonian of reduced dimensionality. In this reduced Hamiltonian, the influence of interactions with states outside the relevant subspace is incorporated through additive energy shifts and modified couplings.

2.1 Schrieffer Wolff transformation and Effective Hamiltonian

Schrieffer Wolff transformation is a basic method to get the low energy dynamics from the full Hamiltonian.This method is broadly used in many body systems. t-J model can be derived from Hubbard model through Schrieffer Wolff transformation, explaining the superconductivity mechanism in high temperature superconductors of copper oxides.In the Fermi-Hubbard model, virtual doublon-hole excitations play a crucial role in mediating antiferromagnetic spin-exchange interactions within the effective model, often referred to as the t-J-3s model..Formally this procedure is described by performing a unitary Schrieffer-Wolff basis transformation.It is the broad use of SW transformation method in strong correlated system.The Schrieffer-Wolff transformation is a perturbative method used to simplify the Hamiltonian of a quantum system by eliminating certain off-resonant couplings, thereby focusing on the essential low-energy dynamics. This technique involves a unitary transformation that decouples the high-energy and low-energy subspaces of the system. The transformed Hamiltonian, known as the effective Hamiltonian, captures the low-energy behavior while disregarding the high-energy excitations.Suppose that the total Hamiltonian contains many-body interactions,perturbation theory allows one to construct a simpler high energy simulator Hamiltonian with only two-body interactions whose low-energy properties (such as the ground state energy) approximate the total Hamiltonian.Therefore,in many papers, SW transformation is used with perturbation theory. Mathematically, the SW transformation is represented as:

where H is the original Hamiltonian, and S is an anti-Hermitian operator chosen such that the transformed H ef f eliminates or reduces the coupling between the low-energy and high-energy states.

2.1.1 Lambda system

Over the past few decades, advancements in atomic, molecular, and optical (AMO) technology have significantly enhanced our ability to study and simulate strongly correlated quantum systems. One of the most common and fundamental models in AMO physics is the lambda system, which consists of three energy levels where two lower levels are connected to an excited state via optical transitions.The scheme is shown in the figure 1 ,which is ubiquitous in atomic physics.

δ k = ω kω (L) k (k = a, b), where hω k = E e -E k and ω (L) a,b are the frequencies of the lasers.In this system, |a⟩ -> |b⟩ is the forbidden transition. We need to use one laser to excite electrons from |a⟩ to |e⟩ and another laser to deexicte electrons from |e⟩ to |b⟩.In order to make more electrons to excite to |b⟩ instead of staying at |e⟩, we can excite |a⟩ into a virtual state (dotted line in the figure 1) and then deexcite electrons to |b⟩. The energy difference between virtual state and |e⟩ is small.We can write the Hamiltonian H = H 0 + V ,where H 0 is the unperturbed Hamiltonian of this system which is independent of time and V is the interaction potential between lasers and atoms. If we choose the middle of |a⟩ and |b⟩ as zero potential. Then H 0 = h(∆ω a ) |a⟩ ⟨a| + h(∆ω b ) |b⟩ ⟨b| + h∆ |e⟩ ⟨e|.In the electromagnetic wave, electric field is much larger than magnetic field.So that we only consider the interaction between electric field and atoms. We see atoms as dipoles ⃗ d which have the energy -⃗ d • ⃗ E in the electric field. We write the electric field as ⃗ E(⃗ r, t) = ⃗ E a cos( ⃗ k a • ⃗ rω (L) a t) + ⃗ E b cos( ⃗ k b • ⃗ rω (L) b t).Because the wavelength of electric field is much larger than the interatomic distance, it is a uniform electric field in the atomic scale. Then the electric field becomes: ⃗ E(t) = ⃗ E a cos(ω (L) a t) + ⃗ E b cos(ω (L)

where E + = 1 2 ( ⃗ E a e -iω (L) a t + ⃗ E b e -iω (L) b t ) and E -= 1 2 ( ⃗ E a e iω (L) a t + ⃗ E b e iω (L) b t ). ⃗ E a only functions on |e⟩ and |a⟩ and ⃗ E b only functions on |e⟩ and |b⟩. The electric dipole can be expressed as:

where |a⟩ , |b⟩ , |c⟩spans a space and I is the Identity matrix. Due to the transition prohibition between |a⟩ and |b⟩ and ⟨i| ⃗ d |i⟩ = 0, where i = a, b, e(the integral of odd function in the symmetric space is 0).The dipole becomes:

where σ 1 = |a⟩ ⟨e| , σ 2 = |b⟩ ⟨e| , ⃗ d ae = ⟨a| ⃗ d |e⟩ , ⃗ d be = ⟨b| ⃗ d |e⟩.The interaction term can be expressed as:

Remembering that E a only acts on |a⟩ and E b only acts on |b⟩,we define Ω i = -⟨e| ⃗ d• ⃗ Ei|i⟩ h

, i = a, b.Then we put V in the interaction picture.

V I = e i h H0t V e -i h H0t (6) After Schrieffer Wolff(SW) transformation and Rotating Wave Approximation(RWA),the interaction term becomes: V I = ( hΩ * a 2 e it(ωa-ω (L) a ) σ † 1 + hΩ * b 2 e it(ω b -ω (L) b ) σ † 2 ) + ( hΩ a 2 e -it(ωa-ω (L) a ) σ 1 + hΩ b 2 e -it(ω b -ω (L) b

We change the interaction term into the Schrodinger picture after throwing away fast rotating term. V = h   0 0 Ω * a 2 e iω (L) a t 0 0 Ω * b 2 e iω (L) b t Ωa 2 e -iω (L) a t Ω b 2 e -iω (L)

Then the total Hamiltonian becomes:

Ĥ = Ĥ0 + V = h   ∆ω a 0 Ω * a 2 e iω (L) a t 0 ∆ω b Ω * b 2 e iω (L) b t Ωa 2 e -iω (L) a t Ω b 2 e -iω (L)

This Hamiltonian is hard to solve and we can do a rotation transformation to apply this Hamiltonian into the rotating coordinate system. Then the final Hamiltonian becomes: H = h   ∆ω a + ω (L) a 0 Ω * a 2 0 ∆ω b + ω (L) b Ω * b 2 Ωa 2 Ω b 2

Although we obtain the low-energy effective Hamiltonian from the Schrieffer-Wolff (SW) transformation, solving this Hamiltonian remains challenging. To simplify the problem, we can exploit symmetries within the Hamiltonian. If we have a symmetry operator PR under which the Hamiltonian is invariant, this operator commutes with the Hamiltonian.Moreover,if PR is a Hermitian operator, then Ĥ and PR can be block diagonalized. It is clear that these operators form a group. If an operator PR belongs to this group, its inverse PR -1 must also be in the group. Additionally, the product of two operators within the group remains an operator of the group, as they can be considered as acting independently on the Hamiltonian. The associative property is evidently satisfied, which ensures the group structure. Regardless of whether the operators PR correspond to rotations, reflections, translations, or permutations, these symmetry operations do not affect the eigenvalues of the Hamiltonian. If the action of PR on an arbitrary vector composed of l eigenfunctions yields an l × l matrix representation of PR in block diagonal form, the matrix is block-diagonalized into smaller blocks,just like showing in the figure 2. We can then isolate each block and calculate its eigenvalues, significantly simplifying the overall calculation. For instance, using a symmetry operator, we can block-diagonalize a 7 × 7 matrix into the direct sum of a 3 × 3 matrix and a 4 × 4 matrix. Consequently, instead of calculating the eigenvalues of a 7 × 7 matrix, we only need to compute the eigenvalues of the smaller 3 × 3 and 4 × 4 matrices.

Figure 2: The hamiltonian is block diagonalized under symmetry operators 2.2.1. SO(4) in hydrogen atom

The hydrogen atom is a prime example of how symmetry operators can be used to analyze a system. The Hamiltonian of the hydrogen atom is invariant under rotational transformations, making it natural to introduce the SO(3) symmetry group into this system. The dimension of the irreducible representation of SO(3) is 2l + 1, which corresponds to the degeneracy of energy levels. However, according to the Schrödinger equation, when the principal quantum number n is determined, l can take values from 0 to n -1. The actual degeneracy of energy levels is n-1 l=0 (2l + 1) = n 2 , which is greater than 2l + 1. In other words, the SO(3) group does not fully describe the degeneracy (or symmetry) of the hydrogen atom system.

To resolve this, we can show that the hydrogen atom actually possesses SO(4) symmetry. According to central field theory, this system has a Runge-Lenz vector defined as

We can prove that [H, M] = 0, indicating that M is conserved. Here, L is the generator of SO(3), and M acts as another generator of SO(3). However, unlike L, M does not satisfy the standard commutation relations of angular momentum.It satisfies:

We now define I = L+N 2 , K = L-N 2 .

Let the eigenvalues of I 2 and K 2 be i(i + 1)h 2 and k(k + 1)h 2 , respectively. From equation (13), we obtain i = k and E = -me 4 2h 2 1 (2k+1) 2 . This result indicates that only the representation of SO(4) satisfying I = K correctly describes the bound energy levels of the hydrogen atom.

Furthermore, we can choose a matrix representation of the SO(4) group. Since M is a symmetry operator of this system, we can apply the numerical method illustrated in Figure 2 to calculate the energy levels of the hydrogen atom in matrix form.

3 Case study of quantum many body systems

Many one-dimensional (1D) spin models are exactly solvable, making them highly valuable for theoretical studies. For example, the 1D Ising model, the Heisenberg model, and the XY model all have known solutions. The two-dimensional Ising model was exactly solved by Lars Onsager in 1944. However, solving models in higher dimensions is generally much more challenging.

One-dimensional spin models, on the other hand, can be analyzed in detail using various methods, such as the Bethe Ansatz, matrix product states, and quantum field theory. Therefore, 1D spin models are an excellent choice for studying strongly correlated systems. In this paper, we focus on two typical 1D models. As discussed in the previous section, we can use the Schrieffer-Wolff (SW) transformation and symmetries to simplify our calculations. These methods allow us to perform numerical calculations for certain Hamiltonians.

3.1 Exact Diagonalization of 1D spin model

In 1D spin model,Heisenberg model is one of the most famous one. First,we consider a six-particle Heisenberg chain which can be written as follows:

This Hamiltonian satisfies the periodic boundary condition and λ 1 , λ 2 are two variables. We can rewrite this Hamiltonian into the form of spin-up operatorS

Naturally, if a particle at a site has spin-up, we label it as 1; otherwise, we label it as 0. For example, if the spins of six particles are up, up, down, down, up, down, we represent this configuration as |110010⟩. We then use these basis states to construct the matrix form of the Hamiltonian. Note the following actions of the spin operators:

..1...⟩, and S iz |...0...⟩ = -1 2 |...0...⟩. We can use a symmetry operator to block diagonalize the Hamiltonian, simplifying our calculations. Specifically, for this system, the Hamiltonian is invariant under rotation, meaning that the rotation operator is a symmetry operator and commutes with the Hamiltonian. When the symmetry operator is diagonalized in a new basis, the Hamiltonian will also be block diagonalized in that same basis. Therefore, we can use the eigenvectors of the symmetry operator to achieve a block diagonalization of the Hamiltonian.

3.1.1. Symmetry operators in 1D spin model

• Spin Flip Operator When the spin flip operator is applied to a basis state, it flips the spin of each particle: spin-up changes to spin-down, and spin-down changes to spin-up (where 0 represents the spin-down state and 1 represents the spin-up Specifically, the particle at site 1 is swapped with the particle at site 5, the particle at site 2 is swapped with the particle at site 4, and the particles at sites 3 and 6 remain unchanged. If we visualize these 6 particles arranged in a circle, the line from particle 3 to particle 6 acts as a reflection mirror. Due to the periodic boundary conditions, this operator is also a symmetry operator. For example, if the initial state is |110010⟩, applying the parity operator P results in the state |100110⟩. Applying P twice returns the system to the initial state. Therefore, the parity operator has two eigenvalues, 1 and -1.

• Rotation Operator The rotation operator can be constructed to cyclically exchange particles between adjacent sites. For instance, the particle at site 1 moves to site 2, the particle at site 2 moves to site 3, and so on. For example, if the initial state is |110010⟩, applying the rotation operator R results in the state |011001⟩. When this rotation operator R is applied six times, the system returns to the initial state, so R 6 = 1. The operator has six eigenvalues: 1, e i π 3 , e i 2π 3 , -1, e i 4π 3 , and e i 5π 3 . These six eigenvalues correspond to rotations by angles of 0, π 3 , 2π 3 , π, 4π 3 , and 5π 3 , respectively.

3.1.2 Block diagonalized Hamiltonian for 6 sites

After applying the spin flip operator, the Hamiltonian separates into two parts. Each part is represented by a 32x32 matrix, corresponding to the eigenvalues of the spin flip operator, -1 and 1, respectively. The detailed values of these matrices can be found in the Appendix. Next, we apply the parity operator to the Hamiltonian in this new basis. Each 32x32 block is further block diagonalized into two smaller parts: a 12x12 block and a 20x20 block, representing the eigenvalues of the parity operator, -1 and 1, respectively. We then proceed by applying the rotation operator. For the case of 6 sites, the eigenvalues of the rotation operator are 1, -1, -1 2 , and 1 2 after the spin flip and parity operators have been applied. Notably, the eigenvalues -1 2 and 1 2 exhibit a two-fold degeneracy. This degeneracy arises because rotations by angles of π 3 and 4π 3 , as well as 2π 3 and 5π 3 , are equivalent in the new basis, leading to this two-fold degeneracy.

The final block-diagonalized Hamiltonian takes the form in the figure 3:

Figure 3: This is 6 particle block diagonalized Hamiltonian after using spin flip operator,parity operator and rotation operator.Each block in this Hamiltonian has physical meaning.The lowest left block means this block has spin flip symmetry(S=1),parity symmetry(P=1),rotation symmetry(R=1).

-2 -1 0 1 2 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 λ 1 Energy λ 2 = -2 λ 2 = -1

Figure 4: Theoretical ground energy for 1D spin model

In the figure 4, we observe that when we select λ 2 and its inverse, their graphs are symmetric about the y-axis. This symmetry is expected. Since we can arbitrarily choose the spin operators S y and S z , λ 1 and λ 2 can be interchanged in the context of the Hamiltonian. Consequently, when we fix λ 2 and plot the energy as a function of λ 1 , we obtain the same result.

3.1.3 Properties of ground energy

To find the ground energy properties, 6 sites are not enough. We gradually increase site numbers just as shown in the figure 5.

10 12 14 16 18 20 L -13 -12 -11 -10 -9 -8 -7 -6 -5 Energy λ1 = 1,λ2 = 2 Figure 5: Ground energy for different sizesAs size goes bigger, which means closer to actual situation, we find that the ground energy becomes lower and lower because of more spins in the system.

-4 -3 -2 -1 0 1 2 3 4 λ 1 -20.0 -17.5 -15.0 -12.5 -10.0 -7.5 -5.0 -2.5 Energy L=5 L=10 L=16 L=21

Figure 6: Ground energy for different λ 1 when λ 2 = 1

In the figure 6, as λ 2 increases and λ 1 becomes larger, the ground energy decreases continuously. This implies that when the spin preferentially aligns in a particular direction, it becomes easier for the spin to rotate into that direction, leading to a lower energy state. Moreover, as λ 2 increases, the ground energy also decreases. This can be demonstrated by considering operators in other directions as perturbations. The resulting Hamiltonian is given by:

Therefore, as λ 2 increases, the lowest eigenstate tends to approach |000000⟩ for 6 sites, which corresponds to the lowest energy of -3λ 2 . Conversely, when λ 2 is sufficiently small, the ground state becomes |111111⟩ for 6 sites. This indicates the existence of a state with a maximum ground energy. As shown in Fig. 4, the maximum ground energy occurs when λ 2 ≈ ±1. Additionally, it is evident that as the number of sites increases (bringing the system closer to realistic conditions), λ 2 is more likely to be around ±1.

3.1.4 Properties of Sz

Then we want to find the mean value of

where ψ is the ground state of the hamiltonian. First, we set λ 1 unchanged and plot S z versus λ 2 .

8 10 12 14 16 18 20 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 S z even L odd L -4 -3 -2 -1 0 1 2 3 4 λ 2 -8 -6 -4 -2 0 2 4 6 8 S z λ 1 = 1 L=14 L=15 -4 -3 -2 -1 0 1 2 3 4 λ 2 -6 -4 -2 0 2 4 6 S z λ 1 = 1 L=14 L=15 Figure 8: S z for different λ 2 when λ 1 has different valuesIn the figure 7, we find that for even numbers of sites, S z is zero. However, for odd numbers of sites, S z oscillates with λ 2 . Moreover, as λ 2 becomes greater than -1, the amplitude of this oscillation decreases. This indicates that for even numbers of sites, there is no net magnetic moment in the system. This generally implies that the spins are antiparallel, which is a more stable configuration. Due to symmetry breaking or local spin asymmetry, a non-zero expectation value of S z may appear in the ground state. This phenomenon reflects a local spin bias or asymmetry in the system. This behavior is observed when λ 1 ̸ = 1.

However, when λ 1 = 1, the situation changes. For λ 1 = 1, S z oscillates with λ 2 regardless of whether the number of sites is odd or even. For even numbers of sites and λ 2 > -1, S z remains zero. For odd numbers of sites and λ 2 > -1, S z oscillates with a smaller amplitude, just as shown in the figure 8.

3.1.5 Degeneracy of the ground energy

As discussed above, -1 and 1 are two significant values. Additionally, parity also influences the system. Therefore, when determining the degeneracy of the ground state, it is important to take these factors into account carefully.

Table 1: Ground energy degeneracy for even sites

Table 2: Ground energy degeneracy for odd sites

the degeneracy increases to 2 N . In this situation, the Hamiltonian takes the following form:

It is obvious that this hamiltonian has larger degeneracy.It has the relation with the symmetry of the system.

3.1.6 Symmetry of the ground state

Degeneracy is closely related to symmetries. Therefore, it is essential to understand the symmetries of the states.

• S (Spin Flip Operator): In the S = 1 block, the basis states are spin flip symmetric. In the S = -1 block, the basis states are spin flip anti-symmetric.

• P (Parity Operator): In the P = 1 block, the basis states are inversion symmetric. In the P = -1 block, the basis states are inversion anti-symmetric.

• R (Rotation Operator): In the R = 0 block, there is no rotation. In the R = 1 block, the chain is translated by one lattice site, and so on.

Table 3: Ground state symmetry for even sites

Table 4: Ground state symmetry for odd sites

states occur only when R=0 or R=N/2, corresponding to no rotation or a half rotation, respectively.Additionally, the degeneracy changes as λ 1 and λ 2 vary. For instance, when λ 1 < -1 and λ 2 = 1, the degeneracy is 2. In this case, two ground states can be interpreted as a combination of the states when λ < 1 and λ 1 > 1. This is analogous to the ground state shifting from |S = -1, P = 1, R = 0⟩ to |S = 1, P = 1, R = 0⟩ as λ 2 changes from less than 1 to more than 1. During this transition, the degeneracy is broken. A similar phenomenon occurs when λ 1 > 1. For odd numbers of sites, the minimum ground state degeneracy is 2. When the ground state degeneracy is 4, it corresponds to the states |S = ±1, P = ±1⟩.When the energy degeneracy is 2, the ground state must must be of the form |S = ±1, P = 1, R = 0⟩.

3.2 Numerical study of 1D Bose-Hubbard model

The earliest investigations into bosons' behavior in a lattice can be traced back to Philip Anderson's theory on localization phenomena in 1963.The critical work for the Bose-Hubbard model is indeed the 1998 paper by Jaksch et al., which provided a theoretical framework for realizing the Bose-Hubbard model with cold atoms in optical lattices.Immanuel Bloch and his colleagues at the University of Munich first experimentally realized the Bose-Hubbard model in optical lattices and observed the superfluid to Mott insulator transition. This achievement marked the experimental confirmation of the Bose-Hubbard model's predictions. The absence of disorder and the low temperatures at which these experiments are conducted make such systems ideal for studying quantum phase transitions. Theoretically, the simplest Hamiltonian that describes a system of bosonic atoms in an optical lattice is the Bose-Hubbard model. This model can be expressed as follows:

• Mott Insulator in the strong interaction regime, which means U >> t,we can set t = 0 and the 1D Hamiltonian reduces to

Different states are now decoupled.The number states are the eigenstates for Ĥ,which means in order to get the ground state,we need to find the smallest value of h i .Then we obtain

it means the number of particle per site is a specific integer number due to the ratio of U and µ.If we add a particle into one site,which means there is a particle excitation,there will be some change in energy.

This implies that there will be energy gaps in the system. Adding a particle to the system requires significant changes. This phase is referred to as a Mott insulator. Unlike insulators in band theory, where the insulating property arises from the band structure, a Mott insulator's insulating behavior is due to strong interactions between particles.

In real world,the total particle in a system L i n i = N b in the strong interaction(t=0). The Hamiltonian becomes

where L is the number of sites and Nb is the total number of particles.We are primarily concerned with the ground state, which corresponds to minimizing the Hamiltonian.Therefore, our goal is to minimize

By applying the Cauchy-Schwarz inequality, we have:

(a) N b ≤ L For the ground state, n i could only be 1 or 0. For N b sites, there are 1 particle in each site. For the else sites(L -N b), there is 0 particle.Then the ground energy becomes:

The average should be n i = [ N b L + 1 2 ]. However, sometimes n i * L ̸ = N b .Therefore,due to the total particle conservation,particle on each site may not perfectly equals to n i .

(ii) U < 0 If we want to maximize L i n 2 i ,we need to place Nb particles at one site and no particle at other sites.Then eq(23) becomes:

We want to plot the energy with different particles. The maximum ground energy occurs when N b = max([ 1 2 + µ U ], 0). In the figure 9 we choose U = -2 and when µ >= 1, the maximum energy occurs at Nb=0.When µ < 1, the maximum energy occurs at N b = [ 1 2 + µ U ].

-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 µ -10

(a) repultsive interaction(U=2) -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 µ -120 -100 -80 -60 -40

(b) attractive interaction(U=-2) Figure 9: Energy for different chemical potential.The dots are numerical resultsand lines are theoretical results.

Ĥ = -t

we can use Fourier transition bk = 1 √ L L j=1 bj e -ikj ,then the Hamiltonian becomes

Therefore, the ground state for this Hamiltonian is when all particles condense into the k = 0 state in the Fourier space

Then the ground energy becomes:

Due to the delocalization of particles, this state characterizes the superfluid phase. According to Landau's theory, a superfluid phase occurs if V c = min( ϵ|q| |q| ) ̸ = 0, where ϵ|q| represents the energy dispersion relation.

6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 N b -15 -10 -5 0 5 10 15 20 25 Energy µ = -2 µ = -1 µ = 0 µ = 1 µ = 2

(a) Mott Insulator 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 N b -60 -50 -40 -30 -20

(b) Superfluid Figure 10: Energy for different particles in different phaseWe start with the Hamiltonian of the Bose-Hubbard model and use exact diagonalization to solve it. In the following discussion, we apply periodic boundary conditions, meaning that the beginning and end of the chain are connected. We utilize the constrained Bose-Hubbard model, where the number of particles and sites are fixed,just like in the figure 10. In a realworld scenario, both the number of particles and the number of sites would approach infinity. However, due to computational limitations, we are constrained to calculations with a finite number of particles and sites.

3.2.1 Exact diagonalization

It is evident that we use the canonical ensemble to calculate the Hamiltonian. Consequently, we set the chemical potential µ to zero. By definition, the chemical potential represents the energy required to add or remove a particle from the system. There are various methods to calculate the chemical potential and assess whether the system is in the Mott phase or the superfluid phase. We can use definition µ, where E 0 is the ground state of this system.The gap ∆ is finite for the Mott insulator.We interactions. In contrast, superfluids have finite compressibility, which indicates their ability to adjust particle density and flow without resistance. Therefore, in exact diagonalization, the difference in compressibility can be used to determine whether the system is in a Mott insulator phase or a superfluid phase. 0.0 0.2 0.4 0.6 0.8 µ/U 0.2 0.4 0.6 0.8 1.0 1.2 ρ 12 sites t=0.05

Figure 11: The number density ρ as a function of chemical potential µ when we set t/U = 0.05

The slope of the figure 11 is compressibility κ.Clearly, a phase transition occurs whenµ ≈ 0.1. Based on the previous discussion, when the slope of the compressibility is zero, the system is in the Mott insulator phase. Conversely, when the slope is finite, the system is in the superfluid phase. (Due to the limitations of lattice size, the slope may not be perfectly smooth.) By varying the hopping parameter t, we can observe different phase transition points, which can then be plotted on a phase diagram,just like in the figure 12. 0.00 0.05 0.10 0.15 0.20 0.25 t/U 0.0 0.2 0.4 0.6 0.8 1.0 µ/U Figure 12: The phase transition when U=1. The part enclosed by two lines is Mott insulator phase and the outside part is superfluid phase 3.2.2. Mean-field theory

Ĥt = -t ⟨ij⟩ ⟨ b † i ⟩ bj + b † i ⟨ bj ⟩ -⟨ b † i ⟩⟨ bj ⟩ = 2tb 2 0 N s -2tb 0

In this equation,we suppose the transition invariant in the system and due to the phase we choose.Then the order parameter must be real constant so that we choose⟨ b † i ⟩ = ⟨ bj ⟩ = b 0 .N s means total number of sites.The total Hamiltonian can be expressed as:

The summation term is the mean field term, which includes the hopping term and onsite energy term.In order to calculate this Hamiltonian,we need to in the |n b ⟩ basis,like |n b ⟩ = 0, |n b ⟩ = 1, |n b ⟩ = 2... Due to the no upper limit of n b , we need to do truncationn b,max = n.Then in this basis, the matrix of bj =

. . . . . . . . . . . . . . .

. We now have an (n+1)×(n+1) matrix. As n increases, the results become more accurate. To build the Hamiltonian for each site, we need to include the order parameter b 0 . The procedure is as follows:

• 1.Initial Guess: Start with an initial guess for b 0 and incorporate it into the Hamiltonian.

• 2.Diagonalization: Diagonalize the Hamiltonian to obtain the ground state and calculate the average value ofb j .

• 3.Use the average value of b j as the new b 0 and substitute it back into the Hamiltonian.

• 4.Repeat this process several times until the difference between the new b 0 and the old b 0 is within an acceptable range.

The average particle density per site can then be used to determine whether the system is in the Mott insulator phase or the superfluid phase.Then we can plot the phase diagram as figure 13.

0.03 0.06 0.09 0.12 t/U 0.0 0.5 1.0 1.5 2.0 µ/U ρ=1 ρ=2 Figure 13: Mean field phase diagram

Moreover, we can also treat the hopping term as a perturbation, while the interaction energy term remains the nonperturbative part. Using second-order perturbation theory, the ground state energy can be expressed as:

Hamiltonian do not contribute when acting on the state |n b ⟩. Additionally, all odd-order perturbation terms are zero. Therefore, the Hamiltonian can be expressed as: E ′ (n b ) = E(n b ) + A(n b )b 2 0 + B(n b )b 4 0 (34) It is the same as Landau free energy.Due to the Landau theory,in the B(n b ) > 0 situation,the transition occurs at A(n b ) = 0.Then we get: t U = (n b -µ U )( µ U + 1n b ) 2( µ U + 1)

We can compare the theoretical phase diagram with numerical phase diagram in the figure 14. 0.03 0.06 0.09 0.12 t/U 0.0 0.5 1.0 1.5 2.0 µ/U Mott Insulator Mott Insulator Superfluid ρ=1 ρ=2

Figure 14: The smooth curve is the theoretical line and the zigzag line is the mean field exact diagonalization.

However, in 1D systems, quantum fluctuations are more pronounced compared to higher dimensions. Mean-field theory approximates the behavior of particles by averaging their interactions with surrounding particles. While this approach can be effective in higher dimensions, where fluctuations are relatively small, it fails in 1D systems where fluctuations play a critical role. Therefore, using mean-field theory to calculate the 1D Bose-Hubbard model is not appropriate.

3.2.3 Gutzwiller variational ansatz

To describe transitions between the superfluid and insulating states we use Gutzwiller variational wavefunction.

Here the coeffiecient is normalized.After minimizing the wavefunction with respect to the normalization constraint, we obtain the following equation:

To numerically calculate the phase transition of the Gutzwiller Hamiltonian, we use the simulated annealing (SA) algorithm. By employing normalized random functions to simulate the Gutzwiller wavefunction, we can construct the Gutzwiller Hamiltonian. We then generate another Gutzwiller Hamiltonian using a different normalized random function. The decision to accept the new 0.0300 0.0867 0.1433 0.2000 t/U 0.0 0.5 1.0 1.5 2.0 µ/U ρ=1 ρ=2 Figure 15: Gutzwiller phase diagram

We can compare phase diagram in the figure 16. 0.03 0.07 0.11 0.15 t/U 0.00 0.25 0.50 0.75 1.00 µ/U Gutzwiller Meanfield Exact Diagonalization

Figure 16: Red plot is Gutzwiller, blue plot is meanfield and green plot is Exact Diagonalization.The area enclosed by a curve is Mott Insulator and outside is Superfluid phase.

3.2.4 Correlation function

Moreover,we can use correlation function < b † 0 b i > to verify the phase transition.In the superfluid phase, the system exhibits long range coherence so that the attenuation behavior of b † 0 b i is usually polynomial.

where r = |i -0|,α is a positive constant.The long range polynomial attenuation reflects the existence of long range quantum coherence and coherent fluctuations in the superfluid phase.In the Mott insulator phase, the motion of particles in the system is † consider half of the chain.From the figure 17, we know that when t is very small(Mott insulator phase), correlation function decays exponentially. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 i 0.0 0.2 0.4 0.6 0.8 1.0 < b † 0 b i > t = 0.05 t = 1 t = 5

Figure 17: Correlation function < b † 0 b i > for different particle sites when we set different tWhen the number of sites t is large (indicating the superfluid phase), the correlation function decays polynomially. To determine the phase transition, we can use exponential fitting. If the fit residual is sufficiently small, we can conclude that the system is in the Mott insulator phase. Conversely, if the residual is large, the system is in the superfluid phase. To calculate the transition point from the Mott insulator phase to the superfluid phase, we require statistical analysis. First, we fit the correlation function with an exponential function. Then, we calculate the root-mean-square error (RMSE) of the fit, as illustrated in Fig. 17.

where a,b and c are fit parameters of the correlation function.

0.00 0.05 0.10 0.15 0.20 0.25 0.30 t 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 RMSE 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.2 0.4 0.6 0.8 1.0 Chem

superluid. The result is similar to the mean field theory results.It is a great verification of our phase transition.

4 Conclusion

In this paper, we have summarized a variety of tools, including the RWA approximation, symmetry operators, and effective Hamiltonians, which are highly useful in the study of strongly correlated systems for simplifying complex Hamiltonians. We aimed to capture the essential technical aspects of these approaches and provide examples of their application in various physical contexts. We gained deep insights into the 1D spin model by employing three symmetry operators-parity operator, spin flip operator, and rotation operator-to simplify the Hamiltonian. This allowed us to determine ground state properties such as degeneracies and symmetries. Additionally, we studied the 1D Bose-Hubbard model using three methods-exact diagonalization, mean-field theory, and Gutzwiller wave approximation. We constructed phase diagrams from these methods and used correlation functions to verify the results.

Appendix

Appendix A: Detailed derivation of lambda system The dipole-electric field interaction term in the interaction picture

The interaction term can be expressed as:

We put V in the interaction picture.

as an example.Here we give a very useful equation.

In the same way, we can derive that e

In the same way, we can derive V I as follows:

unitary transformation

If we have a unitary transformation U and |ψ⟩ is the solution to the original Schrodinger equation. New state becomes

It is easy to prove that the relation between equation(10). U =    e -iω (L) a t 0 0 0 e -iω (L) b t 0 0 0 1    (50) Rotating Wave Approximation According to time-dependent perturbation theory,we have below:

where u I (t, t 0 ) is the time-revolution operator in the interaction frame;V I (t) is the perturbation if the total Hamiltonian H = H 0 + V (t) and V (t) is small.For ûI (t, t 0 ) we have following equation:

tn-1 t0

The transition amplitude from state m to state n is written as follows:

where ψ 0 n and ψ 0 m are unperturbed states which means the eigenstate of H 0 .In our question, we need to know the amplitude from |a⟩ to |b⟩. We expand the equation (52) into second order.

Because we need to calculate the amplitude from a to b, the interaction term must include terms similar to |b⟩ ⟨a|. Otherwise, due to the orthogonal between |a⟩ , |b⟩ and |e⟩, other terms should be 0.Therefore, the first order perturbation becomes 0 because it doesn't include |b⟩ ⟨a| term.After selecting the term includes |b⟩ ⟨a|,the amplitude becomes:

after integrating, the transition probability including ω a + ω (L) a is much smaller than that including |ω aω (L) a |.The term including ω a + ω (L) a can be neglected.This is called Rotating Wave Approximation:the fast rotating term like e it ′′ (ωa+ω (L) a ) can be neglect.

Appendix B: Scalar operator, vector operator and tensor operator Electric field operator and dipole operator are vector operators.

scalar operator

The generator of rotating is angular momentum. If the expectation value of an operator is invariant under rotation, we call this operator the scalar operator.We have a state |ψ⟩ and after rotation, the state becomes ψ ′ = D(R) |ψ⟩. The eigenvalue of operator V can be expressed as follows: 56) is true for all ψ so that V = D(R) † V D(R). In specific case, D(R) = 1 -iϵJ•n h ,where n is the rotating axis;J is the angular momentum of this system;ϵ is the small rotating angle.After neglecting the higher term of ϵ,the equation (56) becomes:

If an operator V is a vector operator, its expectation value must satisfy V i = Σ j R ij V j ,where R i j is the rotating matrix.Then the equation (57) becomes:

For n in z-axis, we get the rotating matrix: R(n; ϵ) =   1 -ϵ 0 ϵ 1 0 0 0 1

We choose i = 1, 2, 3 and then we can get the commutation relation:

Therefore, If the operator V satisfy equation (60), we can say that V is a vector operator.

tensor operator

we can use two vector operator to structure a tensor operator,like T ij = U i V j , and we will get nine components.We can also write T ij as follows:

The first term is the scalar product invariant under rotation. The second term is an anti-symmetric tensor.The last term is a 3x3 traceless tensor with 5 components.The definition of tensor operator:

we take equation (62) into the equation (57) and we get:

Appendix C: SO(3) generator

For the rotation symmetry(SO(3) group),we have the relation:

where J i , J j , J k is the generator of SO(3).If we need SO(4), we will have 6 generators.We can separate six generators into two groups. One is J 1 , J 2 , J 3 which is the same as that in the SO(3) group and another is L 1 , L 2 , L 3 . They have the relation as follows:

We can mix those six generators as:J +,i = Ji+Li 2 ; J -,i = Ji-Li 2 ,then the relation becomes:

We now have two separate SO(3),which means SO(4) = SO(3) ⊗ SO(3).

Table 5: This table shows the nonzero value of 1D spin model Hamiltonian(32*32 matrix)for 5 sites in the original basis as an example.The subscript denotes rows and columns of Hamiltonian.

H 26,26 λ 2 /4 H 26,27 (1 + λ 1 )/4 H 26,32 (1λ 1 )/4 H 27,3 (1 -λ 1 )/4 H 27,12 (1 + λ 1 )/4 H 27,23 (1 + λ 1 )/4 H 27,26 (1 + λ 1 )/4 H 27,27 (-3*λ 2 )/4 H 27,29 (1 + λ 1 )/4 H 28,4 (1 -λ 1 )/4 H 28,11 (1 -λ 1 )/4 H 28,24 (1 + λ 1 )/4 H 28,25 (1 -λ 1 )/4 H 28,28 λ 2 /4 H 28,30 (1 + λ 1 )/4

Table 6: This table shows the Hamiltonian after using spin flip operator

H 20,20 λ 2 /4 H 20,22 (1 + λ 1 )/4 H 20,32 (1λ 1 )/4 H 21,6 (1 + λ 1 )/4 H 21,13 (1 + λ 1 )/4 H 21,19 (1 + λ 1 )/4 H 21,21 (-3*λ 2 )/4 H 21,24 (1 -λ 1 )/4 H 21,25 (1 + λ 1 )/4 H 22,5 (1 -λ 1 )/4 H 22,14 (1 + λ 1 )/4 H 22,20 (1 + λ 1 )/4 H 22,22 (-3*λ 2 )/4 H 22,23 (1 + λ 1 )/4 H 22,26 (1 + λ 1 )/4 H 23,8 (1 + λ 1 )/4 H 23,15 (1 + λ 1 )/4 H 23,17 (1 -λ 1 )/4 H 23,22 (1 + λ 1 )/4 H 23,23 (-3*λ 2 )/4 H 23,27 (1 + λ 1 )/4 H 24,7 (1 -λ 1 )/4 H 24,16 (1 + λ 1 )/4 H 24,18 (1 -λ 1 )/4 H 24,21 (1 -λ 1 )/4 H 24,24 λ 2 /4 H 24,28 (1 + λ 1 )/4

(1 -λ 1 )/4 H 27,12 (1 + λ 1 )/4 H 27,23 (1 + λ 1 )/4 H 27,26 (1 + λ 1 )/4 H 27,27 (-3*λ 2 )/4

(1 -λ 1 )/4 H 32,29 (1 -λ 1 )/4 H 32,32 (5*λ 2 )/4 In the table 5 and 6, we can clear see that after using spin flip operator, the Hamiltonian is block diagonalized because H 1∼16,17∼32 = 0 and H 17∼32,1∼16 = 0.

Appendix E: Ground energy of 1D spin model for 6 sites For 6 sites Heisenberg model,in the paper,we use symmetry operators to get the ground state.Therefore, we want to know the correctness of this method.We plot the ground 3D pictures in the figure19. Energy -4 -2 0 2 (a) theoretical ground energy -3 -2 -1 0 1 2 3 λ 1 -3 -2 -1 0 1 2 3 λ 2 -6 -5 -4 -3 -2 Energy -6 -5 -4 -3 -2 Energy (b) ground energy using

Figure 19: 6 sites Heisenberg model ground energy 3D figure

We also plot energy versus λ in the figure 20. -4 -3 -2 -1 0 1 2 3 4 λ 2 -6 -5 -4 -3 -2 Energy λ 1 = -1 λ 1 = 0 λ 1 = 1 Figure 20: The dot plot is theoretical for 6 sites and the curve is results for 6 sites using package. We can see that curves perfectly link points.

We can conclude that by using block diagonalization, we can still get the right answer.

Appendix F: Basis construction of 1D Bose-Hubbard model

We need to arrange N b bosonic particles into L sites. Since the particles are bosons, there is no restriction on the number of particles per site. The total number of possible arrangements is given by (N b +L-1)! (L-1)!N b ! ,which can be derived using combinatorial methods. For example, to arrange 2 particles in 3 sites, the possible arrangements are: (2,0,0) (1,1,0), (1,0,1),(0,2,0), (0,1,1), and (0,0,2). We use these orthonormal basis states to construct the Hamiltonian in the figure 21.

(a) basis without using package (b) basis using package

Figure 21: We use two different approaches to construct basis when L=3,N b =2

After we get the basis, we can construct the Hamiltonian to get the ground energy of this system.

Appendix G: Detailed derivation of Gutzwiller Hamiltonian

We take the Gutzwiller wavefunction withf n = (1 -2α 2 ) 1/2 , f n-1 = f n+1 = α.We can expand the hopping energy to the second order.

Moreover, the n is large enough so that we can use the approxiamtion.If we apply µ = U (n -1 2 ) into the chemical and potential term, the relation becomes when U = 8nt (69) Mathematically, there is no significant difference between the Gutzwiller variational ansatz and the mean-field approximation when the number of particles n is sufficiently large in three dimensions. However, in one the Gutzwiller variational ansatz provides a more concrete and detailed description compared to the mean-field approach.

Appendix H: Detailed derivation of Mean field Hamiltonian

From eq(33) and eq (34),we only need to calculate A(n b ).By using equation bj |n b ⟩ = √ n b |n b -1⟩ and b † j |n b ⟩ = √ n b + 1 |n b + 1⟩,we can easily get the second perturbation term: E ′ (n b ) = E(n b ) + (-2tb 0 √ n b + 1) 2 E(n b ) -E(n b + 1) + (-2tb 0 √ n b ) 2 E(n b ) -E(n b -1) = E(n b ) + (2tb 0 ) 2 n b + 1 µ -U n b + n b U (n b -1)µ (70)

After we apply the constant energy 2tb 2 0 ,the average energy per site is

As we discuss above, phase transition occurs when the coefficient is 0,which means 2t( n b U (n b -1)µ + n b + 1 µ -U n b ) = -1 (72) If we unitize t and µ:ω = 2t U and µ ′ = µ U ,we can get the function between ω and µ ′ when given n b .The equation(22) becomes:

Solving this quadratic equation,we get

Moreover, in order to plot pictures more easily, we can rewrite equation 23 like

The wave function becomes:

0.00 0.05 0.10 0.15 0.20 0.25 0.30 t 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 RMSE L = 5 L = 9 L = 12

Figure 22: residual vs t for different sizes

As the number of sites increases, the peak values of the transition point become nearly constant. in the figure 22. This suggests that the transition point approaches the value representative of the true physical system.