Extended degenerate perturbation theory for the Floquet–Hilbert space

In the first step, we separate $\bar{Q}^{\prime}$ into the unperturbed part and the perturbation, the latter having a block-diagonal part and a block-off-diagonal one:

where

The dimensionless parameter $\lambda$ has been introduced to track the order of the expansion. Inserting $\bar{Q}^{\prime}$ into Eq. (3) and assuming $\bar{G}=\lambda\bar{G}^{(1)}+\lambda^{2}\bar{G}^{(2)}+\ldots$, one can collect the terms of the same order. In the first order, this yields

We require $\bar{G}^{(1)}$ be block-off-diagonal, so that $[\bar{G}^{(1)},\bar{Q}^{\prime(0)}]$ is block-off-diagonal as well. The remaining block-diagonal terms immediately yield the first-order term of the effective Hamiltonian according to

or, in the Hilbert space,

Next, we require $\bar{W}_{{\rm X}}^{(1)}=0$, obtaining the equation for the block-off-diagonal terms: $[\bar{G}^{(1)},\bar{Q}^{\prime(0)}]=\bar{V}_{{\rm X}}$. This allows us to find $\bar{G}^{(1)}$ and proceed to the equation for the second-order terms. Continuing in the similar fashion, one obtains

One recognizes that the resulting expressions are identical to the van Vleck high-frequency expansion [28], which is expected since $\bar{Q}^{\prime}$ possesses exactly the same structure as in the usual applications of this expansion. The Hilbert-space expressions, however, do different because the matrix elements of $\bar{V}_{{\rm X}}$ and the quantities $\varepsilon_{\alpha m}^{(0)}$ have a different meaning in our case. The second-order term of the effective Hamiltonian results as

Expressions for the third-order terms of the effective Hamiltonian are provided in Appendix A.

The formula (14) could alternatively be obtained by a direct application of the conventional degenerate perturbation theory (DPT) widely used in the context of ordinary Hilbert space [31], provided one assumes that all of the states $|\alpha m\rangle\!\rangle$ of a given block ($m=0$, say) constitute the degenerate subspace. The condition $n\neq m$ in the summation over the intermediate states $|\beta n\rangle\!\rangle$ in Eq. (14) corresponds precisely to the skipping of states belonging to the degenerate subspace. We will refer to the Eqs. (13) and (15) as the results of the extended degenerate perturbation theory (EDPT) since the degenerate subspace is extended to include all levels in a Floquet zone.

For comparison, let us consider the conventional DPT, whereby the degenerate subspace of ${\cal F}$ consists only of the states that are exactly (or nearly) degenerate. We start by diving the Hilbert space ${\cal H}$ into two subspaces, $\mathbb{D}^{0}$ and $\mathbb{D}^{1}$, so that $\mathbb{D}^{1}$ contains the states sharing the same value of reduced energy of interest, $\varepsilon_{*}^{(0)}$, while $\mathbb{D}^{0}$ contains the remaining states. That is, for each state $|\alpha\rangle$,

Consequently, we partition each diagonal block of $\bar{Q}^{\prime}$ into two subblocks — one containing the states belonging to $\mathbb{D}^{0}$, and one containing the states of $\mathbb{D}^{1}$. The couplings between these subblocks are then considered to constitute the off-diagonal blocks of $\bar{Q}^{\prime}$, and are therefore treated as a perturbation. This means that the definition of diagonal and off-diagonal blocks is changed from the one given in Eq. (10) to

Here, $|\alpha\rangle\in\mathbb{D}^{A}$ and $|\alpha^{\prime}\rangle\in\mathbb{D}^{A^{\prime}}$ so that the Kronecker delta $\delta_{A^{\prime}A}$ is unity if both states $|\alpha^{\prime}\rangle$ and $|\alpha\rangle$ belong to the same subspace — either degenerate or not — and zero if they belong to different subspaces. With these definitions, and $\langle\!\langle\alpha^{\prime}m^{\prime}|\bar{V}|\alpha m\rangle\!\rangle=\langle\alpha^{\prime}|\hat{V}_{a-a^{\prime}+m^{\prime}-m}|\alpha\rangle$ as before, one finds

in the first order and

in the second. In the last expression, $B$ is the subspace number which $|\beta\rangle$ belongs to, i.e., $|\beta\rangle\in\mathbb{D}^{B}$. The terms of the effective Hamiltonian constructed using DPT are denoted here by $\hat{w}^{(n)}$, and they are composed of two uncoupled blocks. In fact, when using DPT, one is only interested in the block describing the degenerate subspace, and it can be diagonalized separately from the other block.

