# Dynamical recovery of SU(2) symmetry in the mass-quenched Hubbard model

###### Abstract

We use non-equilibrium dynamical mean-field theory with iterative perturbation theory as an impurity solver to study the recovery of symmetry in real-time following a hopping integral parameter quench from a mass-imbalanced to a mass-balanced single-band Hubbard model at half-filling. A dynamical order parameter is defined to characterize the evolution of the system towards symmetry. By comparing the momentum dependent occupation from an equilibrium calculation (with the symmetric Hamiltonian after the quench at an effective temperature) with the data from our non-equilibrium calculation, we conclude that the symmetry recovered state is a thermalized state. Further evidence from the evolution of the density of states supports this conclusion. At the same time, we find the order parameter in the weak Coulomb interaction regime undergoes an approximate exponential decay. We numerically investigate the interplay of the relevant parameters (initial temperature, Coulomb interaction strength, initial mass-imbalance ratio) and their combined effect on the thermalization behavior. Finally, we study evolution of the order parameter as the hopping parameter is changed with either a linear ramp or a pulse. Our results can be useful in strategies to engineer the relaxation behavior of interacting, quantum many-particle systems.

## I Introduction

Research on non-equilibrium quantum phase transitions have seen dramatic progress and in the past decade.Lindner et al. (2011); Gull et al. (2011); Wang et al. (2013); Aoki et al. (2014) Many experimental and theoretical studies have focused on pump-probe experiments in solid state systems where the system is driven out-of-equilibrium by a pump laser,Manmana et al. (2007); Rigol et al. (2007); Kollath et al. (2007) or in cold atom systems drivenMentink and Eckstein (2014); Mentink et al. (2015); Eckstein and Werner (2016); Bukov et al. (2016) by Coulomb interaction strength change. In a real-time quantum phase transition, symmetry breaking or recovery often plays an important role. For example, the phase transition between paramagnetic and anti-ferromagneticWerner et al. (2012); Tsuji et al. (2013); Tsuji and Werner (2013) behavior as a function of Coulomb interaction quench is associated with the spontaneous symmetry breaking of lattice symmetries and time-reversal; the transition to a Floquet topological insulator in BiSe,Lindner et al. (2011); Wang et al. (2013) Graphene,Oka and Aoki (2009) and bilayer (LaNiO)/(LaAlO) thin filmsDu et al. (2017a); Du and Fiete (2017) is triggered by the time-reversal symmetry breaking induced by shinning with a circularly polarized laser but does not entail any corresponding breaking of lattice symmetries.

If the electronic system driven out-of-equilibrium is strongly correlated, theoretical and numerical techniques to deal with the system are limited, which makes this problem especially challenging. The Hubbard model is widely considered the simplest model to capture the most essential features of strongly correlated systems (either solid state materials or cold atom systems). Non-equilibrium dynamical mean-field theory (DMFT),Georges et al. (1996); Gull et al. (2011); Aoki et al. (2014) in which the original lattice problem is mapped onto a Anderson impurity problem with a self-consistently determined bath, has proved to be a powerful tool in solving the Hubbard model. The non-equilibrium real-time evolution of physical observables can be obtained within the framework of DMFT.

Even though the methods to solve the Anderson impurity model in equilibrium are well developed, solving the model out-of-equilibrium remains difficult and is still under development.Wolf et al. (2014); Balzer et al. (2015); Cohen et al. (2015); Dong et al. (2017) Numerically, the robust impurity solver in equilibrium, hybridization expansion continuous time quantum Monte Carlo, can suffer from the dynamical sign problem and the simulation time is usually rather short. The numerically exact impurity solver, weak coupling continuous time auxiliary-field quantum Monte Carlo, is limited to a single-band model at half-filling and a short time evolution. The dynamical sign problem increase exponentially with evolution time, which limit its applicability. Analytically, the iterative perturbation theory (IPT) impurity solver at weak and the non(one)-crossing approximation (NCA, OCA) at strong Coulomb interactions have been shown to be powerful impurity solvers to capture the physical picture in the weak and strong Coulomb interaction regimes, respectively.Tsuji and Werner (2013); Eckstein and Werner (2010)

