# Supersolid phases of dipolar fermions in a 2D lattices bilayer array

## Abstract

Supersolid phases as a result of a novel coexistence of superfluid and density ordered checkerboard phases are predicted to appear in ultracold Fermi molecules confined in a bilayer array of 2D square optical lattices. We demonstrate the existence of these phases within the inhomogeneous mean field approach. In particular, we show that tuning the interlayer separation distance at a fixed value of the chemical potential produces different fractions of superfluid, density ordered and supersolid phases.

###### pacs:

67.85.-d, 03.75.Ss, 03.75.Gg## I Introduction

Several phases of matter appear in the quantum degenerate regime only, namely, superconductivity, superfluidity (SF) and supersolid (SS) phases. While the two formers have been successfully explained, the supersolid phase (1); (2), consistent with thermodynamic stability criteria (3), remains as a subject of theoretical investigation. On the other hand, not at necessarily low temperatures, but also manifesting many body quantum statistical behavior, the High- (HTc) superconductivity phenomenon still continues as an open question in the context of condensed matter. Experiments with ultracold neutral gases are at the present time the closest candidates to quantum simulate, and thus address the description, of such not quite yet understood quantum phases (4); (5); (6); (7); (8); (9).

Recent experimental studies have shown how quantum phases as SF, SS, Mott insulator and charge density wave, emerge from competing short- and long-range interactions among ultracold Bose atoms confined in an optical lattice coupled to a high finesse optical cavity (10). Those phases arise as a result of exploiting the matter-light coupling since in such a case interactions can be tuned on demand. There is however an alternative way of handling either the range and direction of interactions in ultracold neutral gases, which is by confining dipolar atoms or molecules in optical lattices. As it has been shown from the theoretical perspective, the combination of both, the long range anisotropic character of dipolar interactions and the controllable lattice structure where the atoms/molecules lie, make the many-body physics becomes very rich (11); (12); (13); (14).

In this work we consider a model proposed previously (15); (16); (17) to demonstrate that ordered density wave (DW), SS and SF phases can be accessed by changing the external fields that set the system. In Fig. 1 we show a scheme of the quantum simulator that can be created in the laboratory to explore the referred phases. The model system is composed of dipolar Fermi molecules lying in a bilayer array of square lattices in 2D. Although such a configuration has not been realized yet, the current experimental panorama of ultracold dipolar gases, in particular the potential capacity of loading long-lived Fermi molecules of NaK, KRb and NaLi (18); (41) in optical lattices as well as the recently produced ro-vibrational ground state in molecules of NaK, is promising in setting the array here considered.

Previous mean field analysis on dipolar fermions placed onto a single-layer square lattice, with arbitrary orientation and considering dipoles with fixed orientation, have predicted the melting among SF and DW phases (20); (21); (22); (23) and a variety of DW phases (24) respectively. Also, an extended model including a mixture of Fermi molecules with contact interactions, loaded in a bilayer array predicted density ordered phases as well as superfluids phases (15); (25); (26); (27); (28). The possibility of supersolid phases in these dipolar Fermi gases has also been studied (29); (21); (22). On the other side SF, Mott insulating, DW and SS phases of He have been investigated within a mean field context too (32); (31); (30). In the present study we consider dipolar Fermi molecules situated in a double array of parallel optical lattices having dipole orientation perpendicular to the lattice in the presence of a harmonic trap, to demonstrate that in addition to SF and DW patterns there is a region of coexistence in the phase diagram where SS phases emerge. Working within the Bogoliubov-de Gennes (BdG) approach we show that depending upon carefully controlled parameters these phases can be accessed under current experimental conditions.

The paper is organized as follows. In Section II we introduce the model considered in our study and describe the theoretical approach employed. In Section III we illustrate the coexistence and spatial overlap among superfluid and DW phases for several values of the temperature. We summarize our findings in the phase diagram at finite and zero temperature. Finally, we present our conclusions in Section IV.

