Second-order tunneling of two interacting bosons in a driven triple well

Advances in laser technology have enabled studies of quantum tunneling and its coherent control of a single particle in light-induced quantum wells without dissipation [1, 2]. Research attempting to manipulate quantum states has been underway for a long time [3, 4]. The time-periodic driving field is a powerful tool to control the tunneling dynamics and can lead to important phenomena, such as dynamic localization [5, 6], coherent destruction of tunneling (CDT) [7–10] and photon-assisted tunneling [11–13]. In recent years, the effects of interparticle interaction have attracted much attention. It has been shown that adjusting the interaction can give rise to richer behavior, including many-body selective CDT [14, 15] and the second-order tunneling of two interacting bosons [16]. The two-body interaction model is the simplest model for studying the interacting effects and it has received much attention [17–20], since the seminal experimental result was reported [21]. The tunneling dynamics are related to the interplay between the on-site interparticle interaction and external field, the former can be tuned by the Feshbach resonance technique [22, 23].

In the presence of interaction and a periodic external field, the quantum well system may be nonintegrable, which necessitates the use of a perturbation method for an analytical investigation. The multiple-scale technique is a very useful perturbation method and has been extensively employed for different physical systems [24–30]. It has been demonstrated that in the multiple-scale perturbation method the usual high-frequency approximation corresponds to the first-order perturbation correction [30]. In the far-off-resonant strongly interacting regime with a stronger reduced interaction [31], the high-frequency approximation is no longer valid. In this case, the dominant tunneling of doublons, which corresponds to the paired states in [32], is a second-order process of a long time scale that can be described by the second-order perturbation correction. The correlated tunneling of two strongly interacting atoms corresponding to time-resolved second-order tunneling has been observed directly in an undriven double well system [16]. Very recently, Longhi et al [32] have studied the second-order effect of two far-off-resonant strongly interacting bosons in a periodically driven optical lattice by using a multiple-time-scale asymptotic analysis.

As mentioned above, there have been a number of studies on the tunneling dynamics of two interacting atoms, which have focused on systems with double-well or optical lattice potentials. The triple-well system is a bridge between the double well and the optical lattice systems and it is very important for us to fully understand the coherent control of particle tunneling in the quantum wells [33–40]. From the triple-well (or equivalent three-level) system one can study richer physical phenomena, such as the effects of finite and multiple dimension for a few-particle case, the simplest non-nearest-neighbor couple and nonadjacent dipolar interaction [34, 35], and the adiabatical transport of a quantum particle from the left well to the right well with negligible middle-well occupation [33, 40]. The tunneling dynamics of two-particles in a triple-well system have also attracted extensive attention [41–43]; however, research on the second-order effect of the system has not yet been reported.

In this paper, we investigate the coherent control of second-order tunneling for two triple-well confined bosons driven by a high-frequency ac field. By means of the multiple-time-scale asymptotic analysis, we construct the second-order Floquet solutions, including the Rabi oscillation state and quasi-NOOOON state, and characterize the quantum dynamics of the two bosons with a continuous increase of interaction intensity. A far-off-resonant strongly interacting regime is demonstrated in which two bosons initially occupying the same well would form a stable bound pair because of the CDT from the doublons to the unpaired states. Taking into account the second-order tunneling effect, the selected CDT can also occur between the doublons or between the unpaired states on the six-dimensional Hilbert space. The result is confirmed by the Floquet quasienergy analysis, where the Floquet quasienergy band of the three unpaired states exhibits the avoided level-crossings (or new level-crossings) at (or near) the collapse points and the fine structure of quasienergy band of the three doublons shows the different level-crossings beyond the former collapse points. Good agreements between the analytical and numerical results are shown. The theoretical predictions could be verified further under the current accessible experimental setups [16, 21, 37, 40, 44, 45].

We consider two interacting bosons confined in a triple-well potential and driven by an ac field. The Hamiltonian of the system in the tight-binding approximation is described by the three-site Bose–Hubbard model [41–43]

$$\begin{equation} \fl \skew3\hat{H}(t)=-J(\hat{a}_1^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_1+\hat{a}_2^{\dagger}\hat{a}_{3}+\hat{a}_{3}^{\dagger}\hat{a}_2) +\frac{U_0}{2}\sum\limits_{l=1}^3\hat{a}_l^{\dagger} \hat{a}_l^{\dagger} \hat{a}_l \hat{a}_l +\varepsilon(t)(\hat{a}_3^{\dagger}\hat{a}_3-\hat{a}_1^{\dagger}\hat{a}_1), \end{equation} \tag{ 1 }$$

where the operator $\hat {a}_l^{(\dagger )}$ annihilates (creates) a boson in well l; J denotes the nearest-neighbor hopping matrix element; U₀ is the on-site interaction energy; and ε(t) = ε cos(ωt) is the ac driving of amplitude ε and frequency ω. Throughout this paper, we take ℏ = 1 and use the parameter J to appropriately scale the energy, U₀, ε and ω such that these quantities become dimensionless [14]. Experimentally, the triple-well potential can be generated by different methods [37, 40] and the periodic driving can be applied by shaking the system [44, 45]. Alternatively, the driven two-boson Hubbard model can also be simulated by light transport in coupled optical waveguides [46]. Here we have assumed that the three wells are deep enough such that the Wannier functions of the two interacting bosons belonging to different wells have very small overlap. A Fock basis |N_L, N_M, N_R〉 of six-dimensional Hilbert space is useful to describe the system, where N_L, N_M and N_R are the number of bosons localized in the left, middle and right wells, respectively, with N_L + N_M + N_R = 2. The quantum state |ψ(t)〉 of the system is expanded as the linear superposition of the Fock states

$$\begin{eqnarray} &&\fl |\psi(t)\rangle =c_1(t)|2,0,0\rangle+c_2(t)|0,2,0\rangle+c_3(t)|0,0,2\rangle+ c_4(t)|1,1,0\rangle +c_5(t)|1,0,1\rangle\nonumber\\ &&+\,c_6(t)|0,1,1\rangle, \end{eqnarray} \tag{ 2 }$$

where c_j(t) (j = 1,2,...,6) denotes the time-dependent probability amplitudes of finding the two bosons in the six different Fock states and obey the normalization condition $\sum_{j\!=\!1}^6 |c_{\!j}(t)|^2\!=\!\!1$. Inserting equations (1) and (2) into Schrödinger equation $\mathrm {i}\partial _t|\psi (t)\rangle =\hat {H}(t)|\psi (t)\rangle$, obtains the coupled equations for the amplitudes c_j(t)

$$\begin{eqnarray} \mathrm{i}\dot{c}_1&=&[U_0-2\varepsilon(t)]c_1-\sqrt{2}Jc_4,\nonumber\\ \mathrm{i}\dot{c}_2&=&U_0c_2-\sqrt{2}J(c_4+c_6),\nonumber\\ \mathrm{i}\dot{c}_3&=&[U_0+2\varepsilon(t)]c_3-\sqrt{2}Jc_6,\\ \mathrm{i}\dot{c}_4&=&-\varepsilon(t)c_4-J(\sqrt{2}c_1+\sqrt{2}c_2+c_5),\nonumber\\ \mathrm{i}\dot{c}_5&=&-J(c_4+c_6),\nonumber\\ \mathrm{i}\dot{c}_6&=&\varepsilon(t)c_6-J(\sqrt{2}c_2+\sqrt{2}c_3+c_5).\nonumber \end{eqnarray} \tag{ 3 }$$

Although it is difficult to obtain exact analytical solutions of equation (3), we can approximately study some interesting phenomena in the high-frequency regime with ω ≫ J. To do so, we rewrite the interaction strength as U₀ = mω + u for |u| ⩽ ω/2, m = 0,1,2,... with u being the reduced interaction strength [31], and make the function transformations $c_1(t)=a_1(t)\exp [-\mathrm {i}U_0t+2\mathrm {i}\varphi (t)]$, $c_2(t)=a_2(t)\exp (-\mathrm {i}U_0t)$, $c_3(t)=a_3(t)\exp [-\mathrm {i}U_0t-2\mathrm {i}\varphi (t)]$, $c_4(t)=a_4(t)\exp [\mathrm {i}\varphi (t)]$, c₅(t) = a₅(t) and $c_6(t)=a_6(t)\exp [-\mathrm {i}\varphi (t)]$, with a_j(t) being the slowly varying functions and $\varphi (t)=\int _0^t \varepsilon \cos (\omega \tau )\,{\mathrm { d}}\tau =\frac {\varepsilon }{\omega }\sin (\omega t)$. Then, equation (3) is transformed into the coupled equations in terms of a_j(t). Under the high-frequency approximation, the rapidly oscillating functions included in the equations can be replaced by their time average, such that the equations of a_j(t) become [43]

