Synchronization of micromasers

We investigate synchronization effects in quantum self-sustained oscillators theoretically using the micromaser as a model system. We use the probability distribution for the relative phase as a tool for quantifying the emergence of preferred phases when two micromasers are coupled together. Using perturbation theory, we show that the behavior of the phase distribution is strongly dependent on exactly how the oscillators are coupled. In the quantum regime where photon occupation numbers are low we find that although synchronization effects are rather weak, they are nevertheless significantly stronger than expected from a semiclassical description of the phase dynamics. We also compare the behavior of the phase distribution with the mutual information of the two oscillators and show that they can behave in rather different ways.


I. INTRODUCTION
Self-sustained oscillators do not have a preferred phase, but when two or more of them are weakly coupled together a phase preference can emerge spontaneously, an effect known as synchronization. Self-sustained oscillators are ubiquitous in nature and synchronization effects have been widely studied across the physical and biological sciences [1]. Synchronization has also been studied in quantum optical systems such as the laser, although generally focussing on regimes where approximate semiclassical descriptions work well [2,3]. In the last few years there has been considerable interest in studying the synchronization of oscillators and related systems [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] close to threshold or at low excitation levels where semiclassical approaches break down and fully quantum mechanical calculations are required. Recent theoretical work has explored different ways of quantifying synchronization in quantum oscillators [5,7,14,15,20], as well as investigating the connection between it and measures of correlation such as mutual information and entanglement [5,8,10,15,19]. Detailed comparisons have also been made between the predictions of quantum models and those of related semiclassical or classical descriptions [7,12].
Studies of synchronization effects in the quantum regime have largely concentrated on the behavior of simple model systems such as van der Pol oscillators [7,9,10,12,15,17,19,21] (together with closely related models [22]), though a number of other systems including atomic ensembles [13,16] and optomechanical oscillators [5,6,12,18] have also been investigated. In this article we investigate synchronization in a very different model system consisting of two weakly coupled micromasers.
The micromaser is a self-sustained oscillator consisting of a microwave cavity driven by a steady flow of excited atoms which interact strongly with a particular cavity mode [23][24][25]. The micromaser was used to carry out a range of pioneering experiments in quantum optics [23]. However, it has also become possible to engineer sys-tems with similar behavior in the solid-state using, for example, superconducting [30,31] or optomechanical [32] devices.
The micromaser makes a very interesting model system with which to explore synchronization effects in the quantum regime because it displays a very rich range of dynamical behaviors, including strongly non-classical features which go well beyond those found in simpler systems like the quantum van der Pol oscillator. Furthermore, an exact steady-state solution is available for the density operator of the micromaser [24] and important dynamical properties such as the linewidth [26][27][28][29] have been studied extensively.
Using a probability distribution for the relative phase, we explore how a preference for a particular phase (or phases) emerges when two micromasers are weakly coupled together in different ways. We investigate the behavior of the phase distribution over a wide range of parameters ranging from a semiclassical regime where photon occupation numbers are large to a quantum regime where occupation numbers are small and the steady-state of the system can be strongly non-classical. We derive a simple Fokker-Planck equation for the phase distribution assuming large photon occupation numbers and compare it with numerical calculations using the full master equation of the system. We find that whilst the Fokker-Planck equation provides a good description of the phase distribution in the semiclassical regime, it substantially underestimates the extent to which a preferred phase emerges in the quantum regime. We also compare the behavior of the phase distribution with the mutual information and entanglement of the micromasers and find that the behavior is somewhat different in each case.
This work is organized as follows. We introduce our model of the coupled micromaser system and briefly review the key properties of the uncoupled micromaser in Sec. II. Then we introduce the relative phase distribution in Sec. III. We show how perturbation theory can be used to understand the behavior of the phase distribution in the weak coupling limit in Sec. IV. Then in Sec. V we derive a simple analytic formula for the phase distribution in the semiclassical limit and compare it with numerical calculations. We examine the behavior of other measures of correlation between the micromasers in Sec. VI. Finally, we summarize our findings and discuss possible directions for future work in Sec. VII.