Notably, construction of $\hat{w}_{{\rm D}}^{[2]}$ is not equivalent to simply neglecting in $\hat{W}_{{\rm D}}^{(2)}$ the couplings between the blocks $\mathbb{D}^{0}$ and $\mathbb{D}^{1}$. Considering two states $|\alpha^{\prime}\rangle$ and $|\alpha\rangle$ of the degenerate subspace, one has $\langle\alpha^{\prime}|\hat{w}_{{\rm D}}^{(2)}|\alpha\rangle\neq\langle\alpha^{\prime}|\hat{W}_{{\rm D}}^{(2)}|\alpha\rangle$. This point will be discussed further in Section III.2.

The above formalism also allows one to construct schemes that are intermediate between DPT and EDPT. Instead of including in $\mathbb{D}^{1}$ only the states whose reduced energies are equal to the reduced energy of interest $\varepsilon_{*}^{(0)}$, one can include all states in a certain interval $[\varepsilon_{*}^{(0)}-\Delta\varepsilon/2,\varepsilon_{*}^{(0)}+\Delta\varepsilon/2]$ (with $\Delta\varepsilon<\omega$). The size of $\mathbb{D}^{1}$ will then increase, resulting in larger numerical effort required to calculate the eigenvalues, but these should approximate the exact values more accurately. The choice $\Delta\varepsilon=\omega$ corresponds to EDPT, whereby all states are assigned to $\mathbb{D}^{1}$. We will not, however, consider the accuracy of these intermediate schemes further in this work, instead focusing on the two limiting cases, DPT and EDPT.

It follows that the necessary condition for the convergence of the expansion (3) is

which has to be satisfied for all states $|\alpha\rangle$, $|\alpha^{\prime}\rangle$ and all integers $n$, except those excluded in the summation in Eq. (15) (EDPT case) or Eq. (19) (DPT case). As long as the numerator does not vanish, this condition might be violated for $|\varepsilon_{\alpha^{\prime}}^{(0)}-\varepsilon_{\alpha}^{(0)}|\approx\omega$, which is possible if the reduced energies of the two states $|\alpha^{\prime}\rangle$ and $|\alpha\rangle$ are on the opposite boundaries of the FZ. A simple resolution is to shift the FZ in a such way so that there are no levels near one or both boundaries. However, if the reduced energies $\varepsilon_{\alpha}^{(0)}$ fill the FZ densely, then the problem cannot be circumvented. The EDPT will be applicable in those cases only if the couplings between states on the boundary of the FZ can be disregarded. On the other hand, in case of DPT this issue does not arise for the states of the degenerate subspace if one centers the FZ around $\varepsilon_{*}^{(0)}$. Since $|\alpha\rangle$ and $|\alpha^{\prime}\rangle$ enumerate only the states of the degenerate subspace, the absolute values of denominators in Eq. (19) are no smaller than $\omega/2$. However, the problem appears in the DPT case as well once the third-order term, provided in Appendix A, is included since it contains differences between the reduced energies of the nondegenerate subspace.

The driven Bose–Hubbard model is defined by the Hamiltonian [35]

where $J$, $U$, and $F$ control the strengths of, respectively, the nearest-neighbor hopping, the on-site interaction, and the external driving (we study monochromatic driving of frequency $\omega$). The first sum runs over nearest-neighbor pairs, while the remaining ones run over all lattice sites. In the last sum, $x_{j}$ is the $x$-coordinate of site $j$, in units of the lattice constant. A gauge transformation $\hat{H}(t)=\hat{U}^{\dagger}(t)\hat{H}^{\prime}(t)\hat{U}(t)-\mathrm{i}\hat{U}^{\dagger}(t)\mathrm{d}_{t}\hat{U}(t)$ with $\hat{U}(t)=\exp\!\left(-\mathrm{i}\frac{F}{\omega}\sin\omega t\sum_{j}x_{j}\hat{n}_{j}\right)$ shows that the effect of the driving amounts to a renormalization of the hopping strength [26]:

The Fourier image of the resulting Hamiltonian is given by

where ${\cal J}_{m}(x)$ denotes the Bessel function of the first kind of order $m$.