Due to the limitations of the impurity solvers for non-equilibrium DMFT, the Falicov-Kimball (FK) model provides a first attack on strongly correlated electronic systems because the projected impurity model is exactly solvable and provides important information on the Hubbard model.Freericks et al. (2006); Eckstein and Kollar (2008). The FK model is well studied both in equilibrium and out-of-equilibrium.Freericks and Zlatić (2003); Freericks et al. (2006) This raises the questions, “What is the physical connection between the FK model and the Hubbard model?” and “To what extent does the FK model reveal the physics of the Hubbard model?”. To answer these questions, previous works focused on studying the two models separately, and compared the results of each model with each other.

For cold atom systems, Eckstein et al.Eckstein and Kollar (2008); Eckstein et al. (2009, 2010) studied the evolution of physical quantities as a function of time driven by a Coulomb interaction quench in the FK model and the Hubbard model, respectively. For the FK model driven by a Coulomb interaction quench, a non-thermal steady state (not a thermalized state) results that can be statistically described by the generalized Gibbs ensemble.Eckstein and Kollar (2008) For the Coulomb interaction quenched Hubbard model, a pre-thermalization behavior and dynamical phase transition are observed.Eckstein et al. (2009) For a solid state system driven by a constant electric field, Freericks et al.Freericks et al. (2006); Freericks (2008) studied the FK model and found damped Bloch oscillations. Eckstein et al.Eckstein and Werner (2011) studied the Hubbard model: except for the damped Bloch oscillation observed in FK model, the current decays to zero and remains there. Fotso et al.Fotso et al. (2015) compared the thermalization behavior of the FK model and the Hubbard model driven by a DC electric field. The FK model can have one of two generic evolution behaviors: (1) either monotonic or oscillatory approach to an infinite-temperature steady state or (2) either monotonic or oscillatory approach to a non-thermal steady state. In addition to the above two features, the Hubbard model can exhibit an extra feature by evolving to an oscillatory state.

In contrast to previous studies, in this work we would like to build the connection between the FK model and the Hubbard model by quenching the hopping parameter of the frozen species (the one that is not able to hop on the lattice) in the FK model to the Hubbard model. In order to avoid the singularity of the FK model (bandwidth for one species is zero), we use the mass-imbalanced Hubbard model with large hopping asymmetry and study the time-dependent evolution of observables following a quench between the mass-imbalanced and mass-balanced Hubbard model. Here mass-imbalance (mass-balance) means the spin- and spin- hopping parameter are unequal (equal) to each other. Previous work has studied in the hopping parameter quench of the Hubbard model for equal strength hopping of two spin species, and can be solved as an Coulomb interaction quench problem with scaled time.Tsuji et al. (2011) In this work, we show a quench on only one hopping parameter leads to rather different physics.

One of the central results of our work is that a dynamical phase transition appears. To put our results in context, it useful to summarize related work that also found dynamical phase transitions. By quenching the Coulomb interaction between two different phase regimes in equilibrium, Tsuji et al.Tsuji et al. (2013) studied the dynamical phase transition between an antiferromagnetic and paramagnetic state. Two dynamical transition points are observed with one the thermal transition and the other related to a non-thermal antiferromagnetic phase.

By contrast, we study the evolution of symmetry recovery by quenching from the mass-imbalanced to mass-balanced Hubbard model. This can be experimentally realized in cold atom systems by tuning the lattice potential amplitude and the recoil energy.Nguyen et al. (2015) In the mass imbalanced Hubbard model, the symmetry is broken. By quenching the hopping integral of one spin species to be the same as the other one, the Hamiltonian recovers its symmetry. However, the evolution of physical observables as a function of time remains unclear. Our work fills that gap. We address the following questions: (1) Is the symmetry recovered state the same as the equilibrium thermalized state? (2) What is the dependence of the evolution process on the Coulomb interaction, temperature, and the initial mass imbalance? (3) How does the time-evolution change if we set up the quench process as a linear ramp or pulse shape?

In this work, we show there an order parameter can indeed serve as a criteria for a dynamical phase transition. As the symmetry in observables is recovered following the quench, the system is thermalized at the same time. We show that the evolution of the order parameter has a monotonic dependence on the mass imbalance, temperature, and Coulomb interaction. The pulse shape influences the time evolution.

