Layer Anti-Ferromagnetism on Bilayer Honeycomb Lattice

Bilayer honeycomb lattice, with inter-layer tunneling energy, has a parabolic dispersion relation, and the inter-layer hopping can cause the charge imbalance between two sublattices. Here, we investigate the metal-insulator and magnetic phase transitions on the strongly correlated bilayer honeycomb lattice by cellular dynamical mean-field theory combined with continuous time quantum Monte Carlo method. The procedures of magnetic spontaneous symmetry breaking on dimer and non-dimer sites are different, causing a novel phase transition between normal anti-ferromagnet and layer anti-ferromagnet. The whole phase diagrams about the magnetism, temperature, interaction and inter-layer hopping are obtained. Finally, we propose an experimental protocol to observe these phenomena in future optical lattice experiments.

Bilayer honeycomb lattice, with inter-layer tunneling energy, has a parabolic dispersion relation, and the inter-layer hopping can cause the charge imbalance between two sublattices. Here, we investigate the metal-insulator and magnetic phase transitions on the strongly correlated bilayer honeycomb lattice by cellular dynamical mean-field theory combined with continuous time quantum Monte Carlo method. The procedures of magnetic spontaneous symmetry breaking on dimer and non-dimer sites are different, causing a novel phase transition between normal anti-ferromagnet and layer anti-ferromagnet. The whole phase diagrams about the magnetism, temperature, interaction and inter-layer hopping are obtained. Finally, we propose an experimental protocol to observe these phenomena in future optical lattice experiments.
T he bilayer honeycomb lattice (BHL) has attracted enormous interest in both experimental and theoretical research. Lots of novel phenomena have been found in BHL, for instance, the quantum Hall effect, quantum spin Hall effect, and chiral superconductivity [1][2][3][4][5][6][7][8][9][10][11] . However, the real charge and magnetic order induced by Coulomb interaction are still challenges in the strongly correlated BHL [12][13][14][15][16] . A quadratic dispersion relation, signed by two touching bands in the corners of Brillouin zone, can be found in BHL, which is driven by the interlayer hopping. In addition, the spontaneous symmetry can be broken by the dimers when the inter-layer hopping changes. Some amazing phases may emerge, such as layer anti-ferromagnetic phase and paramagnetic insulator phase. Previous work mainly focus on the electronic properties of BHL [17][18][19][20][21][22][23][24] , and the phase diagram for the magnetic phase transition induced by the interaction and dimers are still absent. Moreover, the progress of optical lattice provides us a useful tool to set a controllable and clearness experimental platform to simulate the strongly correlated BHL, in which the interaction between trapped fermionic cold atoms can be tuned by the Feshbach resonance [25][26][27][28][29][30][31] .
The dynamical mean-field theory (DMFT) has been proved to be a very useful and effective tool [32][33][34][35] , which has made significantly progress as in the study of metal-insulator phase transition. The DMFT is exact when investigating the strongly correlated system in the infinite-dimension, in which the self-energy is independent of momentum. However, in low-dimensional systems, the quantum fluctuation and short range correlations play important roles, which are ignored in DMFT. The cellular DMFT (CDMFT), as a cluster extension of DMFT, effectually incorporates the spatial correlations by mapping the many-body problem into finite clusters embedding in a self-consistent media [36][37][38][39] . In two-dimensional systems CDMFT is more precise than DMFT, and is more effective to investigate the phase transition in low and multi-component systems 31,40,41 .
In this report, we investigate the finite temperature metal-insulator and magnetic phase transition in strongly correlated bilayer honeycomb lattice (BHL) by combining the cellular dynamical mean-field theory (CDMFT) with continue-time quantum Monte Carlo (CTQMC) method 42 . A phase transition from paramagnetic phase to anti-ferromagnetic phase occurs by investigating the magnetization. In a proper value of inter-layer hopping, a novel layer anti-ferromagnetic phase emerges. The layer anti-ferromagnetic phase can be broken by the increasing inter-layer hopping, and a paramagnetic phase can be found. Moreover, the nonlocal inter-layer hopping plays an important role on localizing the free election and modifying the spatial distribution of the electron in lattice sites, especially in dimer sites. For instance, the dimer sites of our system are double occupied, and the nondimer sites are single occupied. We have also presented the DOS, double occupancy, and fermi surface below, which can be directly detected in future experiments.