$$\begin{eqnarray} \mathrm{i}\dot{a}_1&=&-\sqrt{2}JJ_{m}\left(\frac{\varepsilon}{\omega}\right)a_4\,{\rm e}^{{{\rm i}}ut},\nonumber\\ \mathrm{i}\dot{a}_2&=&-\sqrt{2}J\left[(-1)^mJ_{m}\left(\frac{\varepsilon}{\omega}\right)a_4+J_{m}\left(\frac{\varepsilon}{\omega}\right)a_6\right]{\rm e}^{{{\rm i}}ut},\nonumber\\ \mathrm{i}\dot{a}_3&=&-\sqrt{2}J(-1)^mJ_{m}\left(\frac{\varepsilon}{\omega}\right)a_6\,{\rm e}^{{{\rm i}}ut},\nonumber\\\\ \mathrm{i}\dot{a}_4&=&-J\left[\sqrt{2}J_{m}\left(\frac{\varepsilon}{\omega}\right)a_1+(-1)^m\sqrt{2}J_{m}\left(\frac{\varepsilon}{\omega}\right)a_2\right]\,{\rm e}^{-{{\rm i}}ut}-JJ_{0}\left(\frac{\varepsilon}{\omega}\right)a_5,\nonumber\\ \mathrm{i}\dot{a}_5&=&-JJ_{0}\left(\frac{\varepsilon}{\omega}\right)(a_4+a_6),\nonumber\\ \mathrm{i}\dot{a}_6&=&-J\left[\sqrt{2}J_{m}\left(\frac{\varepsilon}{\omega}\right)a_2+(-1)^m\sqrt{2}J_{m}\left(\frac{\varepsilon}{\omega}\right)a_3\right]{\rm e}^{-{{\rm i}}ut}-JJ_{0}\left(\frac{\varepsilon}{\omega}\right)a_5,\nonumber \end{eqnarray} \tag{ 4 }$$

where J_m is the mth-order Bessel function of the first kind and e^±iut are the slowly varying functions for a small u value.

In [27], Frasca studied the localization in a driven two-level system beyond the high-frequency approximation by means of the multiple-scale perturbation method. Recently, Longhi [30] proposed that the well-known high-frequency approximation commonly used to study CDT corresponds to the first-order perturbation approximation of the multiple-time-scale asymptotic analysis. If the first-order correction term vanishes in the perturbation treatment, then the high-order corrected terms become important. Noticing that, for a set of fixed external field parameters, the dynamical behavior of the system (4) is related to the self-interaction intensity. In this work, we are not concerned about the very strong interaction (e.g. U₀ ⩾ 6ω) since for such a interaction we need to consider not only the usual on-site atom-interaction strength but also the interactions between atoms on neighboring lattice sites [17], which is beyond the considered case.

When the condition J₀ = 0 is satisfied, equation (4) shows that CDT occurs for the weakly interacting case (m = 0, |u| ≪ ω), which leads all the first derivatives of the probability amplitudes to zero. This can be further confirmed by calculation of the Floquet quasienergies of the system in section 4. Besides, for the resonant strongly interacting case (u = 0,m = 1,2,...), CDT for doublons is observed when the condition J_m = 0 is satisfied, which leads the first derivatives $\dot a_{\!j}(t)\ (j=1,2,3)$ of doublon-state amplitudes to zero. This is consistent with that of two interacting electrons in quantum dot arrays by numerical computation of the Floquet quasienergies [47]. As an example, we show time evolutions of the probabilities P_j(t) = |c_j|² = |a_j|² (j = 1,2,...,6) for the resonant case with J = 1, U₀ = ω = 80, ε/ω = 2.405 and P₂(0) = 1, P_j≠2(0) = 0, as in figure 1, where the first-order result (the circular points) from equation (4) is confirmed by the direct numerical simulation (the curves) of equation (3). From figure 1, we can see that transitions between the doublons with probabilities P_j, j = 1,2,3 in figure 1(a) and the unpaired states with probabilities P_j, j = 4,6 in figure 1(b) happen periodically for the resonant case. We will come back to this property to make a comparison with the difference from the far-off-resonant case in next section.

Zoom In Zoom Out Reset image size

It is known that equation (4) is a good approximation of equation (3) only for small values of the reduced interaction strength, |u| ≪ ω. When the |u| values tend to their maximum |u| = ω/2, the functions e^±iut vary moderately quickly compared to the rapidly oscillating driving field. Consequently, in equation (4), although e^±iut may be replaced by their average value of zero [43] such that probability amplitudes a₁(t), a₂(t) and a₃(t) of the doublons are frozen approximately, the effectiveness of the high-frequency approximation is partly lost. Particularly, in the case of moderate |u| values, the values are neither very small nor too large and, therefore, equation (4) is no longer a good approximation. For a stronger reduced interaction with a larger |u| value, we need to employ other approximation methods and to explore the second-order tunneling effects.

We will now consider the far-off-resonant case with a stronger reduced interaction to investigate the tunneling dynamics of the system, by means of multiple-time-scale asymptotic analysis. In the high-frequency regime, we set σ = J/ω as a small positive parameter and t' = ωt is the rescaled time. The probability amplitudes a_j(t') (j = 1,2,...,6) are expanded as a power series of σ

$$\begin{equation} a_j(t')=a_j^{(0)}(t')+\sigma a_j^{(1)}(t')+\sigma^2 a_j^{(2)}(t')+\cdots. \end{equation} \tag{ 5 }$$

Since the high-order infinitesimal can be neglected in the high-frequency regime, we approximately rewrite the probability amplitudes as the leading order a_j(t') = a⁽⁰⁾_j(t') = A_j(t'). Thus, |A_j|² = |c_j|² (j = 1,2,...,6) denotes the probabilities of finding the two bosons in the six different Fock states in equation (2). According to the perturbation analysis in the appendix, we readily obtain that such amplitudes are slowly varying functions in time, which satisfies the following linear equations with constant coefficients:

$$\begin{eqnarray} \mathrm{i}\frac{{\rm d}A_1}{{\rm d}t'}&=&2\sigma^2(A_1\rho_1+A_2\rho_2),\nonumber\\ \mathrm{i}\frac{{\rm d}A_2}{{\rm d}t'}&=&2\sigma^2\left[2A_2\rho_1+(A_1+A_3)\rho_2\right],\\ \mathrm{i}\frac{{\rm d}A_3}{{\rm d}t'}&=&2\sigma^2(A_3\rho_1+A_2\rho_2),\nonumber \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \mathrm{i}\frac{{\rm d}A_4}{{\rm d}t'}&=&-\sigma J_{0}\left(\frac{\varepsilon}{\omega}\right)A_5-2\sigma^2(2A_4\rho_1+A_6\rho_2),\nonumber\\ \mathrm{i}\frac{{\rm d}A_5}{{\rm d}t'}&=&-\sigma J_{0}\left(\frac{\varepsilon}{\omega}\right)(A_4+A_6),\\ \mathrm{i}\frac{{\rm d}A_6}{{\rm d}t'}&=&-\sigma J_{0}\left(\frac{\varepsilon}{\omega}\right)A_5-2\sigma^2(2A_6\rho_1+A_4\rho_2),\nonumber \end{eqnarray} \tag{ 7 }$$

where ρ_i (i = 1,2) are set as (see the appendix)

$$\begin{equation} \hspace{14pt}\rho_1=\sum\limits_{n'=-\infty}^{\infty}\frac{J^2_{n'}(\frac{\varepsilon}{\omega})}{\frac{U_0}{\omega}+n'}, \quad \rho_2=\sum\limits_{n'=-\infty}^{\infty}\frac{J_{n'}(\frac{\varepsilon}{\omega})J_{-n'}(\frac{\varepsilon}{\omega})}{\frac{U_0}{\omega}+n'} \end{equation} \tag{ 8 }$$

for U₀/ω + n' ≠ 0. Therefore, equations (6) and (7) are always definable and applicable, except for the resonant case in which ρ_i tends to infinity. It is worth noting that for a stronger reduced interaction obeying at least |u| > J, the value of any term in the summations of equation (8) is less than σ⁻¹, such that σ²ρ_i may be a second-order quantity and equations (6) and (7) could be applicable as a set of second-order equations. In fact, |ρ_i| < σ⁻¹ implies that the inequality |U₀/ω + n'| = |n + n' + u/ω| > σ = J/ω holds for any pair {n∈[0, ∞), n'∈(− ∞, ∞)}, which results in the conditions |u| > J and ω ⩾ 2|u| > 2J validating the second approximation. Particularly, we will numerically illustrate the parameter regions of |u|, ω and U₀, where the second-order perturbation method is more perfect for |u| ⩾ 10 J = ω/8 with U₀ = ω + u and for |u| ⩾ 5 J = ω/16 with U₀ = 2ω + u. Combining equation (6) with equation (7), we note that dynamics of the three doublons (i.e. the two bosons occupying the same site for A_j(t') with j = 1,2,3) is decoupled from that of the three unpaired states (i.e. the two bosons occupying the distinct sites for A_j(t') with j = 4,5,6).