Our paper is organized as follows. In Sec.II, we describe the mass-imbalanced Hubbard model and illustrate how we calculate several physical observables within time-dependent dynamical mean-field theory, such as the momentum dependent occupation. We also define the order parameter we use to characterize the symmetry. In Sec.III, we characterize the symmetry recovered state as a thermalized state. The dependence on the Coulomb interaction, the initial temperature, and the initial mass-imbalance ratio is studied. We compute the evolution of the order parameter for a linear ramp and a pulse change of the hopping parameter in Sec.IV. Finally, in Sec.V we summarize the main conclusions of this work.

## Ii Model and Method

The time-dependent mass-imbalanced single-band Hubbard model at half-filling is given by,Dao et al. (2012); Liu and Wang (2015); Philipp et al. (2017); Sekania et al. (2017); Du et al. (2017b)

(1) |

where () creates (annihilates) an electron at site with spin , and is the corresponding density operator. The notation indicates the hopping is restricted to nearest neighbors, () is the time dependent (independent) hopping integral parameter for spin- (spin-) electron ( is reserved to denote time). The time dependence of the hopping parameter can be introduced in optical lattices through lasers. Here denotes the time independent Coulomb interaction strength between spin- and spin- fermions occupying the same site. Throughout this paper, we fix the spin- hopping integral to be time independent and set () as our unit of energy (time). The mass imbalance is restricted to lie between 0 and 1. The system is initially prepared in the thermal equilibrium state of the mass imbalanced Hubbard model with and finite repulsive Coulomb interaction . Here the symmetry of the system is broken. The quench dynamics are studied by fixing the Coulomb interaction to be finite while quenching the spin- hopping integral to be from an initial state, where is the ramp time of the hopping parameter change.

We consider a Bethe lattice, which has a semi-elliptic density of states,

(2) |

with half bandwidth . The mass-imbalanced Hubbard model (1) can be solved exactly using non-equilibrium dynamical mean field theory (DMFT),Georges et al. (1996); Freericks et al. (2006); Eckstein et al. (2009); Gull et al. (2011); Aoki et al. (2014) which maps the lattice model self-consistently onto a single-site Anderson impurity model. We use non-equilibrium dynamical mean field theory with iterative perturbation theory as an impurity solver to solve the mass imbalanced Hubbard model at finite temperature. We enforce the paramagnetic solution and half-filling of both spin- and spin- electrons. In the Hubbard model, these constraints can be fulfilled by explicitly symmetrizing over the two spin spices and setting the chemical potential to be , respectively. Away from this mass balanced Hubbard model limit, we again enforce half-filling by fixing . However to ensure the paramagnetic solution at half-filling, we symmetrize the Weiss’s functions in the Keldysh time contour using particle-hole symmetry: . The DMFT self-consistent condition for the Bethe latticeEckstein and Kollar (2010) density of state is

(3) |

The expectational value of an observable at time is given by,

(4) |

where is the partition function of the non-interacting Hamiltonian at , is the time evolution operator. The momentum dependent density matrix is written as

(5) |

where is the lesser Green’s function at equal time . The momentum-dependent occupation depends only on because the self-energy is momentum independent. The kinetic energy is given by

(6) |

The Coulomb interaction energy is given by

(7) |

where denotes the Keldysh contour.Eckstein et al. (2010) The total energy is

(8) |

The real frequency represented retarded Green’s is function is

(9) |

The density of states are calculated from the exact relation

(10) |

The dynamical order parameter we use to characterize the symmetry is,

(11) |

where is time of the hopping integral ramp to , and is the semi-elliptic density of states defined in Eq.(2) with . If the quenched system is thermalized after long enough time, the effective temperature is calculated by numerically solving the equation,Eckstein et al. (2010)

(12) |

where is the same as Eq.(8), and is the Hamiltonian after quench.

## Iii thermalization driven by hopping quench

The prethermalization behavior in the paramagnetic case occurs when the momentum integrated quantities (Coulomb interaction and spin-resolved kinetic energy) thermalize faster than the momentum dependent quantities (momentum dependent distribution).Tsuji et al. (2013) Thus, we must study the momentum dependent observables as a function of time to determine if the system is thermalized. Depending on the system details, there may be convenient quantities for studying the thermalization. For example, in the Coulomb interaction quench at zero temperature problem, the Fermi-surface dis-continuity in the momentum-dependent occupations can serve as a good criteria because the Fermi-surface jump disappears at the thermalized finite temperature.Eckstein et al. (2009) In our calculation we use the momentum integrated order parameter in Eq.(11) for symmetry as the specific criterion in our case.