II. MICROMASER MODEL
We consider a model system consisting of two micromasers, assumed for simplicity to be identical, that are coupled together weakly. As is the case for classical oscillators [1,33,34], the behavior can be very sensitive to the form of the coupling as well as its strength. We investigate two specific forms for the coupling involving either additional terms in the Hamiltonian of the system (coherent coupling) or additional dissipative terms in the master equation (dissipative coupling). The starting point for our model is the standard master equation description for the micromaser [24,25], to which we shall simply add additional terms to describe the coupling.
The master equation for the density operator of the two micromasers, ρ, (in the interaction picture) takes the formρ where here L 1 [ρ] and L 2 [ρ] describe the uncoupled dynamics of the two micromasers and the interaction between them is given by L c [ρ]. The dynamics of each individual micromaser is controlled by a balance of interactions between the cavity and the flow of atoms which pass through it, and between the cavity and its electromagnetic environment which gives rise to losses. The atoms can be traced out of the master equation so that the dynamics of the system is captured by the terms [24,25] where a j is the lowering operator for a mode of the cavity j (j = 1, 2), N is the rate at which atoms pass through the cavity and φ is the Rabi angle which quantifies the strength of the atom-cavity interaction. Note that we have adopted units of time such that the cavity loss rate, γ, is unity and we have taken the zero-temperature limit for simplicity. For coherent coupling between the two micromasers the interaction is described by the Hamiltonian with ε the coupling strength (scaled by γ), and hence in this case For dissipative coupling the master equation includes terms which describe an additional loss channel for the system whose properties depend on the state of both modes [10,12,17] In the following we will focus mainly on the regime where the coupling is very weak ε ≪ 1. Before going on to look at how either coherent or dissipative coupling affects the system, we briefly review the most important properties of the uncoupled micromaser (ε = 0). The steady-state density operator for a micromaser is diagonal in the number state representation and the probability of finding the cavity in the n-th number state is [24,25] where K is a constant determined by normalization. Although the density operator can be written in terms of an apparently simple formula, the state of the system changes dramatically as a function of the atomic flow rate, N , and the strength of the atom-cavity interaction, parameterized by the pump parameter: Θ = φN 1/2 . Fig.  1 illustrates the behavior of the average occupation number n and the number fluctuations, measured by the Fano factor F = ( n 2 − n 2 )/ n , as a function of Θ and N ; the inset shows an example of the full number distribution P n . The system has a threshold at Θ = 1, above which a limit-cycle emerges-manifesting in the P n distribution as a peak at non-zero n. After initially growing very rapidly in size, the limit-cycle then gets progressively smaller as Θ is increased (up until Θ ∼ 4). There is a strong peak in F around threshold and it then drops below unity, a signature of number squeezing [35].
For Θ > 4 the behavior becomes more complicated and the micromaser moves between a range of different states. It undergoes dynamical transitions between limit-cycles with different average energies [24] and can co-exist in a mixed state involving two limit-cycle states (seen as two peaks in the P n distribution). At certain specific values of φ such that sin(φ √ m + 1) = 0 with m = 0, 1, 2, ... the system becomes trapped : P (n>m) = 0 because the matrix element that generates transitions between the m and m + 1 number states vanishes. These trapping states [23] have a number state distribution that can be extremely sharply peaked. For example, at φ = π/ √ 2 (m = 1 trapping state) one typically finds P 1 ≫ P 0 and as N is increased the system gets closer and closer to being exactly in the n = 1 number state.
The extremely narrow P n distributions that the micromaser displays above threshold come close to reaching what might be thought of as the most quantum of limitcycle states-pure number states. However, the micromaser's steady-state isn't always strongly non-classical. The average occupation number of the micromaser is roughly proportional to N for fixed Θ and the nonclassical features are strongest for either relatively small N or larger values of Θ. In contrast, when n ≫ 1 and n ≫ φ 2 the dynamics of the number distribution, P n , is well described by a Fokker-Planck equation in which n is treated as a continuous variable [24], which corresponds to the semiclassical limit of the system.

