Non-equilibrium dynamics of non-linear Jaynes-Cummings model in cavity arrays

We analyze in detail an open cavity array using mean-field description, where each cavity field is coupled to a number of three-level atoms. Such system is highly tunable and can be described by a Jaynes-Cummings like Hamiltonian with additional non-linear terms. In the single cavity case we provide simple analytic solutions and show, that the system features a bistable region. The extra non-linear term gives rise to a rich dynamical behaviour including occurrence of limit cycles through Hopf bifurcations. In the limit of large non-linearity, the system exhibits an Ising like phase transition as the coupling between light and matter is varied. We then discuss how these results extend to the two-dimensional case.


Introduction
The use of cavity QED tools is now ubiquitous across different areas of physics ranging from quantum information [1, p 75] to detection of dark matter [2]. Specifically, the atoms held in optical cavities play a vital role in studies of many-body physics [3]. Such systems are natural implementations of many-body Jaynes-Cummings (JC) or Dicke Hamiltonians [4]. Their high tunability and the possibility of achieving strong light matter coupling or probing the dynamics in real time make them very attractive experimental platforms. The prospect of probing phase transitions and the associated critical phenomena with these platforms have been put forward e.g. in [5][6][7] and the non-equilibrium dynamics of the Dicke model has been theoretically investigated in [8]. Specific many-body phenomena that can be studied with cavity QED include for example the physics of spin glasses [9,10], the emergence of gauge fields and the related quantum Hall effect [11] or the selforganization of the atomic motion [12][13][14], to name a few. The self-organization has been subsequently observed in the experiments [15,16].
So far we have mentioned only studies concerning a single cavity-generalizations to multiple cavity arrays implementing the Hubbard physics have been reviewed in [17]. Although appealing in principle, the realization of many efficiently coupled cavities, each hosting a discrete-level quantum system is a challenging task. To make such experiment scalable requires miniaturization of the cavities. One possibility is the use of microcavities in photonic crystals [18,19]. A further option is to use integrated optical circuits, where in principle arbitrary waveguide forms can be created with high precision by laser engraving in the silica substrate [20]. They have been successfully used for the demonstration of a quantum gate operation [21], creation of classical and quantum correlations [22,23], multi-photon entangled state preparation [24], quantum random walk [25], discrete Fourier transform [26] or Bloch oscillations [27].
While it has been demonstrated that cavities can be fabricated by creating the Bragg grating during the laser writing process [28,29], there is now an active experimental effort to combine the waveguides with atomic microtraps on a single device [30,31].
Motivated by these developments and the prospects of studying many-body physics using integrated optical circuits with trapped atoms, we theoretically analyze the non-equilibrium physics of such system, which we take to be a two-dimensional cavity array, where each cavity hosts a number of atoms. We first derive the effective Hamiltonian describing the system in section 2. Using this Hamiltonian we find mean-field (MF) equations of  and Here, w aux is an auxiliary frequency used in the adiabatic elimination which, in principle, can be chosen arbitrarily. Its physical motivation and interpretation will be discussed momentarily in section 3.1.

MF treatment of the nonlinear JC model
The effective Hamiltonian (2.2) is the starting point in this section, where we analyze the effect of cavity loss and pump on the dynamics. Also, it can be seen from the form of the Hamiltonian and the expressions (2.3), that the parameters are highly tunable through varying W D , e and D s . We first perform a MF analysis of a simpler system with only one cavity mode, which already contains rich physics as we show below. We then turn back to the multi-mode two-dimensional setup in section 3.2.

