Localization control of few-photon states in parity-symmetric photonic molecules under balanced pumping

We theoretically investigate the problem of localization control of few-photon states in driven-dissipative parity-symmetric photonic molecules. We show that a quantum feedback loop can utilize the information of the spontaneously-emitted photons from each cavity to induce asymmetric photon population in the system, while maintaining a balanced pump that respects parity symmetry. To better understand the system's behaviour, we characterize the degree of asymmetry as a function of the coupling between the two optical cavities. Contrary to intuitive expectations, we find that in some regimes the coupling can enhance the population asymmetry. We also show that these results are robust against experimental imperfections and limitations such as detection efficiency.


Introduction
Spontaneous symmetry breaking (SSB) in quantum field theories is a mechanism whereby the ground state breaks one of the symmetries of its Hamiltonian. In many other contexts, SSB is defined broadly as the existence of a solution (not necessarily of lowest energy) that does not respect a symmetry associated with the underlying equations of motion; it often arises due to a nonlinear-induced transition between symmetric and symmetry-broken phases. In systems with many degrees of freedom, these symmetry breaking instabilities can lead to several intriguing phenomena such as spontaneous pattern formation (periodic modulations, spirals, vortices, etc) and soliton wave propagation. These effects can be observed in a variety of physical contexts such as fluids dynamics [1], chemical reactions [2], Bose-Einstein condensates (BECs) [3] and exciton-polariton systems [4].
One of the most important platforms for investigating different symmetry paradigms is classical nonlinear optics. There SSB in the form of modulation instabilities, formation of different localized solutions (bright/dark/vortex solitons, breathers, etc.) has been already observed and investigated thoroughly in both continuous [5] and discrete setups [6,7]. A minimal model for classical SSB in optics consists of two interacting modes under nonlinear conditions. Among several possibilities, this can be implemented for example by using Kerr nonlinearity in either: (1) two identical, linearly coupled waveguides (or cavities under external driving) [6]; or (2) single micro-ring resonators supporting two nonlinearly coupled (via cross-phase modulation) clockwise and counter-clockwise optical modes. In the first system, at low power levels, the symmetric solution (having equal intensity in both waveguides or cavities) is stable. As the power input level (or the driving pump in the case of cavities) is increased beyond a certain threshold, this rather symmetric state becomes unstable and the nonlinear self-trapping effects will amplify any small symmetrybreaking perturbation. Thus in this regime, the nonlinear stable eigenmodes exhibit asymmetric intensity distribution (more localization in one waveguide versus the other). In the micro-ring arrangements, the situation is rather similar with the instability breaking the chirality of the states (see [8,9] for recent experimental demonstration).
While SSB and multistabilities in nonlinear optical setups are well-studied, a more subtle situation arises when considering the quantum regime of such systems. Particularly, the linearity of quantum mechanics (interaction terms in the Hamiltonian still lead to linear terms in the Fock space description), together with pronounced quantum noise can blur the SSB effects. To further complicate matters, until recently controlled experiments that can probe this regime [10][11][12] were lacking.
Recently, however, an experimental work aiming at bridging this gap has demonstrated SSB in two coupled photonic crystal laser nanocavities with only 150 intracavity photons [11]. Subsequently, two theoretical studies modeled a coherently driven version of this experimental setup using driven-dissipative two-site Bose-Hubbard Hamiltonians [13,14]. These works showed that weak signatures of SSB can be traced The schematic for the considered setup: two coupled cavities with spontaneous emission measured using photodetectors. The laser jointly drives both cavities with a state-dependent Rabi frequency. The photodetectors D 1 and D 2 are monitoring the losses of cavity 1 and 2, respectively. The input of these detectors is fed into a device, which uses this input to determine the strength of the pumping.
in the quantum regime with few photons (1 ∼ 30). These experimental and theoretical results indicate that one will soon reach a regime where quantum effects will play a role in SSB experiments. In these studies the driving was chosen to be time-independent.
In this paper, we consider a similar system made of identical coupled optical cavities under time-dependent parity-symmetric coherent driving and investigate the following question: can one control the degree of statistical asymmetry (as defined by the timeaveraged photon population imbalance between the two cavities) without breaking the mirror-symmetry of the driving (see figure 1)? Counterintuitively, we show that this can be achieved by using quantum feedback control [15,16] with temporal pump modulation. It is important to note that while we allow the strength of the pump to vary with time, its profile remains mirror symmetric at any given instant.
The manuscript is organized as follows: in section 2 we present the theoretical model of our setup. We consider the driven, dissipative Bose-Hubbard dimer (BHD) as a simple model system, present our feedback scheme and provide a formal definition of our asymmetry measure. The feedback scheme is applied to our system in section 3, where we show that feedback can produce asymmetry, and we discuss the robustness of our setup. We present our conclusions and outlook in section 4.

