Dissipative Josephson effect in coupled nanolasers

Josephson effects are commonly studied in quantum systems in which dissipation or noise can be neglected or do not play a crucial role. In contrast, here we discuss a setup where dissipative interactions do amplify a photonic Josephson current, opening a doorway to dissipation-enhanced sensitivity of quantum-optical interferometry devices. In particular, we study two coupled nanolasers subjected to phase coherent drivings and coupled by a coherent photon tunneling process. We describe this system by means of a Fokker-Planck equation and show that it exhibits an interesting non-equilibrium phase diagram as a function of the coherent coupling between nanolasers. As we increase that coupling, we find a non-equilibrium phase transition between a phase-locked and a non-phase-locked steady-state, in which phase coherence is destroyed by the photon tunneling process. In the coherent, phase-locked regime, an imbalanced photon number population appears if there is a phase difference between the nanolasers, which appears in the steady-state as a result of the competition between competing local dissipative dynamics and the Josephson photo-current. The latter is amplified for large incoherent pumping rates and it is also enchanced close to the lasing phase transition. We show that the Josephson photocurrent can be used to measure optical phase differences. In the quantum limit, the accuracy of the two nanolaser interferometer grows with the square of the photon number and, thus, it can be enhanced by increasing the rate of incoherent pumping of photons into the nanolasers.


Introduction
Among the striking effects of quantum coherence, the Josephson effect is one of widest used in nowadays technologies [1]. To name a few examples, the SQUID (superconducting interference device) is one of the oldest and most sensitive magnetic sensors [2,3], and Josephson junctions are integral building blocks to construct artificial two-level systems in quantum information processing [4,5]. In precision metrology, the Josephson effect has been used as a practical standard of voltage and the elementary charge, e [6,7].
The Josephson effect occurs when two quantum systems having both well-defined quantum phases, φ 1,2 , are weakly coupled so that quantum tunneling is enabled between them. It is manifested as a net current, between the two subsystems depending on the phase difference, ∆φ = φ 1 − φ 2 . This effect was discovered by Josephson, who predicted a macroscopic electric current along a superconducting tunnel junction [8,9]. Since then, extensions of these ideas have been proposed and tested in several platforms, such as Bose-Einstein Condensates [10,11,12,13,14,15], superfluid 3 He [16], photonic [17,18,19] and optomechanical systems [20] and polaritons [21,22]. While the effect of dissipation has been studied in some cases [17], it has generally been considered as detrimental for the observation of the Josephson current and its applications.
In this work we unveil a dissipative Josephson in a fully photonic setup consisting of two coupled nanolasers. We show that this system presents interesting non-equilibrium phases and it also has an exciting outlook for applications in quantum metrology and sensing. In particular, we introduce a model of an interferometer consisting of two single qubit lasers coherently coupled through a photon tunneling term, as depicted schematically in Fig. 1. Nanolasers are probably the most fundamental example of active dissipative systems with a non-trivial non-equilibrium phase diagram. They can be implemented in several platforms such as photonic [23,24,25,26], plasmonic [27], or nano-mechanical systems [28]. A natural extension from the single nanolaser model into the many-body regime arises when we consider networks of local nanolasers coupled by means of photon tunneling terms of the form H t = t(a † 1 a 2 + a 1 a † 2 ), [29,30]. Here, interesting phenomena may appears as a result of the interplay between coherent tunneling and on-site non-linear dissipative terms. In our setup those nonlinear dissipative terms are responsible for sustaining the quantum coherence that generates the Josephson effect and its possible applications in interferometry.
This article presents the following results: (i) We derive a semi-classical description of a dissipative interferometer consisting of two nanolasers coupled through a coherent photon tunneling process and subjected to coherent drivings with phase difference, ∆φ. (ii) We identify two limiting regimes of the steady-state of this system. If the coherent coupling is small, the lasers are phase locked to each individual coherent driving. As we increase the coherent coupling between nanolasers, the system goes through a nonequilibrium phase transition into a steady-state in which phase coherence is lost. (iii) We present analytical and numerical evidence of the existence of a photonic Josephson current between the two nanolasers. This current causes a photon number imbalance in the photonic steady-state proportional to sin(∆φ). (iv) We analyse the performance of this system as an interferometer and we show that there is an amplification effect by which the accuracy grows as √ n in the shot noise limit, with n the number of photons in the nanolasers. Furthermore, the accuracy is optimal close to the critical point of the lasing phase transition.