Clearly, for |u| > J, equations (6) and (7) describe the second-order approximation, where the time evolution of any doublon state amplitude is a second-order long-time-scale process, since its time derivative is proportional to only the second-order constant σ². Similar results have been seen previously in the tight-binding optical lattice [32]. The second-order coupling coefficients of equations (6) and (7) are proportional to the parameter σ²ρ₂, which describes the second-order tunneling rate of the system. The nonzero tunneling coefficient means that the tunneling can occur between the three doublons based on equation (6) and between the three unpaired states based on equation (7). In figure 2, we plot the factor ρ₂ of the second-order tunneling coefficients as a function of the driving parameters ε/ω and self-interaction intensity U₀/ω, where figure 2(b) is the plan view of figure 2(a). Combining equation (8) with figure 2 we can see that the factor ρ₂ tends to infinity for any integer value of U₀/ω and arbitrary value ranges of ε/ω, while its values are small enough for the considered far-off-resonant regime. Note that, in figure 2 the very great ρ values are not shown because we have avoided the integer values of U₀/ω through selecting a rational step such that the multiple-time-scale asymptotic analysis holds. In the second approximation, equations (6) and (7) show that the CDT between the doublons will occur provided that the condition ρ₂ = 0 is satisfied.

Zoom In Zoom Out Reset image size

According to equations (6) and (7), we make an exact comparison of tunneling rates between the three doublons and three unpaired states for two different initial conditions as follows. Firstly, for the two bosons initially occupying the middle well [i.e. P₂(0) = 1, P_j≠2(0) = 0], we seek the analytical solutions P_j(t) (j = 1,2,...,6) from equations (6) and (7). To do so, we make the function transformations $A_1=A'_1\exp (-\mathrm {i}\frac {2J^2\rho _1}{\omega }t)$, $A_2=A'_2\exp (-\mathrm {i}\frac {4J^2\rho _1}{\omega }t)$ and $A_3=A'_3\exp (-\mathrm {i}\frac {2J^2\rho _1}{\omega }t)$. Inserting these expressions into equation (6) yields the coupled equations