Physical model
The driven, dissipative Bose-Hubbard dimer model includes coherent and incoherent evolution for the two cavities.
The coherent evolution is given by the Hamiltonian H = H BH + H drive with the Bose-Hubbard (BH) Hamiltonian [17] H BH and a driving part H drive : Here we have presented the Hamiltonians in the rotating frame with respect to the laser driving frequency. With a j we denote the annihilation operator for cavity j. The first term of equation (1) describes the detuning ∆ = ω c −ω L of the (equal) cavity frequencies ω c from the laser frequency ω L , the second term is an occupation-dependent nonlinear interaction with strength U/2 in each cavity, and the third term describes exchange of excitation with strength J between the cavities. The coherent driving symmetrically addresses both cavities with real Rabi frequency Ω(t), which can be time dependent.
The incoherent evolution in the model represents out-coupling of photons from the cavity mode to other environmental electromagnetic modes. We describe the evolution of the state ρ of the BHD using the following Lindblad master equation: where D[a j ]ρ = a j ρa † j − 1 2 a † j a j ρ + ρa † j a j terms describe the effect of the photon outcoupling on the system. Throughout the work we use γ as the unit of energy.
Since our feedback scheme, depicted in figure 1 and discussed in detail below, leads to stochastic contributions in the evolution of ρ it is convenient to use a stochastictrajectory-based description from the outset. This is, in particular, numerically advantageous as the dimension of the state is the square root of the dimension of the density matrix ‡. Individual trajectories are propagated with the following stochastic Schrödinger equation (SSE) [15]: where the time dependent populations of the two sites j = 1, 2 are given by and the stochastic increments dξ j have the following properties: Here E[· · ·] denotes the average over trajectories. The solution of equation (3) is approached by taking the average of many trajectories. The stochastic equation, ‡ For the calculation used in the present work we typically have a basis with dimension ∼ 400.
equation (4), was solved using the XMDS package [18,19]. Note that equation (6) implies that the stochastic increment is either 0 or 1. Besides being numerically more efficient, the SSE (4) is also a model for certain types of measurements which will be relevant for the feedback scheme presented below. This equation (4), with the given stochastic statistics in equations (6, 7), describes the system evolution with photodetection of the photons emitted from each cavity, assuming perfect detection efficiency and no delay in time between photon emission and detection [15]. A single quantum trajectory corresponds to a single run of an experiment. The detection events correspond to the stochastic noise increment dξ j (t) ∈ {0, 1} at any time t, for the detector at site j. The first term in equation (4) thus represents a photodetection event and corresponding 'jump' in the state when dξ j = 1, and no detection when dξ j = 0.
We will make use of this photodetection model for our feedback scheme. The photodetection record (history of 'clicks' from the detector) provides information about the population in each cavity [20]. In fact, given the photodetection record and knowledge of the Hamiltonian parameters, we have perfect knowledge of the quantum state of the system at any given time. We can then employ properties of the state in our feedback scheme to control the state evolution.

The feedback scheme
In our feedback scheme the value of the driving at time t is determined by the state of the system, i.e, Ω(t) → Ω · f (ψ t ), where Ω determines the overall strength of the driving and f (ψ t ) contains the time-dependence via ψ t ≡ |ψ(t) . Accordingly the driving Hamiltonian reads Through this driving the Hamiltonian H(t) in equation (4) becomes explicitly statedependent. Note that at any given time, the driving addresses both sites with the same driving strength. In the following we set Ω = γ. Thus all relevant information of the (time-dependent) driving strength is encoded in the feedback function f . We are now interested in what effect this state-dependent driving will have on the system, and in particular, how asymmetry can be engineered. To this end we consider a simple functional form for the driving feedback: Here the time-dependent total population and the time-dependent population difference are defined as The argument c ≡ {c 0 , c s , c a } indicates the parametric dependence of the populations (and f ) on the choice of the parameters c 0 , c s , c a of the feedback function. These coefficients can tune the state dependence from the symmetric case c a = 0 to more general asymmetric state-dependence. For c a = c s = 0, constant driving with no feedback is recovered.