Dimer of single-qubit lasers
Our system consists of two single-qubit lasers coupled by coherent photon tunneling. This scheme can be implemented in several setups including circuit QED (for example using the ideas proposed in [31]) and trapped ion phonon lasers [28]. The discussion below is focused on the case of optical nanolasers, however, our results are independent of any particular realization in photonics, vibronics or optomechanics. Figure 1. (a) General scheme: two single-qubit lasers are coupled by a coherent photon hopping term with rate t c . Each nanolaser is and subjected to incoherent qubit pumping with rate γ, photon loss with rate κ, and a periodic driving field with amplitude 1,2 = | |e iφ1,2 . A Josephson photocurrent is generated if there is a finite difference between the optical phases, θ 1 , θ 2 , at each nanolaser. (b) Conventional Josephson effect: an electric current is generated across a junction separating two superconductors with different phases.
Each nanolaser laser consists of a two-level system (qubit) with levels |g and |e , coupled on resonance with a cavity mode through a Jaynes-Cummings type interaction. The two-level system is incoherently pumped with a rate γ and the photonic mode has a decay rate κ. In addition, the cavity interacts with a weak coherent driving on resonance with the photonic mode. The two nanolasers are connected by a photon coherent tunneling term that couples the photonic cavities. In an interaction picture rotating at the mode frequency, the following master equation for the system density matrix, ρ, captures the complete quantum dynamics of the system (we use units such We use the notation The Hamiltonian in Eq. (1) is The first term represents the qubit-field coupling, with strength g. The second one describes external driving terms acting on the nanolasers, with driving amplitudes j = | |e iφ j . In this work we will consider that both driving amplitudes have the same strength, | |, but there may be a phase difference, ∆φ = φ 2 − φ 1 . Finally, the last term represents the coherent photon tunneling term, with amplitude t c , This coherent coupling occurs in systems of single-mode nano-cavity arrays [32,29], superconducting circuits [30] or nano-mechanical systems such as trapped ions (where it takes the form of a phonon-tunneling term, see [33,34]). The minus sign is added to the coupling so that for positive t c the lowest energy mode is the symmetric or center of mass mode. At first sight, we would expect a term like Eq. (2) to induce some synchronization of the phase between nanolasers, however, we will prove later on that this intuition is wrong and, actually, strong tunneling constants, t c tend to destroy the quantum coherence in the system. The steady-state of each nanolaser is governed by the parameter such that, for C p < 1, nanolasers are in a non-lasing steady-state, whereas C p > 1 is the lasing phase, with a number of photons that roughly scales like γ/κ, as we show below.
In this paper, we will assume that the local nanolaser dynamics is much faster than the photon tunneling term, and in particular, t c γ, such that we can safely assume that the lasing transition stays at C p = 1. The mean-field prediction for the number of photons is [35], From Eq. (4), we learn that the parameter that determines the effective size of our nanolaser model is actually the ratio γ/κ, which determines any finite-size scaling effects and plays a role that is analogous to the number of sites in a quantum lattice model (see Appendix A).

Effective non-linear photonic master equation
We expect the most interesting physics occurring in a lasing regime of large photon numbers in which each nanolaser can be approximated by a self-sustained quantum oscillator. This regime may be attained for a strong pumping of the qubit such that All along this paper we will be working in this regime, in which we can simplify our model by adiabatically eliminating the qubit dynamics [36]. This step will ultimately allow us to obtain a semiclassical description in terms of a Fokker-Planck equation. After the qubit's adiabatic elimination we get the following master equation (see Appendix A for details),ρ with the coefficients: ρ f = Tr qubit (ρ) is the reduced density matrix of the photonic subsystem, obtained after tracing out the qubit. The last two terms in Eq. (6), proportional to the coefficient B, account for the non-linear matter-light interaction and they are ultimately responsible for the non-equilibrium phase transition into the lasing phase. Eq. (6) is strictly valid only below the critical point, C p < 1, and slightly above it, C p 1, since it has been obtained under the assumption of total qubit population inversion (see Appendix A). To quantify better the regime of validity of this approximation we can calculate the ground-state population of the qubit in the steady-state, σ + j σ − j ss , and check whether it can really be neglected. Actually, by using the equations for the single-qubit case, we find [37], We confirm that condition σ + j σ − j ss 1 is met close to the phase transition at C p 1, and in the non-lasing phase (C p < 1). It is, however, highly desirable to extend Eq. (6) well into the lasing phase (C p > 1), to fully understand the systems's behavior. This can be done with a proper renormalization of the coefficient that takes into account the neglected terms in the adiabatic elimination. As we show in Appendix A, this procedure amounts to replacing in Eq. (6). The new parameter B r includes the effect of processes of higher order in g, and it ensures the right prediction on the photon number in the steady-state.

Semiclassical Fokker-Planck equation
Although we are dealing with only two nanolasers, the solution of Eq. (6) is computationally demanding because a very high number of photonic Fock states must be included in any exact numerical calculation. However, in the limit of high-photon numbers an analytical approach based on phase-space methods can be used to further simplify Eq. (6). Concretely, we shall employ the Glauber-Sudarshan P representation [36] of the effective master equation. This representation is defined as the pseudoprobability distribution satisfying where |α is the coherent state |α = exp (αa † − α * a)|0 . The function P (α, α * ) plays the role of a classical probability distribution over |α α|, with the normalization condition d 2 αP (α, α * ) = 1. Expectation values of normal ordered operators can be evaluated with the identity, The substitution of Eq. (10) into Eq. (6) transforms the master equation into a Fokker-Planck equation for P (α, α * , t) [38], see Appendix C. We can achieve further simplification by working with polar coordinates, α j = r j e iθ j (j = 1, 2).
The radial components, r j , are related to the photon number observable in each cavity, In the last equation and along the rest of this work, we will understand O as refering to both quantum average, or average with respect to the semiclassical distribution, P , depending on whether O is an operator or a Fokker-Planck variable. The key to simplifying the Fokker-Planck equation is the observation that the fluctuations in the radial components, r j , may be neglected as long as we are deep enough in the lasing regime (C p > 1 and n j 1). Thereby we can assume that radial variables remain close to their steady-state values r j ≈ r 0 j . In this case we can trace out the radial variables and consider a reduced description in terms of phase variables, θ j . Formally, this is accomplished by assuming a factorized P (α 1 , α 2 ) ≈ R(r 1 )R(r 2 )P θ (θ 1 , θ 2 ). Each R(r j ) is a Gaussian distribution properly normalized around r 0 j (average value of r j ), which corresponds to the radial probability distribution of each nanolaser in the lasing regime. This procedure is discussed in details in Appendix C, and leads to and equation that depends on the phases of the nanolasers only, in which n 0 j = (r 0 j ) 2 stands for the steady-state average number of bosons at site j. In the lasing regime and in the absence of tunneling, this quantity is independent of the site and is given by n 0 j = n 0 , with In the homogeneous case we have r 0 j = r 0 , which leads to a homogeneous Fokker-Planck equation for the phases, representing a uniform photon density together with a single phase diffusion rate, We will see later that this picture has to be corrected to account for photon imbalance induced by the Josephson current between nanolasers. Eq. (16) is one of the most important results for our work. We remark that the novel element in this equation is the coupling between phases induced by t c . Actually, this equation is closely related to the dissipative Kuramoto model, however, in our phase model there is a coherent coupling term which differs from the usual dissipative couplings considered in coupled laser or synchronization models. This will have severe consequences in the non-equilibrium phase diagram of the model, as we see below.

Dissipative phase transition induced by coherent photon tunneling
We investigate now the effect of the coherent photon coupling in the steady-state of our Fokker-Planck equation (16) and show that, surprisingly, it does not lead to any synchronization effect between the nanolasers. On the contrary a coherent photon tunneling process leads to the loss of quantum coherence in the system.

Phase-locked and non-phase-locked steady-states
To simplify the discussion we consider first an homogeneous driving with φ 1 = φ 2 = φ, and r 0 1 = r 0 2 = r 0 . In Eq. (16) there are two limiting cases: (i) /r 0 t c . In this limit we expect that the system is well approximated by two independent single qubit lasers phase-locked to the driving fields, as explored in Ref. [37]. In this case the Fokker-Planck equation can be exactly solved and we get Hence, here we find non-zero coherences | a j | = 0. We will refer to this steadystate as a phase-locked (PL) phase.
In the NPL phase the photon tunneling between nanolasers is the dominant process, and it induces a total loss of phase coherence in the system.
(ii) t c | |/r 0 , where we expect that the coherent tunneling dominates the system's dynamics. To understand this limit it is useful to study Eq. (16) with | | = 0, Eq. (19) can be solved by a change of variables to collective coordinates, In the new variables we find the steady-state solution, with vanishing coherences, a j = 0. Furthermore, in this limiting case, the twopoint correlation G = a † 1 a 2 also becomes zero, which implies the absence of any synchronization between the modes [39]. We will refer to a steady-state where t c is dominant, and the system does not retain any coherence, as non-phase-locked (NPL) phase.
We arrive to the somehow counter-intuitive conclusion that the presence of a coherent photon tunneling term does not induce any correlations between nanolasers. This situation differs strongly from the case of a dissipative coupling as induced, e.g. by incoherent tunneling through evanescent modes or an intermediate lossy cavity [40], in which coupling does induce a phase correlation between nanolasers.

Phase transition between PL and NPL steady-states: driving with homogeneous phases
What happens in the intermediate regime between the PL and NPL steady-states identified above? To address this question, we assume that we have a constant driving term, , and we increase the tunneling from t c /r 0 to t c /r 0 . In particular, we are interested to know whether a dissipative phase transition separates the two phases. The coherences | a j | can be used as the order parameter to distinguish between the PL and NPL steady-states. We also need to define a valid thermodynamic limit to establish the existence of critical properties. Even though this is a two-site system, a thermodynamic limit is obtained by letting the number of photons in the steady state, n 0 , play the role of the system size [41,37]. The number of photons is essentially regulated by the ratio γ/κ as shown in Eq. (4), so that the thermodynamic limit will be reached in a limit of strong pumping, γ κ. We have solved numerically Eq. (16) by discretizing the angular variables, θ 1 , θ 2 , in a number of values, n c , running from 0 to 2π. This procedure allows us to express P θ as a vector and ∂P θ /∂t as a non-Hermitian matrix, and to calculate numerically the steady-state solution, ∂P θ /∂t = 0. We expect that this numerical method is accurate as long as n c 2π or, equivalently, that the angular probability distribution function does not change much within a phase interval ∆θ = 2π/n c . We have checked that in all the calculations shown in this work the numerical results have converged for the values of n c used. Numerical calculations allow us to find the reduced probability distributions, where in the last terms we have re-expressed the reduced probability for α 1 in polar coordinates. An analogous definition holds for P 2 (α 2 ). In Fig. 3 we show our numerical results for reduced probability distributions for increasing coherent coupling t c . We observe two main qualitative effects. First, the probability distribution spreads in phase space as we increase t c , in agreement with the expected transition from the PL to the NPL phase. Second, we observe a rotation of the probability distribution: at values t c ≈ 0 the two nanolasers are phase-locked to θ = 0, however as we increase t c , the distribution rotates to an increasing angle around θ = π/2. This effect can be qualitatively understood with the equations of motion for the coherences of the bosonic modes, All non-linear and dissipative effects are included in the single nanolaser contribution to the time-evolution of the coherences, d a j /dt| nl . The second and third terms in the r.h.s. of Eqs. (23,24) are the external drives and coherent coupling terms, respectively. Note that the external drives -proportional to | | -contribute with a term that is out of phase by an angle π/2. The coherent coupling can be understood as an additional external drive with a phase and amplitude determined by the coherence in the nearby cavity mode. If we assume that t c is small, then cavities are phase locked to an angle θ j = φ j − π/2. However, as we increase t c , the photon tunneling process induce an effective external driving, with a phase rotated by an angle π/2. This explains qualitatively the rotation of the probability distribution in Fig. 3.
To gain a quantitative understanding, we have explored numerically the PL -NPL phase transition as a function of t c , see Fig. 4. The coherence in each of the nanolasers is calculated with the distribution P θ by using the expression, We find in Fig. 4 (a) that an abrupt transition between a PL and a UL phase happens at t c = 1, at which the coherences seem to become a non-analytic function of t c . The same behaviour is qualitatively observed in the correlation function between cavity modes in Fig. 4 (b). In order to evaluate the finite-size scaling at this critical point we consider a constant coherent coupling t c . We then scale such that we stay at the critical point /r 0 = / √ n 0 = t c . Increasing the number of photons has thus the effect of decreasing the phase decoherence rate which scales as D p ∝ 1/n 0 . Our numerical results show that the coherence at the critical point follows a power-law dependence, and we find the critical exponent β ≈ 0.63 from our numerical calculations (see Fig.  4 (d)). The critical point is thus a non-analytical point in the thermodynamical limit. Finally, by considering the different contributions to the bosonic coherences we have checked the self-rotation of the phase-locking angle induced by the coherent coupling, which is apparent in an increase in the average sin(θ) relative to cos(θ) as we approach the critical point in Fig. 4 (b).

PL -NPL phase transition with inhomogeneous phases
So far we have considered the case of homogeneous driving. An interesting behaviour is found if we study the case of different driving phases φ 1 = φ 2 . Consider for simplicity that we fix φ 1 = 0 and change φ 2 . We plot numerical results in Fig. 5. We compare three cases by keeping the same scaling criteria, namely, we fix /r 0 = 1 and calculate the coherence in the nanolaser system as a function of both t c and the angle φ 2 . We observe the PL/UL phase transition at t c = 1 for values φ 2 = 0, π, in agreement with our results in Fig. 4. However, as we approach the value φ 2 = π/2, we observe that the critical t c required to enter into the UL phase increases to values t c 1. In other words, close to φ 2 − φ 1 = π/2 the PL phase is more robust to a coherent coupling between nanolasers. This effect is more pronounced for large photon numbers and thus lower values of the phase decoherence rate (like in Fig. 5 (b) and (c)). We can understand this effect by, again, looking at the Eqs. (23,24). If we assume small coherent couplings, then the condition φ 1 = 0 implies, according to Eq. (23), that the coherent driving of the first cavity mode induced phase locking into θ 1 = −π/2. Eq. (24) then becomes ȧ 2 = ȧ 2 nl + | |e iφ 2 −π/2 + t c | a 1 |. We find that the coherent driving on the second nanolaser and the photon tunneling term have the same phase, only if φ 2 = π/2. This qualitative argument explains the trend observed in the numerical calculation that the PL phase is more resilient to the coherent coupling term, t c , when ∆φ = π/2. Figure 5.
We plot the value of the average optical coherence in the nanolaser dimer, C = | a 1 | 2 /n 0 + | a 2 | 2 /n 0 /2. This quantity is calculated by using the Fokker-Planck equation, Eq.

Josephson photo-current
By restricting ourselves to the previously mentioned PL steady-states, interesting phenomena emerge as a result of the interplay between on-site interactions and the coherent tunneling term (t c ). In particular, we will show that a photonic analog of the Josephson current is generated between the two nanolasers in the case of a finite difference between the optical phases into which they are locked.
To get a qualitative understanding of this effect, let us examine the evolution of the average photon number, ṅ j , in each nanolaser (j = 1, 2), which can be calculated by means of the Heisenberg equations obtained with the master equation, Eq. (6). We distinguish two contributions. First, we find a term that is given by the non-linear local dynamics of the nanolasers laser, which, in the limit of large photon numbers, reads Note that, apart from corrections arising from the coherent drive, proportional to j , the steady-state solution in the lasing phase is n j = n 0 , given by Eq. (4). Second, we have a contribution arising from the coherent hopping of photons between sites, (analogously for ṅ 2 | tc ). For weak coherent coupling, t c , and small values of j , the lowest order approximation of Eqs. (27,28) can be found by assuming that the effect of the tunneling and coherent drive is negligible in the calculation of the two-point correlation function, a † 1 a 2 . In that case, we find an approximate expression, which fundamentally indicates a Josephson effect by which a phase difference, ∆φ, of the coherent drivings induces a photo-current between sites. This current arises from the coherent nature of the coupling, and it does not appear in the case of a dissipativemediated coupling [40]. Finally, in the limit of perfect phase-locking, there will be a fixed phase relation between the nanolaser and the coherent drive, such that * We can obtain now an approximate expression for the steady average photon number by adding up (27,28), and imposing the conditions ṅ j | nl + ṅ j | t = 0, leading to Hence, we find an imbalanced average photon number given by, The result in Eq. (32) arises from the balance between the individual nanolaser dynamics and the photon tunneling between nanolasers. Note that we have written the solution so as to make explicit the dependence n J ∝ γ/κ, which implies that the photon number imbalance can be increased by increasing the incoherent pumping rate. This indicates that the system can act as a good sensor to estimate the phase difference, ∆φ, by simply measuring the imbalance of the average photon number in the steady state of the system.

Numerical investigation of the Josephson current between nanolasers
So far, our results apply in a limit of strict phase-locking and small coherent couplings. However, we have seen in previous sections that the effect of t c is to destroy the coherence in the coupled nanolasers system. We have to expect that it is not possible to simply increase the Josephson photocurrent by means of increasing the coupling t c , since at some point the system will enter into the NPL, incoherent, steady-state.
To investigate the interplay between Josephson photocurrent and the PL-NPL phase transition we need to carry out numerical calculations that bring us beyond the approximations considered in the previous subsection, in particular beyond the approximation of small t c values. Unfortunately, exact diagonalization of Eq. 1 is a numerically challenging task, since we need to deal with a very large Hilbert space. We resort to two approximate methods: (a) Self-consistent Fokker-Planck equation.
In our first approach we use the inhomogeneous Fokker-Planck equation, Eq. (14), which is valid deep in the lasing regime, and we complement it with an approach to account for the evolution of the different photon numbers in each of the two nanolasers. In particular, we rewrite Eqs. (27,28) in terms of the radial variable of the Fokker-Planck equation, via the identities The condition ṅ j = 0 reads, where the angular averages, sin(θ 1 − θ 2 ) and sin(φ j − θ j ) have to be evaluated relative to some angular probability distribution, P θ (θ 1 , θ 2 ). We proceed by using the following iterative algorithm: (i) We solve Eq. (16) with initial values r 1 = r 2 = √ n 0 . (ii) We use the probability distribution, P θ (θ 1 , θ 2 ), obtained in step (i) to evaluate all the angular terms in Eqs. (34), and we solve those equations to obtain new values r 1 , r 2 . (iii) We solve again the angular equation Eq. (16) with the new values r 1 , r 2 obtained in step (ii) (iv) The last two steps are repeated until convergence is reached.
(b) Gutzwiller ansatz. To obtain results beyond the semiclassical limit described by the Fokker-Planck equation, we use a Gutzwiller ansatz that is valid in the limit of small photon tunneling rates. The Gutzwiller anstatz approximates the steady state of Liouvillian (1) by assuming a separable state of the form ρ ≈ ρ 1 ρ 2 , where each ρ j follows a local Liouvillian, where α j = Tr (ρ j a j ). We seek numerically a self-consistent solution of the set of equations. Whereas this ansatz neglects quantum and classical correlations between nanolasers, it does allow us to include an exact description of the single nanolaser.
To summarize the two methods: the self-consistent Fokker-Planck equation method allows us to include correlations between phases, but it relies on the validity of the semiclassical approximation. The Gutzwiller ansatz, on the other hand, allows us to describe exactly the quantum dynamics at the single nanolaser level, but it does not allow to include correlations. In the lasing regime and with small t c , the two methods must yield the same results. Our numerical calculations are presented in Fig. 6, where we plot the photon number imbalance caused by the Josephson photocurrent. There is a reasonable agreement between the Gutzwiller ansatz and the self-consistent Fokker-Planck approaches. At large values of the coherent coupling t c , the Gutzwiller ansatz does not converge numerically to a steady-state value, and we must assume that neglecting correlations between cavities is not a valid approximation. The self-consistent Fokker-Planck equation, on the other hand, is more robust and converges up to higher values of the coupling term, t c . We see in Fig. 6 that the Josephson effect decreases for large values of t c , as expected from our discussion on the PL-NPL transition in the previous section. Figs. 6 (a) and (b) shows that the photocurrent imbalance grows with t c up to a certain value at which increasing the coupling is detrimental to the coherence between nanolasers, since the system enters into the NPL steady-state.

The photonic Josephson current as a metrological resource.
The result in Eq. (32) indicates that the nanolaser dimer considered in this work is very sensitive to a phase difference φ = φ 2 − φ 1 . If we calculate the derivative, we find that the sensitivity of the system grows linearly with γ/κ, which physically can be interpreted in terms of the bosonic amplification of the Josephson photocurrent. Furthermore, we also find that the sensitivity is increased as we approach the critical point of the lasing phase transition at C p = 1.
We have compared the prediction of the sensitivity around the critical point by using the self-consistent numerical methods introduced in the previous section. We have confirmed that the critical point is indeed the optimal operating point of view for detecting a phase difference, as shown in Fig. 7. This result can be qualitatively explained from the dependence of the photon number dynamics on the typical rate scale A = C p κ, see Eq. 27: close to C p , the local photon number dynamics slows down, thus leading two a stronger effect from the Josephson photocurrent. Finally, we can also estimate the uncertainty in the value of a phase measurement carried out by the coupled nanolaser system in the quantum limit. The latter assumes that experimental error in photon number measurements is entirely due to the quantum fluctuations of the photon number observable. The latter can be estimated to be, where we have used the calculation of photon number fluctuations obtained with the single nanolaser radial Fokker-Planck equation, Eq. (C.11). Eq. (37), together with Eq. (36) allows us to estimate the uncertainty in the estimation of a small phase difference, δφ, from a measurement of the difference in photon numbers in the coupled nanolasers, The last expression shows that the error grows with C p , and it scales like δφ ∝ 1/ γ/κ ≈ 1/ √ n. This expression shows that the amplification effect survives in the quantum limit, in which the accuracy can be enhanced by increasing the number of photons, for example by increasing the excitation rate of the qubits, γ.

Conclusions & Outlook
We have presented a theoretical study of two nanolasers coupled by a photon tunneling term. We have arrived to two main conclusions. The first is that a photon tunneling term is a source of decoherence which can ultimately destroy the phase locking of each individual nanolaser to an external driving field. The second conclusion is that, in a limit of small photon tunneling rates, a photonic Josephson effect is induced that can be used to measure the phase difference between external fields.
The model of interferometer proposed in this article may be implemented in circuit QED systems. In this experimental platform, single-qubit nanolasers have already been demonstrated [24], and schemes for controlling the properties of single-qubit laser light, including the generation of entangled states of light, have been proposed [31]. Photon tunneling is actually the main mechanism that couples microwave cavities in circuit QED [30], thus making this system an ideal experimental platform for implementing lattices of coupled nanolasers. Our ideas can also be translated to the realm of phononics, for example in trapped ion setups. Here, the coupling between internal electronic states acting as quibts and the vibrational degrees of freedom can be controlled with optical lasers. Single ion phonon lasing has actually been experimentally demonstrated [28]. Vibrational modes of coupled trapped ions can often be described in terms of local phononic modes, with the Coulomb interaction between ions inducing a phonon hopping process [33], an effect that has been experimentally observed [42,43]. Trapped ions actually offer the possibility to study few-mode coupled lasers or include highly controllable qubit-phonon interactions [44] .
Our work could be extended to larger lattice sizes, something that could potentially allow high sensitive estimation of phase gradients. Additionally, the simultaneous implementation of coherent couplings and dissipative-mediated couplings [40] is expected to exhibit interesting dissipative phase transitions of photonic phases, also interesting for studying quantum synchronization. An intriguing possiblity here is the investigation of non-reciprocal couplings and topological effects, as well as topological amplification [45]. Finally, these results invite to study further Josepshon effects or configurations, like the a.c. Josepshon effect or the SQUID [46].

Appendix A. Adiabatic elimination
In this appendix we show how the adiabatic elimination of the qubits leads to Eq. (6). For this we employ a generalization of the procedure discussed in Refs. [37,40]. Since we assume that photon tunneling rates smaller than local energy scales, we consider the single nanolaser case for this derivation.
Firstly, we trace over the qubits from the master equation. Since the effect of the photon decay, κ, is not affected by the adiabatic elimination procedure, we set κ = 0, and reintroduce it later in the calculation, where ρ f is the reduced density matrix of the field after tracing over the qubit states, ρ f = Tr q {ρ} = g|ρ|g + e|ρ|e , and we employ the notation ρ ge = g|ρ|e . To obtain a closed equation for ρ f we need to eliminate the operators ρ ge , ρ eg from Eq. (A.1). The corresponding equations of motion for these operators are, The operators ρ ge and ρ eg may now be adiabatically eliminated (in the limit γ κ, g, | |) from Eq. (A.1) by takingρ ge ≈ 0 in Eq. (A.2) and substituting their steady-state solutions, Now we can use Eq. (A.3) and insert it into Eq. (A.1) and get, We need now equations of motion for ρ gg and ρ ee , which can be derived again from Eq. (1), leading tȯ We may now reach a perturbative solution to the steady-states of Eqs. (A.5, A.6) in terms of the field density matrix ρ f by adiabatically eliminating ρ gg andρ gg . This is done by takingρ gg ≈ 0 in Eq. (A.6), The ground state population of each qubit is expected to be negligible as a result of the fast pumping of the qubits. Thus, in first order we can assume ρ gg ≈ 0 and ρ ee = ρ f − ρ gg ≈ ρ f . A second order approximation is achieved by inserting this first order approximation into Eq. (A.7), hence The desired closed equation for ρ f is accomplished by plugging Eqs. (A.8, A.9) into Eq. (A.4), leading tȯ We can check, by using the commutation relations of a that the equation can be written in Lindbladt form, with definitions Ag 2 /γ, B2g 4 /γ 3 . If we include now the photon decay term, we finally get to Eq. (6).

Appendix B. Renormalization of the adiabatic equation
Eq. (A.10) is a perturbative description in the limit of adiabatic elimination of the fast qubit dynamics. This equation is strictly valid below the critical point, C p < 1, and slightly above it, C p 1 [37], which implies a severe limitation in the range of applicability of our results.
For the single qubit laser the average number of bosons in the steady state predicted by Eq. (6) is, which tends to zero for large values of the pumping parameter. Eq. (B.1) can be obtained by finding the steady-state value of the photon number operator and neglecting photon number fluctuations. This prediction differs from the mean field result [37], The two expressions agree only close to C p = 1. Still, one may think of a renormalization procedure such that allows us to perform a summation of the neglected terms in the perturbative series with respect to (g/γ) 2 . The root of the perturbative nature of Eq. (A.10) is given by the truncation performed in Eq. (A.8). We search for a new adequate parameter β that takes into account the remaining terms of the series such that, We discuss now a procedure to compute β by using an exact relation held in the steady state which is easily inferred from the Eq. (1) of the Letter. Let us calculate the Heisenberg equation of motion for the observable N = a † a + σ z /2, which reads, for the single site case, Figure B1.
In the steady state, d N dt = 0, which implies where n = a † a. Notice that this is an exact relation. On the other hand, by taking traces in both sides of (B.3), we find the following equation We equate now Eqs (B.5) and (B.6), and take the limit n 1, and get, This eventually leads to a renormalization of the parameter B that amounts to B → B/C p ≡ B r . The latter leads to a correction of Eq. (B.1), such that so that the number of photons finally actually agrees with the mean result. Furthermore, in Fig. B1 (a) we present numerical results that show that the photon number calculated with Eq. (6) agrees very well with an exact single qubit laser calculation using Eq. (1). In Fig. B1 (b) we show calculations carried out with the renormalized B r parameter, that show that the number of photons curves converge to the mean-field result for large values of γ/κ. Note that for smaller values of n it is possible to make an n-dependent definition of β. One could then carry out an iteration procedure (e.g. estimate n with β = 1/C p , then use this value to re-calculate β, calculate n again, and so on until reaching convergence). The added benefit of this is, however, not very high, since our renormalization procedure works very well in practice.

Appendix C. Fokker-Planck equation
Appendix C.1. Derivation of the Fokker-Planck equation Let us sum up the derivation of the Fokker-Planck equations used in this work. Recall that we employ the Glauber-Sudarshan P representation of the density matrix [36], defined by where |α is the coherent state |α = exp (αa † − α * a)|0 . We transform the master equation (A.10) with the help of the following relations for coherent states, a|α α| = α|α α|, With Eqs. (C.2) we can write the master equation as an equation of motion for P (α, α * , t), after an integration by parts with the assumption of zero boundary conditions at infinity. This procedure is simplified in the limit |α| 2 1, and g γ in which we drop any contribution smaller than B|α| 2 α, as B happens to be a very small coefficient compared to A, B/A ∝ (g/γ) 2 1. In doing so, we arrive at the Fokker-Planck equation claimed in the Letter, Our equation of motion can be written in polar coordinates with the aid of the equivalences, Eq. (C.3) then reads, This is a very complicated equation, whose analytical or even numerical solution if very challenging. We are going to see below how a simplification can be justified in the lasing regime.

Appendix C.2. Single nanolaser Fokker-Planck equation
Before addressing the case of two coupled nanolasers, let us consider the solution of a single nanolaser Fokker-Planck equation without any external drive ( = 0). In this case, P (r, θ) = R(r)P θ (θ), with P θ = 1/(2π). The radial function, R(r), satisfies the differential equation, The steady-state solution of Eq. (C.6) with dR/dt = 0 is, where N = ∞ 0 drrR(r), is a normalization constant. In the lasing regime (C p > 1) and in the limit of large photon numbers (γ/κ > 1), the steady-state radial probability function takes the form, The radial distribution can be used to calculate the mean and variance of the photon number, n = (r 0 ) 2 , (C.10) ∆ 2 n = n 2 − n 2 ≈ 1 2 γ κ . (C.11) Note that Eq. (C.10) agrees with the mean-field calculation (n mf = (r 0 ) 2 ). Also, Eq. (C.11) agrees with the calculation of the variance from a coherent state in the limit of C p 1, for which we recover the well-known relation ∆ 2 n = n [47].

Appendix C.3. Reduced Fokker-Planck equation for two nanolasers
In the lasing regime, we can simplify Eq. (C.3) by assuming that the radial variables are settled around their steady-state mean values r j ≈ r 0 j , and P (α, α ) can be well approximated by a factorized form, P (r, θ) = R(r 1 )R(r 2 )P θ (θ 1 , θ 2 ), (C.12) where each R(r j ) is a properly normalized Gaussian function around r 0 j , Eq. (C.12) is valid as long as the phase variable evolves on much slower time scales as the radial variables. The typical times scales for phase τ φ and radial τ r variables can be estimated, in the lasing regime and assuming values of r j close to r 0 j , from Eq. (C.5), 1/τ θ = A 2 n = κC p , (C.14) Additionally, the relaxation time associated to the photon number variable can also be estimated by means of the Heisenberg equations, see Eq. (27). We find, thus, that τ r ≈ nτ φ , such that in the large photon number limit, the nanolaser radial coordinate can be considered to relax very fast compared with the phase variable. This behaviour is consistent with other theoretical descriptions of the laser's dynamics [47].
If we now integrate both sides of Eq. (C.5) in the radial variables ∞ 0 rd r, the result is an equation of motion for the angular part. In this derivation it is important to notice that the terms ∂R j /∂t can be cancelled together with other terms by using the identity Eq. (C.6). Also, the integration of the first derivative ∂ r R is eliminated through the relation, After grouping terms, the resulting equation adopts the form claimed in equation (6)  Notice that this a sort of quantum version of the stochastic Kuramoto model, that is, it represents two Fokker-Planck equation coupled by a term that was originated by coherent photon tunneling.