Controllable Finite-Momenta Dynamical Quasicondensation in the Periodically Driven One-Dimensional Fermi-Hubbard Model

TL;DR

This paper investigates how periodic driving influences the formation and controllability of doublon quasicondensates in a one-dimensional Fermi-Hubbard model, revealing mechanisms to engineer non-equilibrium states in cold-atom and solid-state systems.

Contribution

It demonstrates that periodic driving can control the momentum of doublon quasicondensation, breaking SU(2) symmetry and enabling engineered non-equilibrium condensates.

Findings

Doublons dynamically quasicondense at band edges.

Periodic drive breaks eta-SU(2) symmetry.

Driving amplitude controls condensate momentum.

Abstract

In the strongly interacting limit of the Hubbard model localized double-occupancies form effective hard-core bosonic excitations, called a doublons, which are long-lived due to energy conservation. Using time-dependent density-matrix renormalisation group we investigate numerically the dynamics of doublons arising from the sudden expansion of a spatially confined band-insulating state in one spatial dimension. By analysing the occupation scaling of the natural orbitals within the many-body state, we show that doublons dynamically quasicondense at the band edges, consistent with the spontaneous emergence of an eta-quasicondensate. Building on this, we study the effect of periodically driving the system during the expansion. Floquet analysis reveals that doublon-hopping and doublon-repulsion are strongly renormalised by the drive, breaking the eta-SU(2) symmetry of the Hubbard model. Numerical simulation of the driven expansion dynamics demonstrate that the momentum in which doublons quasicondense can be controlled by the driving amplitude. These results point to new pathways for engineering non-equilibrium condensates in fermionic cold-atom experiments and are potentially relevant to driven solid-state systems.

Equations

$$\begin{array}{rl}{H}_{\text{kin}}& =\sum _{i=1}^{L-1}\sum _{\sigma }\left(c_{i\sigma }^{†}{c}_{i+1,\sigma }+\text{h.c.}\right),\\ H_{\text{int}}& =\sum _{i=1}^L{n}_{i,\uparrow }{n}_{i,\downarrow },\end{array}$$

$$|{\psi }_{\text{init}}⟩=\left(\prod _{i\in O}{c}_{i\uparrow }^{†}{c}_{i\downarrow }^{†}\right)|\text{vac}⟩,$$

$$H_{\text{full}}(\tau )=H_{\text{hub}}+\theta (-\tau )H_{\text{con}}+\theta (\tau )H_{\text{drv}}(\tau ).$$

$$H_{\text{drv}}(\tau )=\frac{A}{2}\cos (\Omega \tau )\sum _i(-1{)}^i{n}_i,$$

$$\Gamma \sim Ct\exp [-\alpha (U/t)\log (U/t)],$$

$$H_{\text{eff}}=\text{mathbbm}P\left[\frac{J_0}{2}\left(H_{\text{hop}}+H_{\text{rep}}\right)\right]\text{mathbbm}P,$$

$$\displaystyle\begin{split}H_{\text{hop}}&=\sum_{i=1}^{L-1}\left(d^{\dagger}_{i}d_{i+1}+{\rm h.c.}\right),\\ H_{\text{rep}}&=\sum_{i=1}^{L-1}\Big{(}n_{d,i}(1-n_{d,i+1})+{\rm h.c.}\Big{)}.\end{split}$$

$$⟨{\eta }_N|d_i^{†}{d}_j|{\eta }_N⟩=\frac{N(L-N)}{L(L-1)}{e}^{\mathrm{i}\pi (i-j)},$$

$$\mathcal{U}({\tau }_2,{\tau }_2)=e^{-\mathrm{i}K({\tau }_2)}{e}^{-\mathrm{i}({\tau }_2-{\tau }_1)H_F}{e}^{\mathrm{i}K({\tau }_1)}.$$

$$\displaystyle\begin{split}\tilde{H}_{\text{eff}}(\nu,\Omega)=\mathbbm{P}\left[\frac{J(\nu,\Omega)}{2}\Big{(}H_{\text{hop}}+\Delta\left(\nu,\Omega\right)H_{\text{rep}}\Big{)}\right]\mathbbm{P},\end{split}$$

$$J(\nu ,\Omega )=J_0\sum _{m=-\infty }^{\infty }(-1{)}^m\frac{{\mathcal{J}}_m(\nu {)}^2}{1+m\Omega /U},$$

$$\Delta (\nu ,\Omega )=\frac{J_0}{J(\nu ,\Omega )}\sum _{m=-\infty }^{\infty }\frac{{\mathcal{J}}_m(\nu {)}^2}{1+m\Omega /U},$$

$${\rho }_{ij}(\tau )=⟨\varphi (\tau )|d_i^{†}{d}_j|\varphi (\tau )⟩,$$

$$n_q(\tau )=\frac{1}{L}\sum _{ij}{e}^{-\mathrm{i}q\left(i-j\right)a}{\rho }_{ij}(\tau ),$$

$$q_{\mathrm{c}}a\sim \pm \arccos (-\Delta )=\pm \arccos (\Delta )+\pi .$$

$${\rho }_{ij}(\tau )\to {\rho }_{ij}(\tau )e^{\mathrm{i}(i-j)\arccos (\Delta )},$$

$$⟨d_i^{†}{d}_j{⟩}_{\text{melt}}=⟨d_i^{†}{d}_j{⟩}_{\text{gs}}{e}^{-\mathrm{i}\theta (i-j)},$$

$${\rho }_{ij}\propto \frac{e^{\mathrm{i}{q}_{\mathrm{c}}(i-j)a}}{\sqrt{|i-j|}},$$

$$S_E(i)=-\sum _{\alpha }{\Lambda }_{\alpha }^{[i]}\log \left({\Lambda }_{\alpha }^{[i]}\right),$$

$$\delta D(\tau )=8{\left(\frac{t}{U}\right)}^2{\sin }^2\left(\frac{1}{2}U\tau \right).$$

$$F|{\varphi }_n(\tau )⟩=\left[H(\tau )-\mathrm{i}{\partial }_{\tau }\right]|{\varphi }_n(\tau )⟩={\epsilon }_n|{\varphi }_n(\tau )⟩,$$

$$⟨\text{leavevmode}⟨f,g⟩\text{leavevmode}⟩=\frac{1}{T}{\int }_0^{T}d\tau f^{\ast }(\tau )g(\tau ).$$

$$⟨\text{leavevmode}⟨m|⟨v|F|u⟩|n⟩\text{leavevmode}⟩=⟨\text{leavevmode}⟨m|⟨v|H(\tau )|u⟩|n⟩\text{leavevmode}⟩+n\Omega {\delta }_{n,m}.$$

$$R(\tau )=$$

$$=$$

$$\bar{H}(\tau )=\mathrm{i}\dot{R}(\tau )R^{†}(\tau )+R(\tau )H(\tau )R^{†}(\tau ).$$

$$\begin{array}{rl}R(\tau )& H_{\text{kin}}{R}^{†}(\tau )\\ =& -t\sum _{i\sigma }\exp \left[\mathrm{i}\frac{A}{\Omega }\sin (\Omega \tau ){\left(-1\right)}^i\right]c_{i\sigma }^{†}{c}_{i+1\sigma }+\text{h.c.}\\ =& -t\sum _{i\sigma m}{\mathcal{J}}_m(\nu )\exp \left[\mathrm{i}m\Omega \tau (-1{)}^i)\right]c_{i\sigma }^{†}{c}_{i+1\sigma }+\text{h.c.},\end{array}$$

$$\text{mathbbm}{P}_D=$$

$$\text{mathbbm}{M}_n=$$

$${\bar{F}}^{\prime }=\bar{F}/U={\bar{F}}_0^{\prime }+\lambda {\bar{F}}_1^{\prime },$$

$${\bar{F}}_1^{\prime }=-\sum _{i\sigma mn}{\mathcal{J}}_m(\nu )c_{i\sigma }^{†}{c}_{i+1\sigma }\otimes |n+m(-1{)}^i⟩\text{leavevmode}⟩⟨\text{leavevmode}⟨n|+\text{h.c.}.$$

Full text

I Introduction

Strongly correlated quantum systems are well known to exhibit a wide variety of novel phenomena like antiferromagnetism, the fractional quantum Hall effect and high-$T_{c}$ superconductivity. If such systems are driven out of equilibrium the emerging physics is expected to be richer still. So far only a small portion of this phenomenological landscape has been explored experimentally, and even less is understood theoretically. As a result the non-equilibrium dynamics of quantum many-body systems one of the most challenging branches of modern physics. Yet it is attracting growing attention due to the spectacular experimental advances in numerous complex quantum systems, ranging from cold-atom Greiner et al. (2002); Bloch et al. (2008), photonic Christodoulides et al. (2003); Rechtsman et al. (2013); Noh et al. (2017), optomechanical Ludwig and Marquardt (2013) and condensed matter platforms Fausti et al. (2011); Stojchevska et al. (2014); Wang et al. (2013). In particular the intense interest stems from the ability to implement controllable strong perturbations to a system and subsequently measure its properties in real-time with a resolution commensurate with the intrinsic microscopic time scales Aoki et al. (2014). This capability has opened up new spectroscopies for probing non-equilibrium dynamics as well as new approaches for manipulating them Giannetti et al. (2016a).