Asymmetry measure
We now introduce an asymmetry parameter to characterize the imbalance in the photon population: with Here the summation over indicates the averaging over the time-average of M individual trajectories. The time-average is performed by integrating from t i to t f . We are primarily interested in the relative population difference. Thus, in equation (12), we have normalized by the average total population. Accordingly, the values of A( c) range from −1 to 1. If cavity 1 is higher populated than cavity 2, then A( c) is negative. Note that switching the sign of c a changes the sign of A( c). Thus for c a = 0 we expect that A = 0, and this is indeed what we find in the numerical simulations.
Henceforth when we refer to asymmetry we refer to the definition in equation (12).

Localization results
Now that we have defined the model and its governing equation, we return to our initial question: can one control the degree of asymmetry in a parity-symmetric setup without breaking the mirror-symmetry of the driving? To this end, we now consider a concrete example with J = 0.05γ, U = 0.3γ and ∆ = −1.65 γ (we discuss this choice of parameters in the Conclusions). Interestingly, even with such a crude parameter scan, one sees a clear evidence of localization control as characterized by a non-zero value of A. Furthermore, for each c 0 there is a broad, smoothly-connected region which gives maximal asymmetry A ≈ 0.4. This observation persists even if we repeat the simulations by varying c 0 . We emphasize here that this occurs despite the fact that we employ balanced driving that respects parity symmetry. What makes this result particularly interesting is that for few photons in the system (∼2-8, as seen in figure 4), quantum fluctuations play a significant role. The driving, despite its parity symmetry, makes use of these fluctuations to generate asymmetry in the photon population. This intriguing observation is the central result of this work.
Next, we investigate how the coupling coefficient impacts the degree of localization. From a computational perspective, finding large values for the asymmetry measure A (ideally the maximum value) for a given set of parameters (∆, U and J) requires an expensive search in the three dimensional feedback parameter space spanned by the vector c. In order to overcome this difficulty, we employ the Nelder-Mead optimization scheme [25] (see appendix A for more details). This approach reduces the computational cost considerably while yielding large values for the asymmetry parameter A, which we believe to be very close to the optimal asymmetry that can be achieved (see appendix A).  . Naively, one would expect that A should drop as a function of |J| since coupling provides a pathway for photons to move between the two cavities, making it harder for the control to favor one cavity over the other. Indeed this is the case for positive values of J as can be seen in Fig. 3 (a). Surprisingly however, for negative values of J, we find that the asymmetry increases with |J|, approaching large values (one should keep in mind that a value of 1.0 means perfect localization on one cavity). Our simulations thus indicate that the naive picture is not always correct and that coupling can provide a pathway for increasing asymmetry. For each value of J, the optimized values of the feedback parameters are shown in Fig.  3 (b). One sees that these parameters stay roughly constant over the entire J range,  the largest variation being in c 0 . For completeness, we have also computed the values of A with no feedback (c s = c a = 0), where we found that A = 0. As expected, we also found that A = 0 for c a = 0.

