Emergence of chaos and controlled photon transfer in a cavity-QED network

We develop optimal protocols for efficient photon transfer in a cavity-QED network. This is executed through stimulated Raman adiabatic passage (STIRAP) scheme, where time-varying capacitive couplings (with carefully chosen sweep rate) play a key role. We work in a regime where semiclassical limit is valid and we investigate the dynamical chaos caused by the light-matter coupling. We show that this plays a crucial role in estimating the lower bound on the sweep rate for ensuring efficient photon transfer. We present Hermitian as well as an open quantum system extension of the model. Without loss of generality, we study the three cavity and four cavity case and our results can be adapted to larger networks. Our analysis is also significant in designing transport protocols aimed for nonlinear open quantum systems in general.

Population transport through a nonlinear network [such as multi-mode Bose-Hubbard (BH) systems] results in intricate physics of various types of instabilities 25,26 . Apart from energetic instability 26 due to nonlinear eigenstates (of the problem in the semi-classical limit), chaos can play major role in determining transfer efficiency 25 . Therefore, a judicious control of system parameters is crucial to tackle such sensitive physical processes. Nonlinear STIRAP consisting of interacting atomic Bose-Einstein (BEC) condensates has been analyzed semiclassically 25 as well as in a quantum manybody framework 2 , and the role of various instabilities have been investigated theoretically. It has been shown that the adiabatic conditions for such processes get modified due to emergence of chaos 25 . Another platform to investigate nonlinearity is a c-QED lattice. This platform precisely implements the Jaynes-Cummings nonlinearity which is very different from the BH nonlinearity in atomic BEC. In addition to this, a dispersive regime of a c-QED can mimic the BH nonlinearity (Kerr type). Therefore, a c-QED lattice, being an efficient quantum simulator, demands an extensive analysis of nonlinear transport. Chaotic signature in systems where a single cavity is involved 1,28,30-32 has been investigated and such systems can be considered to be a good testing bed * amit.dey@icts.res.in † manas.kulkarni@icts.res.in for quantum-classical correspondence 1 of chaos. In a linear trimer of cavities 33 , control of non-directed (unlike STIRAP scheme) single photon transfer is proposed by tuning the ratio of inter-cavity tunnelings in ultrastrong light-matter coupling regime. Although nonlinear contribution is studied for adiabatic light passage in terms of excitation power dependence 34 , to the best of our knowledge, role of chaos in these optical processes remained elusive so far. The Jaynes-Cumming (JC) interactioninduced nonlinearity is exploited in coupled c-QED systems and delocalized-localized phases have already been realized 15,16,35 . Furthermore, driven-dissipative preparation of exotic steady states in extended cavity systems paved the avenue of controlling photon propagation in scaled-up architectures 15 . Therefore, a deeper understanding of aspects of nonlinear dynamics (such as efficient photon transfer) of these systems will significantly add to the existing control strategies and it is much needed to open up myriad of technological applications 21 . Developing such protocols warrants a deep understanding of nonlinear systems and subsequently bringing in important notions (for e..g., chaos) can play a paramount role in engineering the systems to ensure efficient transfer.
In this paper, we investigate a c-QED based STIRAP and show that dynamical chaos sets the lower bound for the sweep rate (which quantifies how fast one tunes the coupling strength), resulting in efficient photon transfer. Without loss of generality, we study the case of three and four cavities, and by efficient photon transfer, we mean, a nearly 100%, transfer of photons from the first cavity to the last cavity with almost no occupation of the intermediate cavities during the time evolution. Quantifying chaos by Lyapunov exponent (LE) in the semiclassical limit, we make connection with the sweep rate. This sets the lower bound on the sweep rate for the tuning parameters and thereby helps in achieving a strategy to ensure nearly 100% transfer of photons in an interacting/nonlinear system. Such a successful strategy for interacting systems lays a strong foundation for (i) establishing photon mediated communication by minimizing dissipation in a quantum network (because intermediate cavities are essentially empty in the process) and (ii) qubit-state transfer and its readout at the terminal cavity.
Starting from the model Hamiltonian, we work out the semiclassical equations of motion as we are in a regime where it is valid. We obtain the stationary point (SP) solutions for the Hermitian problem at every stage of sweep (see J 1 , J 2 sweep in Fig. 1). These can be typically multi-valued but we track a special SP (SSP) branch that leads to a near perfect transfer from the first to the last cavity with negligible content in the intermediate cavities. We present the STIRAP time dynamics for various sweep rates and analyze chaotic effects. For SSP branch, we present LE analysis at different sweep stages. This characterises the chaotic aspects of the system. The results after including the inevitable presence of dissipation in experiments has also been discussed in supplementary material 46 . Although much of our analysis relies on a semi-classical approximation, we have successfully demonstrated the consistency with an exact quantum calculation (see supplementary material 46 ). We then summarise our findings and discuss the future outlook. Model and dynamical equations: The c-QED STIRAP given by the schematic Fig. 1 can be described by the time-dependent Hamiltonian given by, whereĤ a = ω aâ †â + Ω aŝ z a + g a (â †ŝ− a + h.c.) describes the Jaynes-Cummings Hamiltonian for the cavity labeled by 'a' (cavity-a).â destroys a photon with frequency ω a in cavity-a, g a denotes light-matter (photon-qubit) interaction in cavity-a. We definen a,b,c as the photon number operator in cavity-a, cavity-b, cavity-c re-spectively. The two-level system with transition frequency Ω a,b,c is described by the spin operatorsŝ α a (where α ≡ {x, y, z}). The time-dependent couplings are defined as Gaussian pulses J 1,2 (t) = Ke − t−t 1,2 τ 2 (with the sequence t 1 > t 2 ), where τ is the pulse width and therefore 1/τ is the sweep rate (measured henceforth in units of K). Therefore, the Hamiltonian (Eq. 1) is explicitly time-dependent in a rescaled timet ≡ t/τ . Throughout the paper we use t 1 /τ = 3.697, t 2 /τ = 2.4242 and K = 1. The evolution under the time-dependent Hamiltonian (Eq. 1) depends on the sweep rate of J 1,2 (t). The STIRAP scheme can be implemented, for e.g., in coupled optical waveguides through variation of the spacings between the central and terminal waveguides 34,[39][40][41][42][43] . Similar variations of couplings in time domain, is also, in principle implementable 44,45 . Although there is no nonlinear photon-photon interaction in the Hamiltonian, perturbative treatment in dispersive regime shows that the lightmatter interaction induces such interaction 37 . However, we deal with a resonant situation (ω a,b,c = Ω a,b,c ) where the effect of anharmonicity is most manifest. For a linear system (with g a,b,c = 0), Eq. (1) is just a standard STIRAP Hamiltonian, which ensures the ex- Here, |A (|C ) is the state vector when the total population resides at cavity-a (cavity-c). This particular eigenstate does not project on to |B and acts as a 'dark state' 38 , i.e., B|Ψ 0 = 0. However, as Θ varies from 0 to π/2 (i.e., a complete sweep of J 1,2 , see Fig. 1), the state vector |Ψ 0 rotates from |A to |C resulting population transfer. We define transfer efficiency as T = nc(t→∞) na(t=0) . A semiclassical analog of an eigenstate is a solution of semi-classical equations of motion that do not evolve (i.e., a SP solution). Since, we will work in a regime where semi-classical approximation is valid, the analog of the dark state mentioned above is a SSP branch that leads to a perfect transfer from cavitya to cavity-c (with zero content in cavity-b throughout the time dynamics). Therefore, for a successful adiabatic passage, we sweep the couplings J 1 , J 2 slow enough (i.e, 1/τ is relatively small), so that the time evolution of the quantities follows the SSP branch.
In this paper, we deal with a more complicated nonlinear situation and seek similar photon transfer mechanism. In the large photon number limit, we exploit the approximation âŝ + a ≈ â ŝ + a and treat the problem in a semiclassical framework 35,36 . For brevity, we use the following notations: â → a, ŝ − a → s a , and ŝ z a → s z a . Using Heisenberg equation of motion and incorporating the above approximation we obtain the dynamical equations given by the following five equations, In obtaining Eqs. (2)-(4) we consider a setup where only third cavity has qubit-cavity coupling, i.e., g a = g b = 0, g c = 0. This is because, we want {n a , n b , n c , s c , s z c } ≡ {N, 0, 0, 0, −0.5} to be a SP at t = 0 which is experimentally more feasible. The more general case, g a,b,c = 0 will have more complicated SP solutions which att = 0 may not have easily experimentally implementable initial conditions. In addition, we consider ω a = ω c = ω b − ∆ without loss of generality and write Eqs. (2)-(4) in the rotating frame of frequency ω a,c . Here, ∆ is the detuning of cavity-b from cavity-a/cavity-c. We work with the Hermitian problem governed by Hamiltonian (Eq. 1), subjected to the conservation n a + n b + n c + s z c + 1/2 = N . We show results pertaining to non-Hermitian processes in the supplementary material 46 . It turns out that, the non-Hermitian results for a complete sweep reflect similar features as the Hermitian case, provided the dissipation rates are considerably less than 1/τ .