One exciting example of this has been the enormous progress over the past decade in ultra-fast THz pump-probe experiments on solid-state systems Orenstein (2012); Nicoletti and Cavalleri (2016); Giannetti et al. (2016b). By strongly driving low-energy structural or electronic degrees of freedom of a solid Rini et al. (2007) the ultra-fast melting of equilibrium long-ranged order, such as charge-density waves Fausti et al. (2011); Schmitt et al. (2008); Miyano et al. (1997); Cavalleri et al. (2001); Perfetti et al. (2008), magnetic order Stanciu et al. (2007); Ehrke et al. (2011), and orbital order Först et al. (2011) has been demonstrated. Even more remarkably, recent experiments have also used strong external modulations to induce superconducting order far from equilibrium in several different materials Hu et al. (2014); Kaiser et al. (2014a); Mankowsky et al. (2014); Mitrano et al. (2016). The long-term goal of this approach is ultimately to design and control quantum materials properties “on demand” by using driving to stabilize order that is otherwise inaccessible thermally Basov et al. (2017); Mankowsky et al. (2016). From a theoretical perspective these experiments raise important fundamental questions about what mechanisms exist for the emergence of order in driven systems, some of which have been explored in a number of recent studies Denny et al. (2015); Sentef et al. (2016); Knap et al. (2016); Subedi et al. (2014); Coulthard et al. (2017); Kennes et al. (2017).

Complementary to real materials, there have been equally spectacular experiments with systems of ultra-cold atomic gases in optical lattices Hofstetter and Qin (2018); Jördens et al. (2008); Schneider et al. (2008). These ‘synthetic’ solids provide near ideal quantum simulations of Hubbard-like Hamiltonians. Consequently they offer unique perspectives on non-equilibrium dynamics of interacting systems, owing to the unprecedented tunability of their time-dependent hopping amplitudes and inter-particle interactions, as well as the ability to engineer novel initial states Bloch et al. (2008). Cold-atom systems have thus opened up many long-standing non-equilibrium problems to exquisite experimental scrutiny, such as quenching across quantum phase transitions Greiner et al. (2002), controlling exchange interactions Anderlini et al. (2007); Trotzky et al. (2008), transport effects with strong interactions Brantut et al. (2012); Stadler et al. (2012); Schneider et al. (2012); Lebrat et al. (2018) as well as the influence of integrability Kinoshita et al. (2006); Langen et al. (2015) and many-body localization Schreiber et al. (2015) on closed system thermalization Eisert et al. (2015). Of particular relevance to our work are recent experiments Görg et al. (2018); Messer et al. (2018) where the effect of periodic driving on the basic interactions in fermionic cold-atom systems have been unravelled.

Motivated by all these developments, here we focus on a particularly intriguing example of the spontaneous emergence of order, namely dynamical quasicondensation, predicted Rigol and Muramatsu (2004) and observed experimentally in cold-atom systems Vidmar et al. (2015). Our aim is to assess whether a similar effect occurs with fermions, with broader implications for electronic systems. Dynamical quasicondensation manifests from a unit-filled Mott insulating region of strongly interacting bosons confined to the centre of a one-dimensional optical lattice with spacing $a$. Upon quenching the confinement this inhomogeneous initial state expands out into the surrounding empty lattice. It is found that the momentum distribution of the hard-core bosons quickly develops sharp peaks at quasimomenta $q_{\rm c}=\pm\pi/2a$ and power-law decaying spatial correlations, signalling unconventional current-carrying quasicondensation Rigol and Muramatsu (2004, 2005a). This unusual phenomenon has since been explained by showing the time-evolved state in this case is always an eigenstate of a time-dependent emergent Hamiltonian, nontrivially related to the underlying Hamiltonian governing the system Vidmar et al. (2017).

In this work we examine whether dynamical quasicondensation occurs with the expansion of a spatially confined band insulating state in the fermionic Hubbard model. In this case the initial state is a cluster of double occupations (doublons), which for strong repulsive interactions is highly energetic. However, when this energy far exceeds the single particle bandwidth energy conservation demands that the decay of doublons occur through multi-particle scattering processes that are exponentially suppressed with increasing energy Sensarma et al. (2010). The stability of such repulsively bound pairs has been confirmed experimentally Winkler et al. (2006); Strohmaier et al. (2010); Sensarma et al. (2010) and their distillation dynamics Heidrich-Meisner et al. (2009) has recently been observed Xia et al. (2015); Scherg et al. (2018). Since doublons are long-lived bosonic quasiparticles the question of whether they Bose (quasi)-condense has been addressed. It was found that they do condense at the band-edge $q_{\rm c}=\pm\pi/a$, leading to adiabatic proposals Rosch et al. (2008); Kantian et al. (2010) for generating the much sought-after $\eta$-condensate Yang (1989). Using time-dependent density matrix renormalization group (td-DMRG) methods Vidal (2003); White and Feiguin (2004); Daley et al. (2004); Al-Assam et al. (2017); Schollwöck (2011) we demonstrate here that the sudden expansion of the band insulator also undergoes dynamical $\eta$-quasicondensation with definitive signatures emerging within 10’s of hopping times.

We significantly expand the relevance of quench induced dynamical quasicondensation by examining its interplay with a simultaneously applied strong periodic driving. The effect of time periodic external fields can be captured by Floquet theory Shirley (1965); Eckardt and Anisimovas (2015); Bukov et al. (2015) and has been successfully used to predict and explain wide ranging phenomena in condensed matter and cold-atoms systems. Seminal examples include induced topological effects for cold-atoms via optical lattice modulation Jotzu et al. (2014); Struck et al. (2012) or exposure to circularly polarised light Grushin et al. (2014); Kitagawa et al. (2011) in solids, dynamical localization induced Mott transitions Dunlap and Kenkre (1986); Eckardt et al. (2005), effective repulsive to attractive interaction conversion by band-flipping Tsuji et al. (2011), and the renormalization of the super-exchange interaction Mentink et al. (2015); Coulthard et al. (2017); Jose et al. (2017). Here we use Floquet theory to show that the driving breaks the $\eta$-SU(2) symmetry of the Hubbard model Kitamura and Aoki (2016) and that by tuning above resonance a wide range of doublon dynamics is realizable, including one where they are non-interacting and directly mimic hard-core bosons Rigol and Muramatsu (2004). We compare this effective theory to td-DMRG numerical calculations to show that even finite frequency driving applied on short time scales can accurately control the momentum at which quasicondensation occurs.

This paper is organised as follows. In Sec. (II), we introduce the driven Hubbard model and initial state that is the focus of this work. Sec. (III) is devoted to deriving and analysing a simpler effective theory. This begins in Sec. (III.1) and Sec. (III.2) by employ a combination of a strong coupling expansion and Floquet theory to reduce the full driven problem into a quench of an effective doublon Hamiltonian with a driving amplitude dependent anisotropy. Building on this Sec. (III.3) and Sec. (III.4) then describe td-DMRG calculations that solve the quench dynamics of this effective Hamiltonian for large systems. Sec. (IV) then returns to analyse numerically the full driven Hubbard model by confirming in Sec. (IV.1) and Sec. (IV.2) that the behaviour seen in the effective model also manifests when the interactions are moderate and when a finite frequency drive is applied abruptly. In Sec. (IV.3) we examine the realisation of the full driven Hubbard model with current cold-atom experiments and analyse how signatures of driving controlled dynamical quasicondensation in the momentum distribution will appear in real-time measurements. Finally we conclude in Sec. (V).

II Model and setup

In this work we investigate dynamics of the fermionic Hubbard model in one dimension (1D) with $L$ sites and open boundaries. The Hamiltonian is $H_{\text{hub}}=-tH_{\text{kin}}+UH_{\text{int}}$ composed of single-band kinetic and on-site interacting contributions (taking $\hbar=1$ throughout)

[TABLE]

with the hopping amplitude $t$, and the on-site repulsion $U\geq 0$. Here, $c^{\dagger}_{i,\sigma}$ operators create a spin-$\sigma=\{\uparrow,\downarrow\}$ fermion localized on lattice site $i$, with the corresponding density for spin-$\sigma$ on the site being $n_{i,\sigma}=c^{\dagger}_{i,\sigma}c_{i,\sigma}$ and $n_{i}=n_{i,\uparrow}+n_{i,\downarrow}$.