$$\begin{equation} \mathrm{i}\frac{{\rm d}A'_1}{{\rm d}t}=\frac{2J^2\rho_2}{\omega}A'_2\,{\rm e}^{-{{\rm i}}\frac{2J^2\rho_1}{\omega}t}, \end{equation} \tag{ 9 }$$

$$\begin{equation} \mathrm{i}\frac{{\rm d}A'_2}{{\rm d}t}=\frac{2J^2\rho_2}{\omega}\,{\rm e}^{{{\rm i}}\frac{2J^2\rho_1}{\omega}t}(A'_1+A'_3), \end{equation} \tag{ 10 }$$

$$\begin{equation} \mathrm{i}\frac{{\rm d}A'_3}{{\rm d}t}=\frac{2J^2\rho_2}{\omega}A'_2\,{\rm e}^{-{{\rm i}}\frac{2J^2\rho_1}{\omega}t}. \end{equation} \tag{ 11 }$$

Combining equations (9) with (11) produces $\frac {{\mathrm { d}}(A'_1+A'_3)}{{\mathrm { d}}t}=\frac {4J^2\rho _2}{\omega }A'_2\,{\mathrm { e}}^{-{\rm{i}}\frac {2J^2\rho _1}{\omega }t}$. Eliminating (A'₁ + A'₃) from this equation and equation (10) yields the decoupled equation

$$\begin{equation} \mathrm{i}\ddot{A}'_2+\frac{2J^2\rho_1}{\omega}\dot{A}'_2+\mathrm{i}\frac{8J^4\rho_2^2}{\omega^2}A'_2=0. \end{equation} \tag{ 12 }$$

This is a second-order linear equation with constant coefficients, whose general solution is well-known, $A'_2=A'_+\exp (\mathrm {i}\chi _+t)+A'_-\exp (\mathrm {i}\chi _-t)$ for the parameters $\chi _{\pm }=\frac {J^2}{\omega }\rho _1\pm \frac {J^2}{\omega }\sqrt {\rho _1^2+8\rho _2^2}$ and the undetermined constants A'_± adjusted by the initial conditions. Under the above initial conditions, the general solution becomes the special solution

$$\begin{equation} A'_2(t)={\rm e}^{{{\rm i}}\frac{J^2\rho_1}{\omega}t}\left[\cos(\omega_1t)-\mathrm{i}\frac{\rho_1}{\rho_1^2+8\rho_2^2}\sin(\omega_1t)\right] \end{equation} \tag{ 13 }$$

with $\omega _1=\frac {J^2}{\omega }\sqrt {\rho _1^2+8\rho _2^2}$. Thus, under the initial conditions P₂(0) = 1, P_j≠2(0) = 0, the analytical probabilities from equations (6) and (7) can be readily obtained as P₄(t) = P₅(t) = P₆(t) = 0, and

$$\begin{equation} P_2(t)=|A'_2(t)|^2=\frac{\rho_1^2}{\rho_1^2+8\rho_2^2}+\frac{8\rho_2^2}{\rho_1^2+8\rho_2^2}\cos^2(\omega_1t), \end{equation} \tag{ 14 }$$

$$\begin{equation} P_1(t)=P_3(t)=\frac{4\rho_2^2}{\rho_1^2+8\rho_2^2}\sin^2(\omega_1t), \end{equation} \tag{ 15 }$$

which corresponds to a Rabi oscillation state of equation (2). Secondly, for the two bosons initially occupying the left and right wells (i.e. P₅(0) = 1, P_j≠5(0) = 0), we easily obtain another Rabi oscillation solution with probabilities P₁(t) = P₂(t) = P₃(t) = 0, P₅(t) = cos²(ω₂t) and $P_4(t)=P_6(t)=\frac {1}{2}\sin ^2(\omega _2t)$ for $\omega _2=\sqrt {2}JJ_{0}(\frac {\varepsilon }{\omega })$ if we neglect the second-order small quantities of equation (7) in the high-frequency regime. For example, selecting the parameters ω = 80 ≫ J = 1, ε/ω = 2 and U₀/ω = 1.5. Therefore, the tunneling period T₁ = π/ω₁ ≈ 89 corresponds to the second-order tunneling effect while the tunneling period T₂ = π/ω₂ ≈ 10 corresponds to the first-order tunneling effect.

In figure 3, we numerically plot the time evolutions of the probabilities P_j (j = 1,2,...,6) based on equation (3) for the above two initial conditions, the circular points correspond to the above analytical results. Obviously, the analytical results are in perfect agreement with the numerical simulations. The zero probability of the unpaired states in figure 3(a) and the zero probability of the doublons in figure 3(b) mean the CDT from the doublons to the unpaired states. This CDT leads to the interesting two-boson correlated tunneling.

Zoom In Zoom Out Reset image size

Furthermore, we can see the selected CDT between the doublons for arbitrary nonzero parameters of equation (6) and between the unpaired states for $J_{0}(\frac {\varepsilon }{\omega })=0$ in equation (7), which displays different second-order effects. Firstly, under the initial conditions $A_1(0)=\frac {1}{\sqrt {2}}=- A_3(0)$ and A_j≠1,3(0) = 0, from equation (6) we obtain the analytical solutions A_j≠1,3(t) = 0 and $A_1(t)=-A_3(t)=\frac {1}{\sqrt {2}} \exp (-\mathrm {i} \frac {2J^2\rho _1}{\omega } t)$ with the second-order small phase. Inserting these into equation (2), and applying the relations among A(t), a(t) and c(t), produces the quantum state $|\psi (t)\rangle =\frac {1}{\sqrt {2}} \exp \{-\mathrm {i}[U_0t+\frac {2J^2\rho _1}{\omega } t-\frac {2\varepsilon }{\omega }\sin (\omega t)]\}\{|2,0,0\rangle -\exp [-\mathrm {i}\frac {4\varepsilon }{\omega }\sin (\omega t)]|0,0,2\rangle \}$. This superposition state between |2,0,0〉 and |0,0,2〉 may be called the quasi-NOOOON state, compared to the well-known NOON state [48]. In such a state, the two bosons occupy the left well and the right well with the same probability $P_1=P_3=\frac 1 2$, and the middle-well occupation probability always vanishes, which is shown in figure 4(a) analytically (circular points) and numerically (curves). This means the selected CDT between the three doublon states. Secondly, if we take the parameter to obey $J_{0}(\frac {\varepsilon }{\omega })=0$, then the dominating dynamics is decided by the second-order correction terms in equation (7). In figure 4(b), we numerically plot the time evolutions of the probabilities P_j (j = 1,2,...,6) based on equation (3) for two bosons initially occupying the left and middle well, respectively (i.e. P₄(0) = 1, P_j≠4(0) = 0), where the parameter is set as ε/ω = 2.405 corresponding to the first zero of $J_{0}(\frac {\varepsilon }{\omega })$ and the other parameters are the same as those of figure 3(b). In this case, we readily calculate the analytical probabilities, because of $\dot {A}_5(t)=0$ in equation (7). Taking the initial conditions A₄(0) = 1 and A_j≠4(0) = 0, we immediately obtain the analytical solutions A₁ = A₂ = A₃ = A₅ = 0, $A_4=\exp (\mathrm {i}\frac {4J^2\rho _1}{\omega }t)\cos (\frac {2J^2\rho _2}{\omega }t)$ and $A_6=\exp (\mathrm {i}\frac {4J^2\rho _1}{\omega }t)\sin (\frac {2J^2\rho _2}{\omega }t)$. The corresponding probabilities read P₁ = P₂ = P₃ = P₅ = 0, $P_4=|A_4(t)|^2=\cos ^2(\frac {2J^2\rho _2}{\omega }t)$ and $P_6=\sin ^2(\frac {2J^2\rho _2}{\omega }t)$, which are plotted by the circular points in figure 4(b). Here the two bosons perform Rabi oscillations between the state |1,1,0〉 and the state |0,1,1〉 with period $T_3=\frac {\pi \omega }{2\,J^2|\rho _2|}\approx 132$ which is in the same order of magnitude as the above second-order tunneling period T₁ ≈ 89. Here the zero probability in the unpaired state |1,0,1〉 means the selected CDT between the three unpaired states. The analytical result agrees with the numerical result in figure 4.

Zoom In Zoom Out Reset image size

In order to further confirm the analytical results in equations (4), (6) and (7), we define the time-averaged total probability of finding the two interacting bosons in the three doublons as $\langle S\rangle =\langle P_1\rangle +\langle P_2\rangle +\langle P_3\rangle =\frac {1}{\tau }\int _0^{\tau }(P_1+P_2+P_3)\,{\mathrm { d}}t$ for τ = 200(J⁻¹). The normalization means the time-averaged total probability in the three unpaired states being 1 − 〈S〉. Taking the initial conditions P₂(0) = 1,P_j≠2(0) = 0 and parameter J = 1 from equation (3) we numerically give 〈S〉 as the function of the self-interaction U₀ for ω = 80 and ε/ω = 2.405, as in figure 5(a). It is shown that the time-averaged total probability in the three doublons possesses different features, for the multiphoton resonant points, the far-off-resonant regions, and the near-resonant regions, respectively. Firstly, at each of the resonant points (i.e. U₀ = mω with m = 1,2,...,5 being integer) 〈S〉 drops to the lowest points which means that the separation probability 1 − 〈S〉 of the two bosons is the largest, as shown in figure 1. Secondly, the two bosons can also be separated in the near-resonant regions; however, the time-averaged total probability 〈S〉 tends to one and the separation probability tends to zero rapidly as increasing the reduced interaction strength |u|. Finally, in the far-off-resonant regions with larger |u| values, 〈S〉 is always equal to 1. Let half-width of the valley centered at mth resonant point of figure 5(a) be |u_m|. The largest two half-widths for fitting 〈S〉 ≈ 1 can be estimated as |u₁| ≈ 10 J and |u₂| ≈ 5 J from the first and second valleys, as shown in the insets of figure 5(a). This indicates that a selected CDT between the doublons and the unpaired states can happen in the parameter regions |u₁| ⩾ 10J = ω/8 for U₀ = ω + u, and |u₂| ⩾ 5J = ω/16 for U₀ = 2ω + u, which belonged to the far-off-resonant strongly interacting regime in this paper. Such a CDT enables the two bosons form a stable bound pair and cannot move independently for a stronger reduced interaction [47]. Comparing with the similar phenomena in [21], it can be understood that any doublon state has a potential energy offset U₀ with respect to the unpaired states, such that breaking up of the stable pair is suppressed due to the band structure and energy conservation of the system. Only for the resonant case can the boson pair be separated, which happens because of the resonant absorbtion of energy from the ac driving field. In the weakly interacting regime, CDT is expected to occur in figure 5(a) because ε/ω = 2.405 is the first root of J₀(ε/ω) = 0, which we will consider further in the next section. To show that the above analysis is generic in the high-frequency regime, we plot the time-averaged total probability 〈S〉 of the doublons as a function of the ratio U₀/ω of the interaction and the driving frequency in figure 5(b) with ε = 160 and U₀ = 200. Figure 5(b) explicitly shows that the lowest points of 〈S〉 appear at the resonant points U₀/ω = n for integer n in the high-frequency regime (n ⩽ 30, i.e. ω ⩾ 6.66). The results agree with the above analysis on figure 5(a). The plots similar to figures 5(a) and (b) can also be given for the parameter regions ω∈[10,80] and U₀∈[10,200], respectively, which are not shown here. It is worth noting that in the plot similar to figure 5(a), lower frequency ω is associated with a smaller |u_m| value.

Zoom In Zoom Out Reset image size

It is well known that CDT can occur at the collapse points of the Floquet quasienergy spectrum [30, 49, 50], so the above-mentioned tunneling properties will be confirmed by the Floquet quasienergy analysis as follows:

The Floquet theory provides a powerful tool to analyze the dynamics of a time-periodic quantum system [51]. According to the Floquet theory, the solutions of the time-dependent Schrödinger equation can be written as |ψ_k(t)〉 = e^−iE_kt|ϕ_k(t)〉, with |ϕ_k(t)〉 being the Floquet states and E_k Floquet quasienergies. In analogy to the Bloch solutions for the spatially periodic system, the quasienergy can only be determined up to a integer multiple of the photon energy ω. For the sake of definiteness, it is usually assumed to vary in the first Brillouin zone −ω/2 < E ⩽ ω/2. The Floquet states inherit the period of the Hamiltonian and are eigenstates of the time evolution operator for one period of the driving

$$\begin{equation} U(T,0)=\mathcal{T}\exp\left[-\mathrm{i}\int^T_0H(t)\,{\rm d}t\right], \end{equation} \tag{ 16 }$$

where $\mathcal {T}$ is the time-ordering operator and T = 2π/ω is the period of the driving. Noticing that eigenvalues of U(T,0) are $\exp (-\mathrm {i}E_kT)$, the quasienergies of this system can be determined directly so long as we diagonalize U(T,0). In figure 6, selecting the parameters as J = 1 and ω = 80, we show the numerical results of the quasienergy spectra as the functions of driving parameters ε/ω for U₀ = mω, m = 0,1 with zero reduced interaction, respectively. In figure 6(a), for two noninteracting bosons, the quasienergy spectrum shows collapses at some fixed values of the driving parameters, for which $J_{0}(\frac {\varepsilon }{\omega })=0$. The inset of figure 6(a) is an enlargement of quasienergies near the collapse point ε/ω = 2.405, corresponding to the first zero of $J_{0}(\frac {\varepsilon }{\omega })$ and shows an exact level-crossing at this collapse point, this is analogous to a single boson in a triple-well system [35]. In the resonant regime with U₀ = ω = 80, the quasienergies of figure 6(b) show that the crossings of some quasienergies and the avoided crossing of the other quasienergies appear at the zero points (ε/ω = 3.832,...) of the first-order Bessel function J₁(ε/ω). From equation (4) with m = 1,u = 0 we know that the first derivative of the probability amplitudes of the three doublons $\dot {a}_j\, (j=1,2,3)$ are equal to zero and the three unpaired states $\dot {a}_j\, (j=4,5,6)$ do not vanish at the zero points of J₁(ε/ω), which indicates that the CDT occurs between the three doublons for the crossing of quasienergy bands and the tunneling could appear between the three unpaired states for the avoided crossing of quasienergy bands. In the considered case u = 0, equations (6) and (7) are no longer valid and any quasienergy may be associated with both the doublons and the unpaired states; however, different quasienergies may have different projections onto the doubly or singly occupied basis vectors. In [47], Creffield and Platero distinguished the quasienergies according to the scheme of whether the Floquet states project mainly onto doubly or singly occupied basis vectors. Based on this scheme, we also numerically compute the different kinds of the quasienergies and distinguish them in figure 6, where the thin lines indicate the quasienergies associated with the states mainly projecting onto doubly occupied basis vectors and the heavy lines correspond to the singly occupied states. In figure 6(a) for U₀ = 0, we show that the curve of zero energy is a coincidence of the thin and heavy lines, while this coincided line is separated into two near curves in figure 6(b) with U₀ = ω. Although the level crossings occur for all states in figure 6(a), they only occur for the doublons in figure 6(b), which shows the different tunneling properties (as in the above-mentioned analytical results).

Zoom In Zoom Out Reset image size

When the reduced interaction strength is sufficiently larger (e.g. |u| > J), the quasienergy spectrum is divided into two energy bands, which correspond to the three doublons and the three unpaired states, respectively (as shown in figure 7), where the quasienergies of doublons (or unpaired states) aperiodically oscillate near u values (or 0 value), hence the width of the energy gap between the two bands is proportional to the |u| value. For a weaker interaction with U₀ = u = 2, CDT between the different doublons and between the different unpaired states can be realized for the same driving parameters, as indicated by the level-crossing points in figure 7(a), when the ratio of the field amplitude ε and the field frequency ω is a root of the equation J₀(ε/ω) = 0. Precise agreements between the numerical results based on equation (3) and the analytical results from equations (19), (20) and (23) are observed in figure 7(a) for a sufficiently larger range of ratio ε/ω. The small deviation between the both results in figure 7(a) indicates that the second-order perturbation method is perfectly applicable only for some suitable parameter regions. We have also investigated the quasienergy spectra for |u| = 5,6,... which are not shown here. The results verify that, in the far-off-resonant strongly interacting regime, ω/2 ⩾ |u| ⩾ |u|₁ ≈ 10 J is just the above suitable parameter regions. Interestingly, the energy band corresponding to the doublons becomes narrower and the energy gap tends to widen with the increase of self-interaction intensity from |u| = 2 < |u|₁ to |u| = 30 > |u|₁. A wider gap means that the quantum transition between the both kinds of states is hard to occur. A narrower band necessitates analysis the Floquet quasienergy spectrum from both cases of the unpaired states and the doublons, respectively.

Zoom In Zoom Out Reset image size

4.1. Avoided level-crossing of unpaired states

In the far-off-resonant strongly interacting regime, selecting the parameters as J = 1, ω = 80 and U₀ = 50 (i.e. m = 1 and u = −30), from equation (3) we numerically plot quasienergy spectrum versus ε/ω in figure 7(b). In this figure, the quasienergies corresponding to the unpaired states shows collapses when ε/ω are the roots of J₀(ε/ω) = 0; however, the energy band corresponding to the doublons has collapsed into an approximate straight line. The inset of figure 7(b) is an enlargement of quasienergies corresponding to the three unpaired states near the first collapse point ε/ω ≈ 2.405. The fine structure of energy spectrum exhibits that the pseudocollapse point is converted to an avoided crossing point at ε/ω ≈ 2.405 and two different crossing points due to the second-order correction terms in equation (7).

To explain the numerical result, from equation (7) we analytically calculate the quasienergies corresponding to the three unpaired states. Note that the period of functions $\exp [-\mathrm {i}\varphi (t)]$ and $\exp [\pm 2\mathrm {i}\varphi (t)]$ is T. Therefore, we can construct the Floquet states by setting [35] $A_j(t)=B_{\!j}\exp (-\mathrm {i}Et)$ (j = 4,5,6) for the three unpaired states with constant B_j, then rewriting equation (7) as the time-independent form

$$\begin{eqnarray} E\!B_4&=&-JJ_{0}\left(\frac{\varepsilon}{\omega}\right)B_5-2\frac{J^2}{\omega}(2B_4\rho_1+B_6\rho_2),\nonumber\\ E\!B_5&=&-JJ_{0}\left(\frac{\varepsilon}{\omega}\right)(B_4+B_6),\\ E\!B_6&=&-JJ_{0}\left(\frac{\varepsilon}{\omega}\right)B_5-2\frac{J^2}{\omega}(2B_6\rho_1+B_4\rho_2).\nonumber \end{eqnarray} \tag{ 17 }$$

The existence condition for the nontrivial solution of equation (17) reads

$$\begin{equation} \left| \begin{array}{@{}clr@{}} E+4\displaystyle\frac{J^2}{\omega}\rho_1 & J{J}_{0}\left(\displaystyle\frac{\varepsilon}{\omega}\right) & 2\displaystyle\frac{J^2}{\omega}\rho_2 \\ J{J}_{0}\left(\displaystyle\frac{\varepsilon}{\omega}\right) & E & J{J}_{0}\left(\displaystyle\frac{\varepsilon}{\omega}\right) \\ 2\displaystyle\frac{J^2}{\omega}\rho_2 & J{J}_{0}\left(\displaystyle\frac{\varepsilon}{\omega}\right) & E+4\displaystyle\frac{J^2}{\omega}\rho_1 \end{array} \right|=0. \end{equation} \tag{ 18 }$$

From equation (18) we obtain three Floquet quasienergies corresponding to the three unpaired states

$$\begin{eqnarray} E_1&=&-\frac{2(2J^2\rho_1-J^2\rho_2)}{\omega},\nonumber\\ E_2&=&\frac{-2J^2\rho_1-J^2\rho_2-\rho_3}{\omega},\nonumber\\ E_3&=&\frac{-2J^2\rho_1-J^2\rho_2+\rho_3}{\omega}, \end{eqnarray} \tag{ 19 }$$

where we have set

$$\begin{equation} \rho_3=\sqrt{(2J^2\rho_1+J^2\rho_2)^2+2J^2\omega^2J_{0}^2\left(\frac{\varepsilon}{\omega}\right)}. \end{equation} \tag{ 20 }$$

We now compare the analytical results of equations (19) and (20) with the numerical computation based on the original equation (3). A typical behavior of quasienegies near the first crossing point is plotted in the inset of figure 7(b). It is clearly shown that the analytical result (the circular points) is in perfect agreement with the direct numerical computation (the curves).

4.2. A fine structure of quasienergy spectrum of doublons

Next, we examine some detailed features of the quasienergies corresponding to the three doublons. In [32], Longhi et al proposed that in a lattice system, CDT can be realized between the doublons and between the unpaired states for the same parameters, namely the field parameters take the second root ε/ω = 5.52 of J₀(ε/ω) = 0 and the interaction intensity obeys U₀/ω = 2.58 corresponding to ρ₂ = 0. Here, for the triple-well system we prove a similar result and exhibit a fine structure of quasienergy spectrum of doublons, based on the analytical Floquet solutions of equations (6)–(8). According to equation (6), CDT occurs between the doublons if the condition ρ₂ = 0 is satisfied. From equation (8), we have ρ₂ = 0 at ε/ω ≈ 0.95 in the region ε/ω∈[0,8] for U₀/ω = 1.6, and have ρ₂ = 0 at ε/ω ≈ 1.20, 2.02, 5.52 and 5.74 in the same region for U₀/ω = 2.58. Selecting two different values of the self-interaction intensity, from equation (3) we numerically plot quasienergy spectrum versus ε/ω in figures 8(a) and (b), respectively, with the insets being enlargements of quasienergies of the three doublons. At the points fitting ρ₂ = 0, the level-crossing of two quasienergies will occur. This indicates that CDT for doublons can be observed at the crossing points of the partial levels. The predictions of the perturbation analysis can be confirmed by direct numerical computation of the temporal evolution of the boson occupation probabilities from equation (3) (not depicted here).

Zoom In Zoom Out Reset image size

We then analytically calculate the quasienergies corresponding to the three doublons. Substituting the Floquet solutions $A_j(t)=B_j\exp [-\mathrm {i}(E-U_0)t]$ (j = 1,2,3) into equation (6), we obtain easily

$$\begin{eqnarray} (E-U_0)B_1&=&2\frac{J^2}{\omega}(B_1\rho_1+B_2\rho_2),\nonumber\\ (E-U_0)B_2&=&2\frac{J^2}{\omega}[2B_2\rho_1+(B_1+B_3)\rho_2],\\ (E-U_0)B_3&=&2\frac{J^2}{\omega}(B_3\rho_1+B_2\rho_2)\nonumber \end{eqnarray} \tag{ 21 }$$

with the existence condition of the nontrivial solution

$$\begin{equation} \fl \left(E-U_0-4\frac{J^2}{\omega}\rho_1\right)\left(E-U_0-2\frac{J^2}{\omega}\rho_1\right)^2-8\frac{J^4}{\omega^2}\rho_2^2\left(E-U_0-2\frac{J^2}{\omega}\rho_1\right)=0. \end{equation} \tag{ 22 }$$

From this equation we obtain three Floquet quasienergies corresponding to the three doublons as

$$\begin{eqnarray} E_4&=&U_0+2\frac{J^2}{\omega}\rho_1,\nonumber\\ E_5&=&U_0+\frac{3J^2\rho_1+\sqrt{J^4\rho_1^2+8J^4\rho_2^2}}{\omega},\\ E_6&=&U_0+\frac{3J^2\rho_1-\sqrt{J^4\rho_1^2+8J^4\rho_2^2}}{\omega}.\nonumber \end{eqnarray} \tag{ 23 }$$

We note that the quasienergies E_j for j = 4,5,6 should be converted to E'_j in the first Brillouin zone (i.e. E'_j = E_j − mω = u + E_j − U₀). Thus, the quasienergies are mainly decided by the reduced interaction strength u in the high-frequency regime and their positions depend on the sign and the size of u. At the same time, the quasienergies are no longer degenerate due to the strong interaction-induced second-order corrections. For example, U₀ = 2.58ω ≈ 206.4, so the quasienergies are rewritten as E'_j = E_j − 3ω for ω = 80 (j = 4,5,6). We plot the analytical quasienergies versus ε/ω for U₀/ω = 1.6 and 2.58 in the insets of figures 8(a) and (b) as circular points, respectively, which are in perfect agreement with the direct numerical computations (the curves) based on the original equation (3).

We have investigated the tunneling dynamics of two bosons in a high-frequency driven triple well for a continuously increasing interaction intensity by means of the multiple-time-scale asymptotic analysis. In the far-off-resonant strongly interacting regime, we have analytically constructed the second-order Floquet solutions, including the Rabi oscillation state and quasi-NOOOON state. It is shown that the dominant tunneling effect of doublons is a second-order process, similar to two bosons in a driven optical lattice [32]. For a stronger reduced interaction, we make an exact comparison of tunneling rates between the doublons and unpaired states, and find that two bosons initially occupying the same well will form a stable bound pair and the stable unpaired states can also be kept for two initial unpaired particles. Therefore, the selected CDT between the doublons or between the unpaired states can occur for different values of interaction intensity. However, in the near-resonant case such initially paired bosons can separate due to the multiphoton resonance. Furthermore, we calculate the Floquet quasienergy spectrum and demonstrate that for the reduced interaction strength obeying |u| > J, the quasienergy is divided into two energy bands corresponding to the three doublons and the three unpaired states. The width of the energy gap between the two bands is proportional to the |u| value. The prediction on the CDT is confirmed by the quasienergy spectra, in which the avoided level-crossings and new level-crossings near the collapse points are exhibited for the three unpaired states due to the second-order corrections. While for the three doublons, a fine structure of quasienergy spectrum up to the second-order is displayed, by which we show the different level-crossings beyond the former collapse points.

These analytical results are consistent with the direct numerical computations from the time-dependent Bose–Hubbard Hamiltonian for some sufficiently large values of the reduced interaction strength |u|, which allows the wide parameter regions ω∈[10,80] and U₀∈[10,200] to fit experimental requirements. The second-order results of the long time scale may be conveniently applied to adiabatic manipulation [4, 52] of the paired-particle correlated tunneling experimentally. In the experiments [16, 21], the second-order tunneling of two strongly interacting bosons has been observed in the undriven optical lattice [21] and undriven double well [16]. The triple-well potential has also been prepared by using different methods [37, 40]. Combining these with the periodical shaking method [44, 45] and the Feshbach resonance technique [22, 23] has enabled us to show that experimental verification of the theoretical prediction in this paper could be expected under the currently accessible setups [16, 21, 37, 40, 44, 45]. Particularly, the results from the triple-well system can be extended to the three-level system [33] and triple-quantum-dot system [9], exhibiting richer new physics.

This work was supported by the NNSF of China under grant numbers 11175064 and 11375059, the Construct Program of the National Key Discipline of China and the Hunan Provincial NSF (grant numbers 11JJ7001 and 14JJ4052).

In the high-frequency regime, σ = J/ω is a small positive parameter. Let t' = ωt be the rescaling dimensionless time variable, then we rewrite equation (4) as

$$\begin{eqnarray} \mathrm{i}\frac{{\rm d}a_1}{{\rm d}t'}&=&-\sqrt{2}\sigma a_4\,{\rm e}^{\left[{\rm i}\frac{U_0}{\omega}t'-\mathrm{i}\varphi(t')\right]},\nonumber\\ \mathrm{i}\frac{{\rm d}a_2}{{\rm d}t'}&=&-\sqrt{2}\sigma\left[a_4 \,{\rm e}^{\left[{\rm i}\frac{U_0}{\omega}t'+\mathrm{i}\varphi(t')\right]}+a_6\,{\rm e}^{\left[{\rm i}\frac{U_0}{\omega}t'-\mathrm{i}\varphi(t')\right]}\right],\nonumber\\ \mathrm{i}\frac{{\rm d}a_3}{{\rm d}t'}&=&-\sqrt{2}\sigma a_6\,{\rm e}^{\left[{\rm i}\frac{U_0}{\omega}t'+\mathrm{i}\varphi(t')\right]},\nonumber\\ \mathrm{i}\frac{{\rm d}a_4}{{\rm d}t'}&=&-\sigma\Big[\sqrt{2}a_1\,{\rm e}^{\left[-{\rm i}\frac{U_0}{\omega}t'+\mathrm{i}\varphi(t')\right]}+\sqrt{2}a_2 \,{\rm e}^{\left[-{\rm i}\frac{U_0}{\omega}t'-\mathrm{i}\varphi(t')\right]}+a_5 \,{\rm e}^{[-{\rm i}\varphi(t')]}\Big],\nonumber\\ \mathrm{i}\frac{{\rm d}a_5}{{\rm d}t'}&=&-\sigma\Big[a_4 \,{\rm e}^{{{\rm i}}\varphi(t')}+a_6 \,{\rm e}^{[-{\rm i}\varphi(t')]}\Big],\nonumber\\ \mathrm{i}\frac{{\rm d}a_6}{{\rm d}t'}&=&-\sigma\Big[\sqrt{2}a_2 \,{\rm e}^{\left[-{\rm i}\frac{U_0}{\omega}t'+\mathrm{i}\varphi(t')\right]}+\sqrt{2}a_3 \,{\rm e}^{\left[-{\rm i}\frac{U_0}{\omega}t'-\mathrm{i}\varphi(t')\right]}+a_5 \,{\rm e}^{{{\rm i}}\varphi(t')}\Big] \end{eqnarray} \tag{ A.1 }$$

with $\varphi (t')=\frac {\varepsilon }{\omega }\sin t'$. At first, we transformed the independent variable t' into the multiple-time-scale variables T_n = σⁿt', n = 0,1,2,..., and then replaced the time derivatives by the expansion

$$\begin{equation} \frac{\rm d}{{\rm d}t'}=\partial_{T_0}+\sigma\partial_{T_1}+\sigma^2\partial_{T_2}+\cdots. \end{equation} \tag{ A.2 }$$

At the same time, we expand a_j(t')  (j = 1,2,...,6) as the power series of σ

$$\begin{equation} a_j(t')=a_j^{(0)}(t')+\sigma a_j^{(1)}(t')+\sigma^2 a_j^{(2)}(t')+\cdots. \end{equation} \tag{ A.3 }$$

Substituting equations (A.2) and (A.3) into equation (A.1), and collecting the terms of the same order, we obtain a hierarchy of approximation equations of different orders in σ. At leading order σ⁰, one has

$$\begin{equation} \frac{\partial a_j^{(0)}}{\partial T_0}=0, \quad a_j^{(0)}=A_j(T_1,T_2,\ldots), \end{equation} \tag{ A.4 }$$

where the amplitudes A_j(T₁,T₂,...) are functions of the slow time variables T₁,T₂,... but are independent of the fast time variable T₀. At order σ, one obtains the coupled equations:

$$\begin{eqnarray} &&\mathrm{i}\frac{\partial a_1^{(1)}}{\partial T_0}=-\mathrm{i}\frac{\partial A_1}{\partial T_1}-\sqrt{2}A_4\,{\rm e}^{\left[{\rm i}\frac{U_0}{\omega}T_0-\mathrm{i}\varphi(T_0)\right]},\nonumber\\ &&\mathrm{i}\frac{\partial a_2^{(1)}}{\partial T_0}=-\mathrm{i}\frac{\partial A_2}{\partial T_1}-\sqrt{2}\Big[A_4\,{\rm e}^{\left[{\rm i}\frac{U_0}{\omega}T_0+\mathrm{i}\varphi(T_0)\right]}+A_6\,{\rm e}^{\left[{\rm i}\frac{U_0}{\omega}T_0-{\rm i}\varphi(T_0)\right]}\Big],\nonumber\\ &&\mathrm{i}\frac{\partial a_3^{(1)}}{\partial T_0}=-\mathrm{i}\frac{\partial A_3}{\partial T_1}-\sqrt{2}A_6\,{\rm e}^{\left[{\rm i}\frac{U_0}{\omega}T_0+\mathrm{i}\varphi(T_0)\right]},\nonumber\\\\ &&\mathrm{i}\frac{\partial a_4^{(1)}}{\partial T_0}=-\mathrm{i}\frac{\partial A_4}{\partial T_1}-\Big[\sqrt{2}A_1\,{\rm e}^{\left[-{\rm i}\frac{U_0}{\omega}T_0+\mathrm{i}\varphi(T_0)\right]}+\sqrt{2}A_2\,{\rm e}^{\left[-{\rm i}\frac{U_0}{\omega}T_0-\mathrm{i}\varphi(T_0)\right]}+A_5 \,{\rm e}^{[-{\rm i}\varphi(T_0)]}\Big],\nonumber\\ &&\mathrm{i}\frac{\partial a_5^{(1)}}{\partial T_0}=-\mathrm{i}\frac{\partial A_5}{\partial T_1}-\Big[A_4\,{\rm e}^{{{\rm i}}\varphi(T_0)}+A_6\,{\rm e}^{[-{\rm i}\varphi(T_0)]}\Big],\nonumber\\ &&\mathrm{i}\frac{\partial a_6^{(1)}}{\partial T_0}=-\mathrm{i}\frac{\partial A_6}{\partial T_1}-\Big[\sqrt{2}A_2\,{\rm e}^{\left[-{\rm i}\frac{U_0}{\omega}T_0+\mathrm{i}\varphi(T_0)\right]}+\sqrt{2}A_3\,{\rm e}^{\left[-{\rm i}\frac{U_0}{\omega}T_0-\mathrm{i}\varphi(T_0)\right]}+A_5\,{\rm e}^{{{\rm i}}\varphi(T_0)}\Big].\nonumber \end{eqnarray} \tag{ A.5 }$$

For the convenience of our discussion, we simplify equation (A.5) as i∂a⁽¹⁾_j/∂T₀ = −i∂A_j/∂T₁ + G⁽¹⁾_j(T₀) for j = 1,2,...,6. To avoid the occurrence of secular growing terms in the solution a⁽¹⁾_j, the solvability condition [30, 32]

$$\begin{equation} \mathrm{i}\frac{\partial A_j}{\partial T_1}=\overline{G_j^{(1)}(T_0)} \end{equation} \tag{ A.6 }$$

must be satisfied, where the overline denotes the time average with respect to the fast time variable T₀ (i.e. the dc component of the driving term G⁽¹⁾_j(T₀)). The amplitudes a_j at order σ are given by

$$\begin{equation} a_j^{(1)}=-\mathrm{i}\int\left[G_j^{(1)}(T_0)-\overline{G_j^{(1)}(T_0)}\,\right]{\rm d} T_0. \end{equation} \tag{ A.7 }$$

Employing equations (A.5) and (A.6), one gives

$$\begin{eqnarray} \mathrm{i}\frac{\partial A_1}{\partial T_1}&=&\mathrm{i}\frac{\partial A_2}{\partial T_1}=\mathrm{i}\frac{\partial A_3}{\partial T_1}=0,\nonumber\\ \mathrm{i}\frac{\partial A_4}{\partial T_1}&=&\mathrm{i}\frac{\partial A_6}{\partial T_1}=-J_{0}\left(\frac{\varepsilon}{\omega}\right)A_5,\\ \mathrm{i}\frac{\partial A_5}{\partial T_1}&=&-J_{0}\left(\frac{\varepsilon}{\omega}\right)(A_4+A_6).\nonumber \end{eqnarray} \tag{ A.8 }$$

So the solutions of order σ read

$$\begin{eqnarray} a_1^{(1)}&=&\sqrt{2}\mathrm{i}A_4F_0(T_0),\nonumber\\ a_2^{(1)}&=&\sqrt{2}\mathrm{i}[A_4F_1(T_0)+A_6F_0(T_0)],\nonumber\\ a_3^{(1)}&=&\sqrt{2}\mathrm{i}A_6F_1(T_0),\nonumber\\\\a_4^{(1)}&=&\mathrm{i}[\sqrt{2}A_1F^*_0(T_0)+\sqrt{2}A_2F^*_1(T_0)+A_5F^*_2(T_0)],\nonumber\\ a_5^{(1)}&=&\mathrm{i}[A_4F_2(T_0)+A_6F^*_2(T_0)],\nonumber\\ a_6^{(1)}&=&\mathrm{i}[\sqrt{2}A_2F^*_0(T_0)+\sqrt{2}A_3F^*_1(T_0)+A_5F_2(T_0)]\nonumber \end{eqnarray} \tag{ A.9 }$$

with

$$\begin{eqnarray} F_0(T_0)&=&\sum\limits_{n'}J_{n'}\left(\frac{\varepsilon}{\omega}\right)\frac{\exp[\mathrm{i}(\frac{U_0}{\omega}-n')T_0]}{\mathrm{i}(\frac{U_0}{\omega}-n')},\nonumber\\ F_1(T_0)&=&\sum\limits_{n'}J_{n'}\left(\frac{\varepsilon}{\omega}\right)\frac{\exp[\mathrm{i}(\frac{U_0}{\omega}+n')T_0]}{\mathrm{i}(\frac{U_0}{\omega}+n')},\\ F_2(T_0)&=&\sum\limits_{n'\neq 0}J_{n'}\left(\frac{\varepsilon}{\omega}\right)\frac{\exp(\mathrm{i}n'T_0)}{\mathrm{i}n'}.\nonumber \end{eqnarray} \tag{ A.10 }$$

We note that the probabilities to find the two strongly interacting bosons in the same wells are constants in time up to the first-order time scale T₁ from equation (A.8). Therefore, we need to consider the asymptotic analysis up to the order σ². Following the same procedure outlined above, one has i∂a⁽²⁾_j/∂T₀ = −i∂A_j/∂T₂ + G⁽²⁾_j(T₀), with

$$\begin{eqnarray} &&\fl G_1^{(2)}=\sqrt{2}\mathrm{i}J_{0}\left(\frac{\varepsilon}{\omega}\right)A_5F_0(T_0)-\sqrt{2}\mathrm{i}(\sqrt{2}A_1F^*_0(T_0)+\sqrt{2}A_2F^*_1(T_0)+A_5F^*_2(T_0))\nonumber\\ &&\times \exp\left(\mathrm{i}\frac{U_0}{\omega}T_0-\mathrm{i}\varphi(T_0)\right),\nonumber\\ &&\fl G_2^{(2)}=\sqrt{2}\mathrm{i}J_{0}\left(\frac{\varepsilon}{\omega}\right)A_5(F_0(T_0)+F_1(T_0))-\sqrt{2}\mathrm{i}\left[\vphantom{A_5F_2^*(T_0))\exp\left(\mathrm{i}\frac{U_0}{\omega}T_0-\mathrm{i}\varphi(T_0)\right)}(\sqrt{2}A_1F^*_0(T_0)+A_5F_2^*(T_0)\right.\nonumber\\ &&+\sqrt{2}A_2F_1^*(T_0)) \exp\left(\mathrm{i}\frac{U_0}{\omega}T_0+\mathrm{i}\varphi(T_0)\right)+(\sqrt{2}A_2F_0^*(T_0)+\sqrt{2}A_3F_1^*(T_0)\nonumber\\ &&+\left.A_5F_2^*(T_0))\exp\left(\mathrm{i}\frac{U_0}{\omega}T_0-\mathrm{i}\varphi(T_0)\right)\right],\nonumber\\ &&\fl G_3^{(2)}=\sqrt{2}\mathrm{i}J_{0}\left(\frac{\varepsilon}{\omega}\right)A_5F_1(T_0)-\sqrt{2}\mathrm{i}(\sqrt{2}A_2F^*_0(T_0) +\sqrt{2}A_3F^*_1(T_0)+A_5F_2(T_0))\nonumber\\ &&\times \exp\left(\mathrm{i}\frac{U_0}{\omega}T_0+\mathrm{i}\varphi(T_0)\right),\nonumber\\ &&\fl G_4^{(2)}=\mathrm{i}J_{0}\left(\frac{\varepsilon}{\omega}\right)(A_4+A_6)F_2^*(T_0)-\mathrm{i}\left[2A_4F_0(T_0) \exp\left(-\mathrm{i}\frac{U_0}{\omega}T_0+\mathrm{i}\varphi(T_0)\right)+2(A_4F^*_1(T_0)\right.\nonumber\\ &&\left.+A_6F_0^*(T_0)) \exp\left(-\mathrm{i}\frac{U_0}{\omega}T_0-\mathrm{i}\varphi(T_0)\right)+(A_4F_2(T_0)+A_6F_2^*(T_0))\exp(-\mathrm{i}\varphi(T_0))\right],\nonumber\\ &&\fl G_5^{(2)}=\mathrm{i}J_{0}\left(\frac{\varepsilon}{\omega}\right)A_5(F_2(T_0)+F_2^*(T_0))-\mathrm{i}\Big[(\sqrt{2}A_1F_0^*(T_0)+\sqrt{2}A_2F_1^*(T_0)+A_5F_2^*(T_0))\nonumber\\ &&\times \exp(\mathrm{i}\varphi(T_0))+(\sqrt{2}A_2F_0^*(T_0)+\sqrt{2}A_3F_1^*(T_0)+A_5F_2(T_0))\exp(-\mathrm{i}\varphi(T_0))\Big],\nonumber\\ &&\fl G_6^{(2)}=\mathrm{i}J_{0}\left(\frac{\varepsilon}{\omega}\right)(A_4+A_6)F_2(T_0)-\mathrm{i}\left[2(A_4F_1(T_0)+A_6F_0(T_0)) \exp\left(-\mathrm{i}\frac{U_0}{\omega}T_0+\mathrm{i}\varphi(T_0)\right)\right.\nonumber\\ &&\left.+(A_4F_2(T_0)+A_6F_2^*(T_0))\exp(\mathrm{i}\varphi(T_0))+2A_6F_1(T_0)\exp\left(-\mathrm{i}\frac{U_0}{\omega}T_0-\mathrm{i}\varphi(T_0)\right)\!\right].\nonumber\\ \end{eqnarray} \tag{ A.11 }$$

Then the solvability condition at order σ² reads

$$\begin{eqnarray} \mathrm{i}\frac{\partial}{\partial T_2}A_1&=&\overline{G_1^{(2)}}=2(A_1\rho_1+A_2\rho_2),\nonumber\\ \mathrm{i}\frac{\partial}{\partial T_2}A_2&=&\overline{G_2^{(2)}}=2(2A_1\rho_1+(A_1+A_3)\rho_2),\nonumber\\ \mathrm{i}\frac{\partial}{\partial T_2}A_3&=&\overline{G_3^{(2)}}=2(A_3\rho_1+A_2\rho_2),\nonumber\\ \mathrm{i}\frac{\partial}{\partial T_2}A_4&=&\overline{G_4^{(2)}}=-2(2A_4\rho_1+A_6\rho_2),\nonumber\\ \mathrm{i}\frac{\partial}{\partial T_2}A_5&=&\overline{G_5^{(2)}}=0,\nonumber\\ \mathrm{i}\frac{\partial}{\partial T_2}A_6&=&\overline{G_6^{(2)}}=-2(2A_6\rho_1+A_4\rho_2), \end{eqnarray} \tag{ A.12 }$$

where we have set

$$\begin{equation} \rho_1=\sum\limits_{n'}\frac{J^2_{n'}(\frac{\varepsilon}{\omega})}{\frac{U_0}{\omega}+n'}, \quad \rho_2=\sum\limits_{n'}\frac{J_{n'}(\frac{\varepsilon}{\omega})J_{-n'}(\frac{\varepsilon}{\omega})}{\frac{U_0}{\omega}+n'} \end{equation} \tag{ A.13 }$$

for U₀/ω + n' ≠ 0. Thus the evolution of the amplitudes A_j up to the second-order long time is given by

$$\begin{equation} \frac{{\rm d}A_j}{{\rm d}t'}=\left(\frac{\partial}{\partial T_0}+\sigma\frac{\partial}{\partial T_1}+\sigma^2\frac{\partial}{\partial T_2}\right)A_j \end{equation} \tag{ A.14 }$$

from equations (A.3) and (A.4). The corresponding probability amplitudes read

$$\begin{equation} a_j(t')=A_j(t')+o(\sigma)+o(\sigma^2)+\cdots, \end{equation} \tag{ A.15 }$$

in the high-frequency regime, where the high-order small terms can be neglected. Substituting equations (A.4), (A.8) and (A.12) into equation (A.14), we obtain the two sets of coupled equations, equations (6) and (7).