Stationary Point Solutions:
The SP solutions are parametrised byt and are therefore independent of τ . To obtain SP solutions at varioust, we define the quantityĥ ≡Ĥ(t) − µ(t)(n a +n b +n c +ŝ z c + 1/2) where µ(t) is a Lagrange multiplier (chemical potential) ensuring conservation. As in previous section, we derive Heisenberg equations of motion with respect toĥ and by setting the time-derivatives to zero, we write below four equations, The above four equations along with the constraint n a + n b + n c + s z c + 1/2 = N and the fact that s 2 c +s * 2 c 2 + s z 2 c = 1/4 (spin length) gives us six equations and six unknowns (a, b, c, s c , s z c , µ). This has been solved numerically to yield SP solutions at respectivet values (Fig. 2). In particular, we look for SSP branch that facilitates cavity-a to cavity-c transfer by sweeping Θ(t) from 0 to π/2.
Numerical Results: In this section, we show the SP solutions for g c = 0.2K and present the time dynamics for various sweep rates 1/τ (Fig. 2). We present SP solutions only for n SP a and n SP c since we are interested in photon number of source (cavity-a) and terminal cavity (cavityc). Among the various SP branches, there is a SSP branch that starts as {n SP a = N = 20, n SP c = 0} att = 0 and ends as {n SP a ≈ 0, n SP c ≈ N } with negligible value of n SP b . This special branch is the nonlinear analog of 'dark state'. This SSP branch helps us to choose the correct initial conditions {a(0), b(0), c(0), s c (0), s z c (0)} that need to be subjected to Eqns. 2-4 for a chosen value of pulse width (τ ) in J 1 (t) and J 2 (t). In addition, we need to make sure that the time dynamics remains on the SSP branch. This can be achieved by choosing optimal (elaborated below) sweep rate 1/τ subsequently leading to efficient photon transfer. Fig. 2 captures two main aspects. One is the SP so-lutions which are dependent only ont. The other aspect is the real time dynamics (which depends on sweep rate 1/τ ) where one wishes it to adiabatically follow the SSP branch that leads to efficient photon transfer. It is worth recapping the adiabatic theorem which states that upon slowly varying the parameters of the Hamiltonian, the system remains in the Hamiltonian's instantaneous eigenstate. We adapt a similar intuition for its semi-classical limit which forms the basis for the standard STIRAP scheme 38 . This implies that sweep rates cannot be too fast irrespective of whether the Hamiltonian is interacting or non-interacting. On the other hand, the sweep rates cannot be too slow if the system is interacting (chaotic to be more precise) as demonstrated in  On the other hand, in (c) and (d), we present two cases both in which we initiate the system att > 0 on the SSP branch. One case is in which the system is initiated immediately after the location of the inset in Fig. 2c Fig. 2c shows the regime where SSP branch (marked by red arrows) shows a sharp change in SP solutions. and caption therein. Fig. 2 (a) and (b) demonstrate near perfect transfer when the couplings are swept relatively faster (red circles). A relatively slow sweep (blue dots) shows smooth following of SSP branch only till the onset of chaos. Our findings therefore demonstrate that the standard notion of adiabaticity is contradicted when the system is chaotic 25 . To ensure a smooth following of the SSP branch even in the chaotic window, the sweep rate should be sufficiently faster. To show that chaos is the only origin of this breakdown, we show in Fig. 2 (c) and (d) that if the system is initiated by avoiding the chaotic region (shown as blue dots, slower sweep), then it follows the SSP branch. Therefore, the standard notion of adiabaticity holds here. In the inset region of Fig. 2(c), one can notice a sharp change in the SP solutions. To follow the SSP branch in this region the sweep rate needs to be slower than the rate of change of of the corresponding energy of the SSP solution w.r.t.,t. This is the reason why we have small deviations from SSP branch in Fig. 2(a) for the faster sweep rate case (red circles). The same faster sweep gives a smooth following of SSP branch [ Fig. 2