We will consider the regime of strong interactions $U\gg t$ large enough that strongly correlated effects become readily apparent both in and out of equilibrium Scherg et al. (2018). We take the initial state of the system to be the ground state of $H_{\text{hub}}+H_{\text{con}}$, where $H_{\text{con}}=\sum_{i=1}^{L}v_{i}n_{i}$ is a box confinement potential, with $v_{i}=0$ for sites $i\in O$ inside a contiguous patch of the chain $O$ of size $N=|O|$, and $v_{i}\gg U$ otherwise. We assume the system has a total of $N_{\uparrow}=N_{\downarrow}=N$ electrons making the region $O$ doubly-filled and the ground state

[TABLE]

which is a spatially confined band insulator, a portion of which is depicted in Fig. 1(a).

We consider the dynamics of this system in real time $\tau$ simultaneously subject to two different time-dependent perturbations which together give a Hamiltonian

[TABLE]

The first perturbation describes the sudden switch off at $\tau=0$ of the confining potential, resulting in a quench whose expansion dynamics melts $\left|\,\psi_{\text{init}}\right\rangle$, as depicted in as depicted in Fig. 1(b). The second perturbation describes the abrupt switch on at $\tau=0$ of an external periodic driving

[TABLE]

describing an oscillating on-site energy alternating between sub-lattices with a frequency of $\Omega$ and amplitude $A$, as illustrated in Fig. 1(c). We will be considering high frequency cases where $\Omega\gtrsim U$ is the largest energy scale in the system. The driving term $H_{\text{drv}}(\tau)$ can be realizable in cold-atom experiments by laterally modulating in time the position of a zig-zag optical lattice configuration Zhang and Jo (2015).

We calculate the time-evolution of this system using td-DMRG Vidal (2003); White and Feiguin (2004); Daley et al. (2004); Al-Assam et al. (2017); Schollwöck (2011) with a time-step $\delta\tau<T/50$ where $T=2\pi/\Omega$ is the drive period, and a matrix product bond dimension $\chi>1000$ sufficient to ensure the discarded weight $\epsilon_{\rm disc}<10^{-6}$. In our numerical calculations we take the region $O$ to be the right half of a chain with $L\geq 2N$ giving an initial state of the form $\left|\,\psi_{\text{init}}\right\rangle=\left|\,\Updownarrow\Updownarrow\Updownarrow...\Updownarrow 0...000\right\rangle$. This asymmetric setup, where particles can only expand in one direction, allows access to longer timescales compared to the more experimentally motivated case where $O$ is located at the centre of the system that is overall twice as large. Once $L$ is large enough that no reflections occur at either open boundary in the simulated time we find the same results as those of the symmetric setup Heidrich-Meisner et al. (2008).

III Quenched effective model

Accurately computing the dynamics of the full driven Hubbard model $H_{\text{full}}(\tau)$ on long time-scales is extremely challenging. For this reason, and to give a deeper understanding of the physics, we begin our analysis by deriving an equivalent quench problem for a simpler effective model, applicable in the regime of the strong interactions $U\gg t$ and high-frequency $\Omega\gtrsim U$ drives.

III.1 Undriven system

For the Hubbard model $H_{\text{hub}}$ in the strongly interacting limit, doublons are known to be repulsively bound long-lived excitations, owing to their binding energy far exceeding the single-particle bandwidth Winkler et al. (2006); Strohmaier et al. (2010). The decay of doublons within the Hubbard model, in the presence of background holes, is dominated by multi-particle processes that generate many particle-hole pairs. A diagrammatic perturbative argument Sensarma et al. (2010) indicates that the rate of doublon decay $\Gamma$ is exponentially suppressed with $U/t$ as

[TABLE]

where the constants are $\alpha\approx 0.82$ and $C\approx 1.6$. This motivates examining the expansion dynamics of $\left|\,\psi_{\text{init}}\right\rangle$ using an effective model that explicitly conserves doublons.

Given $\left|\,\psi_{\text{init}}\right\rangle$ contains a maximally localized doublon domain for its filling, it is a highly excited state of $H_{\text{hub}}$ with an energy $E\approx UN$. This places it predominantly in the highest well-isolated band of Hubbard eigenstates when $U\gg t$. An effective Hamiltonian $H_{\text{eff}}$ describing just this highest set of eigenstates is derived in a standard perturbative approach accounting for virtual transitions to lower eigenstates. This is entirely analogous to the well-known derivation of the isotropic antiferromagnetic Heisenberg spin model describing the lowest-lying eigenstates of the half-filled Hubbard model Anderson (1959). The result reveals that the system is governed to second order in $t/U$ (up to a constant) by

[TABLE]

where $\mathbbm{P}$ is the projector onto real-space configurations containing no singly occupied sites and

[TABLE]

Here the allowable empty $\left|\,0\right\rangle$ or doubly occupied $\left|\,\Updownarrow\right\rangle$ local charge states are equivalent to the absence or presence of a hard-core boson described by the doublon creation operator $d^{\dagger}_{i}=c^{\dagger}_{i\downarrow}c^{\dagger}_{i\uparrow}$ obeying $(d^{\dagger}_{i})^{2}=0$ and associated number operator $n_{d,i}=d^{\dagger}_{i}d_{i}$. For $U>0$, the doublon hopping term $H_{\text{hop}}$ gives a single-doublon band with its minimum at the zone boundary $q=\pm\pi/a$, while the interaction term $H_{\text{rep}}$ gives repulsion between neighbouring doublons and holes. In Eq. (6) these terms have an identical coupling given by super-exchange interaction $J_{0}=4t^{2}/U$.

The isotropy between $H_{\text{hop}}$ and $H_{\text{int}}$ makes the ground state of $H_{\text{eff}}$ with $N$ doublons a so-called $\eta$-pair state $\left|\,\eta_{N}\right\rangle\sim(\eta^{+})^{N}\left|\,\textrm{vac}\right\rangle$, where $\eta^{+}\sim\sum_{i}(-1)^{i}d^{\dagger}_{i}$ creates a doublon at $q=\pm\pi/a$ momentum. This is a direct manifestation in $H_{\text{eff}}$ of the celebrated $\eta$-SU(2) symmetry of the Hubbard model, $\left[H_{\text{hub}},\eta^{+}\right]=U{\eta^{+}}$, first introduced by Yang Yang (1989). Consequently, $\left|\,\eta_{N}\right\rangle$ is an exact eigenstate of $H_{\text{hub}}$ with an energy $E=UN$ for any $U$. Furthermore, since

[TABLE]

the $\eta$-pair state displays staggered off-diagonal long-range order consistent with $N$ doublons Bose condensing at the zone edge $q_{\rm c}=\pm\pi/a$. Various proposals for generating in the Hubbard model an $\eta$-condensate, or states with $\eta$-like correlations, have been put forward. These include adiabatic switching of an optical lattice confinement and superlattice potentials Kantian et al. (2010), flipping of the band structure induced by driving the attractive Hubbard model Kitamura and Aoki (2016), and as an eigenstate of a “dark” Hamiltonian created from a Hubbard model with spin dephasing Buča et al. (2019).

III.2 Generalizing to driven system

For $\tau\geq 0$ the driven Hubbard model in Eq. (3) is time-periodic $H_{\text{full}}(\tau)=H_{\text{full}}(\tau+T)$. Analogous to Bloch’s theorem for discrete spatial translational symmetry, discrete time translational symmetry constrains the solutions of the time-dependent Schrödinger equation via Floquet’s theorem Eckardt and Anisimovas (2015). The key result is that the time-evolution operator $\mathcal{U}(\tau_{2},\tau_{1})$ of the driven system between two times $\tau_{1}$ and $\tau_{2}$ can be decomposed into a product of three unitaries

[TABLE]

Here $H_{F}$ is the time-independent Floquet effective Hamiltonian which generates continuous evolution between $\tau_{1}$ and $\tau_{2}$. This is sandwiched between two unitaries generated by the hermitian kick operator $K(\tau)$, with parametric dependence on the application time $\tau$ such that $K(\tau)=K(\tau+T)$. The kicks describe micromotion within each drive period. We will discuss its impact on the observables in Sec. (IV.3). In particular we find that the initial kick has no effect on the density or momentum distributions, meaning the effective Hamiltonian $H_{F}$ is sufficient to fully characterise the dynamics at times $\tau$ that are integer multiples of $T$, namely the stroboscopic time evolution of the driven system.