## Ii Model

We consider Fermi molecules of dipole moment and mass lying in two parallel square lattices of lattice constant separated by a distance , and a harmonic trap with frequency (see Fig. 1). In the presence of an electric field perpendicular to the layers, the dipoles align along the same direction. Fermions in the same layer repel each other always, however, dipoles in different layers attract each other at short range, while also repelling each other at large distances. Thus, interaction between fermions in the same and in different layers is given respectively by

(1) |

where is the in-plane distance between two fermions. Greek indices label the layer where the molecule is placed. Thus, superscripts () indicate that interaction occurs between fermions in same (different) layers. For clarity, we denote the intralayer interaction by , and the interlayer interaction by . The system is described by the Hubbard model with the Hamiltonian given by , with

where are the standard creation and annihilation operators and , is the in-plane energy dispersion of the ideal Fermi gas within the tight binding approximation, being the hopping among nearest neighbors. and are the Fourier transforms of and , respectively and the number of sites. The terms containing the harmonic confinement are written in the Fock basis of sites, where the vector position in the lattice is denoted by , and the frequency of the harmonic trap that confines the molecules. In what follows, all the energies will be scaled with respect to . We also introduce two relevant physical quantities: the dipolar interaction strength , with the effective mass, and the dimensionless parameters and .

The proposed model can be mapped into a system of fermions in two different hyperfine spin states () confined in a 2D lattice. The terms and describe repulsive and attractive interactions among fermions in the same and different hyperfine states respectively. We should stress that is an interaction that is attractive at short distances while becoming repulsive at long distances . Thus, by controlling the separation between the layers , the proposed model represents a promising candidate to study the quantum phases from competing short- and long-range interactions as well as attractive versus repulsive interactions. The inclusion of the harmonic potential plays also a crucial role since, as it is well known, the global thermodynamics of the phase transition is qualitatively different from that of the homogeneous case (33). A remarkable signature of this fact is that the magnitude of the coherence length can be as large as the typical confinement distance.

To investigate the physics of the model described above, we use mean-field theory including the usual BCS pairing terms and the Hartree contributions(34). We expect this approximation to be reasonably accurate in the weakly interacting regime. The mean-field Hamiltonian is diagonalized by solving the Bogoliubov-de Gennes equations (BdG) (35). This mean-field approach is commonly used for studying competing magnetic and superconducting phases in the context of high superconductivity (13); (36) and in the context of ultra cold fermions (37). It also has been recently employed for describing strongly correlated systems, like effective -wave interaction and topological superfluids in quantum gases in lower dimensions systems (1D and 2D)(8). The equation to be solved is,

where the matrix elements incorporate the tunneling among nearest neighbors , the effect of the harmonic confinement and the inter site interaction on the Hartree level that is expected to dominate (34), that is, with the Kronecker delta for nearest neighbors. is the superfluid parameter. The eigenvalues denoted by are self-consistency obtained through the usual relations and with the local Bogoliubov quasiparticle amplitudes, being the Fermi distribution and the expectation value of . For simplicity we shall assume that both layers are equally populated.

To include the effects of the harmonic trap, in addition to the usual order parameters that globally describe the system, we introduce local order parameters. The global density order parameter is given by , where the vector identifies either, a checkerboard pattern , or a stripe density order or . The local density order parameter is given by , where the index denotes that the sum runs over the first and second nearest neighbors. This local order parameter describes a checkerboard density pattern in a sublattice centered at site . The superfluid local order parameter is given by , being the average superfluid behavior studied through .

## Iii Supersolid: coexistence of Superfluid and Density Ordered phases

To determine the density profile and the behavior of the gap across the lattice, we solve BdG equations for lattices of size , maintaining fixed the value of chemical potential (38). This restriction causes that the total number of fermions is increased with temperature, that is, at zero temperature the number of fermions is while for this value is increased to . We also keep fixed the values of the interaction strength and the harmonic frequency at and , respectively.