(c) and (d)] if we initiate the system after the inset region.
We have also done analysis for higher g c values and our findings still holds (see supplementary material 46 ).
Lyapunov exponent analysis: In this section, we will  do a LE analysis that will quantify chaos and this will turn out to be an important ingredient in executing the STIRAP scheme explained earlier. One starts with infinitesimally separated initial conditions and extracts LE from diverging copies of trajectories. Due to the fact that n a + n b + n c + s z c + 1/2 = N , our phase space is bounded. Consequently, the distance between the trajectories does not grow monotonically and at long times this may produce false estimate of LE. To circumvent this issue we exploit the prescription of resetting the phase space distance between trajectories in Refs. 47-51. The method is described as follows: Two phase-space trajectories are chosen, so that they initially differ by phase-space distance δ 0 . The trajectories are allowed to diverge for a time step ξ and the new distance between the trajectories δ 1 is reset to δ 0 . This procedure is repeated M times where M is large. The LE is then defined as, λ M = lim δ0→0 where λ M means that we have computed the Lyapunov exponent upto the time t = M ξ. Here, δ j denotes the deviation between the trajectories before the j th reset. The maximum LE (λ max ) can be obtained by taking the limit M → ∞ and a positive (zero) λ max indicates a chaotic (non-chaotic) behavior. Resetting at every step ensures that the deviation of trajectories is well within the phase space boundary.
In Fig. 3, in the left coloumn, we plot λ M as a function of dimensionless time KM ξ. Each of the figure is for a representativet value. It is to be noted that specifying at, fixes the parameters of the Hamiltonian J 1 , J 2 therefore making the Hamiltonian explicitly time independent. The procedure we employ to generate Fig. 3a -Fig. 3c is the following: We create two infinitesimally seperated copies, A, B such that P (A) = P SP (t) and P (B) = P SP (t) + δP (t) where P denotes the set {a, b, c, s c , s z c }. It is worth re-emphasizing that H(t) is explicitly time-independent for a particulart. As can be seen in Fig. 3 (left coloumn), we see non-chaotic (Fig.  3a), chaotic (Fig. 3b), and again a non-chaotic (Fig. 3c) behaviour. This is also reflected in the spreading features of phase-space points (right coloumn of Fig. 3). Keeping in mind, that we are interested in a transfer of photons from cavity-a to cavity-c, as a section of our phase space, we choose n a − n c on the y-axis and the conjugate variable φ a − φ c on the x-axis for the right coloumn of Fig.  3. Fig. 3 (e) suggests that such chaotic stages should be crossed quickly (by choosing a sufficiently fast sweep rate) to minimise the effect of chaos. For higher g c , the analog of Fig. 3 shows higher Lyapunov exponents and more prominent phase-space spreading (see supplementary material 46 ). In Fig. 4, we show the maximum LE (λ max ) as a function oft for three values of g c . As can be seen, for each g c , there is a windowt left <t <t right where λ max > 0. Therefore, LE analysis plays a paramount role in obtaining a window oft where the system is chaotic. In particular, for the case of g c = 0.2K, the vertical lines in Fig. 2 are obtained by the above analysis.
Transfer efficiency: In Fig. 2 we showed that the dynamics and the transfer efficiency have strong dependence on the sweep rate. In Fig. 5 we plot the transfer efficiency T w.r.t 1/τ for various g c . Fig. 5 (a) demonstrates that the slow sweep boundary does not exist for noninteracting case, implying no presence of chaos. However, the sweep must not be too fast which will result in violation of adiabaticity. This feature is reflected in the low-transfer sudden region beyond the black dashed lines in Fig. 5. The black dashed lines (1/τ fast ) in Fig. 5 are constructed such that T becomes less than 0.95 beyond the line. It is to be noted that 95% is generally regarded as satisfactory high efficiency 38 . In Fig. 5(b,c,d), we show the interacting case (g c = 0). Compared to the noninteracting case [ Fig. 5 (a)] additional chaos-dominated region emerges for the interacting case (g c = 0). The dotted red vertical line (1/τ slow ) sets a lower bound on the sweep rate below which T < 0.95. In other words, the sweep rate for high efficiency needs to satisfy 1/τ fast > 1/τ > 1/τ slow . It is to be noted that during real time dynamics, the time spent within the chaotic window (t right −t left )τ is finite. This means that the relevant LE for chaotic spreading is the finite time LE (see Fig. 3). This automatically implies, 1/τ slow < λ peak max where λ peak max is the peak of λ max for a given g c (see Fig. 4). Therefore, this finding relates LE to the lower bound on sweep rate. For higher g c , 1/τ slow increases thereby shrinking the efficient region and widening the chaotic regime. Four-cavity STIRAP scheme: To demonstrate scalability, we extend the three-cavity network to a fourcavity network (cavity-a, cavity-b, cavity-c, cavity-d) where cavity-d houses a qubit. We couple the cavities by counter-intuitive tunneling sequence shown in the inset of Fig. 6 (a). Fig. 6 (a) shows near-unity transfer from cavity-a to cavity-d for a linear closed-system case. Fig. 6 (b) demonstrates near perfect efficiency for the nonlinear case (g d = 0) when the parameters are swept at a faster rate. For a slower sweep rate (not shown here), similar to the three cavity case [blue dots in Fig. 2 (a) and Fig. 2  (b)], we find that the efficiency is hindered by chaos.
Conclusions and Outlook: We have demonstrated a protocol for achieving high transfer efficiency in an interacting c-QED STIRAP network. Such protocols are far from obvious given the fact that we are dealing with a scalable interacting system. To the best of our knowledge, for the first time, we have found and exploited the deep connection between LE and nonlinear STIRAP schemes. While, our protocols are developed on a semiclassical platform, we show that the resulting optimal choice of parameters successfully achieve our target for the fully quantum case (see supplementary material 46 ). Our findings are immensely useful for quantum communication and state transfer in cavity-based quantum networks 56,57 and for nonlinear waveguide optics 38 .
Future outlook includes adapting these protocols in different fields where variety of engineered Hamiltonians are achieved (for e.g., Optomechanics 52,53 ). It is interesting to generalize our scheme to higher dimensional systems 3 and complex geometries 1 . An open fundamental question is connecting Out-of-time-Ordered Correlator (OTOC) 5,55 and LE in our STIRAP setup especially because STIRAP is a unique platform to access both chaotic and non-chaotic regimes.
Supplementary Material for "Emergence of chaos and controlled photon transfer in a cavity-QED network"