In general the computation of $H_{F}$ is a highly non-trivial problem. As detailed in Appendix A, we use an approach based on diagonalizing the so-called Floquet quasi-energy operator $F$ in an extended Floquet-Hilbert space after transforming to the frame rotating with the driving via the unitary $R(\tau)=\exp[-\mathrm{i}\int d\tau^{\prime}H_{\text{drv}}(\tau^{\prime})]$. Specifically, for the driven Hubbard model, $F$ is composed of diagonal blocks that are replicas of the Hubbard Hamiltonian $H_{\text{hub}}$, each shifted in energy by $\Omega$ relative to its neighbouring block, and with $H_{\text{kin}}$ renormalized by the drive by $\mathcal{J}_{0}(\nu)$, where $\nu=A/\Omega$ and $\mathcal{J}_{m}(\nu)$ is the $m$th Bessel function of the first kind. The off-diagonal couplings between blocks $m$ and $m^{\prime}$ is given by $H_{\text{kin}}$ renormalized by $\mathcal{J}_{m^{\prime}-m}(\nu)$, corresponding to fermion hopping accompanied by an exchange of $m^{\prime}-m$ quanta with the drive.

The diagonalization of $F$ is accomplished approximately by again employing standard degenerate perturbation theory in which we isolate simultaneously the band of highest-energy eigenstates of $H_{\text{hub}}$ within one block and account for corrections due to virtual transitions to eigenstates both in this block and all others. To second order in $t/U$ we obtain a driving dependent effective Hamiltonian

[TABLE]

where the modified super-exchange and anisotropy are

[TABLE]

obtained by summing over all blocks $m$.

The $(-1)^{m}$ factor in $J(\nu,\Omega)$ arises because the hopping of a doublon consists of two single particle hops in the same direction, as in Fig. 2(a). If the first hop is Fourier shifted by $m\Omega$ with an amplitude $\mathcal{J}_{m}(\nu)$ then the second hop returning to the same initial band would be a Fourier shift of $-m\Omega$ with an amplitude $\mathcal{J}_{-m}(\nu)=(-1)^{m}\mathcal{J}_{m}(\nu)$. In contrast the doublon-hole repulsion strength arises from the hopping of a single fermion forwards and back, with the same amplitude $\mathcal{J}_{m}(\nu)$, as in Fig. 2(b). Consequently, the effective model has a driving induced anisotropy that breaks the $\eta$-SU(2) symmetry of the undriven model Kitamura and Aoki (2016); Nocera et al. (2017). The reason for this is that the effective model describes charge degrees of freedom and the driving couples directly to charge. Had we considered instead the more conventional lowest-lying eigenstates of the half-filled Hubbard model, where the effective model describes spin degrees of freedom, no such anisotropy would be induced since the driving does not break the spin SU(2) symmetry.

In the undriven limit we have $\lim_{\nu\rightarrow 0}~{}J(\nu,\Omega)=J_{0}$ and $\lim_{\nu\rightarrow 0}~{}\Delta(\nu,\Omega)=1$, recovering Eq. (6). For any finite $\nu$ in the high-frequency driving limit, $\Omega\gg U\gg t$, only the $m=0$ contribution to Eq. (11) and Eq. (12) survives giving $\lim_{\Omega/U\rightarrow\infty}~{}J(\nu,\Omega)=J_{0}\mathcal{J}_{0}(\nu)^{2}$ and $\lim_{\Omega/U\rightarrow\infty}~{}\Delta(\nu,\Omega)=1$, so isotropy is preserved with a renormalized super-exchange. Thus $\eta$-SU(2) symmetry breaking arises when $\Omega/U$ is finite. Here we will study driving frequencies $\Omega\gtrsim U$, while avoiding direct resonant couplings between Hubbard excitations to first order in $t$ Bukov et al. (2016).

In Fig. 3 we plot $J(\nu,\Omega)$ and $\Delta(\nu,\Omega)$ for $U=20t$ for two moderate driving frequencies. In both cases the anisotropy $\Delta(\nu,\Omega)$ is suppressed for moderate non-zero amplitudes $\nu$. The doublon-hole repulsion can even be completely removed at $\nu_{0}$, as highlighted in Fig. 3, where the doublons then behave as non-interacting hard-core bosons. For both $\Omega$’s stronger driving can induce doublon-hole attraction $\Delta(\nu,\Omega)<0$.

III.3 Doublon domain melting

We have established that in the regime $U\gg t$ and $\Omega\gtrsim U$ the stroboscopic dynamics of the periodically driven Hubbard model introduced in Sec. (II) reduces to an effective model. Specifically, the driven dynamics is approximated as a sudden quench at $\tau\geq 0$ of interacting hard-core bosons governed by the effective model in Eq. (10) specified by $J(\nu,\Omega)$ and $\Delta(\nu,\Omega)$ with an initial state $\left|\,11\dots 1100\dots 00\right\rangle$. Analysing this effective model brings several benefits. First, it is computationally much simpler than the full Hubbard model, allowing much larger system sizes to be accessed numerically. Second, it is isomorphic to magnetic domain-wall melting in the XXZ spin model, so insight can be gleaned from the extensive studies on this spin model Sabetta and Misguich (2013); Collura et al. (2018); Misguich et al. (2017) .

For $0\leq\Delta\leq 1$ the quench dynamics of the initially sharp domain generically displays an expansion of bosons outwards from the boundary of the domain, and correspondingly holes inwards into the domain. For a selection of $\Delta$’s, Fig. 4 shows the evolution of the doublon density profile $\langle n_{d,i}(\tau)\rangle$ up to a time $\tau=24/J$. Centred on the domain boundary we see that a melt lobe forms with a size that grows linearly in time for $0\leq\Delta<1$, grows sub-linearly for $\Delta=1$ and stops growing after finite time for $\Delta>1$. For $0\leq\Delta\leq 1$ the melt profile displays a universal form. Formally this is revealed by rescaling the site coordinates in some appropriate way $\bar{x}(\tau)$ such that the melt profile at all times collapses onto itself. To illustrate this Fig. 5(a) shows the melt profile for each site $i$ for times between $\tau=5/J$ and $\tau=24/J$ with $\Delta=1$. Owing to the anomalous super-diffusive transport properties at the isotropic point a power-law rescaling $\bar{x}(\tau)=a(i-N)/(J\tau)^{p}$ with $p=3/5$ Misguich et al. (2017) is found to induce the profile collapse, as shown in Fig. 5(b). For all other interaction strengths $0\leq\Delta<1$, a linear rescaling $\bar{x}(\tau)=a(i-N)/v(\Delta)\tau$ with a $\Delta$ dependent speed $v(\Delta)=J\sqrt{1-\Delta^{2}}$ establishes universality, as depicted in Fig. 5(c)-(f). For $\Delta=0$ the density profile in the universal region is known Antal et al. (1999) to be $\langle n_{d,i}(\tau)\rangle=\arccos(\bar{x}(\tau))/\pi$.

Here our central focus is on the properties of the system, beyond the density profile, within the melt lobe universality region $-1\leq\bar{x}(\tau)\leq 1$ as a function of the interaction $0\leq\Delta\leq 1$. In particular, as first discovered by Rigol and Muramatsu Rigol and Muramatsu (2004), for quenches in the non-interacting $\Delta=0$ limit dynamical quasicondensation occurs. This novel effect is only revealed by examining the full doublon one-particle density matrix (OPDM)

[TABLE]

computed from the time-evolving state $\left|\,\phi(\tau)\right\rangle$ of the effective model. In Fig. 6 a colormap of the OPDM for some representative interaction strengths are shown for three time slices. The expanding black square signifies the universal region. For $\Delta=0$ the universal region expands at the maximum speed $J$ of the single-doublon band and so spans the light-cone of the system dynamics. Even with interactions there remains small contributions to the OPDM throughout the light-cone due a tiny fraction of doublons escaping the domain at a speed $J$ from short-time perturbative adjustments to the sudden quench. However, for $\Delta>0$ the expansion speed $v(\Delta)$ defining the melt lobe becomes increasingly slow in comparison. As we shall show this separation of speeds makes the analysis of the interacting system numerically challenging because the full light-cone must still be captured such that $L\sim N\sim J\tau$ to avoid spurious boundary effects.

By exploiting the equivalence of the $\Delta=0$ limit to non-interacting spinless fermions Rigol and Muramatsu Rigol and Muramatsu (2004) solved numerically exactly the dynamics of large systems over long times. Surprisingly, they found that the evolution of the momentum distribution of the hard-core bosons

[TABLE]

quickly changes from a featureless constant, characterising the localized state at $\tau=0$, to a distribution with a pronounced peak at finite-momentum $q_{\rm c}a=\pi/2$, indicative of a condensate forming dynamically. This result is reproduced here in Fig. 7(a) using td-DMRG.