We will use the Fock basis $|\alpha\rangle$ to refer to the elements of $\hat{H}(t)$, denoting the diagonal ones by $E_{\alpha}=\langle\alpha|\hat{H}_{0}|\alpha\rangle$. An example of the distribution of diagonal elements $E_{\alpha}$ is displayed in Fig. 2(a). The presence of degenerate elements does not require extra care — these degeneracies will directly translate into degeneracies in the first-order term (14) of the effective Hamiltonian. The effective Hamiltonian, once obtained perturbatively, will be diagonalized exactly (numerically). On the other hand, the possibility of resonant transitions requires the application of degenerate perturbation theory in ${\cal F}$. Since the diagonal elements $E_{\alpha}$ are given by integer multiples of $U$, resonant transitions become possible when the condition $pU=q\omega$, where $\{p,q\}\in\mathbb{Z}$, is satisfied. This translates into degeneracies in ${\cal F}$, as depicted in Fig. 2(b) showing the diagonal elements $\varepsilon_{\alpha m}^{(0)}$ [see Eq. (6)] for $3U=2\omega$. In the figure, black horizontal lines delimit the FZ chosen as $[-\frac{\omega}{2},\frac{\omega}{2})$, and the values of $a$ are displayed.

The analysis of accuracy of DPT and EDPT will be performed by choosing the driving strength $F$ and the frequency $\omega$, and calculating the quasienergies as the parameter $U$ is varied. We will focus our attention on the quasienergy of the “driven Mott-insulator (MI) state”, a name we will use for the Floquet mode having the largest overlap with MI state (denoted as $|\text{MI}\rangle$). For $U\gg J$, $|\text{MI}\rangle$ is the ground state of the undriven system, corresponding to $E=0$. Therefore, we center the FZ around $\varepsilon_{*}^{(0)}=0$ by choosing the interval $[-\frac{\omega}{2},\frac{\omega}{2})$.

We note that the interesting features of the quasienergy spectra — the anticrossings indicatory of resonant processes — will appear at certain values of the ratio $U/\omega$. From the definition (5) with $E_{\alpha}=k_{\alpha}U$ where $k_{\alpha}$ is integer, it follows that the denominator in Eq. (20), $\omega|(\varepsilon_{\alpha^{\prime}}^{(0)}-\varepsilon_{\alpha}^{(0)})/\omega-n|$, grows linearly in $\omega$ for fixed $U/\omega$. Therefore, the accuracy of the perturbation theory in the vicinity of resonances is expected to increase with increasing $\omega$, although this reasoning does not take into account the dependence on $\omega$ of the coupling strength [i.e. the numerator in Eq. (20)].

Finally, let us clarify that DPT can be applied not only exactly on resonance (when the degenerate energies of interest exactly coincide), but also in its vicinity. For example, we can calculate the quasienergies for $U$ in the vicinity of $\frac{2}{3}\omega$ by including in $\mathbb{D}^{1}$ the same states as those constituting $\mathbb{D}^{1}$ when $U$ exactly equals $\frac{2}{3}\omega$. Similarly, the same reduction numbers [see Eq. (5)] as those obtained exactly on resonance are used even when $U$ does not exactly equal $\frac{2}{3}\omega$.

We begin the assessment of accuracy of the perturbation theories with the study of a BH system defined on a periodic $1\times 8$ lattice containing 8 bosons; we set $F/\omega=2$ and $\omega/J=20$. The plots in Fig. 3(a) display the quasienergies of ten Floquet modes having the largest overlap with the MI state for $U\in[0,1.5\omega]$. The quasienergies of the driven MI state are highlighted in red. The left plot displays the exact result, while the right one shows the results obtained using third-order EDPT. As explained above, the DPT approach is only applicable in the vicinity of resonances, and cannot be directly used to calculate the quasienergies for such a wide range of $U$. The most pronounced anticrossing seen at $U\approx\omega$ corresponds to first-order creation of particle–hole excitations in the MI state, whereby a particle is annihilated at a certain site and created at its neighboring site [32]. At $U=\omega$, transitions between all levels of the system become resonant, therefore, all diagonal elements of $\bar{Q}^{\prime}$ in the given diagonal block share the same value. Consequently, the actual distribution of quasienergies is almost entirely captured by the first-order effective Hamiltonian (13). Additionally, DPT becomes equivalent to EDPT since all levels of the system are included in the degenerate subspace. The second- and third-order terms of $\hat{W}_{{\rm D}}$ provide a slight improvement and yield results in agreement with the exact ones, as shown in Fig. 3(a). Notably, the EDPT results remain sufficiently accurate away from resonances as well. Exactly on the resonance $U=\omega$, the largest value of the coupling ratio (20) is $r_{\max}=0.12$, which is one of the largest values that can be considered to satisfy the condition (20). Therefore, for smaller values of $\omega$ (and the same value of $F/\omega$), the EDPT is not expected to yield reliable results near the resonance $U=\omega$.