Amit Dey and Manas Kulkarni
International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru -560089, India

I. SEMICLASSICAL-QUANTUM CORRESPONDENCE
Here, we make a comparison between the semiclassical and exact quantum many-body treatments of our model subject to a STIRAP protocol. The linear case (non-interacting) in Fig. S1 (a) and the nonlinear (interacting) fastersweep case in Fig. S1 (b) show remarkably good quantum-classical agreement. In Fig. S1 (c), both semiclassical and quantum treatments predict decreased efficiency due to chaos for a slower sweep case. In this case [ Fig. S1 (c)], we note that the quantum expectation values n a,b,c do not exactly follow the semiclassical oscillations that begin in the chaotic window (in the slower sweep case).
In a chaotic region, for quantum and classical counterparts, both the mechanism and diagnostics are different and this is indeed a subject of active interest. This stems from the fact that for the quantum Hamiltonian there are a large number of many body eigenstates (that depends on dimensionality of Hilbert space). However, in a semi-classical approximated description, the number of SP branches are small and this number is solely dependent on the number of solutions of the semi-classical equations of motion. Quantum measures of chaos include participation number of eigentates 1,2 , level spacing statistics 3,4 , and Out-of-time-ordered correlator (OTOC) measures 5 . The possible analogy between quantum measures, LE analysis of the semi-classical equations in the context of STIRAP scheme is an interesting open question.