Our analysis of the interacting system similarly begins by examining the momentum distribution. A simple energetic argument anticipates that a peak in the momentum will shift towards $qa=\pi$ as $\Delta$ approaches the isotropic point, and that no peak will develop for $\Delta>1$. Initially, the sharp domain wall has an average energy $\langle H_{\text{eff}}\rangle=J\Delta/2$ coming exclusively from the doublon-hole repulsion at the interface $...1\ 1\ 1\ 1\ 0\ 0\ 0\ 0\ ...$. If we now consider a doublon at the interface unbinding and propagating into the empty region as $...1\ 1\ 1\ 0\ 0\ 1\ 0\ 0\ ...$ then the interfacial interaction energy remains unchanged, but the isolated doublon now contributes an additional interaction energy of $J\Delta$. If this melting is to occur conservation of energy demands that the increased interaction is compensated by a reduction in the isolated doublon’s kinetic energy. Given that the single-doublon dispersion is $\epsilon_{\text{hop}}(k)=J\cos(qa)$ this implies that the propagating doublon will possess a momentum

[TABLE]

Assuming the melt lobe is dilute in the long-time limit we thus anticipate the accumulation of melting doublons into this momentum state. The relation Eq. (15) is in agreement with $q_{\rm c}a=\pi/2$ for $\Delta=0$, and also reveals that for $\Delta>1$ the finite bandwidth is insufficient to compensate the interaction precluding condensation at any $q_{\rm c}$.

To confirm these expectations we have computed the momentum distribution from the interacting OPDM. Crucially, to isolate the contributions from the melting, Eq. (14) was not applied directly. Instead, we first restricted the Fourier transform of the OPDM to the universal melt lobe region defining a subsystem centred on the domain boundary with a time-dependent size $\ell(\tau)=2v(\Delta)\tau/a$ for $0\leq\Delta<1$, or $\ell(\tau)=2(J\tau)^{p}/a$ for $\Delta=1$ 111To smooth the integer jumps in the melt lobe size we performed weighted averages between adjacent sizes.. Next, to allow a clean comparison to Eq. (15) given the limited resolution in momentum space available for small systems, we applied a phase-shift to the OPDM as

[TABLE]

before taking the Fourier transform. The purpose of this was to shift any contribution at $q_{\rm c}$ to precisely $q=\pi/a$, which is guaranteed to coincide with a discrete momenta available since $\ell(\tau)$ is always even. After taking the Fourier transform the momentum grid was then shifted back. The resulting $n_{q}$ are shown in Fig. 7 for several interaction strengths and times. The number of doublons in the melt lobe is $N(\tau)=\ell(\tau)/2$.

The evolution of the momentum distributions in Fig. 7(b)-(d) reveal an intriguing parallel to the non-interacting behaviour in Fig. 7(a) – as time progresses we see the emergence of an increasingly sharp peak at finite momentum close to $q_{\rm c}$ predicted by Eq. (15). These peaks are a smoking gun of dynamical quasi-condensation, and indicate that this effect is not simply an anomalous feature of the non-interacting system. Also plotted are the momentum distributions of the half-filled ground state calculated with DMRG using the same interaction strength $\Delta$ and system size $\ell(\tau)$ as the corresponding melt state. All ground states display a peak at $qa=\pi$ that becomes sharper with increasing $\Delta$ and whose magnitude grows with $N$ sub-linearly for $0\leq\Delta<1,$ consistent with quasi-condensation, and linearly for $\Delta=1$ consistent with $\eta$-condensation, as shown in Fig. 7(e).

The behaviour of the melt state peak at $q_{\rm c}$ in Fig. 7(f) displays the same trends. While the ground state peak growth with $N$ settles quickly into it a universal asymptotic form the melt dynamics with interactions lags behind and the form of its growth form has yet to fully stabilise. Nonetheless for $\Delta<0.5$ there is good correspondence between the melt state and ground state peak growth. For $\Delta>0.5$, even though there are severe melt lobe size (or melt time) limitations, the growth appears suppressed. For $\Delta=1$ the peak in the melt state distribution shown in Fig. 7(d) does not occur at exactly at $qa=\pi$, but is instead located at the discrete momentum state directly adjacent to it, suggesting it may approach $qa=\pi$ only in the thermodynamic limit. For $\Delta>1$ (not shown) a small peak emerges, caused by the sharp initial state transients, whose magnitude rapidly saturates with $N(\tau)$.

For $\Delta=0$ it is known that the melt state asymptotically converges to a boosted ground state Vidmar et al. (2017), apparent here already for relatively short times with the similarity of the peaks in Fig. 7(a). A natural question is whether melting gives a boosted ground state more generally in the interacting case. An affirmative answer was obtained in previous work which focused on a initial “soft” domain wall Sabetta and Misguich (2013); Lancaster and Mitra (2010). In this case the initial state was the ground state of a weak confinement potential $v_{i}\propto\tanh(\beta i)$, where the constant $\beta$ controls the wall width. As a result the initial state has a small difference in the density $\delta\ll 0.5$ far to the left $\langle n_{d,i\ll 0}\rangle\sim 0.5+\delta$ and far to the right $\langle n_{d,i\gg 0}\rangle\sim 0.5-\delta$ of the interface. Using both numerical Sabetta and Misguich (2013) and analytic hydrodynamic arguments Lancaster and Mitra (2010) it was shown that the correlations of the melt state in the long-time limit generated around the domain boundary (far from the edges of the chain) have a simple relation to the corresponding half-filled ground state as

[TABLE]

where $\theta=2\delta\arccos(\Delta)$. For the setup we consider here, where the ground states displays a peak at $qa=\pi$ and a “hard” domain wall is used so $\delta=0.5$, we see that this phase relation agrees with Eq. (15) used above. Yet there are some differences in the behaviour of the melt states and ground states shown in Fig. 7. To better assess these difference we now move on to compute other properties of the melt state expected of a quasicondensate.

III.4 Further signatures of quasicondensation

While the emergence of a sharp peak in the momentum distribution is suggestive of a condensation phenomenon a more careful examination of the scaling properties of the OPDM is needed to fully identify the nature of the melt state. Following the analysis for $\Delta=0$ performed by Rigol and Muramatsu Rigol and Muramatsu (2004) we diagonalized the OPDM to find the natural orbitals and their occupations as a function of time after the quench. The emergence of quasicondensation is revealed by the behaviour of the so-called lowest natural orbital (LNO) $\phi_{0}$ with the largest occupation $\lambda_{0}$. For simple unfragmented dynamical quasicondensation the LNO should display, after transients have subsided, a time-invariant universal form with the rescaled of coordinate $\bar{x}(\tau)$. In Fig. 8 the rescaled LNO’s for a variety of $\Delta$’s are shown for times $15\leq J\tau\leq 25$. Identically to the density profile we find that within the melt lobe $|\bar{x}(\tau)|\leq 1$ the form of LNO for all $\Delta$’s is effectively constant in time confirming its universality. Decaying contributions outside the melt lobe are seen in all cases, but are especially prominent for larger interactions due to the slower expansion speed.

The next crucial property for quasicondensation is the scaling of the LNO occupation $\lambda_{0}$ with the number of doublons $N(\tau)$ in the melt lobe. True condensation occurs if $\lambda_{0}(N)=\alpha N$, implying a finite O(1) fraction of particles $\sigma(N)=\lambda_{0}(N)/N$ occupy the LNO in the thermodynamic limit. Quasicondensation is instead a sub-linear power law growth of $\lambda_{0}(N)\propto N^{p}$ with $0<p<1$ so the occupancy fraction vanishes as $\sigma(N)\sim 1/N^{1-p}$. Analogous to our analysis of the momentum distribution we determined LNO’s by diagonalizing the OPDM only inside the melt lobe region $-1<\bar{x}(\tau)<1$ for each $\Delta$. The behaviour of $\sigma(N)$ for the corresponding half-filled ground state for the same $\Delta$ and particle number $N(\tau)$ is shown in Fig. 9(a). For $0\leq\Delta<1$ this shows a power-law decrease in $\sigma(N)$ with $N$, but the decay slowing down with increasing $\Delta$. At the point $\Delta=1$ the LNO density saturates to a constant $\sigma(N)\rightarrow 0.5$ signalling $\eta$-condensation.

In Fig. 9(b) we show $\sigma(N)$ for the melt state. For $\Delta=0$ a power-law decay identical to that of the ground state is seen. For interacting systems $\Delta>0$ the basic trend is similar to the ground state with the decay in $\sigma(N)$ slower than $\Delta=0$ and softening with increasing $\Delta$. This is indicative of stronger dynamical quasicondensation with increasing interactions, consistent with the sharpening momentum distributions with $\Delta$ observed already in Fig. 7(a)-(d). However, for all interacting cases the behaviour of $\sigma(N)$ displays more discernible differences to the ground state than was apparent in the momentum distribution alone. None of the decay curves have reliably converged to a power-law form, even for the weakest interactions, and there are even signs of saturation. However, due to the limitation in reachable $N(\tau)$’s, which is most severe as $\Delta\rightarrow 1$, it is currently inconclusive whether genuine dynamical condensation is occurring.