To assess the accuracy of the methods on a finer scale, we inspect the quasienergy of the driven MI state in the vicinity of $U=\frac{2}{3}\omega$, as shown in Fig. 3(d). The obtained anticrossing is a result of the second-order process whereby two particles of an MI state hop to the site of their common neighbor, producing a triply occupied site and resulting in a state of energy $E=3U$. This process has been analyzed in Ref. [32] using the DPT approach. Plugging the expressions (23) into Eq. (19) one finds that for the resonance condition $3U=q\omega$, the element $c_{3}\equiv\langle 3|\hat{w}_{{\rm D}}^{(2)}|{\rm MI}\rangle$ vanishes for $q$ odd ($|3\rangle$ denotes any one of the states with a triply occupied site reachable starting from $|\text{MI}\rangle$ in two hops). Meanwhile, for $q=2$, the strongest coupling is observed for $F/\omega\approx 2$. The quasienergy of the driven MI state is indeed calculated correctly using DPT, as shown in Fig. 3(d). However, as $U$ is tuned away from resonance, the EDPT yields more accurate results.

As noted in Section II.2, the EDPT and DPT approaches are not equivalent, therefore, the value of $c_{3}$ depends on which method is used to construct the effective Hamiltonian. These values are compared in Fig. 3(c) for $U=\frac{2}{3}\omega$ and a range of coupling strengths $F/\omega$. The curves are quite different in nature: for example, at $F/\omega\approx 3.4$ the DPT curve crosses the zero, while the EDPT curve approaches a local maximum. Vanishing coupling indicates disappearance of the corresponding anticrossing in the quasienergy spectrum [32], which is indeed the case for $F/\omega\approx 3.4$, as confirmed by an exact calculation (not shown). Thus, even though EDPT yields quasienergies with higher accuracy than DPT, the matrix elements of $\hat{W}_{{\rm D}}$ constructed using EDPT do not have such a straightforward interpretation compared to the case when DPT is used. Another difference between EDPT and DPT concerns the higher-order terms of the effective Hamiltonian. In the DPT case, the third-order term $\hat{w}_{{\rm D}}^{(3)}$ gives an insignificant correction to the second-order theory at $U=\frac{2}{3}\omega$. The additional processes appearing in the third-order are the creation of a state featuring three doubly occupied sites (the corresponding matrix element is $\eta=9.2$ times smaller than $c_{3}$) and the process of creating a state with a triply occupied site in three hops ($\eta=25$). In the EDPT case, on the other hand, inclusion of the third-order term $\hat{W}_{{\rm D}}^{(3)}$ gives a substantial improvement in terms of accuracy.

It is instructive to consider the values of the coupling ratios (20) for EDPT. Exactly on resonance $U=\frac{2}{3}\omega$, we find $r_{\max}=0.30$, which is quite large. However, this value comes from the coupling of a state corresponding to $\varepsilon_{\alpha m}^{(0)}=-\frac{1}{3}\omega$ with a state corresponding to $\varepsilon_{\alpha^{\prime}m^{\prime}}^{(0)}=-\frac{2}{3}\omega$, therefore, this coupling does not directly influence the driven MI state, whose reduced energy is $\varepsilon_{*}^{(0)}=0$. The said coupling is indicated in red in Fig. 3(b) depicting the diagonal elements of $\bar{Q}^{\prime}$ at $U=\frac{2}{3}\omega$. On the other hand, all of the degenerate states sharing the value $\varepsilon_{*}^{(0)}=0$ are coupled weakly to states outside the central FZ, as indicated by the blue numbers in Fig. 3(b) (see figure caption for details). This explains the fact that accurate results have been obtained using EDPT despite the condition (20) not being satisfied. We remind that the couplings between the states in the central FZ [see green arrows with numbers in Fig. 3(b)] are taken into account exactly in the EDPT framework. Meanwhile, the couplings of degenerate states with the nondegenerate ones are treated perturbatively in the DPT, which explains why the reported DPT results are less accurate even exactly on resonance.