In Fig.1, we plot the spin-resolved momentum dependent occupation number as a function of energy at different times with fixed Coulomb interaction and initial mass-imbalance ratio . Here, for the purpose of better vision, we only plot half of energy axis (). The other part of momentum distribution () is constrained by which is hold by particle-hole symmetry. The only parameter difference is for Fig.1(a) and for Fig.1(c). The corresponding effective temperature for the final thermalized state is calculated using Eq.(12). We find for case (a) and for case (c), respectively.

At time , the momentum distribution for spin- and spin- electrons are apparently separated. The black dashed line is the thermalized value of the momentum distribution with symmetry. The area encapsulated by spin- and spin- distribution is defined as the order parameters breaking the symmetry. As time evolves, the area is diminished monotonically. (See and 4 in Fig.1, for example.) Finally, at , the area vanishes which indicates the symmetry of the Hubbard model is fully recovered in the time evolution of the states. By comparing the recovered distribution with the thermalized state at , one sees they match each other, indicating that the symmetry recovered state is just the thermalized state. Therefore, the order parameter defined in Eq.(11) can serve as a measure of whether the state is thermalized. In Fig.1 (b), we plot the order parameter as a function of time. We realize the evolution of the order parameter fits an exponential decay reasonably well: . With and , the fitting function is a good approximation of the original data. Figs.1(c-d) give very similar information except some quantitative difference, mainly the small initial order parameter and small decay rate . A systematic discussion of these difference is deferred to future sections of this paper.

To confirm our conclusion that the symmetry recovered state is a thermalized state, we plot the spin-resolved density of states at different times in Fig.2. Here the density of states is calculated by Fourier transforming the two-time retarded Green’s function to the real frequency axis using Eq.(9). In the non-interacting limit, the density of states for the two spin species are identical after the hopping parameter quench. We used a small Coulomb interaction . The density of states for spin- and spin- exhibit very small differences at time . We checked numerically that larger Coulomb interaction will induce a larger difference in the spin-resolved density of states. As time evolves, the density of states for the two spin species move toward each other and finally meet at . By comparing the density of states for the thermalized state with effective temperature , we confirmed our conclusion in the previous paragraph.

The exponential decay with time of the order parameter depends on the Coulomb interaction, the initial mass imbalance ratio, and the initial temperature. We will study the effect of one of the three factors by fixing the other two. In Fig.3, we plot the order parameter of symmetry (deviation) as a function of time for different Coulomb interactions with fixed initial inverse temperature and initial mass imbalance ratio . In the top panels, the order parameter (log scale -axis) at different Coulomb interactions are plotted as a function of time. An approximate exponential decay is observed. The larger the Coulomb interaction, the faster the decay rate of the order parameter. In the lower panels, the decay rate as a function of Coulomb interaction strength is plotted with open dots.

The qualitative behaviors above can be understood from two limits. In the non-interacting limit, the momentum distribution after a hopping quench for different spins are and , respectively, where is the initial mass imbalance ratio. As time evolves, the momentum distribution does not change because the momentum is a good quantum number (no Coulomb scattering since ). This limit has an infinitely long approach time to the thermalized state ( symmetry recovered). As a result, the decay rate is zero in the non-interacting limit. In the infinite Coulomb interaction limit (atomic limit), the mass imbalance quench can be ignored (since the kinetic energy in either case is negligible) and so has no effect on the thermal distribution. In the infinite Coulomb interaction limit we can take the decay rate . Note, in the large Coulomb interaction region with relative low temperature, an antiferromagnetic state will appear. In this paper we restrict ourself to the weak interaction limit and relative high temperature to ensure we have a non-equilibrium paramagnetic solution. Further, based on second order perturbation theory (the first order terms will cancel due to particle-hole symmetry) a quadratic function of Coulomb interaction is plotted in the lower panels as a solid line, where depends on the specific system parameters. The larger the Coulomb interaction, the greater the deviation from quadratic dependence of (higher order perturbation terms are needed).