Robustness of the feedback
While the main objective of the present work is to demonstrate that via symmetric feedback one can obtain quite strong asymmetry in the populations of the cavities, it is nevertheless interesting/important to discuss the robustness of this scheme with respect to imperfections that would always be present in an experimental realization.
Variation in the control parameters c: One important aspect is how sensitively the asymmetry A( c) depends on c. More precisely, will small deviations from the values found from optimization still yield similar asymmetry? To gain insight into this question, it is instructive to look at figure 2. There we see that the region around the maximum is quite broad and smooth. Thus small variations around the maximum will not change the value of the asymmetry significantly. For various J we have checked that small deviations from our numerically optimized values c indeed do not change the asymmetry significantly.
Changes in the driving field are not infinitely fast: In experiments the change in the driving cannot be arbitrarily fast. To take that into account we define the maximum The ideal driving control varies on a timescale similar to that of the mean populations of the two sites, as can be seen from equation (9) and in figure B1. One expects that a reduction of this rate of change will lower the achievable asymmetry. The effect of limiting the rate of change r max in the applied driving is demonstrated in figure 4(a). One sees that even for a slow maximal rate compared to dn 1,2 /dt, one can still obtain pronounced asymmetry. In the inset one nicely sees the finite rate of change in the driving function f (t).
Detection efficiency: When a cavity spontaneously emits a photon, only a fraction of these spontaneous emission events is collected by the respective photodetector. In the following, we use η to denote the fraction of detected photons. The ability to detect the photons emitted from the cavities plays a crucial role in our feedback scheme. Thus one expects that imperfect detection efficiency (η < 1), which is always present in experiments, will strongly influence our results. However, numerically we find that even for quite small efficiencies (η ∼ 0.1) the feedback scheme works surprisingly well, as can be seen in figure 4(b). In this figure we present numerical calculations for the same parameters as in figure 3 for different values of η. Remarkably, we also used the same c values which are optimized for the perfect detection case. We suspect that even higher asymmetry could be achieved if we optimize taking the known detector efficiency into account. It is also instructive to take a look at the inset where a single trajectory with the corresponding driving is shown for η = 0.1. Comparing to the trajectory and feedback function with η = 1 (see figure B1) one sees that for η = 0.1 the feedback is adjusted less often, as expected.
Our numerical implementation is based on the method discussed in pages 190-191 of Ref. [15].
The Hamiltonian and the coupling to an environment are not exactly known. Our feedback scheme relies on the fact that we are able to infer the population of the cavities from measurements. The situation becomes complicated if the parameters entering the Hamiltonian are not known exactly. One also has to keep in mind that our model Hamiltonian is always an approximation to the real experimental system. We suspect that our feedback scheme will still work if the deviations from our ideal situation are not too large. However, such a study (as in [21][22][23]) is beyond the scope of the present work. For each experimental setup this has to be individually investigated.

Conclusions
In this work, we have investigated the problem of few-photon localization in a system of two mirror-symmetric coupled optical cavities modeled by a driven-dissipative Bose-Hubbard Hamiltonian with many-body interactions. We have shown that, in this regime, state-dependent quantum feedback control can be utilized to promote an asymmetric state with more photons in one cavity than the other, even when the driving itself is kept symmetric. Intuitively, this result can be better understood by reference to Maxwell's demon. In the famous problem of gas molecules confined in a container divided symmetrically by a wall having a small hole, the demon can measure the position of the molecules on both sides of the wall and use their random motion to control the opening of the hole to force the molecules to move from one side to the other. Analogously, our scheme is based on a similar principle, except that here randomness originates from quantum fluctuations rather than thermal noise.
In the manuscript, we have demonstrated localization using an exemplary set of parameters. This parameter set is not unique. We have found that the feedback generates localization for the following regions about these parameter values: U/γ = 0.3 ± 0.15 and ∆/γ = −1.65 ± 0.6.
We now comment on the experimental feasibility of our results using optical setups. The main challenge here would be the required fast feedback timescales. Particularly, our scheme relies on feedback in the form of continuous, time-dependent driving whose strength is modulated in time as a function of the cavity populations. Practically, this temporal modulation could be introduced by electro-optic modulators. Typical time scales of these devices are in the order of 1 ns. If we consider a system similar to that studied in [11] and having a damping rate of ∼ 7 ps, we find that a time scale of 0.7 ns for variation in the driving (to increase from its minimum (f min = 0.01) to its maximum value (f max = 4.0)) achieves reasonable results (with r max = 0.04 as discussed in figure 4). Clearly, the required time scales are not far from those provided by current on-shelf electronic components. It is thus foreseeable that near-future electronic technologies will be adequate for use in such experiments.
An attractive alternative would be polariton-based setups in cavity QED [24], which can have a damping timescale on the order of 1 µs. This would render a feedback delay timescale of 1 ns negligible. Here N T is the total number of trajectories and N a is the number of 'sub-ensembles' in which we group these trajectories. For each sub-ensemble j the asymmetry A j is calculated in accordance with the description in the main text (see equation (12) and the discussion afterwards). That means for each trajectory k from sub-ensemble j we integrate n j,k diff and n j,k tot from time t i to t f which are the same for each trajectory. In the simulations shown we have used t i = 0.4 and t f = 200 (in units of γ). The choice t i = 0.4 removes the short initial transients.