III. RELATIVE PHASE DISTRIBUTION
Phase distributions provide a convenient way of characterizing the emergence of a preferred relative phase in systems of coupled oscillators. For a single oscillator the quantum mechanical phase distribution is given by [35], n|ρ|m e i(m−n)ϕ , (7) where |ϕ = ∞ n=0 e inϕ |n is an eigenstate of the Susskind-Glogower operator ∞ n=0 |n n+ 1|. This phase distribution, P (ϕ), also emerges naturally from the Pegg-Barnett description of the phase operator [36] or indeed when one seeks to define a distribution which is the canonical conjugate of the number distribution [37]. In the steady-state only the diagonal components of the uncoupled micromaser density operator are non-zero in the number distribution, so we see immediately that there is no preferred phase and the phase distribution is simply uniform: P (ϕ) = 1/2π.
When two micromasers are coupled, either coherently or dissipatively, a preference emerges for certain values of the relative phase, ϕ − = ϕ 1 − ϕ 2 , but not the total phase ϕ + = ϕ 1 + ϕ 2 . The relative phase distribution takes the form [14,38,39] which can also be rewritten explicitly as a Fourier series where we adopt the notation ρ (p) n,m = n, m+p|ρ|n+p, m . It is also helpful to be able to characterize the emergence of a phase preference using a single number. Starting from the relative phase distribution, one can simply extract the size of the peak relative to the uniform distribution [14], This is something that we will make extensive use of here, though it should be noted that the choice is by no means unique [18,22]. The relative phase distribution for the micromaser system is obtained by solving for the steady-state of the master equation (1) using standard numerical methods [40]. Coherent and dissipative couplings with the same strength give rise to markedly different behavior in the relative phase distribution, as is illustrated in Fig.  2. Coherent coupling leads to much weaker phase locking than dissipative coupling and generates a relative phase distribution which is π-periodic rather than 2π-periodic. This matches the well-known differences in relative phase dynamics of reactively and dissipatively coupled classical oscillators which are usually understood by deriving approximate equations of motion for the relative phase of the oscillator assuming weak coupling [1,34].

IV. PERTURBATION THEORY
Perturbation theory provides a straightforward way of understanding the differences between the quantum mechanical relative phase distributions generated by coherent and dissipative coupling. We begin by considering the case of coherent coupling before moving on to dissipative coupling.

A. Coherent Coupling
Writing the master equation in terms of the number state basis, we finḋ where the terms arising from the coupling are given by and the coefficients which describe the uncoupled evolution are given by [26] In the steady-stateρ n,m = 0, except in the diagonal case (p = 0) where the components must also obey the normalization condition.
Working to first order in the coupling means that we replace the terms in Eq. (12) by their unperturbed values which are all zero-apart from the diagonal ones, ρ n,m = P n P m . Therefore only the equations for the p = 1 components are affected (we only need consider the components with p > 0 which appear in the expression for the relative phase distribution (9)) for which we find Thus at first order, the components ρ (1) n,m are in general non-zero, and pure imaginary (since P n(m) are probabilities), whilst those with p > 1 remain zero. However, since ∆ n,m is no longer zero (since ρ (1) n,m is of order ε) and one finds that the components ρ (2) n,m are all real and proportional to ε 2 . In this case there is no cancelling of the terms in the sum and hence to second order the relative phase distribution takes the form where C 0 is a constant which depends on the parameters of the uncoupled micromasers (φ and N ). This of course is just what we see for the case of coherent coupling in Fig. 2.

B. Dissipative coupling
We now look at what happens for dissipative coupling described by Eq. (5). In this case we can simplify the master equation by taking into account the fact that some of the terms in (5) simply act to increase the effective damping of the oscillators and can be absorbed into the terms which are already present for the uncoupled system by rescaling the parameters:Ñ = N/(1+ε),ε = ε/(1+ε).
Working in the number basis, the master equation with dissipative coupling takes the form given by (11) (though withε andÑ replacing ε and N ) with ∆ (p) n,m = −ε (n + 1)(m + 1)ρ (p−1) n+1,m+1 At first order in the coupling, only the p = 1 term is non-zero, This generates non-zero components ρ m,n so there is no cancellation when they are summed up. Hence to lowest order in the coupling the relative phase distribution for dissipative coupling takes the form where C 1 is a constant (that again depends on the parameters of the uncoupled micromasers) which matches the behavior in Fig. 2. The insights into the general form of the relative phase distribution provided by perturbation calculations are quite general, in the sense that the overall form of the relative phase distribution is entirely determined by the coupling terms (only the constants C 0 and C 1 depend on the details of the micromaser systems). However, perturbation theory is also very useful for exploring the behavior where the state space of the system becomes so large that a direct numerical solution becomes impracticable.