By comparing Fig.3 (b) with and (f) with , we realized the decay rate in (f) is smaller than in (b) for each fixed Coulomb interaction. We conclude from this comparison that the decay rate depends on the initial temperature such that the lower temperature, the longer the time needed to relax to the thermalized state. A systematic study and discussion on the initial temperature dependence is illustrated in Fig.4. By comparing Fig.3 (b) with and (d) with , one sees the decay rate in (d) is smaller than in (b) for each fixed Coulomb interaction. We conclude from this comparison that the decay rate depends on the initial mass imbalance ratio such that the larger mass imbalance ratio (closer to the final mass-balance Hubbard model Hamiltonian), longer the time that is needed to evolve to the thermalized state. A systematic study of the evolution dependence on the initial mass imbalance ratio is illustrated in Fig.5.

By fixing the Coulomb interaction to be and the initial mass-imbalance ratio , we plot the order parameter as a function of time at different temperatures in Fig.4(a). At time , one sees the order parameter is larger if the initial temperature is smaller. Since the Coulomb interaction is weak here, this behavior can be understood from the non-interacting limit. In the non-interacting limit, the order parameter after the hopping parameter quench from to is given by,

(13) |

where is the initial mass imbalance ratio.

The order parameter as a function of temperature is plotted in Fig.4(c) with a solid line for . This plot can be understood physically from the zero temperature and infinite temperature limit. In the zero temperature limit the states with energy () are fully occupied (vacant). As one quenches the hopping integral of spin- electrons, the occupation at each is the same as before (occupied up to ) and the order parameter will be zero. In the infinite temperature limit every state is equally populated and independent of its energy. The order parameter will be zero after hopping parameter quench. In between these two limits, there exist a critical temperature where order parameter increases (decreases) monotonically with the temperature below (above) the critical temperature. The initial order parameters for and are plotted with filled squares and open circles, respectively. One sees the biggest deviation from occurs at the critical temperature. Increasing the temperature beyond the critical one will decrease the deviation from the non-interacting limit. This can be explained as a competition between the kinetic energy and the Coulomb interaction. In high temperature region, the kinetic terms overcome the Coulomb terms and dominate the behavior of the order parameter. This picture in the high temperature region can be further confirmed by plotting the decay rate of the order parameter as a function of temperature in Fig.4(b).

First, we confirm the conclusion that in the low temperature regime, the decay rate increases monotonically with initial temperature. Physically, the relaxation process to the thermalized state is driven by the Coulomb scattering. Increasing the initial temperature will enhance the thermal fluctuations of electrons and enhance the collision probability leading to a larger decay rate. However, the decay rate tends to increase slower and finally saturates in the high temperature regime. When the temperate is high (kinetic energy overcomes the Coulomb energy), the states of the initial equilibrium state tend to converge, and the decay rates tend to saturate. Since the temperature is only well-defined in the initial equilibrium states and the final thermalized states, we limit ourself to the qualitative analysis above.

Finally, we plot the thermalized temperature as a function of initial temperature at in Fig.4(d) with filled squares. A comparison with and are plotted with a solid line and open circles. A critical initial temperature is observed at , where Coulomb interaction tend to increase (decrease) the thermalized temperature comparing to the infinitesimal Coulomb interaction limit.

By fixing the Coulomb interaction to be and initial temperature to be , we plot the order parameter as a function of time for different initial mass imbalance ratios in Fig.5. In Fig.5(a), the order parameters is plotted as a function of time. At time , one sees the order parameter is larger if the mass-imbalance ratio is smaller. The initial order parameter as a function of is plotted in Fig.5(c) with filled squares. This can be understood again by considering the non-interacting limit in Eq.(13). The order parameter at and are plotted in Fig.5(c) with a solid line and open circles, respectively. The order parameter decreases monotonically as is increased until the limit (order parameter is zero). The deviation from the non-interacting limit is larger as the initial mass imbalance ratio decreases. Further, in Fig.5(b), the approximate decay rate is plotted as a function of the initial mass imbalance ratio . Our results indicate that the decay rate decreases monotonically with increasing initial mass imbalance ratio. Finally, the thermalized temperature as a function of for is shown with filled squares. To illustrate the effect of Coulomb interaction, the data for and are plotted with a solid line and open circles, respectively. The thermalized temperature decreases as one increases the initial imbalance ratio. In the limit , the initial temperature equals the final thermalized temperature.