Results
The strongly correlated bilayer honeycomb lattice. Fig. 1(a) shows the lattice structure of BHL. A 1 and B 1 denotes the sites on top-layer, while A 2 , B 2 signs the bottom-layer sites. The inter-layer band connects A 1 (top) and A 2 (bottom) sites. The first Brillouin zone can be found in Fig. 1(c), in which the C, K, M and K9 shows the points with different asymmetry in the k-space. The non-interaction density of states (DOS) for different sites when t 1 5 t is shown in Fig. 1(d), which is different from the mono-layer honeycomb lattice 43 . The low-energy dispersion is quadratic, which is linear in the monolayer case. There is no band gap between the conduction and valence bands in BHL. In this report, we investigate the phase transition on the half-filling BHL. We consider the standard Hubbard model: where t is the intra-layer nearest-neighbor hopping energy, U is the Coulomb interaction, t 1 denotes the inter-layer hopping energy, m is the chemical potential, which is adjusted to keep the system at halffilling, c { isa and c isa denote the creation and annihilation operators of fermionic atoms on site i with spin s respectively, and a shows the layer parameter (a 5 0 means the top-layer, and a 5 1 denotes the bottom-layer), n isa~c { isa c isa corresponds to the density operator. Here we set t 5 1, which is also used as the energy unit in this report.
The metal-insulator phase transition. The t 1 2 U phase diagram is shown in Fig. 2. The phase boundary of the metal and insulator is shown by the red solid line with square plots. When t 1 5 0 the critical point of metal-insulator phase transition is U c 5 4.7, meaning that BHL degenerates into a mono-layer one. This result is the same as Ref. 16. In order to study the magnetic phase transition, we define the staggered magnetization as M 5 S i sgn(i)(, n i" . 2 , n i# .), where i denotes the lattice index, sgn(i) 5 1 for A 1 and A 2 sites and sgn(i) 5 21 for B 1 and B 2 sites. The phase boundary of non-magnetic state and magnetic state is shown by blue solid lines with circle points. As shown in Fig. 2, when 2 , t 1 , 4.8, anti-ferromagnet can be found. When U is increased, a phase transition from gapless antiferromagnetic metal (AFM), where DE 5 0, M ? 0 and DE here is single particle gap, to gapped anti-ferromagnetic insulator (AFI), where DE ? 0, M ? 0, occurs. When t 1 , 2, a paramagnetic metal (PM), where DE 5 0, M 5 0, is formed, and a PM-AFM-AFI transition happens when U is increased. In the value of 0 , t 1 , 0.2, when U . 4.7, a small region, which is named as paramagnetic insulator (PI) is found with DE ? 0 and M 5 0.
Double occupancy is always used to measure the localization of the electrons directly, which is an important parameter to observe the metal-insulator phase transition 44 . In this report, we investigate the double occupancy D occ~L F=LU~1 4 S i vn i: n i; w as a function of   inter-layer hopping t 1 for various interaction U [see in Fig. 3], where F is the free energy. We use D AA to describe the D occ for A 1 and A 2 sites (hollow circles in Fig. 3), and D BB the D occ for B 1 and B 2 sites (solid circles in Fig. 3). When t 1 5 0, D occ decrease when U is increased showing the localization is enforced by the interaction. When t 1 ? t, D AA is different with D BB . It is found that when t 1 is increased, D AA increases while D BB decreases. This suggests the itinerancy of electrons in dimer sites is enhanced due to the increasing inter-layer hopping. However, be different with D AA , D BB decreases while t 1 is increased. These results suggest the intra-hopping between dimer and non-dimer sites is weaken due to the forming of spinpolarized electrons in dimer sites.
We derive DOS from the imaginary time Green's function G(t) to observe the single particle spectral for different interaction U and t 1 by the maximum entropy method 45,46 . Fig. 4(a) and Fig. 4(e) shows the DOS for different inter-layer hopping when U 5 2.5, T 5 0.1 for A sites and B sites. It is found that, in weak interaction, the inter-layer hopping t 1 does not affect the metallic properties of BHL. In Fig. 4(b) and Fig. 4(f), we can find that the system keeps at a metallic state when t 1 5 1.8, and an obvious pseudo-gap is found when t 1 5 2.2. A metal-insulator transition happens when t 1 5 3.2 for U 5 3.5 and T 5 0.1 for both A and B sites. When interaction being large, such as U 5 6.0 [see Fig. 4(c) and Fig. 4(g)], the system stays at an insulating state, which is insensitive of t 1 . Fig. 4(d) and Fig. 4(h) show the DOS for different U when t 1 5 1.0 and T 5 0.1. An opened gap is found when the interaction is increased, which indicates a phase transition from metal to insulator. These results mean that the inter-layer hopping and Coulomb interaction play same roles for the metal-insulator phase transition.
In this report, the Fig. 5 shows the spectral function A(k, v 5 0) near the Fermi surface for various U and t 1 . A(k, v) is defined as: where iv is the Fermionic Matsubara frequency, E k denotes the dispersion relation, S k (iv) corresponds to the k-dependent self-energy, k is the wave vector in the first Brillouin zone. Figs The magnetic phase transition and a novel layer anti-ferromagnetic phase. In strongly correlated BHL, charge imbalance between the two sublattices sites causes different kinds of magnetic spontaneous symmetry breaking, dividing the sites into dimer sites and nondimer sites. In order to study the magnetic order, we use a magnetic order parameter defined as m a~1 N c where a denotes the sublattice A 1 , A 2 , B 1 and B 2 , i a is the lattice index for sublattice a, and N a means the site number for sublattice a [see in Fig. 1(a) and Fig. 1(b)]. Parameter , n is . indicates the electron density in lattice site i with spin index s (we set magnetization of A 1 as positive sign). Fig. 6 shows the evolution of m a as a function of t 1 for U 5 3.8. It is found that when U 5 3.8, m is zero for t 1 , 1.0, denoting that no magnetic order is formed for weak t 1 . A magnetic state with anti-ferromagnetic order is found when t 1 . 1.0. When t 1 . 4.4, the m A 1 and m A 2 decrease to zero, while m B 1 , m B 2 are still finite, meaning a phase with non-magnetic order and magnetic order coexisting. This novel phase is called layer anti-ferromagnetic state, in which the dimer sites are non-magnetic and non-dimer sites are magnetic. The sketches of the possible magnetic order existing in BHL is shown in Figs. 7(a) -(c). Finally, the phase diagram of magnetization about t 1 and U for T 5 0.1 is shown in Fig. 7. A phase transition from paramagnet to antiferromagnet is found when U is increased. A layer anti-ferromagnetic phase is found for large t 1 , in which the magnetic of dimer sites keeps zero while non-dimer sites is nonzero. When U . 3.7, a phase transition from layer anti-ferromagnetic state to paramagnetic state is found. Experimental protocol. We propose an experiment setup to investigate the phase transition in strongly correlated bilayer honeycomb lattice (BHL). The 40 K atoms can be produced as a pure fermion condensate by evaporative cooling 47 , which provides two hyperfine states jF,m F ae 5 j9/2, 29/2ae ; j "ae and jF,m F ae 5 j9/ 2, 27/2ae ; j "ae 48 . Three standing-wave laser beams are used to form the honeycomb lattice, and two extra laser beams along the z direction suppress the tunneling between layers 28 . The potential of optical lattice is given by V h (x, y) 5 V 0 S j51,2,3 sin 2 [k(x cos h j 1y sin h j )1p/2], where h 1 5 p/3, h 2 5 2p/3, h 3 5 0. Then, we use another three standing-wave laser beams with a 2p/3 angle between each other to form triangular lattice. The potential is given by and k y are the two components of the wave vector k 5 2p/l in these two types of lattices, where l 5 738 nm is wavelength of the laser, and V 0 is given in recoil energy E r 5 h 2 k 2 /2m. Inserting the triangular lattice between two layers of honeycomb lattice, the Bernal stacking BHL with trapped 40 K atoms can be formed 49,50 . In BHL the intra-layer is adjusted by the periodic potential of laser beam and t 1 can be tuned by changing the wavelength of laser beam in z direction. The on-site interaction U~ffi ffiffiffiffiffiffi ffi 8=p p ka s E r V 0 =E r ð Þ 3=4 is determined by the s-wave scattering length a s , which can be tuned by Feshbach resonance. The temperature can be extracted from the time-of-flight images 51 . It should be mentioned that, we can detect the numbers of double occupied sites to confirm whether the metal-Mott insulator transition happens. Firstly, we have to increase the depth of the optical lattice to prevent further tunneling of atoms. Next, the energy of the atoms on doubly occupied sites is shifted by approaching a Feshbach resonance. Then the one spin component of atoms on double occupied sites is transferred to a new magnetic sublevel by radiofrequency pulse method. Finally, the double occupancy can be deduced by the absorption images 52,53 .
To get Fermi surface in experimental, we ramp down the optical lattice slowly enough first, and the atoms stay adiabatically in the lowest band while quasi-momentum is approximately conserved. Then the lattice potential is lowered to zero rapidly, by switching off the confining potential and the atoms can expand ballistically for several milliseconds. The Fermi surface can be obtained by a absorption image 54,55 .

Discussion
In this work, we have investigated the metal-insulator transition and magnetic phase transition in strongly correlated bilayer honeycomb lattice using cellular dynamical mean-field theory (CDMFT) combining with continue-time quantum Monte Carlo (CTQMC) method. In lower temperature case we map the phase diagram as a function of interaction U, inter-layer hopping t 1 and magnetization m. It shows that the inter-layer hopping affects the electrons to form spin-polarized electrons, and an insulating state is induced. A layer anti-ferromagnetic phase is found at large t 1 , in which the magnetization of dimer sites is zero while non-dimer keeps finite value. Therefore, the inter-layer hopping t 1 plays an important role to form a singular magnetic spontaneous symmetry breaking phase. Our study may provide a helpful step for understanding the interaction and inter-layer hopping driven metal-insulator transition, the exotic magnetic order with asymmetry breaking and the possible magneticnonmagnetic order coexisting state.

Methods
The cellular dynamical mean-field theory. We combine the cellular dynamical mean-field theory (CDMFT) with continuous time quantum Monte Carlo (CTQMC) method to determine the metal-insulator transition and magnetic phase transition in the strongly correlated bilayer honeycomb lattice. In low-dimensional systems, quantum fluctuations are much stronger than the higher dimensions. The nonlocal effect is much important in this case. Dynamical mean-field theory ignoring the nonlocal correlations leads lots of errors in calculation. Therefore, we use CDMFT, as the advanced method in our work. We map the original lattice onto a 12-site effective cluster embedded in a self-consistent bath field [see Fig. 1(b)]. Starting with a guessing self-energy S(iv) (which is independent of momentum 56 ), we can get the Weiss field G 0 (iv) obtained by the coarse-grained Dyson equation: where iv is Fermionic Matsubara frequency, m is the chemical potential, k is in the reduced Brillouin zone of the super-lattice, and t(k) is hopping matrix for the superlattice. The form of t(k) is: where d 1~e ik : a1 , d 2~e ik : a1{a2 ð Þ , d 3~e {ik : a2 and a 1 , a 2 are real lattice vectors as shown in Fig. 1(a). The cluster Green's function G(iv) can be gotten by the impurity solver. In our work, we use the numerically exact CTQMC simulation as impurity solver and take 5 3 10 6 QMC sweeps for each CDMFT loop 42 . The new self-energy S(iv) is recalculated by the Dyson equation: This iterative loop repeated until self-energy is converged. The CTQMC method as impurity solver can be taken as follows. We start the procedure at partition function, which can be written as:  interaction U and weak inter-layer hopping t 1 , the system is paramagnet. When we increase U, the system undergoes a magnetic phase transition to anti-ferromagnet. When we increase t 1 , the magnetization of A 1 /A 2 sites decrease to zero while B 1 /B 2 sites stay finite. The system becomes layer antiferromagnet. In the region where t 1 . 5.0, the system returns to paramagnetic phase. where T t is time-ordering operator, H 1 t ð Þ~e tH0 H 1 e {tH0 is H 1 in the interaction picture, and Z 0~Tr e {bH0 is a partition function for the unperturbed term. Putting H 1 5 U S i n i# n i" in Eq. 4, the partition function is Here AEae 0 indicates a theromdynamic average with respect to e {bH0 . Using Wick's theorem, for each order in k, T t n s r 1 ð ÞÁÁÁn s r k ð Þ h i 0 s~:,; ð Þcan be written as determinant detD(k): : : : : : : : where G 0 is non-interacting Green's function. There is no spin index in D(k) for the determinants of spin-un and -down being equivalent. Like classical Monte Carlo, by integrand of Eq. 5, we can get the weight of order k where dt 5 b/L is slice of imaginary time. We can get the standard Metropolis acceptance ratio R of adding vertex by the detailed balance condition: Here P kRk11 is the probability to increase the order from k to k 1 1 (P k11Rk the probability to decrease the order from k 1 1 to k), 1 L : N is probability to choose a position in time and space for vertex you intend to add while 1 kz1 is the probability to choose one vertex you intend to remove of from the existing k 1 1 noes. To calculate the ratio R, we have to deal with the function detD(k 1 1)/detD(k).
M s k ð Þ~D {1 s k ð Þ, we can easily get the value of l in matrix form: : Á Á Á : : : : Á Á Á : : Then it is easy to obtain the update M for the order k 1 1 by numerical method: where the factor of the matrix is M' i,j~M k ð Þ i,j zL i,kz1 l {1 R kz1,j , R i,j 5 G 0 (i, l)M(k) l,j and L i,j 5 M(k) i,l G 0 (l, j). For the step k 2 1, we can also get the radio R and update formulas of M(k 2 1):