V. SEMICLASSICAL LIMIT
To understand the behavior in the semiclassical limit we start from the expression for the relative phase distribution in the form of Eq. (9) and focus on the case of dissipative coupling, which is the simplest. The equation of motion for the relative phase is given bẏ n,m .
Using the master equation (11), with the dissipative coupling terms (18), we find that [28] ∞ p=1 e ipϕ− ∞ n,m=0ρ and So far the analysis has been exact. We now introduce some approximations to simplify things [26,28]. We will make use of the fact that the average occupation is large in the semiclassical regime so that n ≫ 1 and n ≫ φ 2 . We will also make two additional assumptions: that the coupling is weakε ≪ 1 and that the number distribution is strongly peaked about n , this will only be the case when the system is above threshold with Θ > 1, but not so large that more than one peak appears in the distribution [see the inset in Fig. 1]. These last two assumptions have nothing to do with semiclassicality (indeed a strongly peaked number distribution can be highly non-classical), but rather they are closely related to those usually made in classical analyses of synchronization [1,41] which rely on the existence of a well defined (single) limit-cycle state.
Our assumptions mean that we need only consider the components ρ (p) n,m for which n, m ≫ p. We proceed by expanding the coefficients in (22) treating p/n, p/m, 1/n and 1/m as small quantities and keeping the lowest order (non-zero) contributions in each case so that we have and In the last line we also expanded the sine assuming (pφ) 2 /n ≪ 1.
Finally, we make use of the assumption that the distributions are strongly peaked about a common average [26,28] and simply replace n and m (together with n + 1 and m + 1) with n so that now Eq. (22) takes the simplified form where we have defined The parameter∆ matches a simple approximate expression for the micromaser linewidth [26][27][28][29] which is valid in the semiclassical regime. Combining Eq. (28) with its complex conjugate leads to a Fokker-Planck equation for the relative phase distribution: The corresponding steady-state distribution is [41] P (ϕ − ) = 1 2πI 0 (ε/∆) eε cos ϕ−/∆ .
The Fokker-Planck equation for the phase distribution of the coupled micromasers is exactly what we would expect for coupled classical oscillators in the presence of noise [1,41]: it describes a competition between phase diffusion and the effects of the coupling which tends to drive the system towards a particular relative phase. However, the origin of the noise which drives the diffusion is nevertheless ultimately quantum mechanical rather than classical and hence it makes sense to see (31) as a semiclassical equation in this context. It is worth noting that a Fokker-Planck equation with the same form emerges in the analysis of coupled lasers far above threshold [3]. Now that we have obtained an expression for the phase distribution in the semiclassical limit we can look in detail at when and how its predictions differ from the full (quantum) dynamics predicted by the master equation. The phase distribution obtained using the semiclassical approximation (31) is compared with the results from a full numerical solution of the master equation for a relatively small pumping rate, N = 5, in Fig. 3. The first thing to note is that the standard expectation of classical synchronization theory is fulfilled: for the weak couplings used here a rather strong change in the phase distribution is combined with a relatively weak change in the average occupation number of the system. Furthermore, the semiclassical phase distribution does a reasonable job of describing the strength of the peak in the phase distribution for Θ values that are not far above threshold. However, the semiclassical calculation systematically underestimates the strength of the phase locking in the quantum regime where n ∼ 1 (i.e. Θ > 3) as shown in the inset to Fig. 3a. Although the phase locking is pretty weak for larger values of Θ, it is nevertheless about twice as strong as predicted by the simple semiclassical calculation. Although a full numerical solution of the master equation becomes very difficult for larger N values we can use perturbation theory to calculate the relative phase distribution provided that we choose a small enough value for the coupling (since it turns out that the range of coupling strengths over which the second order perturbation calculation is a good description varies with N ).  compares the value of S obtained using the perturbation and semiclassical calculations for a range of N values. It shows clearly that the semiclassical solution (31) provides an increasingly accurate description of the strength of the phase locking as N is increased, just as one would expect.