To illustrate the competition among different phases we have selected the cases and . We plot the density profile and the gap parameter profile for several values of the temperature. In particular, we chose and .

First we start considering . From Fig.2 we observe an homogeneous distribution at the center of the trap at zero temperature for both, density and gap profiles. However, at finite temperature, the distribution of the density profile remains homogeneous at the center of the trap, while the gap structure shows a decreasing ratio as temperature is increased until they vanishes at a temperature of . In previous studies (16); (17) reported a BCS superfluid phase in the weakly interacting regime while occurring formation of dimers in the strong interaction regime, leading those dimers to a Bose superfluid. In the present work we focus on the BCS superfluid to DW-Supersolid phase transition. Therefore, the values of the parameters are restricted to those in which the weakly interacting regime is ensured.

When the interlayer spacing is increased, the competition between attractive and repulsive dipole interactions becomes evident. As plotted in Fig. 3, for at low temperatures, there is a large region in the center of the trap having a superfluid order parameter coexisting with a checkerboard density order. This is the signature of a supersolid phase. When the temperature is increased, the radius of the superfluid order parameter shrinks and completely vanishes at a critical temperature of , while the checkerboard phase melts at a temperature of . That is, the supersolid phase exists, in this system, for certain temperatures, as further shown below in the corresponding phase diagram. Cross sections of the local order parameters are shown in Fig.4, where it can be appreciated the presence of a supersolid phase at the center of the trap, as the spatial overlap of both superfluid and DW.

For the repulsive intralayer interaction starts to dominate at the center of the trap. From Fig. 5 one can see that there is a wide region at the center of the trap exhibiting a checkerboard DW pattern. Such patterns persist below temperatures of . We also observe that the superfluid order parameter appears to surround the checkerboard pattern and that such a superfluid disk completely vanishes when the temperature reaches a value of . The cross sections shown in Fig. 6 exhibit a small region where both phases spatially overlap. Other studies (39); (40) with different systems in the presence of a harmonic trap have shown coexistence of phases without spatial overlapping, for instance, in the extended Bose-Hubbard model the Mott insulator and superfluid phases tend to form rings and disks. In 2D those studies agree quantitatively with Quantum Monte Carlo calculations and more sophisticated methods.

Finally, when the interlayer spacing is large enough, the intralayer repulsive interaction dominates over the interlayer attraction and the superfluid order parameter almost vanishes. This behavior is found for , where pairing is inhibited. For larger values of no pairs can be formed but a DW checkerboard pattern still persists at the center of the trap (see Fig. 7). This value of signals the limit from which each layer can be studied separately. The single layer model has been studied previously considering arbitrary dipole moment orientations(21); (22). We found that our predictions are in good agreement with those results, that is, at a given critical temperature that depends on the interaction strength, checkerboard phases for perpendicular orientation of the dipole moment emerge.

In Fig. 8 we plot the two global order parameters and as a function of the temperature for values of in the region of coexistence. As can be appreciated from this figure, Bogoliubov-de Gennes diagonalization predicts continuous phase transitions for the considered model. One can observe that, for a given value of , the critical temperature at which the superfluid phase emerges coincides with that at which the derivative of the DW order parameter shows a discontinuity. Numerical calculations performed for lattices of larger size () lead us to observe how the discontinuity of the derivative in and at the critical temperature becomes more evident as a function of . Namely, the global order parameters and change more abruptly at the critical temperature as the lattice size is increased. It is also important to stress that the referred discontinuity in the derivative becomes less evident when the interlayer spacing is increased.