Another important property of quasicondensation is the decay of long-ranged correlations off-diagonal correlations in the OPDM within the spatial support of the LNO. For the non-interacting case a distinctive power-law decay is found

[TABLE]

with a phase difference of $q_{\rm c}a=\pi/2$ between neighbouring sites signalling quasicondensation at finite momentum Rigol and Muramatsu (2005b). The observed $1/2$ exponent for the algebraic decay of correlations is identical to the ground state of the same non-interacting system, as expected from Eq. (17). In Fig. 10 the absolute value of doublon correlations $|\rho_{N,j}|$ from the domain wall boundary at a time $J\tau=30$ are shown for various $\Delta$’s. Here for all non-zero interactions we see even more visible differences from the power-law decay of the corresponding ground state that are also shown. For $\Delta\leq 0.5$ in Fig. 10(a)-(c) the melt state displays stronger correlations over the melt lobe than the ground state, while in contrast for $\Delta>0.5$ in Fig. 10(d)-(f) they are increasingly weaker. None of the interacting cases considered over the times accessible in our calculations agree with the soft wall hydrodynamic prediction 222We cannot preclude the possibility that calculations on much longer timescales, not currently feasible with td-DMRG for this hard domain wall setup, might recover the behaviour given in Eq. (17) as the interface softens. in Eq. (17).

Overall, we have found distinctive signatures of dynamical quasicondensation within the effective model of doublons arsing from the driven Hubbard model. However, due to the separation of expansion speeds the system sizes reached in our calculations are insufficient to fully quantitatively diagnose the nature of this quasicondensation.

IV Driven Hubbard Model

Having observed that a form of quasicondensation emerges within the effective model we now return to the full driven Hubbard model in Eq. (3). In particular we will now demonstrate that this novel effect is robust beyond the strongly interacting and high-frequency approximations underlying the validity of the effective model. Moreover, we will show that it continues to occur for large but finite interactions $U>t$ and for finite driving frequencies $\Omega>U$ not too close to resonance.

IV.1 Zero driving case

As a baseline we consider first the undriven Hubbard model, corresponding to $\Delta=1$ in the effective model, and examine agreement as $U/t$ is decreased. A readily accessible indicator within td-DMRG of the increased complexity of the full Hubbard model is the entanglement entropy of the time-evolved state $\left|\,\psi(\tau)\right\rangle$

[TABLE]

where $\Lambda^{[i]}_{\alpha}$ are the squared Schmidt coefficients of $\left|\,\psi(\tau)\right\rangle$ for a bipartition of the system between sites $i$ and $i+1$. In Fig. 11(a) $S_{E}(i)$ is plotted for various $U/t$’s for sites $i$ close to the domain wall after a time $J\tau=10$ with $J=2t^{2}/U$. For this short time the slow expansion for $\Delta=1$ gives a melt lobe extending around $\sqrt{10}\sim 3$ sites either side of the wall interface. As expected the strongest interaction $U/t=20$ closely follows the $S_{E}(i)$ of the effective model up to the speed $J$ light-cone of 10 sites, and so captures the small contribution of the ballistically escaping doublons. However, the full Hubbard model solution also has a non-negligible entropy beyond the effective model’s light-cone arising from the partial disassociation of doublons into fast-moving fermions with a speed $2t>J$. The full Hubbard model thus has yet another speed associated to the quench dynamics. As expected this fermion contribution becomes more pronounced with decreasing $U/t$ as doublons become less stable.

Despite having a finite $U/t$ the decay of doublons quickly saturates with time. This is confirmed by computing the deviation in the total number of doublon number $\delta D(\tau)=\langle\hat{D}(0)\rangle-\langle\hat{D}(\tau)\rangle$ with $\hat{D}=\sum_{i}\hat{n}_{d,i}$ in the time-evolved state. Second order time-dependent perturbation theory predicts this observable will behave as

[TABLE]

In Fig. 11(b) we plot $\delta D(\tau)$ for several $U/t$’s showing that it saturates to a constant given by the time-average of Eq. (20) and is thus suppressed as $(t/U)^{2}$.

Given the bounded fraction of fermions generated by finite $U/t$ we next examined the key characteristics of quasicondensation in the full Hubbard model. In Fig. 11(c) the LNO occupation of the doublon OPDM computed from the full Hubbard model is compared to the effective model with the appropriate superexchange $J=2t^{2}/U$. We find good agreement for the LNO growth with time, even for very moderate interaction strengths $U/t=8$. In Fig. 11(d) we verify that a peak in the momentum distribution of the melt-lobe continues to manifest close to $q_{\rm c}a=\pi$ in agreement with the effective model at $\Delta=1$. These results together confirm that the quasicondensation seen in the effective model is extremely robust to the presence of a small amount of initial doublon disassociation. The reason for this is that the fermions rapidly expand beyond the melt lobe region leaving it essentially undisturbed. Such a distillation of doublons and fermions was recently demonstrated experimentally Scherg et al. (2018).

IV.2 Driven case

We now address how well the effective model quench describes the full Fermi-Hubbard model dynamics with an abrupt application of finite frequency driving. Since the Floquet effective model is formally only a stroboscopic description of the driven Fermi-Hubbard model we compare the models in this section at times $\tau$ that are integer multiples of the driving period $T$. This ensures that both the kick operators $K(\tau)$ and transformation $R(\tau)$ from the lab to the rotating frame are equal to identity. We focus on two representative cases with the driving amplitude $\nu$ tuned to a moderate strength, so $\Delta(\nu,\Omega)=0.5$, and also stronger so $\Delta(\nu,\Omega)=0$. In Fig. 12(a)-(d) we plot the LNO occupation in time of the driven Hubbard model for these two $\Delta(\nu,\Omega)$’s for range of $U$’s and two different driving frequencies $\Omega/U=1.5$ and $\Omega/U=1.8$. We observe good agreement with the growth predicted by the effective model over the times examines, despite the finite interactions and driving frequency. As expected the agreement improves when increasing $U/t$ and/or $\Omega/U$. In Fig. 12(e)-(h) we plot the corresponding $\delta D$. For the driven Hubbard model we observe an order of magnitude increase in the doublon loss compared to the undriven case in Fig. 11(a). Furthermore the doublon loss does not saturate and instead displays a roughly exponential increase in time with a rate constant that is suppressed with increasing $\Omega$. This behaviour is a remnant of Floquet heating. The doublon losses also explain the more noticeable lag in the LNO population growth for the lowest values of $U/t$, e.g. seen in Fig. 12(a) and (c).

Given the exponential doublon losses we estimate that in the worse case considered, namely $U/t=8$, $\Omega/U=1.5$ and $\Delta(\nu,\Omega)=0$ in Fig. 12(g), that the quasicondensate predicted by the effective model will be completely depleted by $t\tau\sim 50$. For more favourable parameters, such as $U/t=20$ and $\Omega/U=1.8$ in Fig. 12(h), this time is significantly extended to $t\tau\sim 200$. In this case, examining $t\tau\approx 42$ below this time where doublon loss is negligible, Fig. 13(a) shows that the doublon momentum distribution displays a sharp peak positioned at $q_{\rm c}a=\pi/2$. Remarkably this is very close what is expected from the free expansion of non-interacting hard core bosons in the effective model (also shown), yet is realised here in a driven, strongly interacting fermionic Hubbard model with a moderate drive frequency $\Omega\gtrsim U$. In Fig. 13(b) at weaker driving where $\Delta(\nu,\Omega)=0.5$ we find a peak at $q_{\rm c}a>\pi/2$. Together with the zero driving peak in Fig. 11(d) approaching $q_{\rm c}a=\pi$, we see that a controllable range of quasicondensing momenta are indeed accessible in the driven Hubbard model.

IV.3 Cold-atom implementation

The driven dynamical quasicondensation outlined is realizable in a standard optical lattice experiment with feasible lattice parameters Köhl et al. (2005); Tarruell and Sanchez-Palencia (2018); Scherg et al. (2018); Esslinger (2010). To illustrate this we consider the well studied case of fermionic K40 in an undriven 3D optical lattice potential. One dimensional systems are realized by choosing anisotropic depths $V_{x}=8E_{R}$ along the axis of the chains and $V_{y}=V_{z}=33E_{R}$ for the transverse confinement, where $E_{R}=\hbar^{2}k^{2}_{L}/2m_{a}$ is the recoil energy with $\lambda_{L}=2\pi/k_{L}=532$nm being the laser wavelength and $m_{a}$ is the atomic mass. For this system the nearest-neighbour hopping amplitude of the Hubbard model is $t\sim 0.54$kHz, with small next-nearest-neighbour $t^{\prime}\sim t/50$. Given typical fermionic cold-atoms systems can remain coherent for up to 1s this in principle allows for a total experiment time of $t\tau\sim 500$. This is consistent with a central domain of order 100 hundred sites surrounded by similarly sized empty regions. The band-insulating initial state can be generated using additional magnetic trapping and tuning the system to have attractive interactions via a Feshbach resonance.

