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

We investigate quantum tunneling of two repulsive bosons in a triple-well potential subject to a high-frequency driving field. By means of the multiple-time-scale asymptotic analysis, we evidence a far-resonant strongly-interacting regime in which the selected coherent destruction of tunneling can occur between the paired states and unpaired states, and the dominant tunneling of the paired states is a second order process. Two Floquet quasienergy bands of the both kinds of states are given analytically, where a fine structure up to the second order corrections is displayed. The analytical results are confirmed numerically based on the exact model, and may be particularly relevant to controlling correlated tunneling in experiments.


Introduction
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 timeperiodic 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][8][9][10] and photon-assisted tunneling [11][12][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][18][19][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][25][26][27][28][29][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 secondorder 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][34][35][36][37][38][39][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 triplewell system have also attracted extensive attention [41][42][43]; however, research on the secondorder effect of the system has not yet been reported.
In this paper, we investigate the coherent control of second-order tunneling for two triplewell confined bosons driven by a high-frequency ac field. By means of the multiple-timescale 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].

The model and high-frequency approximation
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][42][43] where the operatorâ ( †) l annihilates (creates) a boson in well l; J denotes the nearest-neighbor hopping matrix element; U 0 is the on-site interaction energy; and ε(t) = ε cos(ωt) is the ac driving of amplitude ε and frequency ω. Throughout this paper, we takeh = 1 and use the parameter J to appropriately scale the energy, U 0 , ε 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 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 6 j=1 |c j (t)| 2 =1. Inserting equations (1) and (2) into Schrödinger equation i∂ t |ψ(t) =Ĥ (t)|ψ(t) , obtains the coupled equations for the amplitudes c j (t) 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 0 = mω + u for |u| ω/2, m = 0, 1, 2, . . . with u being the reduced interaction strength [31], and make the function transformations c 1 , c 5 (t) = a 5 (t) and c 6 (t) = a 6 (t) exp[−iϕ(t)], with a j (t) being the slowly varying functions and ϕ(t) = t 0 ε cos(ωτ ) dτ = ε ω sin(ω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] 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 multipletime-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 selfinteraction intensity. In this work, we are not concerned about the very strong interaction (e.g. U 0 6ω) since for such a interaction we need to consider not only the usual on-site atominteraction strength but also the interactions between atoms on neighboring lattice sites [17], which is beyond the considered case. When the condition J 0 = 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ȧ 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 | 2 = |a j | 2 ( j = 1, 2, . . . , 6) for the resonant case with J = 1, U 0 = ω = 80, ε/ω = 2.405 and P 2 (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.
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 1 (t), a 2 (t) and a 3 (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.

Second-order tunneling in the far-off-resonant strongly interacting regime
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 σ 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 (0) . . , 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: where ρ i (i = 1, 2) are set as (see the appendix) for U 0 /ω + 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 σ −1 , such that σ 2 ρ 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 | < σ −1 implies that the inequality 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 0 , where the second-order perturbation method is more perfect for |u| 10 J = ω/8 with U 0 = ω + u and for |u| 5 J = ω/16 with U 0 = 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 σ 2 . 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 σ 2 ρ 2 , 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 ρ 2 of the second-order tunneling coefficients as a function of the driving parameters ε/ω and self-interaction intensity U 0 /ω, where figure 2(b) is the plan view of figure 2(a). Combining equation (8) with figure 2 we can see that the factor ρ 2 tends to infinity for any integer value of U 0 /ω 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 0 /ω 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 ρ 2 = 0 is satisfied.
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 2 (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 Inserting these expressions into equation (6) yields the coupled equations Combining equations (9) with (11) produces (10) yields the decoupled equation
In (a), the dashed line corresponds to the probabilities P 1,3 , and the thin and the thick solid lines indicate the probabilities P 2 and P 4,5,6 , respectively. In (b) the dashed line corresponds to the probabilities P 4,6 , and the thin and the thick solid lines indicate the probabilities P 5 and P 1,2,3 , respectively. In this figure and in the following figures, all of the circular points indicate the analytical solutions and the curves represent the numerical results.
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.
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 S = P 1 + P 2 + P 3 = 1 τ τ 0 (P 1 + P 2 + P 3 ) dt for τ = 200(J −1 ). The normalization means the time-averaged total probability in the three unpaired states being 1 − S . Taking the initial conditions P 2 (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 0 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 0 = 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|. for U 0 = ω + u, and |u 2 | 5J = ω/16 for U 0 = 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 0 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 (ε/ω) = 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 0 /ω of the interaction and the driving frequency in figure 5(b) with ε = 160 and U 0 = 200. Figure 5(b) explicitly shows that the lowest points of S appear at the resonant points U 0 /ω = 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 0 ∈ [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. 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:

Floquet quasienergy analysis
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 k t |φ 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 where T is the time-ordering operator and T = 2π/ω is the period of the driving. Noticing that eigenvalues of U (T, 0) are exp(−iE k T ), 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 0 = 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 ( ε ω ) = 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 ( ε ω ) 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 0 = ω = 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 firstorder Bessel function J 1 (ε/ω). From equation (4) with m = 1, u = 0 we know that the first derivative of the probability amplitudes of the three doublonsȧ j ( j = 1, 2, 3) are equal to zero and the three unpaired statesȧ j ( j = 4, 5, 6) do not vanish at the zero points of J 1 (ε/ω), 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 = 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 0 = ω. 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).
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 0 = 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 (ε/ω) = 0. Precise agreements between the numerical results based on equation (3)  (19), (20) and (23), the thin dotted lines indicate the zero points of J 0 (ε/ω) and the arrow in 7(b) indicates the amplification position of the inset. 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| 1 ≈ 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| 1 to |u| = 30 > |u| 1 .
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.

Avoided level-crossing of unpaired states
In the far-off-resonant strongly interacting regime, selecting the parameters as J = 1, ω = 80 and U 0 = 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 (ε/ω) = 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[−iϕ(t)] and exp[±2iϕ(t)] is T . Therefore, we can construct the Floquet states by setting [35] A j (t) = B j exp(−iEt) ( j = 4, 5, 6) for the three unpaired states with constant B j , then rewriting equation (7) as the time-independent form The existence condition for the nontrivial solution of equation (17) reads From equation (18) we obtain three Floquet quasienergies corresponding to the three unpaired states where we have set 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).

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 (ε/ω) = 0 and the interaction intensity obeys U 0 /ω = 2.58 corresponding to ρ 2 = 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 ρ 2 = 0 is satisfied. From equation (8), we have ρ 2 = 0 at ε/ω ≈ 0.95 in the region ε/ω ∈ [0, 8] for U 0 /ω = 1.6, and have ρ 2 = 0 at ε/ω ≈ 1.20, 2.02, 5.52 and 5.74 in the same region for U 0 /ω = 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 ρ 2 = 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).
We then analytically calculate the quasienergies corresponding to the three doublons. Substituting the Floquet solutions A j (t) = B j exp[−i(E − U 0 )t] ( j = 1, 2, 3) into equation (6), we obtain easily with the existence condition of the nontrivial solution From this equation we obtain three Floquet quasienergies corresponding to the three doublons as 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 0 ). 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 0 = 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 0 /ω = 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).

Conclusion and discussion
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 0 ∈ [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.
For the convenience of our discussion, we simplify equation (A.5) as i∂a (1) j /∂ T 0 = −i∂ A j /∂ T 1 + G (1) j (T 0 ) for j = 1, 2, . . . , 6. To avoid the occurrence of secular growing terms in the solution a (1) j , the solvability condition [30,32] i ∂ A j ∂ T 1 = G (1) j (T 0 ) (A.6) must be satisfied, where the overline denotes the time average with respect to the fast time variable T 0 (i.e. the dc component of the driving term G (1) j (T 0 )). The amplitudes a j at order σ are given by Employing equations (A.5) and (A.6), one gives So the solutions of order σ read , (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 1 from equation (A.8). Therefore, we need to consider the asymptotic analysis up to the order σ 2 . Following the same procedure outlined where we have set for U 0 /ω + n = 0. Thus the evolution of the amplitudes A j up to the second-order long time is given by , we obtain the two sets of coupled equations, equations (6) and (7).