VI. MUTUAL INFORMATION AND ENTANGLEMENT
We now turn to look at how other measures of correlations, mutual information and entanglement, behave as a function of the type and strength of the coupling between micromasers. Recent work on coupled quantum van der Pol oscillators has suggested that there may be an intimate connection between the emergence of synchronization in coupled quantum oscillators and the behavior of the mutual information [8,15,19] and it has even been suggested that mutual information could serve as an order parameter for synchronization [15].
The mutual information, I, of the coupled micromasers is defined as where ρ is the full density operator, ρ 1 (2) , is the reduced density operator of micromaser 1(2) and S is the von Neumann entropy, S(ρ) = −Tr[ρ ln ρ].
Perturbation theory tells us immediately that the mutual information of micromasers will grow at least quadratically with the strength of the coupling, for both dissipative and coherent coupling. To see this, we can write the density operator ρ, its eigenvalues λ j , and eigenkets |j , as expansions in [42] To first order, the von Neumann entropy is with λ (1) j = j (0) |ρ (1) |j (0) , the first order correction to the eigenvalues. The density operators for uncoupled micromasers are diagonal in the number state basis, and as we have seen in Sec. IV the first order correction terms which form ρ (1) are all off-diagonal in the number state basis; consequently the first order corrections to the eigenvalues are all zero. Hence the first order contributions to I vanish for both coherent and dissipative couplings. Fig. 5 compares the behavior of the mutual information for dissipative and coherent couplings. Not only does I grow quadratically with ε in both cases, but the magnitudes are very similar. This is in sharp contrast to the behavior of the relative phase distributions where the dissipative coupling leads to much stronger features in P (ϕ − ) than the coherent coupling [see Fig. 2], with a linear rather than a quadratic dependence on ε.
These results suggest that the mutual information and the relative phase distribution characterize rather different aspects of the state of the coupled system. Whilst it is certainly to be expected that different measures of correlation between the micromasers will increase with the strength of the coupling, there is no reason why they should increase at the same rate.
Finally, we comment briefly on the extent to which entanglement is generated between coupled micromasers. We use the logarithmic negativity of the system to measure the entanglement [43], where the negativity is N (ρ) = 1 2 i (|λ i | − λ i ), with λ i the eigenvalues of ρ TA , the partial transpose of the density operator of the coupled micromaser system. Figure 6 shows how the entanglement behaves for both types of coupling.
We see in Fig. 6(a) that very little entanglement is in fact generated in this system. For both forms of coupling there is no entanglement except, interestingly, for the values of Θ that correspond to trapping states which occur at integer values, n, such that sin(φ √ n + 1) = 0. Figure  6(b) shows how the entanglement at the n = 1 trapping state (Θ = 4.97 for N = 5) changes with the coupling strength ε, and we see that there is very different behavior depending on how the micromasers are coupled. However, for both cases the entanglement remains very weak, even in the strongly quantum regime of very low photon numbers.

VII. CONCLUSIONS AND DISCUSSION
We have investigated synchronization effects in coupled micromasers using the relative phase distribution as a tool to quantify the emergence of preferred relative phases. We used perturbation theory to show that dissipative coupling between the micromasers leads to a 2π-periodic relative phase distribution whose peak grows linearly with the coupling. In contrast, for coherent coupling the phase distribution is π-periodic and there is a quadratic dependence on the coupling.
The relative phase distribution seems to be a very useful tool for describing synchronization effects in the quantum regime. We were able to derive a Fokker-Planck equation to describe its dynamics in the semiclassical regime, where photon numbers are large, leading to a steady-state distribution which depended on just the coupling strength and the linewidth of the uncoupled micromaser. This derivation also showed quite clearly that the phase dynamics would be much more complicated in the quantum regime where photon occupation numbers are low. Indeed, comparisons of numerical calculations using the full master equation showed that the Fokker-Planck equation substantially underestimated the strength of the features which emerge in the relative phase distribution in the quantum regime. Interestingly, a very similar underestimate of synchronization effects was obtained using a semiclassical model for the case of two van der Pol oscillators [7]. In our case, it seems that low photon occupation number is the key factor that leads to differences between quantum and semiclassical predictions. We also investigated the behavior of the mutual information of coupled micromasers and found that the behavior was rather different to that of the relative phase distribution. This is perhaps not surprising as the relative phase distribution depends on a very specific combination of a sub-set of the elements of the full density matrix of the system. Whilst all forms of correlation can be expected to increase with coupling, at least initially, it seems unlikely that other measures of correlation between the oscillators will increase in precisely the same way as the relative phase distribution.
Our work could serve as a starting point for a number of future studies. For example, it would be interesting to explore synchronization in the bistable regime where the micromasers exist in a mixed state consisting of limitcycles with two different amplitudes. Another possibility would be to explore synchronization effects in systems with more than two micromasers. The perturbation approach that we used here could prove a useful tool in analysing systems with a handful of coupled oscillators where a numerical solution of the full master equation already becomes very challenging because of the potentially very large state space involved. Finally, it would be interesting to investigate in detail the range of couplings which could be achieved in practice with micromasers, as well as solid-state analogs, and the best way to measure features in their relative phase distribution.