In Fig. 9 we present the phase diagram of this model, obtained from Bogoliubov-de Gennes equations for finite temperatures. In the inset we show the phase diagram at zero temperature. The region of coexistence between superfluid and DW phases in both diagrams is the supersolid phase of our system. As expected, in the attractive interaction regime the superfluid phase destroys any density order pattern. When the interlayer spacing becomes comparable with the lattice constant , the superfluid phase and DW start to compete. For there is no formation of density order pattern, while the critical temperature of the BCS superfluid phase decrease monotonously. Close to a density order pattern is formed, then the critical temperature of this phase jumps quickly to a constant value. For values of larger than the critical temperature of the DW phase becomes almost independent of the interlayer spacing. On the other hand, the superfluid parameter suddenly decrease at and then again starts to decrease monotonously. In contrast with the predictions obtained for the single layer system where the critical temperatures may be calculated using the value of the parameters at the center of the trap, this may not be completely true for the system here studied due the possible formation of disks and rings.

The maximum value of the critical temperature for the supersolid phase predicted by our model, considering the parameters of the current experimental systems, is . Although such temperature is one order of magnitude smaller than those measured recently in experiments with fermonic KRb (41) and NaK (18), current efforts in controlling and lowering the temperature of molecules are promising to reach such critical temperatures in the near future.

## Iv Conclusion

We have studied the thermodynamic phases that exhibited dipolar Fermi molecules placed at the sites of a bilayer array of square optical lattices in 2D in the presence of a harmonic confinement. Due to the nature of the dipolar interaction, where attractive and repulsive interactions are present, several phases are shown to form. While at- tractive interaction between molecules in different layers leads to predict superfluid phases, density order phases like checkerboard patterns result from the repulsive interaction. The competition between these phases gives rise to the formation of supersolid phases where both SF and DW phases coexist and spatially overlap. An exhaustive exploration of the space of parameters are summarized in the phase diagrams at zero and finite temperatures, see Fig. 9. Our predictions allowed us to identify clearly the influence of the harmonic potential in the occurrence of the transitions with respect to thermodynamics in the homogeneous case reported in previous literature. The system here studied in combination with the capability of trapping Fermi molecules in optical lattices as well as the recently reported production of rovibrational and hyperfine ground state of NaK molecules constitute a promising candidate to study the competition between BEC and BCS superfluid phases in coexistence with an ordered structure, and thus, offering the opportunity to quantum simulate a supersolid phase in ultracold experiments.

## V Acknowledgments

This work was partially funded by grants IN107014 DGAPA (UNAM) and LN-232652 (CONACYT). A.C.G. acknowledges a scholarship from CONACYT.

### References