Let us also discuss the origin of the discontinuity at $U=\frac{5}{8}\omega$ in the EDPT results seen in Fig. 3(d). For definiteness, consider the group of levels characterized by $\varepsilon_{\alpha}^{(0)}=-\frac{1}{3}\omega$ and reduction number $a=3$, appearing in red in Fig. 3(b), where $U=\frac{2}{3}\omega$ is assumed. These are the levels arising from the unperturbed states of energy $E=4U$. As $U$ decreases, these levels shift downwards, reaching the lower FZ boundary when $U$ attains the value of $\frac{5}{8}\omega$. Decreasing $U$ still further makes this particular group of levels leave the FZ, and another one (appearing in green in the figure) enters from above. However, the reduction number for the levels of the latter group is $a=2$. Since the reduction number directly influences the strength of coupling with other levels, the coupling with the states under consideration changes abruptly as $U$ crosses the value of $\frac{5}{8}\omega$. Although processes such as this one reduce the accuracy of the EDPT, their impact is not that noticeable if the states of actual interest are not strongly coupled to the states crossing the boundaries of FZ. This is certainly the case presently: Our main focus is on the driven MI state, which is not directly coupled to the $4U$ states. Consequently, the artifactual discontinuity in Fig. 3(d) is an order of magnitude narrower than the width of the anticrossing that is of actual interest.

The last case for this parameter set is studied in Fig. 3(e) that displays the quasienergy of the driven MI state in the vicinity of $U=\frac{4}{3}\omega$. The anticrossing is again a manifestation of the second-order process producing a triply occupied site. Despite the width of the anticrossing making up only $\sim 2\%$ of the width of the FZ, it can be calculated accurately using perturbation theory. Again, EDPT is more accurate than DPT, although the former one yields a number of false anticrossings.

We now turn to a BH system described by the same parameter values ($F/\omega=2$, $\omega/J=20$), but this time on a $2\times 4$ lattice periodic in the $x$-direction (we direct the $x$-axis along the longer dimension of the lattice). The quasienergy spectrum calculated for $U\in[0,1.5\omega]$ is shown in Fig. 4(a). We notice that the curve of the driven MI state undergoes many more anticrossings compared to the above case of a one-dimensional lattice, indicating richer system dynamics. The EDPT results are less accurate here, featuring erroneous and noticeable discontinuities. Nevertheless, the quasienergies can be considered calculated qualitatively correctly in the vicinity of the main anticrossing at $U=\omega$.

According to Eq. (22), driving renormalizes the hopping strength only for hopping along the $x$-axis. As a result, the DPT theory no longer predicts that $c_{3}$ vanishes for all $F$ when $U=\frac{m}{3}\omega$ with $m$ odd. Indeed, Fig. 4(b) confirms that numerous anticrossings appear in the vicinity of $U=\frac{1}{3}\omega$. Their presence is predicted by both EDPT and DPT, albeit the exact results are matched only qualitatively. Meanwhile, in the cases $U=\frac{2}{3}\omega$ and $U=\frac{4}{3}\omega$ the EDPT ensures higher accuracy, providing a considerable improvement over DPT [see Figs. 4(c) and (d)].

Finally, we consider the results obtained for the case of a higher driving frequency: $\omega/J=30$. As expected, the accuracy of the perturbation theory is higher in this case, which is confirmed by the plots in Fig. 5. Figure 5(a) shows that EDPT yields accurate results on a large scale, capturing all the essential features of the exact quasienergy spectrum. The close-up views of the spectrum in the vicinity of $U=\frac{1}{3}\omega$, $U=\frac{2}{3}\omega$ and $U=\frac{4}{3}\omega$, shown in Figs. 5(b), (c), and (d), respectively, display that EDPT is capable of providing quantitatively correct results in this regime. The DPT is certainly applicable as well, although the accuracy is lower.

Concluding the discussion of accuracy of the considered methods, let us compare the associated computational costs. Calculation of quasienergies using both DPT and EDPT comes down to a Hermitian matrix diagonalization. In the EDPT case, one is required to diagonalize the matrix whose size is given by the total number of states of the unperturbed systems. In the studied eight-particle systems, there are $N=6435$ such states. In the DPT case, the size is given by the number of states sharing the same value of reduced energy $\varepsilon_{*}^{(0)}$. For example, on a $2\times 4$ lattice and at $U=\frac{2}{3}\omega$, there are 2017 states sharing the value of $\varepsilon_{*}^{(0)}=0$ [cf. Fig. 3(b)]. Both methods thus suffer from exponential growth of the required computation resources, which can only be alleviated by limiting the number of states taken into consideration.

The numerically exact calculation of the quasienergies via the single-period evolution operator $\hat{P}$ requires performing a (unitary) $N\times N$ matrix diagonalization, similarly to the the EDPT case. However, the construction of $\hat{P}$ additionally requires propagating the Schrödinger equation for one period of the drive $N$ times (for each basis state as the initial condition). Although the specific time required to solve the differential equations depends on multiple factors (such as the solver, required tolerance, and hardware), in our experience [36] the exact calculation of the quasienergies for a single value of $U$ took 2.5–6 times longer than the EDPT calculation and, importantly, required $\sim 15$ times more memory (see Appendix B for details).