The implementation of the alternating potential driving term included in the Hamiltonian Eq. (4) is slightly non-standard, but is nonetheless realisable with state of the art optical lattice experiments. Conventionally an optical lattice is driven by shaking the entire trapping potential by a length $\delta(\tau)=\delta_{0}\cos(\Omega\tau)$ via kHz piezoelectric modulation of lattice laser’s phase. Shaking along the $x$-axis then induces an effective linear potential in the non-inertial reference frame $V(x)=-\Omega^{2}\delta_{0}x$. To create an alternating driving potential we propose a scheme based on a zig-zag chain geometry Zhang and Jo (2015); Anisimovas et al. (2016) where odd and even sites are laterally displaced. Shaking perpendicular to the chain but parallel to the zig-zag then induces an effective modulated potential difference between odd and even sites. More details and parameter estimations for this setup are discussed in Appendix 16.

Neither time-of-flight nor in-situ measurements of cold-atom will implement perfect stroboscopic sampling of the driven system. As such we now examine how the momentum distribution of the driven Fermi-Hubbard model deviates from the effective model at times inside a single driving period. Specifically, we time-evolved the full driven Hubbard model in the lab frame given in Eq. (3) with $U/t=20$ until $\tau=80T\approx 16.7/t$ and then frequently measured the momentum distribution over the next driving period at small increments $\Delta t\approx 0.05T$. The contour plots of the resulting distributions for various driving strengths are displayed in Fig. 14. The contour plots begin at a commensurate time $\tau=80T$ so as expected a peak centred at the $q_{c}$ of the corresponding effective model’s $\Delta$ is seen. As time progresses within the drive period the momentum distribution oscillates at a frequency $2\Omega$ by exchanging weight from the primary peak to another secondary peak shifted by a momentum $\pi/a$. The relative amplitude of the secondary peak depends on the driving strength, starting at zero for the undriven $\Delta=1$ case and reaching unity for the strongly driven $\Delta=0$ case.

We confirmed that the transformation from the rotating frame $R(\tau)$ (given by Eq. (24) in Appendix A) is the origin of this oscillatory behaviour. In Fig. 15 the amplitude of the primary and secondary peaks for strong and moderate driving strengths from the full driven Hubbard model’s momentum distribution are directly compared to those of the effective model’s after transforming back to the lab frame. The excellent agreement between the two models indicates that the kick operators have a negligible effect for this observable. Furthermore, the appearance of two distinguishable momentum peaks for measurements at general times is a desirable experimental signature of the alternating driving scheme proposed here since it is robust to time-averaging.

V Conclusions

We have demonstrated that dynamical quasicondensation at a controllable finite momentum can emerge from a band-insulating domain initial state whose expansion is governed by a one-dimensional driven Hubbard model in the strongly-interacting regime. To establish this we first examined an Floquet effective model. We showed how the known dynamical quasicondensation of non-interacting $\Delta=0$ hard-core bosons at $q_{\rm c}a=\pi/2$ not only persists in the interacting systems with doublon-hole repulsion $0<\Delta\leq 1$, but does so with a stronger LNO growth and at a momentum that shifts towards $q_{\rm c}a=\pi$. However, our numerical calculations were not fully conclusive on the nature of this quasicondensation owing to the slowing down of the doublon propagation with increasing $\Delta$ which limited the system sizes accessible. We found significant differences between the melt state and a boosted ground state for the times examined, in contrast to hydrodynamical calculations for soft domain walls. We cannot rule out that these differences are transient and that such a correspondence may yet emerge at much longer times when the initial hard domain has been softened by the expansion. Next, we went beyond the effective model and simulated the full driven Hubbard model. We demonstrated good agreement with the predictions of the effective model for key observables over short times, even at very moderate finite interactions and driving frequencies. The presence of Floquet heating induced doublon loss, not captured by the effective model, was shown to be present but does not preclude quasicondensation for the times examined.

Given that all the ingredients of the setup proposed have been implemented in current state-of-the-art cold-atom setups the effects outlined are in principle within reach of experimental observation Bloch et al. (2008). Indeed a cold-atom quantum simulator might be the most definitive means of answering the open question as to whether dynamical quasicondensation of doublons also occurs in higher spatial dimensions. Moreover, our quench setup presents a potentially fast and robust scheme for creating a much sought-after $\eta$-quasicondensate with cold-atoms Rosch et al. (2008); Kantian et al. (2010).

As an outlook it is interesting to speculate whether dynamical quasicondensation can be relevant to experiments on driven solid-state systems. Several features make this plausible. First, owing to the beam spot-size and limited penetration depth of all pump-probe experiments so far, only a small excitation volume of sample is driven Giannetti et al. (2016b). Consequently the induced excited state is necessarily inhomogeneous and it is an open question whether the subsequent expansion dynamics of charge-carriers out of this excitation volume into the rest of the material plays a crucial role in the non-equilibrium superconducting coherence observed Kaiser et al. (2014b); Mitrano et al. (2016). Second, a sharper domain, more similar to that considered here, can be engineered in a solid by heterostructuring. For example, the time-resolved scrambling of magnetic order in a thin-film due to the expansion of interfacial shockwaves of mobile carriers from a substrate has been observed Först et al. (2015). The possibility of observing instead the dynamical ultra-fast emergence of order from a shockwave is intriguing.

VI Acknowledgements

MWC thanks the University of Bath for URSA funding. SRC gratefully acknowledges support from the UK’s Engineering and Physical Sciences Research Council (EPSRC) under grant No. EP/P025110/1.

Appendix A Floquet expansion

Floquet’s theorem dictates that any Hamiltonian with periodic time dependence will have solutions to the time dependent Schrödinger equation $H(\tau)\Psi(\tau)=\mathrm{i}\partial_{\tau}\Psi(\tau)$ of the form $\left|\,\psi_{n}(\tau)\right\rangle=\exp(-\mathrm{i}\varepsilon_{n}\tau)\left|\,\phi_{n}(\tau)\right\rangle$ where $\varepsilon_{n}$ is the quasi-energy associated with the $T$-periodic Floquet state $\left|\,\phi_{n}(\tau+T)\right\rangle=\left|\,\phi_{n}(\tau)\right\rangle$. In the basis $\{\left|\,\phi_{n}(\tau)\right\rangle\}$ the time dependent Schrödinger equation is recast as

[TABLE]

where $F=H(\tau)-\mathrm{i}\partial_{\tau}$ is the Floquet quasienergy operator.

The task of finding the quasi-energies $\varepsilon_{n}$ and the Floquet modes $\left|\,\phi_{n}(\tau)\right\rangle$ can be achieved by diagonalising $F$ over an extended Hilbert space $\mathcal{F}$. Specifically, $\mathcal{F}=\mathcal{H}\otimes\mathcal{T}$ is the tensor product space of the original Hilbert space $\mathcal{H}$ and the space $\mathcal{T}$ of functions $f(\tau)=f(\tau+T)$ with periodicity $T=\frac{2\pi}{\Omega}$. We denote states in $\mathcal{T}$ with a double angled bracket |f\mathclose{\hbox{{\rangle}}\kern-1.94444pt\leavevmode\hbox{{\rangle}}} and use a scalar product

[TABLE]

In this derivation it is convenient to choose the Fourier basis of time periodic functions |n\mathclose{\hbox{{\rangle}}\kern-1.94444pt\leavevmode\hbox{{\rangle}}}=\exp(\mathrm{i}n\Omega\tau) that are eigenstates of $\partial_{\tau}$ with \partial_{\tau}|n\mathclose{\hbox{{\rangle}}\kern-1.94444pt\leavevmode\hbox{{\rangle}}}=\mathrm{i}n\Omega|n\mathclose{\hbox{{\rangle}}\kern-1.94444pt\leavevmode\hbox{{\rangle}}} obeying \mathopen{\hbox{{\langle}}\kern-1.94444pt\leavevmode\hbox{{\langle}}}n,m\mathclose{\hbox{{\rangle}}\kern-1.94444pt\leavevmode\hbox{{\rangle}}}=\delta_{n,m}. The index $n\in\mathbbm{Z}$ is often called the “photon number” of a Floquet sector. The matrix elements of the quasi-energy operator $F$ using this basis for $\mathcal{T}$ and any basis for $\mathcal{H}$ are then

[TABLE]