II. DISSIPATIVE EFFECTS ON PHOTON TRANSFER
In main text, we deal with Hermitian dynamics governed by Eq. 1 (of the main text) that neglects various decay channels. In this section, we include photon decay and spontaneous qubit decay quantified by rates κ and γ, respectively. Incorporating dissipation, the modified equations of motion can be written as, The dissipative dynamics leads to a steady state with no photon and the qubit at its ground state. If the dissipation rates are comparable to the sweep rate (i. e., κ, γ ∼ 1/τ ) the target goal (of efficient transfer) through STIRAP profile is not achieved. This is because the dissipative effects wash away the possibility of any non-trivial steady state. Therefore, for a practical scenario it is needed that the sweep be completed before the dissipative effects become considerable, i. e., γ, κ ≪ 1/τ .   (b). Interestingly, the rapid oscillation for slower sweep starts exactly at the place where chaos appears for the Hermitian case in Fig. 1 of main text. As long as the dissipation does not become considerable, the system can still be approximated as a nearly closed system and the Hermitian analysis is valid. Therefore, if 1/τ ≫ κ, γ, then for a complete STIRAP cycle, one can make the Hermitian approximation.

III. RESULTS FOR HIGHER gc VALUE
Here we present results for g a,b = 0, g c = 0.4K. Compared to the g c = 0.2K case in Fig. 2 (of main manuscript), chaos is more intense (for e.g., as quantified by λ peak max in Fig. 4 of the main text) and the resulting breakdown appears at comparatively faster sweep rate (than the sweep rate that produces chaotic breakdown in Fig. 2 of main text) in Fig. S3. This observation is also consistent with Figs. 4 and 5 of main text. In Fig. S4 we plot LE for same light-matter couplings as in Fig. S3.