Acknowledgments
For N a = 1 we have A opt = A from equation (12). If the ratio N T /N a is large then also A opt ≈ A. For the values we have used (N T = 1000 and N a = 10) we have found that A opt ( c) is almost identical to A( c).
Since the different A j ( c) are calculated from independent trajectories, their values will be different. We use the spread of these values to estimate the accuracy of A opt ( c). To this end we consider the standard deviation Clearly, by increasing the total number of trajectories N T the accuracy will increase. For our choice of t i , t f , N T , and N a we find σ Aopt values on the order of 0.003.
All results presented are calculated using N S = 30 iterations of the optimization algorithms; that means A opt ( c) is evaluated for N S different choices of c. The c which provides the largest value is then used to calculate A( c). The optimized c parameters used in figure 3 are shown in table A1. The initial parameter ranges for c roughly correspond to the ranges shown in figure 2 of the manuscript.
We believe that our optimization routine produces results close to the global maximum for the asymmetry. Firstly, the parameter scan in figure 2 for J = 0.05γ shows a maximal asymmetry region. Our optimization routine produces feedback parameters corresponding to this region, with a matching asymmetry value, within the given accuracy for figure 2. The feedback parameters are given in table A1 and plotted in figure 3. Secondly, repeated or additional optimization runs have produced consistent results.

Appendix B. Quantum trajectories
Here, we outline the calculation of the system dynamics, and present an illustrative example.  Figure B1.
Trajectories of the mean cavity populations for a particular photodetection record are shown for no feedback (a) and with feedback (b). The 'feedback functions' f (t) are also shown (green curves). In (a) it has a time constant value f (t) = 1.5. The population difference n diff (t) is shown in the middle panel with the probability to find a certain value of n diff in the given trajectory. In the right panel, the histogram p(n diff ) is shown, which is obtained by averaging over many trajectories (we obtain the same result if we consider a very long trajectory). Here, U = 0.3γ, ∆ = −1.65γ and J = 0.0. As an example of the asymmetry measure A (equation (12)) we find A = 0.00 and A = 0.46 for (a) and (b) respectively.
respectively. In figure B1, we illustrate the difference between quantum trajectories with constant driving (a) and with feedback (b). For the constant driving case (a) we have used f (t) = 1.5 = const. In (b) we have used c 0 = 8.61, c s = 2.07 and c a = 1.77. This choice of parameters resulted in a feedback function f (t) that has the same mean T 0 dtf (t)/T = 1.5 as the constant driving case (a).
In the left panel results of single trajectories are displayed and n 1 and n 2 are shown together with the feedback function. The initial state was both cavities empty, i.e. n 1 (0) = n 2 (0) = 0. In the middle panel the population difference n diff (t) is shown together with the probability to find a certain value of n diff in the respective trajectory. Finally in the right panel we show the corresponding histogram p(n diff ), obtained by averaging over many trajectories (we obtain the same result if we consider a longer trajectory). Here we have also removed the first 0.4 time-units during which the initial transient dynamics takes place §.
For the case without feedback (a), most of the time the cavities have equal populations, resulting in a strong peak at n diff = 0 in the histogram. Weak shoulders can be seen symmetrically around n diff = 0, indicating that the system sometimes does not have the same population in the two cavities. Note that the symmetry means that the probability that cavity 2 has more photons than cavity 1 (n diff > 0) is the same as the probability that cavity 1 has more photons than cavity 2 (n diff > 0).
In the case with feedback (b), one clearly sees in the trajectory shown that now cavity 1 consistently has a much higher population than cavity 2. This imbalance is still present when sufficiently many long trajectories are considered (we average 1000 trajectories with duration 200/γ). We see that the histogram has a large broad peak around n diff ≈ 5 and is strongly asymmetric around n diff = 0. The mean value is 2.65.