## Iv Dependence on the ramp shape and pulse form

Experimentally, the change of parameters in the Hamiltonian takes a finite amount of time. To model this, we suppose there exist a linear ramp to achieve the final parameter,

(14) |

where is the time used to achieve the final recovered Hamiltonian. In Fig.6(a-b), we plot the evolution of the order parameter as a function of time for different quench times at fixed Coulomb interaction and inverse temperature . The decay rates for different linear ramps are approximately the same while a longer quench time leads to a longer relaxation time to a thermalized state. This is consistent with our previous study of the decay rate dependence on the initial mass imbalance: In the linear ramp time (), the ratio difference is smaller, so the decay rate in that time region is smaller. At time , the order parameter will be larger than the quenched case ().

Finally, we studied the case in which we have the mass balanced Hubbard model at time and apply a pulse change to the spin- hopping parameter change,

(15) |

with different quench time (width of pulse) . We plot the order parameter as a function of time in normal scale in Fig.6(c) and log scale in Fig.6(d). As the quench time increases, the order parameter is larger. Taken together, we see that the pulse shape can be used as a way to engineer the relaxation behavior of interacting, quantum many-particle systems.

## V discussion and conclusion

In this work, we theoretically studied the dynamical evolution towards symmetry of a system that is quenched from an broken one (mass-imbalanced Hubbard model) to an symmetry recovered one (mass-balanced Hubbard model). This model can be experimentally implemented in cold atom systems. We define the time dependent order parameter (total momentum-integrated difference between spin- and spin- momentum distribution) to characterize the symmetry. By comparing the spin-resolved momentum distribution of the symmetry recovered state (obtained for times such that ) with a thermalized state (an equilibrium state with effective temperature), we conclude the symmetry recovered state is a thermalized state. This conclusion is further confirmed by computing the spin- and time-resolved density of states.

Further, we observe the order parameter undergoes a nearly exponential decay towards the symmetry recovered states. We studied the approximate decay rate and its relation to the initial temperature, Coulomb interaction strength, and the initial mass-imbalance ratio. These dependences are studied by varying one parameter while fixing the other two. We found the order parameter in the weak Coulomb interaction region exhibits a nearly quadratic dependence on , which can be interpreted with second order perturbation theory. For larger Coulomb interaction values, the deviation from quadratic dependence shows higher order terms must be taken into account.

We studied the dependence of approximate decay rate on the temperature. The decay rate increases rapidly with temperate in the low-temperature regime. By contrast, it saturates at higher temperatures (when the Coulomb interaction energy is overwhelmed relative to the kinetic energy). We studied the initial order parameter after the quench and found a critical temperature where the order parameter increases (decreases) for temperatures below (above) the critical temperature. The decay rate towards the thermalized state decreases as the initial imbalance ratio increases. Finally, we studied the dependence on the ramp shape and the pulse shape. Taken together, our results provide a guide to engineer the relaxation behavior of interacting, quantum many-particle systems.

## Acknowledgements

We acknowledge fruitful discussions with Qi Chen, Liang Dong, Quansheng Wu and Shengnan Zhang. We are grateful to Li Huang for a collaboration on a related project. L.D and G.A.F gratefully acknowledge funding from Army Research Office Grant No. W911NF- 14-1-0579, NSF Grant No. DMR-1507621, and NSF Materials Research Science and Engineering Center Grant No. DMR-1720595. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.

## References