Since \mathopen{\hbox{{\langle}}\kern-1.94444pt\leavevmode\hbox{{\langle}}}m|\left<v\right|H(\tau)\left|u\right>|n\mathclose{\hbox{{\rangle}}\kern-1.94444pt\leavevmode\hbox{{\rangle}}} only depends on the difference $(n-m)$, thus the quasi-energy operator $F$ contains replicas of $H$ for each photon number $n$, each shifted in energy by $n\Omega$. Despite the infinite duplication of the many-body Hilbert space, for very high-frequency driving, where $\Omega$ is larger than all matrix elements, standard degenerate perturbation theory can be used to derive an approximate description of a single-band due to corrections from the neighbouring bands.

Here we describe a strong coupling expansion onto both a single Floquet sector and onto a specific band of Hubbard eigenstates with a given doublon number. In the lab frame, the coupling between the Floquet sectors is only non zero for $n-m=\pm 1$, and the coupling strength is comparable to $\Omega$. This is an unsuitable starting point for perturbation theory. However, the coupling can be reduced by transforming into a frame rotating with the driving, where the time dependence becomes a Peierls phase with periodic time dependence. This is achieved with the unitary

[TABLE]

which transforms the Hamiltonian into the rotating frame as

[TABLE]

The first term removes the driving potential $H_{\text{drv}}(\tau)$, and while the second term puts a time-dependent momentum shift onto the kinetic energy term as

[TABLE]

where $\mathcal{J}_{m}(\cdot)$ is the Bessel function of the first kind and $\nu=A/\Omega$ is the dimensionless driving strength.

A perturbative expansion around the relevant set of eigenstates requires two projectors. First, $\mathbbm{P}_{n}$ projects onto the eigenstates with $D$ double occupations

[TABLE]

and $\mathbbm{M}_{n}$ projects onto the $n$ photon sector

[TABLE]

We define the small parameter $\lambda=t/U$ and let the frequency $\Omega$ be of comparable magnitude to the interaction strength $U\sim\Omega\gg t$. We then define the dimensionless Floquet quasi-energy operator in the rotating frame as

[TABLE]

where $\bar{F}^{\prime}_{0}=\sum_{D,n}\mathbbm{P}_{D}\otimes\mathbbm{M}_{n}\,E_{D,n}$ is the operator which counts doublon number and photon index, with a highly degenerate set of eigenvalues $E_{D,n}=D+\left(\Omega/U\right)n$. The perturbation $\bar{F}^{\prime}_{1}$ contains the kinetic coupling of eigenstates of different doublon and photon as

[TABLE]

We denote the perturbative expansion of the quasi-energy operator $\bar{F}^{\prime}$ about the set of eigenstates with exactly $D$ doublons and $n$ photons as $\bar{F}^{\prime}_{D,n}$. For the relevant case of the maximum doublon number $D=N$ and the DC Floquet sector $n=0$ we have

[TABLE]

To second order the only contributions in this effective model are from consecutive hops which break up and reform a doublon. This reformation can occur on the same site the doublon started on, or on a neighbouring site. Since the sign of the Floquet transitions depends on the direction of the hopping, the processes involving two hops in a single direction will pick up different amplitudes to processes which are two hops in opposite directions.

To compute these amplitudes we define two operators. The first is

[TABLE]

which creates a doublon on site $i$ when a fermion of spin $\sigma$ occupied site $j$ and a fermion of spin $\bar{\sigma}=-\sigma$ occupied site $i$. The second is

[TABLE]

which describes the hopping of singly occupied sites over either a doublon or a hole. Using these operators the term $c^{\dagger}_{i\sigma}c_{j\sigma}$ describing hopping of a spin $\sigma$ fermion from $j$ to $i$ is decomposed into

[TABLE]

The energy change to second order is $E_{N,0}-E_{N-1,b}=U-b\Omega$, and the numerator in the expansion of $\bar{F}^{\prime}_{N,0}$ is made of two parts. We automatically drop the $g$ and $h^{\dagger}$ contributions from $\bar{F}^{\prime}_{1}\left(\mathbbm{P}_{N}\otimes\mathbbm{M}_{0}\right)$, since there is no single fermion hop which preserves doublon number and no free fermions available to create a doublon. We similarily drop all $h$ and $g$ type from $\left(\mathbbm{P}_{N}\otimes\mathbbm{M}_{0}\right)\bar{F}^{\prime}_{1}$ for the same reason. Finally we remove the projectors from the middle of the expression since we are automatically guaranteed to step down by one doublon. The expression for $\bar{F}_{N,0}^{\prime}$ then becomes

[TABLE]

Multiplying out the two hopping processes gives

[TABLE]

The sum over lattice site $j$ can be removed by noting that unless the doublon is created by $h^{\dagger}$ on one of the same two sites as $h$, we will change doulbon number. Therefore the only allowed combinations are of the form $h^{\dagger}_{i,j}h_{i,j}$ or $h^{\dagger}_{j,i}h_{i,j}$. We also replace the inner products in $\mathcal{T}$ with the constraints \mathopen{\hbox{{\langle}}\kern-1.94444pt\leavevmode\hbox{{\langle}}}n|m\mathclose{\hbox{{\rangle}}\kern-1.94444pt\leavevmode\hbox{{\rangle}}}=\delta_{n,m} to give

[TABLE]

We omit the left and right projectors $\mathbbm{P}_{N}\otimes\mathbbm{M}_{0}$ from here on since their presence is implied. The three remaining constraints on the four indices denoting the Floquet sectors leaves one independent photon index to sum over as

[TABLE]

The two parts to the Hamiltonian can be interpreted as nearest-neighbour repulsion of a doublon and a hole,

[TABLE]

and a nearest-neighbour hopping of doublons in one direction and holes in the opposite direction,

[TABLE]

The two amplitudes given in Eq. (12) and Eq. (11) are recovered since $\mathcal{J}_{-m}(\nu)=(-1)^{m}\mathcal{J}_{m}(\nu)$.

Appendix B Implementing a zig-zag lattice

One possible way of creating a driving term of the form Eq. (4) is using a one dimensional chain with a lateral displacement of nearest-neighbour sites, i.e. a zig-zag type geometry. Proposals for creating a zig-zag type lattice (aswell as rhombic, sawtooth, etc.) are detailed in Ref. Zhang and Jo (2015). With a specific choices of lattice laser phases they derive a potential of the form

[TABLE]

which is capable of creating a “zig-zag” type lattice under certain parameter choices. A suitable example is as follows, we set $V_{\text{sq}}=0$ to ensure the onsite energy of each Wannier orbital is equal, while $V_{x}=-10E_{R}$, $V_{y}=V_{x}$, $V_{\bar{x}}=2V_{x}$. Here $E_{R}$ is the recoil energy, which for potassium-40 is $E_{R}\approx 17$kHz. An additional deep confinement potential in the $z$ direction of $-20E_{R}\cos^{2}(k_{L}z)$ is used to achieve 1D tubes. The resulting potential is shown in Fig. 16.

For this choice of parameters we have computed the dominant Hubbard Hamiltonian matrix elements using the Wannier MATLAB package as detailed in Ref. Walters et al. (2013). We find that the hopping amplitudes displayed in Fig. 17 are $t=0.03E_{R}\approx 0.6$kHz for nearest-neighbour hopping between A and B sites, $t^{\prime}=7\times 10^{-4}E_{R}\approx 0.01$kHz for next-nearest-neighbour hopping, and $t^{\prime\prime}=1.2\times 10^{-4}E_{R}\approx 0.002$kHz for the hopping across neighbouring one dimensional channels of between A and B sites. The on-site Hubbard interaction is $U\approx 48g/\lambda_{L}^{3}$ where $g/(h\lambda_{L}^{3})=4\pi a_{s}\hbar^{2}/m_{a}h\lambda_{L}^{3}\approx 0.1$kHz using the $s$-wave scattering cross section $a_{s}$ which is 118 Bohr radii for K40 Tarruell and Sanchez-Palencia (2018); Esslinger (2010). The interaction strength $U$ can be independently tuned via a Feshbach resonance. Nearest-neighbour density interactions between A and B sites in the same one dimensional channel are approximately $3$Hz. Finally, the average separation between the two Wannier orbitals is on the order of $\lambda/2$. To obtain a potential difference of $V_{0}/(\hbar\Omega)\approx 1$, sufficient for implementing the strongly driven non-interacting effective model with $\Omega\approx 2U/\hbar=40t/\hbar$, we have $V_{0}=-\Omega^{2}(\lambda_{L}/2)m_{a}\delta_{0}\rightarrow 1=|\Omega(\lambda_{L}/2)m_{a}\delta_{0}/\hbar|\rightarrow\delta_{0}\approx 2\hbar/(\Omega\lambda_{L}m_{a})\approx 200$nm, and so requires a modulation close to the lattice spacing.