The symmetry in the model of two coupled Kerr oscillators leads to simultaneous multi-photon transitions

We consider the model of two coupled oscillators with Kerr nonlinearities in the rotating-wave approximation. We demonstrate that for a certain set of parameters of the model, the multi-photon transitions occur between many pairs of the oscillator states simultaneously. Also, the position of the multi-photon resonances does not depend on the coupling strength between two oscillators. We prove rigorously that this is a consequence of a certain symmetry of the perturbation theory series for the model. In addition, we analyse the model in the quasi-classical limit by considering the dynamics of the pseudo-angular momentum. We identify the multi-photon transitions with the tunnelling transitions between the degenerate classical trajectories on the Bloch sphere.


INTRODUCTION
For decades, the models consisting of interacting nonlinear oscillator modes attract much attention because of their importance for various fundamental concepts of quantum physics, quantum information and nonlinear dynamics.The phenomena present in such models include quantum chaos [1], multi-stability [2], dissipative phase transitions [3], and dynamical tunnelling [4].Nonlinear oscillator networks were suggested as a framework for universal quantum computation [5], and also quantum nonlinear oscillators were suggested as a tool for non-classical states creation such as squeezed states [6], entangled states [7,8] and cat states [9].Also, such models are a fundamental tool to study quantum-classical correspondence [10].
The models of interacting nonlinear oscillators can be realised in various experimental setups including optomechanical systems [11], trapped ions [12] and superconducting circuits [13].The presence of a dynamical bifurcation point was experimentally demonstrated in RF-driven Josephson junctions [14] and transitions between two basins of attraction were observed in a nanomechanical resonator.Also, non-degenerate parametric amplifiers were created based on the Josephson junction arrays [13].
Lots of efforts are devoted to studying the nonlinear oscillator systems in the mesoscopic regime [3,10,15].In this regime, it is possible to use the quasi-classical approximation, but the quantum effects are important as well.The interplay of complex classical dynamics and quantum effects such as quantum tunnelling leads to new interesting phenomena.For example, it was shown that in the model of a single oscillator mode with Kerr nonlinearity, tunnelling affects the switching rate between the stable states [16].
The specifics of tunnelling in such systems are determined by the complex structure of the phase space.For the models exhibiting bi-or multi-stability, there exist several different classical trajectories with the same en-ergies.This opens way for tunnelling transitions between the corresponding quantum states.However, as the tunnelling transition amplitudes are exponentially small, the certain resonant condition should be usually satisfied in order to achieve a non-vanishing tunnelling probability [4].Also, the interconnection between tunnelling and multi-photon transitions was established in many systems [17][18][19][20].In particular, for a single nonlinear oscillator with Kerr nonlinearity, tunnelling between different regions of the phase space is equivalent to simultaneous absorption or emission of many oscillator quanta in the quasi-classical limit.
In this work, we consider the model of two coupled nonlinear oscillators in the rotating-wave approximation (RWA) [21].The classical limit of this model can be described as the dynamics of the pseudo-angular momentum on a two-dimensional sphere.Among the classical trajectories, there exist ones with equal energies, which makes tunnelling transitions between them possible.We show that tunnelling between such trajectories can be described with the perturbation theory in the coupling constant, and tunnelling transition is equivalent to the exchange by many quanta between the oscillators.We use high-order perturbation theory to study these processes and account for both resonant and non-resonant contributions.We identify the resonant condition for tunnelling between different classical trajectories and we prove rigorously that this resonant condition is satisfied simultaneously for many pairs of the oscillator states.Also, it is independent of the value of the coupling constant between oscillators, in contrast with many systems demonstrating similar behaviour (such as [22]).This fact is a consequence of an internal symmetry of the system Hamiltonian which we establish in all orders of the perturbation theory series.Also, we discuss how the presence of such symmetry modifies the energy spectrum and the structure of the eigenstates.These results could have applications for quantum information processing and quantum state manipulation, in particular, the generation of the entangled states of two arXiv:2206.05537v3[quant-ph] 21 Feb 2023 modes.We believe that the predicted features can be observed in systems with high-quality oscillator modes with a pronounced Kerr nonlinearity, such as the plasmon modes of Josephson junction arrays or the phonon modes of trapped ion ensembles.Both mentioned systems represent a natural realization of a system of coupled quantum nonlinear oscillators well isolated from the environment Also, the considered model is closely related to some models [23,24] of dissipative time crystals (DTC) [25,26].The results of this work could be useful to study the quantum effects in the DTC.In addition, the discovered symmetry extends the range of known symmetries in quantum-optical systems [27][28][29][30].