- Lindner et al. (2011) N. H. Lindner, G. Refael, and V. Galitski, Nat. Phys. 7, 490 (2011).
- Gull et al. (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
- Wang et al. (2013) Y. H. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik, Science 342, 453 (2013).
- Aoki et al. (2014) H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Rev. Mod. Phys. 86, 779 (2014).
- Manmana et al. (2007) S. R. Manmana, S. Wessel, R. M. Noack, and A. Muramatsu, Phys. Rev. Lett. 98, 210405 (2007).
- Rigol et al. (2007) M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. Lett. 98, 050405 (2007).
- Kollath et al. (2007) C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007).
- Mentink and Eckstein (2014) J. H. Mentink and M. Eckstein, Phys. Rev. Lett. 113, 057201 (2014).
- Mentink et al. (2015) J. H. Mentink, K. Balzer, and M. Eckstein, Nat. Commun. 6, 6708 (2015).
- Eckstein and Werner (2016) M. Eckstein and P. Werner, Sci. Rep. 6, 21235 (2016).
- Bukov et al. (2016) M. Bukov, M. Kolodrubetz, and A. Polkovnikov, Phys. Rev. Lett. 116, 125301 (2016).
- Werner et al. (2012) P. Werner, N. Tsuji, and M. Eckstein, Phys. Rev. B 86, 205101 (2012).
- Tsuji et al. (2013) N. Tsuji, M. Eckstein, and P. Werner, Phys. Rev. Lett. 110, 136404 (2013).
- Tsuji and Werner (2013) N. Tsuji and P. Werner, Phys. Rev. B 88, 165115 (2013).
- Oka and Aoki (2009) T. Oka and H. Aoki, Phys. Rev. B 79, 081406 (2009).
- Du et al. (2017a) L. Du, X. Zhou, and G. A. Fiete, Phys. Rev. B 95, 035136 (2017a).
- Du and Fiete (2017) L. Du and G. A. Fiete, Phys. Rev. B 95, 235309 (2017).
- Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- Wolf et al. (2014) F. A. Wolf, I. P. McCulloch, and U. Schollwöck, Phys. Rev. B 90, 235131 (2014).
- Balzer et al. (2015) K. Balzer, Z. Li, O. Vendrell, and M. Eckstein, Phys. Rev. B 91, 045136 (2015).
- Cohen et al. (2015) G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, Phys. Rev. Lett. 115, 266802 (2015).
- Dong et al. (2017) Q. Dong, I. Krivenko, J. Kleinhenz, A. E. Antipov, G. Cohen, and E. Gull, Phys. Rev. B 96, 155126 (2017), arXiv:1706.02975 .
- Eckstein and Werner (2010) M. Eckstein and P. Werner, Phys. Rev. B 82, 115115 (2010).
- Freericks et al. (2006) J. K. Freericks, V. M. Turkowski, and V. Zlatić, Phys. Rev. Lett. 97, 266408 (2006).
- Eckstein and Kollar (2008) M. Eckstein and M. Kollar, Phys. Rev. Lett. 100, 120404 (2008).
- Freericks and Zlatić (2003) J. Freericks and V. Zlatić, Rev. Mod. Phys. 75, 1333 (2003).
- Eckstein et al. (2009) M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. Lett. 103, 056403 (2009).
- Eckstein et al. (2010) M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. B 81, 115131 (2010).
- Freericks (2008) J. Freericks, Phys. Rev. B 77, 075109 (2008).
- Eckstein and Werner (2011) M. Eckstein and P. Werner, Phys. Rev. Lett. 107, 186406 (2011).
- Fotso et al. (2015) H. Fotso, K. Mikelsons, and J. K. Freericks, Sci. Rep. 4, 4699 (2015).
- Tsuji et al. (2011) N. Tsuji, T. Oka, P. Werner, and H. Aoki, Phys. Rev. Lett. 106, 236401 (2011).
- Nguyen et al. (2015) D.-B. Nguyen, D.-K. Phung, V.-N. Phan, and M.-T. Tran, Phys. Rev. B 91, 115140 (2015).
- Dao et al. (2012) T.-L. Dao, M. Ferrero, P. S. Cornaglia, and M. Capone, Phys. Rev. A 85, 013606 (2012).
- Liu and Wang (2015) Y.-H. Liu and L. Wang, Phys. Rev. B 92, 235129 (2015).
- Philipp et al. (2017) M.-T. Philipp, M. Wallerberger, P. Gunacker, and K. Held, Eur. Phys. J. B 90, 114 (2017).
- Sekania et al. (2017) M. Sekania, D. Baeriswyl, L. Jibuti, and G. I. Japaridze, Phys. Rev. B 96, 035116 (2017).
- Du et al. (2017b) L. Du, L. Huang, and G. A. Fiete, Phys. Rev. B 96, 165151 (2017b).
- Eckstein and Kollar (2010) M. Eckstein and M. Kollar, New J. Phys. 12, 055012 (2010).