- E. Kim and M. H. W. Chan, Nature 427, 225 (2004).
- H. Choi, S. Kwon, D. Y. Kim and E. Kim, Nat. Phys. 6, 424 (2010).
- M. Mendoza-López and V. Romero-Rochín, arXiv:1604.08078 (2016).
- W. Hofstetter, J. I. Cirac, P. Zoller, E. Demler and M. D. Lukin, Phys. Rev. Lett. 89, 220407 (2002).
- C.-K. Chan, C. Wu, W.-C. Lee, and S. Das Sarma, Phys. Rev. A 81, 023602 (2010).
- B. Liu, X. Li, B. Wu and W. V. Liu, Nat. Comm. 5, 5064 (2014).
- A. Bühler, N. Lang, C. V. Kraus, G. Möller, S .D. Huber and H.P. Büchler, Nat. Comm. 5, 4504 (2014).
- B. Wang, Z. Zheng, H. Pu, X. Zou and Guangcan Guo Phys. Rev. A 93, 031602(R) (2016).
- Y. Fujihara, A. Koga and N. Kawakami, Phys. Rev. A, 81, 063627 (2010).
- R. Landig, L. Hruby, N. Dogra, M. Landini, R. Mottl, T. Donner and T. Esslinger, Nature 532, 476 (2016).
- M. A. Baranov, Phys. Rep. 464, 71 (2008).
- T. Lahaye, T. Koch, B. Frölich, M. Fattori, J. Metz, A. Griesmaier, S. Giovanazzi and T. Pfau, Nature 448, 672 (2007); T. Lahaye, C. Menotti, L. Santos, M. Lewenstein and T. Pfau, Rep. Prog. Phys. 72, 126401 (2009).
- Y. Chen, Z. D. Wang, F. C. Zhang and C. S. Ting, Phys. Rev. B 79, 054512 (2009).
- N. T. Zinner, B. Wunsch, D. Pekker, and D. -W. Wang, Phys. Rev. A 85, 013603 (2012).
- T. I. Vanhala, J. E. Baarsma, M. O. J. Heikkinen, M. Troyer, A. Harju, and P. Törmä, Phys. Rev. B 91 ,144510 (2015).
- Camacho-Guardian, A. and Paredes, R. (2016), Ann. Phys.. doi: 10.1002/andp.201600101.
- F. Ancilotto, Phys Rev. A 93, 053627 (2016).
- Jee Woo Park, Sebastian Will, and Martin W. Zwierlein, Phys. Rev. Lett. 114, 205302 (2015).
- K.K. Ni, S. Ospelkaus, M.H.G. de Miranda, A. Peer, B. Neyenhuis, J.J. Zirbel, S. Kotochigova, P.S. Julienne, D.D, Jin, and J. Ye, Science 322, 231 (2008).
- N. T. Zinner and G. M. Bruun, Eur. Phys. J. D 65, 133 (2011).
- Anne-Louise Gadsbolle, and G. M. Bruun, Phys. Rev. A 85, 021604(R) (2012).
- Anne-Louise Gadsbolle, and G. M. Bruun, Phys. Rev. A 86, 033623 (2012).
- S. G. Bhongale, L. Mathey, S. -W. Tsai, C. W. Clark, and E. Zhao, Phys. Rev. Lett. 108, 145301 (2012).
- K. Mikelsons and J. K. Freericks, Phys. Rev. A, 83 043609 (2011).
- Y. Prasad, A. Medhi, and V. B. Shenoy, Phys. Rev. A 89, 043605 (2014).
- M. A. Baranov, A. Micheli, S. Ronen and P. Zoller, Phys. Rev. A 83, 043602 (2011).
- A. Pikovski, M. Klawunn, G. V. Shlyapnikov and L. Santos, Phys. Rev. Lett. 105, 215302 (2010).
- A. C. Potter, E. Berg, D. -W. Wang, B. I. Halperin and E. Demler, Phys. Rev. Lett. 105, 220406 (2010).
- L. He and W. Hofstetter, Phys. Rev. A 83, 053629 (2011).
- P. Nozières, J. Low Temp. Phys. 156, 9 (2009).
- J. Ye, Phys. Rev. Lett. 97, 125302 (2006).
- S. Rica, Phys. Rev. B 84, 184535 (2011).
- I. Reyes-Ayala, F. J. Poveda-Cuevas, J. A. Seman and V. Romero-Rochín, arXiv:1607.02389v1 (2016).
- E. Müller-Hartmann, Z. Phys. B 74, 507 (1989).
- Self-Consistent Bogoliubov-de Gennes diagonalization was iterated until convergence criteria was reached. The convergence criteria was reached when using the Frobenius norm for the matrices norm.
- Y. Chen, Z. D. Wang, J.-X. Zhu, and C. S. Ting, Phys. Rev. Lett. 89, 217001 (2002).
- B. M. Andersen and G. M. Bruun, Phys. Rev. A 76, 041602(R) (2007).
- The chemical potential remains as a global thermodynamic quantity, even within the local density approximation. V. Romero-Rochín, Phys. Rev. Lett. 94, 130601 (2005).
- J. M. Kurdestany, R. V. Pai, and R. Pandit, Ann. Phys. (Berlin) 524, 234 (2012).
- R. V. Pai, J. M. Kurdestany, K. Sheshadri and R. Pandit, Phys. Rev. B 85, 214524 (2012).
- K.-K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Peer, B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne, D. S. Jin, and J. Ye, Science 322, 231 (2008).