I. MULTI-PHOTON TRANSITIONS BETWEEN TWO COUPLED NONLINEAR OSCILLATORS
We consider the model of two oscillators with linear coupling and Kerr nonlinearities.Let the frequencies of the oscillators be ω 1,2 , and the Kerr shifts per oscillators quanta be α 1,2 .When the detuning between the oscillators ∆ = ω 2 − ω 1 is much smaller than the oscillators frequencies, one can neglect the counter-rotating terms in coupling (RWA) [21], and the model Hamiltonian reads where g is the coupling constant between the oscillators.As the counter-rotating terms are not present in (1), this Hamiltonian commutes with the total number of quanta operator N = â † â + b †b (we will denote its eigenvalues as N ).Therefore, the total Hilbert space H of the considered model ( 1) splits into the direct sum of invariant Hilbert subspaces each corresponding to N quanta in the oscillators: H N , and also Ĥ = N ĤN .The Hamiltonians ĤN which act on the subspace H N can be written with help of bra-ket notation as where |n a , n b = |n a ⊗ |n b are the oscillators Fock states (n a and n b are the numbers of quanta in modes a and b respectively), and The model ( 1) is closely related to the model of a single Kerr oscillator [31] in the classical external field.Namely, in the limit N → ∞, g √ N = const, the Hamiltonian ĤN approaches the Hamiltonian of a single Kerr oscillator in the classical external field with nonlinearity α 1 + α 2 , detuning αµ N /2, and the external field amplitude g √ N .This can be seen from Eq. (2).
In the following, we will focus on the case of relatively weak coupling.Moreover, in Sections II and III, we will use the perturbation theory in g as a theoretical tool.Because of that, let us first consider the model without perturbation.In this case, the Hamiltonian (1) commutes both with a † a and b † b, and its eigenstates are the Fock states |n a , n b of the oscillators.The energy of each state |n a , n b can be found easily from Eq. ( 1) The perturbation operator allows the oscillator to exchange quanta between each other.Namely, the perturbation directly couples |n a , n b and |n a ± 1, n b ∓ 1 .Because of that, the states |n a , n b and |n a + k, n a − k also become coupled in the order k by perturbation, and the transition amplitude between them is proportional to g k .It turns out that at certain values of the detuning between the oscillators, such transitions (we will call them multi-photon) become resonant, and even at small couplings, the oscillators can exchange by k quanta simultaneously between each other.
This occurs because the unperturbed energies of the oscillators Fock states |n, N −n in each subspace corresponding to N quanta have a parabolic dependence on n.The way how a sequence of the values  One can easily see that at integer µ N = m ∈ Z, many energies split into pairs with equal values: m−n,N −m+n for n = 0, . . .m (see Fig. 1(a)).In contrast, for non-integer µ N the unperturbed energies are nondegenerate (see Fig. 1(b)).Also, for the case α 1 = α 2 , the resonance condition is simultaneously satisfied for all N .
When µ N is close to an integer, the degeneracy makes multi-photon transitions possible between the states |n, N −n and |m − n, N − m + n when the perturbation is turned on.In other words, the integer values of µ N at small g correspond to multi-photon resonances between |n, N − n and |m − n, N − m + n .
However, the perturbation not only leads to multi-photon transitions between |n, N − n and |m − n, N − m + n , but also induces the energy shifts of these levels due to the non-resonant coupling with other levels.In general, these shifts could lift the position of the multi-photon resonances from exactly µ N ∈ Z, and to make them dependent on g.The surprising feature of the considered model is that these energy shifts are equal for every degenerate pair |n, N − n and |m − n, N − m + n (see Fig. 2).Because of that, simultaneous multi-photon resonances between many level pairs occur exactly at integer values of µ N even at relatively large values of g.
In the absence of degeneracy in the unperturbed model (for the case of non-integer µ N ), weak coupling between the oscillators leads only to a small O(g) perturbation of the Fock states of the system: the eigenstates remain close to Fock states.This is not the case when the multi-photon resonance condition is satisfied.Due to the reasons indicated above, even at small couplings the structure of the Hamiltonian eigenstates and the temporal dynamics of the system wave function for integer µ N qualitatively differ from the case of non-interacting oscillators.
In this case, the leading-order contributions to the wave functions are symmetric and antisymmetric superpositions of the Fock states |n, N − n and |m − n, N − m + n .For every n = 0, . . ., m, This can be seen from the results of the numerical diagonalization of the Hamiltonian (1).We find the eigenstates of the Hamiltonian (1) in the Fock basis in the form |ψ = n c ,n |n, N − n .The eigenstates are sorted in the ascending order according to their eigenenergies.Then, we plot the matrices of the expansion coefficients c ,n in Fig. 3 for the cases of integer (Fig. 3a) and non-integer µ N (Fig. 3b).For each eigenstates shown in Fig. 3a, there is a single Fock state which has dominant contribution, whereas the eigenstates in Fig. 3b contain equal contributions from two Fock states as in Eq. (5).
The energy difference between the states |ψ ± n,m−n exhibiting the anti-crossing is determined by the multi-photon transition amplitude ω R n,m−n ∝ g m−2n (see Fig. 4)   In addition, from Eq. ( 5) and ( 6), for the initial state: |ψ(0) = n c n |n, N − n , one can find (at integer µ N = m) the following approximate solutions of the non-stationary Schrodinger equation  As can be seen, the system exhibits multi-photon Rabi oscillations between many pairs of the Fock states simultaneously.
In Fig. 5, we show the numerical solution of the Schrodinger equation for different Fock states taken as initial conditions.We plotted the squares of modulus of the overlappings between the wave function |ψ(t) and the bra-states n, N − n|, m − n, N − m + n| for different n.This result is in agreement with Eq. ( 7) for integer µ N .Also, for the case of non-integer µ N = m+δµ N with sufficiently small δµ N , the amplitude of multi-photon Rabi oscillations between |n, N −n and |m−n, N −m+n decreases with increasing δµ N .They completely vanish when (α (see left and right panels of Fig. 5).In addition, there are corrections to Eq. ( 7) which come from the non-resonant contributions of the adjacent Fock states.They lead to the additional modulations of the Rabi oscillations with the amplitude ∼ g/∆ and become more pronounced with increasing n due to the dependence of the transition matrix elements and the energy differences on n.They can be seen on the lower panels of Fig. 5.For other panels, they are also present but not resolved on the plots.Also, let us briefly discuss the effect of weak dissipation on the considered multi-photon Rabi oscillations.In the presence of dissipation, each eigenstate of the system of coupled oscillators obtains a finite lifetime having the same order of magnitude as the decay rates γ a , γ b of the individual bosonic modes [32].The condition necessary to observe the Rabi oscillations between the states |n, N − n and |m − n, N − m + n is that the lifetime of the eigenstates |ψ ± n,m−n is much longer than the period of Rabi oscillations ω R n,m−n , which leads to conditions γ a,b ω R n,m−n .The requirements to observe the multi-photon transitions between two oscillator states become increasingly strict with the increasing order of the multi-photon transition.Thus, in the presence of dissipation, it is realistic to observe several transitions between the states |n, N − n and |m − n, N − m + n with the lowest values of |m − 2n|.
The independence of the anti-crossings position on the value of g can have practical applications for the measurements of the Kerr coefficients of the oscillators.More precisely, assume the detuning ∆ between the oscillators is scanned to measure some observable.Often, the multi-photon resonances lead to peaks/dips in the corresponding dependencies on ∆.According to the above considerations, the positions of these peaks correspond to the integer multiples of the total nonlinearity independently of the coupling constant.This could be used for the measurement of α 1 + α 2 .
In addition, the multi-photon transitions provide a way to prepare the entangled states of two oscillators.For example, according to Eq. ( 7), the initial Fock state |n, N − n evolves into an entangled state at the proper choice of the interaction time (when ω R n,m−n t ∼ π/4).

II. THE PROOF OF THE SYMMETRY PROPERTIES
In this section, we provide a proof (on the physical level of rigour) of the presence of many simultaneous anticrossings at µ N ∈ Z for the model (1).
Let us consider the coupling term between the oscillators as a perturbation.As commented above, at integer µ N = m, the unperturbed energies split into m/2 pairs of degenerate levels.The perturbation leads to two effects: the shifts of the energy levels due to the non-resonant couplings, and the multi-photon transitions between the degenerate levels.The consistent analysis of both of the effects requires the usage of high-order perturbation theory.
Because of the degeneracies, one can not use directly the non-degenerate perturbation theory series in g for the system energies.However, at non-integer µ N one can define the energies n (g) which correspond to the eigenstates evolving from the n-th Fock state |n, N − n after the adiabatic switching of the perturbation.These energies (at non-integer µ N ) can be found from non-degenerate perturbation theory as power series (The expressions for the second and fourth order corrections could be found in Appendix A).
To some extent, this expansion is valid for µ N ∈ Z as well.Although at µ N ∈ Z the states |n, N − n and |m − n, N − m + n are degenerate, the perturbation couples them only in the order |m − 2n|.Therefore, up to this order of the perturbation theory, they can be treated as non-degenerate ones.Also, the perturbation theory coefficients n are rational functions of n and µ N , therefore, they are well-behaved at µ N ∈ Z. So, one can use the expansion (8) up to the order |m − 2n| in the degenerate case.Surprisingly, the following identity is valid for the perturbation theory corrections at µ N = m ∈ Z: Therefore, the degeneracy between the levels |n, N − n and |m − n, N − n + m is not lifted in several lowest orders of the perturbation theory until the order |m − 2n|.This explains the Eq. ( 6) and the results of the numerical diagonalization.
Strictly speaking, the identity (9) makes sense only for the case of the integer µ N = m.However, as the corrections n are rational functions of n and µ N , they can be analytically continued to the case of arbitrary real n and µ N .We will prove (9) for these analytical continuations.
The analytical continuation of n for the case of non-integer n can be made in a natural way for all k.For that, the whole energy spectrum n (g) should be continued to the case of non-integer n.To do this, let us assume that the Hamiltonian (1) acts on the space of all possible real «numbers of quanta» ν with formally defined matrix elements as and analogously for b and b † .As one can see, in the case of an integer ν, the definition of (10) reduces to the usual action of the bosonic creation and annihilation operators in the Fock space.We shall note that â and â † are no longer Hermitian conjugate with each other for the case of non-integer ν, so the Hamiltonian becomes non-Hermitian.However, one can check that the wave functions remain normalizable even for the Hamiltonian at the non-integer ν.
The extended Hamiltonian acts invariantly on each of the subspaces V ν = {|σ, N − σ : σ − ν ∈ Z}.(The subspace V ν consists of the ket vectors |σ such as σ − ν ∈ Z.) Also, the subspaces of states |ν, N − ν with negative integer ν, ν ∈ [0, N ] and integer ν > N are decoupled from each other.So, the extended Hamiltonian for all real ν acts exactly as the original Hamiltonian (1) on the subspace of |n, N − n , n ∈ [0, N ].We emphasize that the extension to the non-integer numbers of quanta should be taken as an auxiliary tool for the proof.Because, as we mentioned above, we are interested in the case of an integer ν which corresponds to the initial Hermitian Hamiltonian.
Until the end of this section, we will work in the subspace of a fixed number of quanta N .We will denote |ν, N − ν as |ν for brevity.
Let us define the action of the extended Hamiltonian on V ν as H ν,N .The operator H ν,N can be written in terms of the states |ν as To prove the symmetry of the perturbation theory corrections to the energy, we will show the exact symmetry of the eigenstates of the extended Hamiltonian with respect to the replacement ν → µ N − ν.Namely, we will prove that H ν,N and H µ N −ν,N have the same eigenenergies ν (g) and µ N −ν (g) respectively as a functions of g.For that, we will show that there exists a linear operator T such where I is the isomorphism between the vector spaces spanned by vectors |ν and |µ N − ν respectively.The existence of the operator T proves that the operators H ν,N and H µ N −ν,N have identical spectra at any g, which means that In Appendix B, we explicitly construct the operator T and show that it can be expressed in the following form: T = U T V −1 , where After we proved the equality of the energies (13), let us turn back to the physically meaningful case of the integer ν and µ N = m.As we mentioned before, the energy ν (g) can be decomposed in perturbation theory series of g (compare with Eq. ( 8)): From the equality ( 13) valid for all orders in g, one concludes that this holds for each term of the perturbation series: From Eq. ( 16), the equality ( 9) can be derived easily.The k-th order corrections are rational functions of µ N and ν and have no singularities at 2ν − µ N ∈ Z unless k 2|µ N − 2ν| (we prove it rigorously in the next sections).Therefore, the equality ( 16) holds even for the case of an integer 2ν − µ N , and also for ν, µ N ∈ Z.This concludes the proof of the equality (9).
We should note that the analogous properties hold for the model of a single nonlinear mode with Kerr nonlinearity driven by the classical external field (recovered as a limit of the considered model at N → ∞).For this model, the equality of the perturbation theory corrections was checked [33] for several lowest orders, and the sketch of the proof was given previously [34].

III. HIGH-ORDER PERTURBATION THEORY FOR THE DEGENERATE ENERGY LEVELS
We have stated that non-degenerate perturbation theory corrections are symmetric with respect to replacement n → µ N − n in all orders.However, additional arguments are needed to relate this result to the physical case of µ N ∈ Z.On one hand, due to the presence of degeneracy in the energy spectrum, one should use the degenerate perturbation theory.However, on the other hand, the off-diagonal matrix elements in the secular equation occur only in the |m − 2n|-th order of perturbation theory.To demonstrate that until the |m − 2n|-th order one can use the non-degenerate perturbation theory and to calculate the multi-photon amplitude via degenerate perturbation theory for higher orders than |m−2n|, it is convenient to apply Green's function formalism.Namely, we consider the operator of Green's function defined as If the spectrum of the problem σ( Ĥ) = { n , |ψ n } is known, Eq. ( 17) can be rewritten as follows Thus, eigenenergies are poles of Green's function and eigenfunctions could be calculated from residues of Ĝ.The matrix element G n,n for each Fock state |n can be calculated from the Dyson equation and reads where Σ n is the self-energy term, defined as the sum of all diagrams which start and finish at n and do not contain the Green's function G n,n .When Σ n is the regular function in the vicinity of ω = (0) n , the position of the pole of G n,n can be found from the equation G −1 n,n (ω) = 0 as power series in g which coincides with the result of usual non-degenerate perturbation theory expansion.However, if there exists another level with the same or close energy (this is |m − n, N − m + n in our case), the self-energy itself acquires a pole in the vicinity of (0) n .The lowest-order perturbation theory term with a pole has the order 2|m − 2n| in g because the corresponding diagram must contain at least one G m−n,m−n Green's function.As the perturbation operator changes the number of quanta by one, at least 2|m − 2n| perturbation vertices are required.
Because of the pole in the self-energy term in the vicinity of (0) n , the Green's function poles positions cannot be found as simple perturbation series.For this case, it is convenient to consider the matrix Green's function In terms of matrix Green's function, Dyson equation could also be written with self-energy matrix The solution of Eq. ( 21) reads The oscillator eigenenergies are the poles of the Green's function.They could be found from the equation det G −1 = 0, and all their dependence on the perturbation is contained in the self-energy matrix Σ.The diagonal terms Σ n,n and Σ m−n,m−n correspond to non-resonant energy shifts and contain contributions of all even powers of g.In contrast, the off-diagonal term Σ n,m−n is proportional to g m−2n and is responsible for resonant multi-photon transition between the states |n, N − n and |m − n, N − m + n .The leading term of its expansion in powers of g reads (see Appendix C) Eqs. ( 21)-( 23) explain the possibility to use the non-degenerate perturbation theory up to the |m − 2n|-th order.According to Eq. ( 23), the off-diagonal terms in (22) (which correspond to the multi-photon resonance), have the order of g m−2n .Therefore, they do not contribute to the k-th order of the perturbation theory when k < |m − 2n|.For k < |m − 2n|, the energies obtained from the secular equation will coincide with those obtained from the nondegenerate perturbation theory.To account for multi-photon transitions, one should consider the perturbation theory of order k |m − 2n|.n,m−n .Further, we will demonstrate that it can be treated as the tunneling amplitude in the quasiclassical limit.

IV. CLASSICAL LIMIT
Multi-photon transitions described above with help of the formalism of the perturbation theory also have a quasiclassical interpretation as tunnelling transitions.For that, one should consider the classical limit of the studied model which is valid for the large number of quanta: N 1, µ N 1, N/µ N = const.To obtain the classical Hamiltonian of the system, one should replace the ladder operators in (1) with classical complex amplitudes [21,35].This results in the following complex Hamilton function Due to the conservation of the total number of quanta N = |a| 2 + |b| 2 , the classical dynamics governed by this Hamiltonian can be described as the dynamics on the surface of the two-dimensional sphere.To show that, one should rewrite the classical Hamiltonian with help of the new pseudo angular momentum variables L x , L y , L z defined as In terms of the components of L, the classical Hamiltonian (24) takes the form The conservation of the total angular momentum L 2 = L 2 x + L 2 y + L 2 z = N (N + 2)/4 (which is equivalent to the conservation of the total number of quanta in the oscillators) allows describing the classical dynamics with help of the classical phase portrait on the Bloch sphere (see Fig. 6): the trajectories in the L i space are defined by the conservation of the number of quanta and the Hamilton function (26).
Due to the quantum-classical correspondence, the eigenstates of the quantum Hamiltonian (1) correspond to a discrete set of classical trajectories on the Bloch sphere.We demonstrate this by comparing the period-averaged values of the classical momenta L x , L z with the quantum-mechanical averages over the eigenstates of (1) for noninteger µ N in Fig. 7 (the average of L y equals zero).As one can see from Fig. 7, the averages calculated from the classical model are close to the quantum averages everywhere except the vicinity of the separatrix.
Let us discuss the structure of the phase portrait in detail.At g = 0, the Hamiltonian is a function of L z only, therefore, the classical trajectories are the circles in the L x , L y plane.At non-zero g, the trajectories are no longer concentric circles.Also, a self-intersecting trajectory (separatrix) emerges which divides the phase portrait into three regions with a stable point inside each one (denoted as «1», «2» and «3», see Appendix D).At a certain critical value of g = g crit , a saddle-node bifurcation occurs, and one of the stable point merges with the unstable stationary state «S».At larger g, only two stable points remain.Depending on the value of the ratio N/µ N , the unstable point «S» can merge with the point «1» or with the point «3».If N/µ N > 1, the point «3» merges with the point «S», otherwise, the point «1» merges with point «S».The value of g crit can be found by analysing the positions of the extrema of the Hamiltonian (26) as a function of L i .After casting the coupling constant to a dimensionless form β = 8g 2 N /(α 1 + α 2 ) 2 µ 3 N (the dimensionless coupling strength), we obtained the following expression for β crit (see derivation in Supplementary Information): where γ = N/µ N .Also, in the limit of γ → ∞, one recovers the result for the model of a single Kerr oscillator with classical driving: at γ → ∞, β crit → 4/27 [36,37].The dependence of the positions of the stable states on the dimensionless coupling strength β is presented in Fig. 8.For each of the equilibrium points i ∈ {1, 2, 3, S}, the value of the polar angle is plotted (see Supplementary Information for more details).At β = β crit (vertical black dotted line), the states «S» (unstable) and «3» (stable) merge.Also, both of the angles ϑ 1 and ϑ 2 approach π/2 at β β crit , which corresponds to the diametrically opposite stationary points in the X-Y plane.
as a function of √ β for the i-th equilibrium point, i ∈ {1, 2, 3, S}.Black dashed line corresponds to β = βcrit.The points numbers correspond to the ones in Fig. 6.Here the parameters are as follows: N = 40, α2/α1 = 0.5, ∆/α1 = 2.5, µN = 30, and g = β(∆ + α2N ) 3 /N (α1 + α2).Now let us discuss how tunnelling modifies the classical picture and establish the relation between tunnelling and multi-photon transitions.Tunnelling transitions are possible providing that there exist classical trajectories with the same energies.Because of the structure of the classical phase portrait, this holds for g < g crit .In this case, tunnelling transitions are possible between the classical trajectories from regions 1 and 3.At small g, the classical trajectories are close to the circles in the X-Y plane (see Fig. 6, left panel) and can be identified with the oscillators Fock states: the Fock state |n, N − n corresponds to the circle with L 2 = N 2 /4, L z = N/2 − n.For this case, one can directly apply the results of Sections I and II and deduce that at the integer µ N , resonant tunnelling transitions are possible between the trajectories with L z = N/2 − n and L z = N/2 − µ N + n.For larger coupling values, the eigenstates are no longer the Fock states, and the classical trajectories are no longer circles in the X-Y plane.However, we state that the condition for resonant tunnelling remains the same for all values of g ∈ [0, g crit ].This is supported by the fact that the calculation of section II is performed in all orders of the perturbation theory.Because of that, it fully takes into account the modification of the classical trajectories at larger g.In addition, the numerical diagonalization of the Hamiltonian (1) demonstrates that the behaviour shown in Fig. 2 (simultaneous anticrossings at the integer µ N ) persists for the whole range g ∈ [0, g crit ].

CONCLUSIONS
For the model of two coupled quantum oscillators with Kerr nonlinearities and linear coupling, we studied the multi-photon transitions between the oscillators.We showed that for certain parameters of the model, the resonant condition for multi-photon transitions is simultaneously satisfied for many pairs of the oscillator Fock states.This holds even for the moderate coupling strength between oscillators, and this is the consequence of a special symmetry of the oscillators Hamiltonian.
The latter is related to the structure of the perturbation series of the model eigenenergies and was proven with help of analytical continuation of the Hamiltonian for non-integer numbers of quanta.
Also, in the quasi-classical limit, the phase space of the two coupled oscillators can be mapped on a sphere, and the multi-photon transitions can be interpreted as tunnelling transitions between the trajectories lying in different regions of the classical phase space.Thus, when the resonant condition is satisfied, tunnel transitions affect the whole region of the classical phase space.
We believe that the results obtained in this work could be relevant for the experiments involving high-quality oscillator modes with low occupation numbers, such as the plasmon modes of the Josephson junction arrays or phonon modes of trapped ions ensembles.In particular, the independence of the multi-photon resonances positions could be used for the measurements of the Kerr coefficients of the modes.Also, the multi-photon transitions provide a way to create the entangled states of two oscillators.In addition, the obtained results could be used for certain models of dissipative time crystals in the quantum regime, and the symmetry discovered in the considered model may allow obtaining new exact results in quantum-optical systems.indeed obey the identity (B7).Also, let us note that T . (B10) are arranged on a parabola depends on the value of the parameter µ N .

Figure 3 .
Figure 3.The expansion coefficients c ,n of the eigenstates of the Hamiltonian (1) at N = 12, α2/α1 = 1, g/α1 ≈ 0.03, (a) ∆/α1 = 6.25 and (b) ∆/α1 = 6.5 are shown.For each pair n, , a unit square centered at the corresponding point of the figure is drawn.The color of the square indicates the magnitude of c n, according to the color bar.
The term Σ n,m−n ω = (0) n can be interpreted as the multi-photon transition amplitude between |n, N − n and |m − n, N − m + n , and Σ n,m−n (0) n equals the frequency of the multi-photon Rabi transitions ω R

Figure 7 .
Figure 7.The average values of the pseudo angular momenta components Lx and Lz over the classical trajectories with the dimensionless energy ˜ (orange dashed and green solid line) and over the eigenstates of the quantum Hamiltonian with the eigenenergies ˜ (red and blue circles).The averages are calculated at the following values of the Hamiltonian parameters: N = 40, α2/α1 = 0.5, ∆/α1 = 0.25, µN = 27, and g/α1 ≈ 1.816.