Single cavity Cavity without pump
Lets first consider the Hamiltonian (2.2) in the single cavity limit with the single mode a, i.e. we drop the indices i and ν (see figure 1(c)). We rewrite the Hamiltonian as˜( with ω being the cavity frequency. Note that the model (3.1) without the λ term is the usual JC model [32]. Below we show, that the λ term is indeed at the origin of intriguing system dynamics (see also [33][34][35][36][37][38][39][40] for various other nonlinear extensions of the JC model).
One can now derive the equation of motion for the operator o according to˙[ where H is given by (3.1). At the same time, any realistic cavity is subject to a decay of the electromagnetic field into the environment. The dissipation process is typically described by means of a master equation (see e.g. [41]), which corresponds to an extra term in the equation of motion for the cavity mode operator,˙k µa a, where κ denotes the cavity loss rate (see figure 1(a)). Introducing the expectation values of the operators a = á ñ = áS ñ = áS ñ a s w , , z , it is now straightforward to derive the the MF equations of motion, which reaḋ In the derivation we have used the MF decoupling † * a á S ñ = a wand a áS ñ = a w z . We have also neglected the spin decay on the transition | | ñ -ñ g s 2 First, we wish to find a steady state solution of the equations of motion (3.3). From (3.3a) and (3.3c) we havẽ Substituting the equation (3.4) to (3.5) yields When C is complex, which occurs only for non-zero cavity decay k ¹ 0, the condition equation (3.6) can be satisfied if s=0 (and consequently a = 0), which results in a trivial solution with empty cavity and all spins 2 This is a justified assumption in implementations with real atoms as the levels | | ñ ñ g s , typically belong to some ground state manifold, where only magnetic dipole transitions are allowed between the states of that manifold. In turn, the spin decay is negligible compared to the cavity decay. down (up) in the steady state. This indicates that in order to obtain some non-trivial physics in the steady state limit and to counteract the cavity losses, one needs to input energy into the system 3 . In our case a natural choice is either through the pumping of the cavity mode or driving directly the | | ñ -ñ g s transition. We focus here on the former case only.

Cavity with pump
The cavity pump can be described by the Hamiltonian where N is the total number of the spins. This can be easily verified as where we have parametrized the total spin as x y . The set of equations (3.8) will serve as the starting point for most of the analysis in this section. Ultimately, we are interested in the effect of the nonlinear term proportional to λ, which makes the system described by (3.1) qualitatively different from the usual JC model. Due to the complexity of the problem, we start by studiyng a simpler situation in section 3.1.1 by putting l = 0 (and as we will see also D = 0 at ). Building upon the l = 0 solution, we then analyze the l ¹ 0 and D ¹ 0 at situations in sections 3.1.2 and 3.1.3 respectively.

λ=0 regime
In order to investigate the steady state solutions, we first put l = 0 to further simplify the problem (see e.g. [8,42,43] for related studies of the phase diagram of the JC and Dicke models). We then turn back to the situation with l ¹ 0 in the next section. With these simplifications, the real equations for the steady state read Solving the set (3.11a) and (3.11b) for a a , R I and substituting to (3.11c)-(3.11e) we obtain for the spin conservation (3.9)˜( This is a 4th order polynomial for w and its solutions in terms of radicals can in principle be found yielding rather complicated expressions, which are not of much practical use. Instead, we will analyze the properties of (3.12) as follows. Since . At the same time the nominator of the rhs is clearly non-positive, so that the non- 3 Note that this is in sharp contrast with the full Dicke Hamiltonian, where the presence of the counterrotating terms guarantees non-trivial solutions even in the absence of the pumping [4,8]. In the model we study, the absence of the counterrotating terms is a direct consequence of applying the rotating wave approximation, as the cavity modes and the | | ñ -ñ g e transition are taken to be at optical frequencies.
positivity of the rhs requires the quadratic polynomial in the denominator to be non-negative. We thus have The roots of this polynomial read˜˜( It can be immediately seen that for k ¹ 0, any D ¹ 0 at would yield imaginary w. This however does not mean that a solution of the original constraint (3.12) is not possible for both k D ¹ , 0 at . Rather, it means that the point D = 0 at has some particular properties which we will investigate further. For D = 0 at , the steady state equations read˜( 3 These solutions are valid only for˜*  h º g N g 2 1 . In the case where w=0, we can use the spin conservation (3.9) to parametrize the spin as 3.17 I R Next we can express a a , R I from (3.15a) and (3.15b) as and corresponds to the value ofg , where q q = -+ cos cos , which can be seen as the rightmost point in figure 2. Note, that the region˜[ ] * * Î g g g , 1 2 admits both solutions (3.16) and ( (3.18) and (3.21)) indicating a bistable behavior. In what follows we refer to the points * * g g , 1 2 as transition points.

Stability study
The stability analysis of the steady state solutions , , , where ( ) p y 4 is some polynomial which is 4th order in y. We thus always have one eigenvalue y=0, which is simply the consequence of the spin conservation law (3.9). Examining numerically the negativity of real part of the roots of ( ) p y 4 we identify stable and unstable solutions. These are depicted by solid and dashed lines respectively in the steady state phase diagram figure 3, where we plot | | a 2 and w as functions of the couplingg . The region II exhibits bistable behavior, whereas regions I and III admit only a single stable solution. We have also verified, by numerically solving the dynamical equations (3.8), that they indeed evolve into the steady state solutions (3.16), (3.18) and (3.21) (not shown).
A comment on the validity of the MF approach is in order. At first sight, the presence of bistable regions in the MF solutions seems to be incompatible with the expected existence of a unique stationary state as it suggests that there are two different states for a given value of the couplingg . However, such situation commonly occurs in the MF treatment of many-body systems that feature a first order phase transition. It can be encountered, for example, in the context of the classical van der Waals gas, described by a MF equation of state [44] or, more recently, in the context of Rydberg gases, where optical bistability was observed in qualitative agreement with the MF predictions [45].
In the next two sections we investigate how the inclusion of the non-zero atomic detuning D at and non-zero coupling λ terms modifies the l D = = 0 at solution.

l ¹ 0 regime
In order to shed light on the effect of the λ term independently of the D at term, we keep D = 0 at . From (3.8) the steady state equations read . We provide the details of this expansion in appendix C. In order to simplify the analytic expressions we take the imaginary part of the pump to be zero, h = 0 I . In this case, to leading order in λ, the solutions read ( ) 3 . 2 7 . This is in contrast to the l = 0 case, where two distinct transition points * * g g , 1 2 were identified. Here, the solution . This indicates that the transition point * g 2 approaches * g 1 as λ is increased until they become identical in the l  ¥ limit. In other words, the region II of the phase diagram figure 3 is shrinking to zero width as λ is increased. For the intermediate values of λ, the values of * * g g , 1 2 can be found numerically (selecting only the physically meaningful solutions of the 6th order polynomial, corresponding to the real values of w). We plot the shrinking of the [ ] * * g g , 1 2 region in figure 4(a). In figure 5 we show the phase diagrams for | | a 2 and w for increasing values of λ. It can be seen from figures 5(a) and (b), that for some values of λ, there is a region with no apparent stable solution. In nonlinear systems, this is typically a signature of the appearance of limit cycles which occur through a Hopf bifurcation as the system leaves the stable fixed point by changing the couplingg [46]. We represent a limit cycle corresponding to the point A in figure 5(b) as the time evolution of the global spin in figure 6.
One can carry the analysis further in the large λ limit and look for asymptotic solutions in the vicinity of the unique transition point * g 1 . The first non-trivial contribution to α is of order l -1 and given by (see appendix C) which is valid for˜*  g g 1 as the solution for˜* > g g 1 is trivial (a = 0). The scaling of the spin observables, say w, is simply obtained from the expansion of (3.27), which read 3.29 for , 3.29 where the sign±depends on the branch of w considered. Note that the scaling (˜) * g g 1 1 2 is characteristic of the Ising type model with  2 symmetry and * g 1 is the phase transition critical point. We compare the asymptotic solutions (3.28) and (3.29) with the solutions obtained by numerically solving (3.25) in figures 5(e) and (f).

λ=0 regime with D ¹ 0 at
Our next aim is to explore the effect of the D at term on the solution in the absence of the nonlinear λ term. In this case, the structure of the solutions is dictated by the 4th order w polynomial (3.12) and it is straightforward to obtain the solutions numerically. For small D at , there is a range ofg -values, which admits four solutions (corresponding to four real solutions of (3.12)), figures 7(a) and (b). As D at is increased, this region eventually disappears leaving us with only two solutions for all values ofg . The latter limit can be simply understood from the MF equations (3.11). In the large D at limit, the leading contribution comes from (3.11c) and (3.11d) with the trivial solutions = = s s 0 x y implying =  w N 1. The evolution of the solutions towards this large D at limit can be seen rather clearly from figure 7(d). Next, we have determined numerically the size of the region ofg admitting four solutions as a function of D at . This is shown in figure 4(b). Clearly, there is some maximal D at after which there are only two possible solutions, as discussed. This limiting value is represented by the red circle in figure 4(b). Note, that the inclusion of the D at term also lifts the transition point * g 1 , i.e. that all solutions are smooth in the vicinity of * g 1 .

Multiple cavities
In this section we seek the generalization of the single cavity case to higher dimensional geometries. For concreteness, we consider a 2D square geometry depicted schematically in figure 1(a). Before diving in the details of the analysis, we motivate this section by asking whether a 2D square geometry offers MF solutions which are qualitatively different from the 1D case studied above. For example, it was shown in [47] in the context of laser driven and interacting Rydberg gases on a square lattice, that a homogeneous system admits a MF solution that breaks the lattice symmetry and exhibits antiferromagnetic (AF) order.
We start our discussion by deriving the 2D MF equations, which follow from the operator equations of motion (A.10) given by the Hamiltonian (2.2). Repeating the argument yielding (3.6) leads to the same conclusion, namely that in the absence of the cavity pump the system posses only a trivial solution, where the cavity modes are empty and all the spins are down (up). We first focus on the situation without the nonlinear term, l = 0. In the following, we consider the cavity pump η, the cavity decay κ and the spin decay γ to be the same for all cavity modes and all spins respectively (the motivation for adding the spin decay will become clear shortly). We obtain the set of MF equationṡ ˙( Here, we allowed for the couplingsg and the photon detunings D ph to be different for rows and columns while taking the remaining parameters to be the same for all cavities. Clearly, the 2D system provides higher tunability by enlarging the parameter space. In the following, we study the problem considering two different perspectives, namely a homogeneous and cluster-MF ansatz.

Square array with homogeneous MF
In order to further simplify the MF equations (3.30), it is reasonable to use the typical MF ansatz, namely that, due to translational symmetry in either direction of the 2D array (along rows or columns), the corresponding cavity fields are the same (ā a = Here N R (N C ) is the number of rows (columns) respectively. Since the fields β and α share the same spin s, β can be given directly in terms of the field α. The situation is then essentially equivalent to the single cavity case analyzed in section 3.1, though with some extra tunability provided by larger number of parameters. For example, one can find the transition point * g 1 by combining (3.31a) and (3.31b) with˜ã b , which is one possible solution of (3.31c). The transition point then corresponds to the situation where w=0, i.e. One can then identify the critical value of e.g.g a for all other parameters fixed. As a consistency check, it is easy to verify that one recovers the single cavity expression by omitting all the 'b' variables and setting N R =1 and = N N C , the number of spins.

Square array with cluster-MF
In the previous section, we have used the homogeneous MF ansatz and found that the solutions correspond effectively to a single cavity case with no intriguing spin configurations. However, it is known and was shown e.g. in [47], that a suitable MF ansatz can lead to a non-trivial configuration (such as antiferromagnetic ordering) even if the steady state MF equations are completely symmetric under exchange of the spin variables. Here we will consider the simplest possible case of such ansatz where a cluster is formed by two adjacent inequivalent spins s w , 1 1 and s w , 2 2 . On the other hand, since the fields α and β couple to both s 1 and s 2 in the same way, we take a b , to be the same along rows (columns) due to translational symmetry. In order to simplify the equations, we now take all the parameters to be the same along all rows and columns (i.e. we set˜˜= = where¯= j 1, 2 labels the different spins of the cluster and we assumed that each field a b , couples to the same number of spins 1 and 2, hence the factor 1/2 in the second term of (3.33a) and (3.33b). Here N denotes the total number of spins, i.e. there are N spins along rows and along columns.
Until now we did not comment on the spin conservation in the 2D case. Going back to the most general situation, where every individual spin and cavity mode is described in terms of the corresponding MF variables¯¯n provided g = 0. In (3.35) = n W w i and S = n s i are the global spin components. One should appreciate that (3.35) actually prevents any spin configuration incompatible with it, including the AF order (which corresponds to =w w 1 2 when using the here considered cluster-MF). In order to see e.g. the AF order, such as in [47], one needs to break the global spin conservation. This is achieved by the inclusion of the spin decay γ, which breaks both the local and global spin conservations (3.34) and (3.35).
Allowing for the spin decay, we can now ask, whether the solution of the cluster-MF equation (3.33) features any non-trivial spin configuration ( ¹ w w 1 2 ). First, it follows from (3.33a) and (3.33b) that a b = . Next, expressing fors j from (3.33c) and substituting to (3.33d), it can be shown that Once again we find that the structure of the equations (3.33) reduces the problem effectively to the single cavity situation described by a simple set of variables a s w , , and this even when allowing for inequivalent spin configurations and the spin decay.
One might argue, that the equivalence to the single cavity case conjectured in sections 3.2.1 and 3.2.2 is an artefact of taking l D = = 0 at . Indeed, when comparing (3.3) and (3.30), one can note that there is a qualitative difference between the 1D and 2D situation. Specifically, there is an extra coupling between the a and b modes, the last term in (3.30a) and (3.30b). We have numerically verified that, in a general situation with variablesā b s w , , , In summary, the MF equations (3.30) with the ansatz considered in sections 3.2.1 and 3.2.2 on a square lattice reduce to effectively one-dimensional description with no intriguing spin configurations. It would be desirable to perform a beyond MF study of the nonlinear two-dimensional model in order to asses the true nature of the steady state and the corresponding spin and field configurations, including their mutual correlations. Also we did not fully exploit the possibilities offered by the proposed implementation, such as taking different geometries of the array or allowing for disordered coupling strengths, which we leave for further investigations.

Conclusion
Motivated by the progress in integrated optical circuits, we have proposed a possible realization of a twodimensional cavity array with trapped atoms, which is a promising scalable quantum architecture. We derived an effective description of the system in terms of JC like Hamiltonian with highly tunable parameters and extra nonlinear terms. We then analyzed the dynamics of the system using a MF approach. We have found a rich behavior including bistable regions, Ising like phase transition or occurrence of limit cycles through Hopf bifurcations. In the present setup, we have not found conceptual differences between the one and twodimensional cases at the level of the MF description and with the geometry considered. We hope that the present work lays down grounds for future studies of the cavity arrays realized with integrated optical circuits. The problems which might be addressed in the future are e.g. going beyond MF description, accounting for more exotic geometries or studying effective spin physics as a low energy limit of the presented cavity Hamiltonian.
where the summation over ν and i in the last two equations is emphasized.  .
The set of equations of order l 1 read C . 3 x y 0 2 0 2 0 2 2 In order to proceed, we realize that the equation  i.e. is unphysical. We are thus left with the solution (˜) . It is easy to find, that * * = g g 1 . This completes the leading order solutions and yields the expressions (3.27a) and (3.27b).
The equations (C.2), (C.3) and (C.6) yield a closed set for spin variables up to the order l -1 and for α up to l -2 . It is straightforward to find the solutions up to the respective order explicitly, giving however rather lengthy algebraic expressions. These can be simplified in specific situations. For example, as we discuss in the main text, when looking for solutions in the vicinity of the transition point * h = g N 2 1 , the leading order contribution to | | a 2 is of order l -2 , namely which gives the relation (3.28). Similarly, the leading contribution to w is of order l 0 yielding the expressions (3.29).