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

The driven, dissipative BHD 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 [27] H_BH and a driving part H_drive:

$$\begin{eqnarray}&&{H}_{\mathrm{BH}}=\displaystyle \sum _{j=1}^{2}({\rm{\Delta }}{a}_{j}^{\dagger }{a}_{j}+\displaystyle \frac{U}{2}{a}_{j}^{\dagger }{a}_{j}^{\dagger }{a}_{j}{a}_{j})-J({a}_{1}^{\dagger }{a}_{2}+{a}_{2}^{\dagger }{a}_{1}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{H}_{\mathrm{drive}}(t)={\rm{\Omega }}(t)\displaystyle \sum _{j=1}^{2}({a}_{j}+{a}_{j}^{\dagger }).\end{eqnarray} \tag{ 2 }$$

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:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-{\rm{i}}[H(t),\rho ]+\gamma \displaystyle \sum _{j}{ \mathcal D }[{a}_{j}]\rho ,\end{eqnarray} \tag{ 3 }$$

where ${ \mathcal D }[{a}_{j}]\rho ={a}_{j}\rho {a}_{j}^{\dagger }-\tfrac{1}{2}({a}_{j}^{\dagger }{a}_{j}\rho +\rho {a}_{j}^{\dagger }{a}_{j})$ terms describe the effect of the photon out-coupling 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 stochastic-trajectory-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) [25]:

$$\begin{eqnarray}{\rm{d}}| \psi (t)\rangle & = & \displaystyle \sum _{j}{\rm{d}}{\xi }_{j}(t)\left(\displaystyle \frac{{a}_{j}}{\sqrt{{n}_{j}(t)}}-1\right)| \psi (t)\rangle \\ & & +{\rm{d}}t\left(\displaystyle \frac{\gamma }{2}\displaystyle \sum _{j}({n}_{j}(t)-{a}_{j}^{\dagger }{a}_{j})-{\rm{i}}{H}(t)\right)| \psi (t)\rangle ,\end{eqnarray} \tag{ 4 }$$

where the time-dependent populations of the two sites j = 1, 2 are given by

$$\begin{eqnarray}&&{n}_{j}(t)=\langle \psi (t)| {a}_{j}^{\dagger }{a}_{j}| \psi (t)\rangle \end{eqnarray} \tag{ 5 }$$

and the stochastic increments ${\rm{d}}{\xi }_{j}$ have the following properties:

$$\begin{eqnarray}&&{\rm{d}}{\xi }_{k}(t){\rm{d}}{\xi }_{j}(t)={\rm{d}}{\xi }_{k}(t){\delta }_{{kj}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{ \mathcal E }[{\rm{d}}{\xi }_{j}(t)]=\gamma {n}_{j}(t){\rm{d}}t.\end{eqnarray} \tag{ 7 }$$

Here ${ \mathcal E }[\cdots ]$ denotes the average over trajectories. The solution of equation (3) is approached by taking the average of many trajectories. The stochastic equation, equation (4), was solved using the XMDS package [28]⁵ . 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 [25]. A single quantum trajectory corresponds to a single run of an experiment. The detection events correspond to the stochastic noise increment ${\rm{d}}{\xi }_{j}(t)\in \{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 ${\rm{d}}{\xi }_{j}=1$, and no detection when ${\rm{d}}{\xi }_{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 [29]. 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.

In our feedback scheme the value of the driving at time t is determined by the state of the system, i.e, ${\rm{\Omega }}(t)\to {\rm{\Omega }}\cdot f({\psi }_{t})$, where Ω determines the overall strength of the driving and f(ψ_t) contains the time-dependence via ${\psi }_{t}\equiv | \psi (t)\rangle$. Accordingly the driving Hamiltonian reads

$$\begin{eqnarray}&&{H}_{{\rm{drive}}}(t)\to {H}_{{\rm{drive}}}({\psi }_{t})={\rm{\Omega }}\cdot f({\psi }_{t})\displaystyle \sum _{j=1}^{2}({a}_{j}+{a}_{j}^{\dagger }).\end{eqnarray} \tag{ 8 }$$

Through this driving the Hamiltonian H(t) in equation (4) becomes explicitly state-dependent. 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:

$$\begin{eqnarray}&&f({\psi }_{t};\vec{c})={c}_{0}-{c}_{{s}}\cdot {n}_{\mathrm{tot}}(t;\vec{c})-{c}_{{a}}\cdot {n}_{\mathrm{diff}}(t;\vec{c}).\end{eqnarray} \tag{ 9 }$$

Here the time-dependent total population and the time-dependent population difference are defined as

$$\begin{eqnarray}&&{n}_{\mathrm{tot}}(t)={n}_{1}(t)+{n}_{2}(t),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{n}_{\mathrm{diff}}(t)={n}_{2}(t)-{n}_{1}(t).\end{eqnarray} \tag{ 11 }$$

The argument $\vec{c}\equiv \{{c}_{0},{c}_{s},{c}_{a}\}$ indicates the parametric dependence of the populations (and f) on the choice of the parameters c₀, 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.

We now introduce an asymmetry parameter to characterize the imbalance in the photon population:

$$\begin{eqnarray}&&A(\vec{c})=\displaystyle \frac{{N}_{\mathrm{diff}}(\vec{c})}{{N}_{\mathrm{tot}}(\vec{c})},\end{eqnarray} \tag{ 12 }$$

with

$$\begin{eqnarray}&&{N}_{{\rm{diff}}}(\vec{c})=\displaystyle \frac{1}{M}\displaystyle \sum _{{\ell }=1}^{M}\displaystyle \frac{1}{T}\left[{\int }_{{t}_{i}}^{{t}_{f}}{n}_{{\rm{diff}}}^{({\ell })}(t;\vec{c}){\rm{d}}{t}\right],\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{N}_{{\rm{tot}}}(\vec{c})=\displaystyle \frac{1}{M}\displaystyle \sum _{{\ell }=1}^{M}\displaystyle \frac{1}{T}\left[{\int }_{{t}_{i}}^{{t}_{f}}{n}_{{\rm{tot}}}^{({\ell })}(t;\vec{c}){\rm{d}}{t}\right].\end{eqnarray} \tag{ 14 }$$

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, and T = t_f − t_i.

Equations (13) and (14) provide the unconditional population difference and total population, respectively. Unconditional means that these values do not depend on a particular photodetection record and the corresponding (conditional) feedback; this independence is desirable for our measure. The independence from the stochastic evolution of individual trajectories is achieved through the averages over trajectories and over time. It can similarly be achieved by averaging over a single long trajectory, neglecting the initial transient behavior.

We are primarily interested in the unconditional, relative population difference. We consider that the averaged (unconditional) difference between the cavity populations is meaningful when scaled by the total cavity populations, as implemented in equation (12). Accordingly, the values of $A(\vec{c})$ range from −1 to 1. If cavity 1 is higher populated than cavity 2, then $A(\vec{c})$ is negative. Note that switching the sign of c_a changes the sign of $A(\vec{c})$. Thus for c_a = 0 we expect that A = 0, and this is indeed what we find in the numerical simulations. This includes the case of no feedback.

Henceforth when we refer to asymmetry we refer to the definition in equation (12).

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).

Figures 2(a)–(c) show the asymmetry measure A as function of the feedback parameters c_a and c_s for three different values of c₀. Each point on these figures has been generated by using 100 trajectories with a time span of 200/γ. Note that c_a = 0 gives A = 0, as expected.

Zoom In Zoom Out Reset image size

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₀ 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₀. We emphasize here that this occurs despite the fact that we employ balanced driving that respects parity symmetry. Quantum fluctuations cause (transient) imbalanced photon populations. 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. Furthermore, we have only few photons in the system (∼2–8, as seen in figure 4), such that the quantum fluctuations cause dramatic changes in the cavity populations. Fluctuations that are large compared to the cavity populations can quickly lead to population asymmetry, but this asymmetry is also unstable. Nonetheless, in our setup the driving is still effective in generating asymmetry. Note that relatively small fluctuations (such as for higher photon populations) can also result in a large and stable asymmetry in the populations (using feedback), but the timescale for this asymmetry to arise is longer than with larger fluctuations.

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 $\vec{c}$. In order to overcome this difficulty, we employ the Nelder–Mead optimization scheme [30] (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).

Figure 3(a) plots the asymmetry measure A as a function of the coupling J when Δ = −1.65 and U = 0.3 (same parameters as above). 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 figure 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 figure 3(b). One sees that these parameters stay roughly constant over the entire J range, the largest variation being in c₀. 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.

Zoom In Zoom Out Reset image size

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 $\vec{c}$: One important aspect is how sensitively the asymmetry $A(\vec{c})$ depends on $\vec{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 $\vec{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 rate of change by

$$\begin{eqnarray}&&{r}_{\max }=| {\rm{d}}f/{\rm{d}}t| .\end{eqnarray} \tag{ 15 }$$

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).

Zoom In Zoom Out Reset image size

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 $\vec{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 pp 190–191 of [25].

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. In particular, given a photodetection record, we propagate the wave function according to equation (4) (or for imperfect detection, we propagate the stochastic master equation as in pp 190–191 of [25]), imposing jumps whenever a photon is detected. From the known state, we calculate the cavity populations. 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 [31–33]) is beyond the scope of the present work. For each experimental setup this has to be individually investigated.

Individual trajectories are calculated by propagating the state vector $| \psi (t)\rangle$ according to the SSE (4). Both cavities are initially empty, before pumping populates them with photons. The trajectories become different because of the stochastic events corresponding to photon emission. Through the feedback $f({\psi }_{t},\vec{c})$ (from equation (9)) the evolution of the state also depends parametrically on the choice of $\vec{c}$.

For numerical reasons, we impose the bounds ${f}_{\min }\leqslant \,f({\psi }_{t},\vec{c})\leqslant {f}_{\max }$, with f_min = 0.01 and f_max = 4.0. If equation (9) gives a value of f(t) larger (smaller) than f_max (f_min) then f(t) is set to be f_max (f_min). The lower bound prevents the cavities from remaining empty. The upper bound prevents occupation of very high photon-number states in either cavity.

For the parameters that we consider, we can truncate the Hilbert space of the state vector to less than 20 photons in each cavity. Thus the dimension of the Hilbert space is 20 × 20 = 400.

To demonstrate the effect of the feedback we consider the case of uncoupled cavities, i.e., J = 0. For the detuning and the nonlinearity we choose Δ = −1.65γ, and U = 0.3γ, 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₀ = 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 ${\int }_{0}^{T}{\rm{d}}{tf}(t)/T\,=\,1.5$ as the constant driving case (a).

Zoom In Zoom Out Reset image size

In the left panel results of single trajectories are displayed and n₁ and n₂ are shown together with the feedback function. The initial state was both cavities empty, i.e. n₁(0) = n₂(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.