The temporal coherence of a photon condensate: A quantum trajectory description

TL;DR

This paper models the temporal coherence of a photon condensate using a quantum trajectory approach, analyzing effects of interactions and dissipation on coherence properties and validating a laser physics-inspired model.

Contribution

It introduces a quantum-mechanical model for photon condensate dynamics, including effects of Kerr interactions and dissipation, and validates a laser physics-based heuristic in this context.

Findings

Observation of long-$\tau$ antibunching effect.

Identification of exponential and Gaussian decay regimes of phase correlations.

The ratio of coherence times inversely relates to number fluctuations.

Abstract

In order to study the temporal coherence of a single-mode dye-cavity photon condensate, a model is developed for the dynamics which treats the condensate mode on a quantum-mechanical level. The effects of driving-dissipation and Kerr interactions on the number fluctuations are studied analytically and numerically, including the finding of a long-$\tau$ antibunching effect. Depending on the interaction strength, we quantitatively observe an exponential Schawlow-Townes-like decay or Gaussian Henry-like decay of phase correlations. The adequacy of a heuristic phasor model originating from laser physics in describing number and phase dynamics is validated within the experimentally relevant parameter regime. The ratio of the first and second order coherence times is shown to be inversely proportional to the number fluctuations, with a prefactor that varies smoothly throughout the crossover between canonical and grandcanonical statistics.

Figures

Tables

Table 1. Table 1: Experimental parameters corresponding to Rhodamin 6G at 560 nm. Schmitt ( 2018b ) . Further, we have taken for the simulations n ¯ = 1000 ¯ 𝑛 1000 \overline{n}=1000 , which somewhat smaller than in typical experimental setups to reduce the relevant timescales. We have chosen the finite value of interaction strength such that the phase and number statistics are significantly altered ( σ = U ​ n ¯ 2 T 𝜎 𝑈 superscript ¯ 𝑛 2 𝑇 \sigma=\frac{U\overline{n}^{2}}{T} of order one), while the average interaction energy per particle U ​ n ¯ 𝑈 ¯ 𝑛 U\overline{n} remains small compared to other energy scales. The notion of ‘without interactions’ refers to U = 0 𝑈 0 U=0 , whereas ‘without losses’ means κ = 0 𝜅 0 \kappa=0

Parameter	value
$B_{12}$	2.5 kHz
$M_{tot}$	${10}^9$
$\Delta$	$-2.4kT/\hslash$
$T$	$300K$
$U$	${10}^{-5}{B}_{12}{M}_{\text{tot}}=6.4\times {10}^{-7}{k}_{B}T/\hslash$
$w_{\uparrow }/w_{\downarrow }$	1
$\kappa$	$2.2GHz=8.3\times {10}^{-4}{B}_{12}{M}_{\text{tot}}$

Equations

$$dn=-dN+dM.$$

$$dn=-dN+dM-dK.$$

$$d{M}_{\uparrow }=dN-dM+dP$$

$$\text{partialderivative}nt=B_{21}(X-n)(1+n)-B_{12}n(M-X+n)-\kappa n,$$

$$X(\overline{n})=\frac{\overline{B_{12}n(M+n)+B_{21}n(n+1)+\kappa n}}{\overline{B_{12}n+B_{21}(n+1)}}.$$

$$\frac{B_{21}(\omega )}{B_{12}(\omega )}=\frac{w_{\downarrow }}{w_{\uparrow }}{e}^{-\beta (\omega -{\omega }_0)}.$$

$$\frac{B_{21}(n)}{B_{12}}=\frac{w_{\downarrow }}{w_{\uparrow }}{e}^{-\beta (\Delta +Un)}.$$

$$\frac{{\mathcal{P}}_{n+1}}{{\mathcal{P}}_n}=\frac{X-n}{M-X+n+1}\frac{B_{21}(n)}{B_{12}}$$

$$\frac{{\mathcal{P}}_n}{{\mathcal{P}}_0}=\frac{\left(\genfrac{}{}{0ex}{}{X}{n}\right)}{\left(\genfrac{}{}{0ex}{}{M_{\text{tot}}-X+n}{n}\right)}{e}^{-\beta \left[\Delta n+(U/2)n^2\right]},$$

$${\mathcal{P}}_n\propto {\left(\frac{\overline{M_{\uparrow }}}{\overline{M_{\downarrow }}}\right)}^n{e}^{-\beta E[n]}=\exp \left(-\frac{(\Delta -\mu )n+(U/2)n^2}{T}\right),$$

$$M_{\text{eff}}=\frac{M_{\text{tot}}}{2}{\left[1+\cosh \left(\beta (\Delta +U\overline{n})\right)\right]}^{-1},$$

$$\text{partialderivative}t$$

$$-\left[\frac{\overline{B_{21}}X}{\overline{n}}+(B_{12}+\overline{B_{21}})\overline{n}+\beta U\overline{B_{21}}{M}_{\uparrow }(1+\overline{n})\right]\delta n(t)$$

$$+\mathcal{O}(\delta n(t{)}^2),$$

$${\Gamma }_2$$

$$=\left[1+\sigma +\frac{{\overline{n}}^2}{M_{\text{eff}}}\right]\frac{B_{12}\overline{M_{\downarrow }}}{\overline{n}},$$

$$d\left[\delta n(t{)}^2\right]$$

$$=-2{\Gamma }_2\delta n(t{)}^{2}dt+R(\overline{n}+\delta n+1)dt+\gamma (\overline{n}+\delta n),$$

$$g^{(2)}(0)-1={\eta }^2=\frac{\overline{\delta n^2}}{{\overline{n}}^2}=\frac{1}{1+\sigma +\frac{{\overline{n}}^2}{M_{\text{eff}}}}.$$

$$M_{\uparrow }=X(\overline{n})-\overline{n}=\frac{(B_{12}M+\kappa )\overline{n}+(B_{12}+B_{21}){\overline{n}}^2{\eta }^2}{B_{21}(\overline{n}+1-\beta U{\overline{n}}^2{\eta }^2)+B_{12}\overline{n}}.$$

$$\frac{\overline{n{n}^{\prime }}}{{\overline{n}}^2}-\frac{\overline{nn}}{{\overline{n}}^2}\approx \frac{{\overline{n}}^2(B_{12}+B_{21}{)}^2}{B_{12}{B}_{21}{M}_{tot}}\kappa {\tau }_c^{(2)},$$

$$E(n)=\Delta n+\frac{U}{2}{n}^2,$$

$$v_p=\text{derivative}En=\Delta +Un.$$

$$dn=(R-\gamma -\kappa )ndt,$$

$$d\theta =\frac{1}{\sqrt{n}}\sin (\varphi )d{M}_{\text{spont}}.$$

$$d\mathrm{Var}[\theta ]=\overline{\frac{R}{2n}}dt.$$

$$dn=(R-\gamma -\kappa )ndt+(1+2\sqrt{n}\cos \varphi )d{M}_s,$$

$$d\mathrm{Mean}[n]=d\overline{n}=-(\gamma +\kappa )\overline{n}dt+R(\overline{n}+1)dt,$$

$$d\mathrm{Var}[n]=2(R-\gamma -\kappa )\mathrm{Var}[n]dt+R(2\overline{n}+1)dt.$$

$$\hat{\mathrm{H}}=\Delta {\hat{\mathrm{a}}}^{†}\hat{\mathrm{a}}+\frac{U}{2}{\hat{\mathrm{a}}}^{†}{\hat{\mathrm{a}}}^{†}\hat{\mathrm{a}}\hat{\mathrm{a}},$$

$$|\tilde{\psi }⟩=$$

$$+\sqrt{\kappa }\hat{\mathrm{a}}d{Z}^{\ast }+\left(\frac{\hat{\mathrm{a}}}{\text{expectationvalue}\hat{\mathrm{a}}}-1\right)dN+\left(\frac{{\hat{\mathrm{a}}}^{†}}{{\text{expectationvalue}\hat{\mathrm{a}}}^{\ast }}-1\right)dM]|\psi ⟩.$$

$${\text{expectationvalue}*{\hat{\delta }}_{\theta }{\hat{\delta }}_{\theta }}_{\text{stat}}=\frac{1+\sqrt{2(1+\frac{\gamma +R}{\kappa })}}{4\text{expectationvalue}\hat{\mathrm{n}}}.$$

Full text

††preprint: APS/123-QED

I Introduction

The last few decades, a vast amount of research has been done to the related phenomena of Bose-Einstein Condensates (BECs) and quantum fluids Pitaevski and Stringari (2016); Bramati and Modugno (2015), for reasons of fundamental interest as well as their exploitation for quantum simulation.

The traditional platform for the realization of quantum fluids are ultracold atoms, but another fascinating platform is a BEC of photons. Because their effective mass is many orders of magnitude lower, condensation can take place at room temperature. In order to achieve thermalization, a successful approach has been to use effective photon-photon interactions in a nonlinear material, resulting in the condensation of polaritons: coherent superpositions of photons and excitons Carusotto and Ciuti (2013), of which the temporal coherence has been studied in Whittaker and Eastham (2009); *Whittakerpolaritons2; Haug et al. (2012).

An alternative approach, where there is no need for the photons to be hybridized, consists of letting the photons interact with their environment Klaers et al. (2010a, b); Rajan et al. (2016); Nyman and Walker (2018); Kirton and Keeling (2015); Greveling et al. (2018); Hesten et al. (2018). This is realized experimentally by frequent absorptions and emissions of the photons by dye-molecules (the gain medium), that are themselves thermalized by collisions with solvent molecules. It has been shown  Kirton and Keeling (2013); Schmitt et al. (2015) that there is a smooth crossover between laser physics and photon BEC. The distinction is that in the latter, almost full thermalization is reached with minor leakage, whereas in a laser, which is driven further out-of-equilibrium, gain and losses dominate the dynamics. A consequence is that lasers generally need population inversion for the gain medium, while a photon BEC does not.

An interesting, presently realized, feature of photon BECs in a dye microcavity system, as opposed to other BECs, is that the number of photons does not remain fixed through time, even in the absence of losses. If the gain medium is sufficiently large such that saturation effects are small, the condensate will exhibit grand-canonical statistics Klaers et al. (2012); Schmitt et al. (2014); Weiss and Tempere (2016); Proukakis et al. (2017). This is remarkable, because the non-interacting Bose gas is one of the few systems where ensemble equivalence is not satisfied Huang (2014).

Along with these statistical fluctuations of the particle number, one observes fluctuations of the phase Schmitt et al. (2016); Schmitt (2015); *thesisJulianEng; Marelic et al. (2016), as first predicted by de Leeuw et al. (2014). These can be understood from a heuristic phasor model (HPM), originally proposed in laser physics Henry (1982); Scully and Zubairy (1997), because the phase fluctuations at sufficient particle number correspond to the standard Schawlow-Townes broadening Scully and Zubairy (1997). The dynamics seems richer though, because in the grand-canonical regime the particle number can become zero so that the phase of a subsequent new photon is entirely random, known as a ‘phase jump’. Although predictions from the HPM seem to match experiment so far, the theoretical understanding is still limited.

In particular, the phase jump picture suggests that a dramatic increase of temporal coherence occurs when the second order coherence falls below roughly $g^{(2)}(0)=1.5$, where the probability of having zero photons in the condensate starts to become exponentially small Schmitt et al. (2016). We will show here that throughout the crossover between canonical and grandcanonical statistics, the behaviour of the phase coherence time behaves much smoother than expected by the phase jump analysis.

We also study the influence of a weak Kerr nonlinearity on the coherence. For ultracold atoms, the temporal coherence properties in the presence of interactions have been studied in Sinatra et al. (2009) for bosons and Kurkjian et al. (2016) for fermions. Nonlinearities are not necessary for condensation to occur and are thus often disregarded Schmitt (2015); *thesisJulianEng. Nevertheless, they naturally emerge in experimental setups, where a not entirely decoherent dye induces a natural Kerr effect Radonjić et al. (2018). Furthermore, it is possible to engineer these interactions on purpose, which makes the photons behave more like polaritons and might open the possibility for effects such as superfluidity Berman et al. (2017). It has also been proposed to realize similar behavior with $\chi^{(2)}$-nonlinear materials Majumdar and Gerace (2013). Although the precise intrinsic value of the Kerr nonlinearity in the commonly used dye molecules is still subject of debate (a recent discussion is given in Radonjić et al. (2018)), we show that already small values of the Kerr nonlinearity can subsantially affect the coherence properties. Because the interaction strength can vary many orders of magnitude, we here focus on the qualitative influence on interaction scales where the number and phase dynamics are altered significantly, while the associated energy scale per particle remains small compared to other energy scales. Apart from the instantaneous Kerr effect, delayed interaction may also arise in photonic condensates because of thermo-optical effects Alaeian et al. (2017). The latter are out of the scope here because they are expected to contribute significantly only on timescales much longer than our simulations, although in principle our model could be extended to include them.

Throughout this work, we will assume for simplicity the presence of only a single photon mode. In many experimental setups, this is valid because of the lower population of excited modes, the fact that these do not interact and their small overlap with the condensate Klaers et al. (2011). Analytical predictions considering only a single mode have also been successful in describing in phenomena as the decay of second-order coherence Schmitt (2015); *thesisJulianEng. As confinement is improving with novel experiments, restriction to single mode becomes more realistic in practice. Nevertheless, beyond this single-mode description, photonic condensates in some parameter regimes can also exhibit effects arising from multiple modes that lead to interesting physics: if the occupation of other modes is large and the finite size of the reservoir is important, mode-competition can become important Hesten et al. (2018). Large occupation of other modes in a substantially interacting system, will additionally result in Beliaev-Landau scattering which may give corrections to the single-mode results Pitaevski and Stringari (2016).

In the next section, Sec. II, we describe the photon condensate system in more detail, and study the number statistics. We revisit standard results with inclusion of a Kerr nonlinearity and present thermodynamic and dynamical estimates for the number fluctuations, which are verified by solving stochastic rate equations. We also show that weak interactions can alter significantly the density distribution and intensity correlations; and further that finite losses can cause long-$\tau$ antibunching.

In Sec III, we tackle the semiclassical HPM in presence of Kerr interactions, and point out some conceptual issues stemming from its the heuristic nature.

To verify the validity of the HPM, we provide in Sec. IV a model where the condensate mode is treated on a fully quantum-mechanical level, while the dye molecules are described by classical rate equations. As we aim to describe single-shot experimental realizations, this naturally leads to the quantum trajectory (wave-function Monte Carlo) formalism yielding a stochastic Schrödinger equation that describes an open system under continuous measurement Breuer and Petruccioni (2002); Carmichael (2008). Because photon numbers in the cavity can easily reach order $10^{4}$, it is numerically very demanding to solve this Schrödinger equation exactly on a truncated Fock space. Therefore, at high densities, we make the variational ansatz that the state is Gaussian in density and phase, which has shown to be good in describing dephasing properties Verstraelen and Wouters (2018). In addition to reducing computational complexity, the variational equations also yield theoretical insight in the density-phase dynamics. Quantum trajectories have previously also been used for the study of the related polariton condensates Wouters (2012).

Using this quantum trajectory model, we show that, despite its heuristic character, the HPM is valid for experimentally relevant quantities as their predictions from both methods match. Finally, In Sec. V we study the properties of first order coherence in photon condensates, including the influence of Kerr interactions and reservoir size, and relate the first and second order coherences to each other . Both Schawlow-Townes and Henry phase diffusion mechanisms are observed, and we show that the presence of phase jumps only a limited quantitative effect on the phase decay. We formulate our conclusions in Sec. VI.

II Number statistics

II.1 Model for a driven-dissipative, interacting photon condensate

We assume a reservoir containing $M_{\text{tot}}$ two-level dye molecules, of which $M_{\uparrow}(M_{\downarrow})$ are in the excited state (ground state). Because of the rapid decoherence due to collisions with solvent molecules Kirton and Keeling (2013), these numbers can be treated as classical integers. We consider this reservoir to be coupled with a single, quantum mechanical, photon mode with number operator $\operatorname{\hat{n}}$. Because we restrict ourselves in this section to number statistics, it is sufficient to consider only number eigenstates (Fock states) with $n$ photons. The corresponding rate equation for the photons is then given by

[TABLE]

Here $dN$ and $dM$ are Poisson processes describing molecular absorption at rate $B_{12}M_{\downarrow}n$ and molecular emission at rate ${B_{21}M_{\uparrow}(n+1)}$, respectively, where $B_{12}$ and $B_{21}$ are the modified Einstein coefficients for absorption and emission Klaers et al. (2012). For convenience, we define $\gamma:=B_{12}M_{\downarrow}$ and $R:=B_{21}M_{\uparrow}$. In the absence of losses, the total number of excitations $X=M_{\uparrow}+n$ is conserved such that $dM_{\uparrow}=-dn$ and $dM_{\downarrow}=dn$. For this system, many previous analytical results regarding the number distribution and correlation functions have been summarized in Schmitt (2015); *thesisJulianEng. In practice, the quality factors of cavities are restricted so that an additional loss at rate $\kappa n(t)$ Klaers et al. (2010a) takes place, which is compensated by an additional pumping of molecules. In a typical experimental setup, $\gamma\gg\kappa$, so that we expect losses not to affect the number statistics significantly. As we are concerned with ensemble statistics only, we have the freedom to choose the photon-counting unraveling for external losses Gardiner and Zoller (2004), such that they are also modelled by a Poisson process $dK$. This results in a rate equation

[TABLE]

In order to keep the long-time expectation value of $X$ constant, the molecule reservoir is pumped through a process $dP$ with constant rate $\kappa\overline{n}$, where $\overline{n}$ is the average particle number such that

[TABLE]

and $M_{\downarrow}=M_{tot}-M_{\uparrow}$. From equation (2), the average evolution is Schmitt (2015); *thesisJulianEng

[TABLE]

from which one finds that in the steady state

[TABLE]

In general, the Einstein-coefficients $B_{12},B_{21}$ of the absorption and emission processes are related by the Kennard-Stepanov law (see the supplementary material of Klaers et al. (2012)):

[TABLE]

Here, $\omega$ is the frequency of the photon, $\omega_{0}$ the frequency of the atomic transition between ground and excited state of the dye, $\beta=\frac{1}{T}$ the inverse temperature (we set $k_{B}=\hbar=1$ throughout this work) and $w_{\uparrow,\downarrow}$ are weight coefficients taking into account the internal molecular density of (rovibrational) states in the electronic ground and excited level. For noninteracting photons with frequency $\omega_{c}$, the frequency dependence in (6) depends only on the detuning $\omega_{c}-\omega_{0}=\Delta$.

A Kerr effect adds by definition an interaction energy $E_{\text{int }}(n)=\frac{U}{2}n^{2}$ to the photons. The $n$th photon now carries a frequency $\omega=\omega_{c}+[E_{\text{int}}(n)-E_{\text{int}}(n-1)]\approx\omega_{c}+Un$. As long as the interaction strength and particle number are not excessively large, only the relative difference between $B_{21}(\omega)$ and $B_{12}(\omega)$ is important, so that we can treat $B_{12}$ as constant. With this, we rewrite (6) as

[TABLE]

II.2 Number statistics

As in Ref. Klaers et al. (2012) the steady state number distribution can be found by assuming detailed balance

[TABLE]

such that

[TABLE]

where the notation between brackets refers to binomial coefficients. In the limit of an infinite reservoir, (9) reduces to a Bose-Einstein distribution

[TABLE]

where the chemical potential is $\mu=T\log{\left(\overline{M_{\uparrow}}/\overline{M_{\downarrow}}\right)}$.

In Fig. 1, we see that the predictions in number distributions (9), closely match numerical results from rate equations (2), for a set of parameters corresponding to recent experiments Schmitt (2018b), as summarized in Tab. 1. As for the interaction strength, we have taken the value $U=10^{-5}\hbar M_{\text{tot}}B_{12}$, corresponding to a dimensionless interaction parameter $\tilde{g}=\frac{U}{\Omega\hbar}\approx 9.95\times 10^{-5}$ for trapping frequency $\Omega\approx 8\pi\times 10^{10}~{}Hz$, which is larger than the most common estimates of the natural Kerr effect Klaers et al. (2011); Radonjić et al. (2018). Importantly, by defining the effective reservoir size

[TABLE]

one can distinguish in the noninteracting case a canonical regime with Poissonian number statistics ($\overline{n}^{2}\gg M_{\text{eff}}$), a grandcanonical regime with Bose-Einstein statistics ($\overline{n}^{2}\ll M_{\text{eff}}$) and a transition region ($\overline{n}^{2}\approx M_{\text{eff}}$).

For our parameters in Tab. 1, we have an effective reservoir size $M_{\text{eff}}=7.6\times 10^{7}\gg\overline{n}^{2}=10^{6}$ such that the system is rather on the grandcanonical side of this crossover. Definition (11) only weakly depends on $U$ to the extent that it changes the average energy per photon. However, also for fixed $M_{\text{eff}}$, the number distribution is significantly altered by interactions: the number fluctuations are reduced, as predicted by Weiss and Tempere (2016). This can be understood thermodynamically: by approximating $n$ to be continuous and integrating (10) over positive $n$, we obtain the grand-canonical partition function $\mathcal{Z}$. From the associated free energy $F=-T\log\mathcal{Z}$, the equation of state $\overline{n}(\mu)$ can be obtained, as wel as the amount of number fluctuations $\operatorname{Var}{[n]}(\mu)$. Eliminating $\mu$, the ratio $\eta=\frac{\sqrt{\overline{n^{2}}-\overline{n}^{2}}}{\overline{n}}$ can be expressed as a function of the interaction parameter $\sigma=\frac{U\overline{n}^{2}}{T}$, as shown on Fig. 2. Asymptotically, The relative amount of fluctuations $\eta$ decreases as $1-\sigma$ for small $\sigma$ and as $\sigma^{-1/2}$ for large interactions $\sigma$. For our parameters, we find $\sigma=0.64$, corresponding to $\eta=0.75$, whereas our noninteracting condensate has $\eta=0.99$.

For a finite reservoir, $\mu$ becomes dependent on $n$, so that the above thermodynamic analysis is no longer valid. $\eta$ can also be computed self-consistently from the number distribution (9), or analytically from the dynamical argument below.

II.3 Second-order coherence

The second order coherence time $\tau_{c}^{(2)}$, the decay time of $g^{(2)}(\tau)=\frac{\overline{n(t)n(t+\tau)}}{\overline{n}^{2}}$, can be obtained by an extension of the approach in Schmitt (2015); *thesisJulianEng. In linear approximation, $B_{21}=\overline{B_{21}}(1-\beta U\delta n(t))$, such that the average fluctuation evolves as

[TABLE]

where the driven-dissipative nature has been disregarded. Accordingly, the number correlations decay at rate

[TABLE]

As long as $U\overline{n}\ll\Delta$, the change in $\Gamma_{2}=(\tau_{c}^{(2)})^{-1}$ from a noninteracting condensate is directly proportional to the interaction parameter $\sigma$. Note that the same expression (13) is also obtained if $B_{21}$ is treated as constant while considering frequency dependence in $B_{12}$. We verify the validity of decay rate (13) by comparison to numerical simulations by rate equations in Fig. 3.

For generic parameter values, $g^{(2)}(0)=1+\eta^{2}$ can also be estimated by the stationary solutions of the expectation value of

[TABLE]

from which

[TABLE]

For an infinite reservoir ($M_{\text{eff}}/\overline{n}^{2}\rightarrow\infty$), we can compare this dynamical result directly with the thermodynamic result above. As we see in Fig. 2, the relative difference in predictions for $\eta$ remains smaller than 5% for all values of $\sigma$, and (17) becomes exact in the limit of large $\sigma$. We can attribute the deviations at smaller $\sigma$ to higher order contributions in $d\delta n$.

Such an estimate of $\eta$ is also useful to initiate the number of excited molecules in the numerical simulations. As can be seen from Eq. (5), both the average number of photons and its fluctuations determine the number of excited molecules:

[TABLE]

In previous analytic discussions, we have disregarded the driven-dissipative nature of the system. Intuitively, one might expect that losses will reduce the second order coherence time. Such an effect is not observed in the number distribution, as seen in Fig. 1.

However, as we see in Fig.3, $g^{(2)}(\tau)$ is altered on timescales long compared to the photon correlation time $\tau_{c}^{(2)}$, where a long-$\tau$ antibunching effect appears. This can be explained in the following way: if at $\tau=0$ the stochastic particle number is $n=\overline{n}+\delta n$, the corresponding fluctuation modifies the excitation number by $\Delta X=-\kappa\,\delta n\tau_{c}^{(2)}$. Because $\overline{n}=\overline{n}(X)$, the particle number expectation value is altered to $\overline{n}^{\prime}=\overline{n}+\partialderivative{\overline{n}}{X}\Delta X$ where from (5), $\partialderivative{X}{\overline{n}}\approx\frac{B_{12}B_{21}M_{tot}}{\overline{n}^{2}(B_{12}+B_{21})^{2}}$, so that $g^{(2)}(\tau)$ changes in the order of

[TABLE]

where we used that in the grandcanonical regime, $\delta n\sim\overline{n}$. Expression (19) predicts a decrease in $g^{(2)}(\tau)$ of about $0.012$. We verify this numerically by averaging $g^{(2)}(\tau)$ over a time-interval $10^{4}B_{12}^{-1}M_{\text{tot}}^{-1}-5\times 10^{4}B_{12}^{-1}M_{\text{tot}}^{-1}$, which is sufficiently larger than $\tau_{c}^{(2)}$ but smaller than the timescale of fluctuations in total excitation number $\tau_{X}$ and obtain a value $g^{(2)}(\tau_{c}^{(2)}<\tau<\tau_{X})-1=-0.012$ for the experimental $\kappa$ and $g^{(2)}(\tau_{c}^{(2)}<\tau<\tau_{X})-1=0.003$ without losses, in agreement with our results. Of course, one always has $g^{(2)}(\infty)=1$ as initial and final state become entirely independent, but the relaxation from deviation (19) only takes place over the timescale of the dynamics of $X$, namely $\tau_{X}=(\kappa\partialderivative{\overline{n}}{X})^{-1}$, which is of order $10^{5}~{}B_{12}^{-1}M_{tot}^{-1}$ in our simulations. Because of the different values for $\tau_{c}^{(2)}$ and $\eta$, this effect is weaker in our simulations of the interacting condensate.

III The Heuristic Phasor Model for phase evolution

III.1 A semiclassical model

Next, we proceed to the evolution of the phase. Some descriptions of the temporal coherence of a photon condensate rely on a heuristic phasor model (HPM) Schmitt (2015); *thesisJulianEng, originally developed in laser physics Henry (1982); Scully and Zubairy (1997), that we will repeat here with the addition of Kerr interactions. According to the HPM, one considers the field to be classical such that the state is defined by a single phasor (corresponding to a coherent state). The energy of a state with n photons, relative to a situation where all excitations are in the dye, is given by

[TABLE]

where $\Delta$ is the detuning between the cavity and dye transition. As we treat the phasor as an order parameter, the phase oscillates at speed Pitaevski and Stringari (2016)

[TABLE]

As statistical properties remain the same in a rotating frame, we may replace $\Delta$ in numerical simulations by $\Delta^{\prime}=0$ (no rotation for the vacuum) or $\Delta^{\prime}=-U\overline{n}$ (no rotation on average) for convergence. Absorption, stimulated emission and external losses are here treated as deterministic currents that retain the coherence entirely: they are modeled by an evolution of the particle number

[TABLE]

where $Rn$ is the rate of stimulated emission, $\gamma n$ the rate of absorption and $\kappa n$ the rate of external losses through the mirrors. Spontaneous emissions into the condensate mode are then taken into account as an additional Poisson process $dM_{s}$ with expectation $\overline{dM_{s}}=Rdt$. For each spontaneous emission, a vector $e^{i\phi}$ with unit magnitude and random phase is added to the field $\sqrt{n}e^{i\theta}$, with the motivation that this represents an additional photon that is fully incoherent Abraham and Smith (1977). In fact , such a spontaneous emission does not deterministically change the photon number with one and may even reduce it. As before, $R$ and $\gamma$ depend on $U\overline{n}$ and $M_{\uparrow}$, which is evolved simultaneously with the photon number. Since the photon number now varies continuously, the number of molecules in the ground ($M_{\downarrow}$) and excited ($M_{\uparrow}$) states are no longer integers.

III.2 Predictions

If $M_{\uparrow},M_{\downarrow}$ are sufficiently large such that $R$ and $\gamma$ remain approximately constant, the following analytic discussion applies: from geometric reasons, the phase evolves as Henry (1982)

[TABLE]

From this, we obtain $d\operatorname{Mean}(\theta)=0$ and

[TABLE]

Meanwhile, the photon number evolves as

[TABLE]

such that

[TABLE]

and

[TABLE]

Despite its simplicity, the HPM is able to give a simple explanation for Schawlow-Townes phase diffusion Scully and Zubairy (1997).

A special case occurs when the particle number vanishes entirely. Here, Eq. (24) becomes singular. This means that, in absence of other photons, the phase of the spontaneous emission is entirely random over the interval $[0,2\pi[$ and a so-called phase jump occurs Schmitt (2015); *thesisJulianEng.

In Figs. 1 and 3, we see that the predictions of the HPM match the exact values for the number statistics, as described in the previous section. In Fig. 4 (a,b,c) we show the evolution of phase and particle number of a typical HPM sample.

Despite these successes, there are some conceptual difficulties with the HPM. First of all, by treating the field as classical, any squeezing effects are disregarded. Secondly, the continuous variation of the field is inconsistent with the fact that the emission and absorption are discrete processes at the level of the dye. Thirdly, it is physically dubious that only spontaneous emission and neither absorption nor stimulated emission cause shot noise. In the phasor model, the overly large noise from the spontaneous emission actually mimics the shot noise from all loss and gain processes. Note that even though a single spontaneous emission event adds one photon on average, this number has an uncertainty of $\sqrt{n}$, such that a spontaneous emission can even decrease the photon number. Finally, regarding the phase, one may wonder what is so special about these spontaneous emission events that these influence $\theta$ whereas other processes do not.

IV The Quantum Trajectory Description

IV.1 A quantum model

In order to address the questions posed at the end of the previous section, we will study the photon field on a fully quantum mechanical level. Because of its driven-dissipative nature, the photon condensate is an open quantum system. The study of an open system starts with the distinction between system and environment. Here, we will treat the condensate mode as the system to be modelled by a quantum-mechanical stochastic wavefunction $\ket{\psi(t)}$. We will use the notation $\langle\cdot\rangle$ to denote quantum expectation values with respect to this wavefunction. The gain medium on the other hand is modeled classically: of the $M_{\text{tot}}$ dye molecules, the integer amounts $M_{\uparrow}(t)$($M_{\downarrow}(t)$) are in the excited(ground-) state as described by stochastic rate equations given in Sec. II. This is justified because of the frequent collisions with solvent molecules that lead to thermalization Kirton and Keeling (2015).

The coherent evolution of the photons, in the frame rotating at the dye transition frequency, is governed by the Hamiltonian

[TABLE]

where operators $\operatorname{\hat{a}}(\operatorname{\hat{a}^{\dagger}})$ annihilate (create) a photon, defining also the number operator $\operatorname{\hat{n}}=\operatorname{\hat{a}^{\dagger}}\operatorname{\hat{a}}$. Again, for numerical purposes we can go to a rotating frame and replace $\Delta$ by arbitrary $\Delta^{\prime}$. In absence of Hamiltonian dynamics, the evolution of the photon field is governed by three processes: gain (corresponding to emission of the dye molecules) occurring at rate $R(\operatorname{\expectationvalue{\operatorname{\hat{n}}}}+1)$, absorption at rate $\gamma\operatorname{\expectationvalue{\operatorname{\hat{n}}}}$ and external losses through the mirrors at rate $\kappa\operatorname{\expectationvalue{\operatorname{\hat{n}}}}$. Because of the discrete nature of the dye excitations, the first two processes are naturally described through a photon-counting unraveling as if the photon number is ‘measured’ by the dye. External losses are modelled by a heterodyne unraveling, firstly because heterodyne detection is typically performed in experiments on this leaking current, and secondly because it keeps the wave function localized in phase space, which will be helpful. By combining these processes, we readily obtain a stochastic Schrödinger equation Breuer and Petruccioni (2002)

[TABLE]

Here, the tilde denotes that the left hand side describes the unnormalised wavefunction, to be renormalised after every timestep. $dZ=\frac{1}{\sqrt{2}}(dW_{x}+idW_{p})$ is a complex (Itô) Wiener noise process such that $|dZ|^{2}=dt$ and $dN$, $dM$ describe Poisson processes as defined in section II, where the role of $n$ is now replaced by $\expectationvalue{\operatorname{\hat{n}}}$. Because the photon number typically becomes mesoscopic, exact evolution of (IV.1) in a truncated Fock basis rapidly becomes computationally unfeasible. A number of variational approaches have been proposed for efficient simulation to this extent Daley (2014); Casteels et al. (2018). Here, the fact that the system is a single bosonic mode with large occupation, combined with a visual inspection of the Wigner function in Fig. 5 leads us to model the field as being Gaussian in particle number and phase: an ‘$N\Theta$-Gaussian’ state as was also used in Verstraelen and Wouters (2018). Under this assumption, the field $\ket{\psi}$ is entirely characterized by the expectation values $\operatorname{\expectationvalue{\operatorname{\hat{n}}}},\operatorname{\expectationvalue*{\operatorname{\hat{\operatorname{\theta}}}}},\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{n}}\operatorname{\hat{\delta}_{n}}}},\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{\operatorname{\theta}}}\operatorname{\hat{\delta}_{\operatorname{\theta}}}}}$ and $\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{n}}\operatorname{\hat{\delta}_{\operatorname{\theta}}}}_{\text{sym}}}$, where $\operatorname{\hat{\delta}}_{O}:=\operatorname{\hat{O}}-\expectationvalue{\operatorname{\hat{O}}}$ is defined as the fluctuation of operator $\operatorname{\hat{O}}$. The full variational equations for these expectation values are given in appendix A. Note that if external decay vanishes, or if it is described as a photon-counting process, we find that $\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{n}}\operatorname{\hat{\delta}_{n}}}}\rightarrow 0,\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{\operatorname{\theta}}}\operatorname{\hat{\delta}_{\operatorname{\theta}}}}}\rightarrow\infty$ and the equations reduce reduce to stochastic rate equation as in section II. Interestingly, it are thus the (small) losses that allow us to define a phase, and avoid the condensate to become a particle number eigenstate Mølmer (1997).

Because in the grand-canonical regime the particle number fluctuates also to vanishingly small densities we will use a combination of both the exact (IV.1) and variational (50)-(54) evolution. That is, we define treshold particle numbers $n_{\text{trans},\searrow}$, $n_{\text{trans},\nearrow}$. When $\operatorname{\expectationvalue{\operatorname{\hat{n}}}}<n_{\text{trans},\searrow}=200$ we perform the exact evolution and when $\operatorname{\expectationvalue{\operatorname{\hat{n}}}}>n_{\text{trans},\nearrow}=240$ we perform the variational evolution. In the intermediate window $n_{\text{trans},\searrow}<\operatorname{\expectationvalue{\operatorname{\hat{n}}}}<n_{\text{trans},\nearrow}$, the method of the previous regime remains in use. The transitions between methods are discussed in more detail in Appendix B. Whereas the variational equations are solved by a straightforward Euler method with direct addition of the Poisson increments, for proper numerical convergence of the exact trajectories, an approach where the deterministic and diffusive evolution is separated from the individual jumps Bruti-Liberati and Platen (2007) is used.

IV.2 Predictions

Similar to the HPM, we calculate the evolution of the moments of the field while treating $\gamma,R$ as constant. We will use the variational equations from Ap. A as a starting point, considering highest order $\mathcal{O}(\operatorname{\expectationvalue{\operatorname{\hat{n}}}})=\mathcal{O}(\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{n}}\operatorname{\hat{\delta}_{n}}}})=\mathcal{O}(\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{\operatorname{\theta}}}\operatorname{\hat{\delta}_{\operatorname{\theta}}}}}^{-1})$ (this relation between orders can be seen from a coherent state where $\operatorname{\expectationvalue{\operatorname{\hat{n}}}}=\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{n}}\operatorname{\hat{\delta}_{n}}}}$ and relation (57)). Under the assumption $\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{n}}\operatorname{\hat{\delta}_{\operatorname{\theta}}}}_{\text{sym}}}\equiv 0$ and neglecting noise, one obtains as stationary solution for the phase variance (53),

[TABLE]

Substituting in (52) results in

[TABLE]

where $dW_{\theta}$ denotes a real Wiener process from which

[TABLE]

and

[TABLE]

In trajectory simulations, the total variance of an observable $\operatorname{Var}[\operatorname{\hat{O}}]=\operatorname{\operatorname{Var}_{1}}{\operatorname{\hat{O}}}+\operatorname{\operatorname{Var}_{2}}{\operatorname{\hat{O}}}$ can be decomposed into the intra-trajectory variance $\operatorname{\operatorname{Var}_{1}}{\operatorname{\hat{O}}}=\overline{\expectationvalue{\delta_{\operatorname{\hat{O}}}\delta_{\operatorname{\hat{O}}}}}$ and the inter-trajectory variance $\operatorname{\operatorname{Var}_{2}}{\operatorname{\hat{O}}}=\operatorname{Var}\expectationvalue{\operatorname{\hat{O}}}$ Breuer and Petruccioni (2002). Here, from our stationary assumption, $d\operatorname{\operatorname{Var}_{1}}{\operatorname{\hat{\operatorname{\theta}}}}=0$ such that $d\operatorname{Var}[\operatorname{\hat{\operatorname{\theta}}}]=d\operatorname{\operatorname{Var}_{2}}{\operatorname{\hat{\operatorname{\theta}}}}$.

Similarly, we obtain the stationary solution for the particle number from (51):

[TABLE]

where we have used $\kappa\ll\gamma,R$. Substituting in (50) results in

[TABLE]

and

[TABLE]

Comparing Eqs. (32),(33),(35),(IV.2) with the predictions from the HPM (24),(26),(27) the similarity is evident. In line with the expectations, the HPM predicts a variance increase proportional to the gain $R$, whereas this is symmetrical in $R,\gamma,\kappa$ according to our trajectory result. In the steady-state however, $R\approx\gamma+\kappa$ at least to highest order in $\operatorname{\expectationvalue{\operatorname{\hat{n}}}}$, such that the time-averaged result is the same. This correspondence is also reminiscent to the result for a laser Scully and Zubairy (1997): also there, the phasor model predicts a phase diffusion proportional to the gain coefficient $R$ whereas a more detailed quantum-mechanical derivation makes clear that $\frac{R+\kappa}{2}$ is the correct quantity, although both become equivalent at threshold where $R\approx\kappa$.

In Fig. 4 (d,e,f) a representative evolution of phase and particle number during a single trajectory simulation is shown. On short times there is a qualitative difference with the HPM: discrete steps corresponding to emission events are replaced by scale-invariant noise, induced by the heterodyne detection and fueled by all emission and absorption processes. However, current experiments cannot resolve these small fluctuations on short timescales, as long as they do not affect correlation functions. It is clear from fig. 4 that the predictions of both methods become qualitatively indistinguishable on longer timescales. We see also that for both methods, phase diffuses more rapidly at lower photon numbers, as predicted by Eqs. (24),(33).

V First order coherence

In Fig. 6 the first order correlation function $g^{(1)}(\tau)=\frac{\overline{\alpha^{*}(t)\alpha(t+\tau)}}{\overline{\absolutevalue{\alpha(t)}^{2}}}$ is shown, and again predictions of the HPM agree well with the exact result, obtained by trajectories. In the noninteracting case, the decay of correlations is clearly exponential, whereas it is Gaussian in presence of interactions. The decay can be understood from multiple points of view, let us start with the noninteracting case.

V.1 Schawlow-Townes effect from the canonical to the grandcanonical regime: the influence of ‘phase jumps’

According to the HPM picture, absorption and stimulated emission only affect the phasor radially, as $d\sqrt{n}=\frac{dn}{2\sqrt{n}}=\frac{R-\gamma}{2}\sqrt{n}\,dt$, where we used eq. (22), and neglected the external losses. The stochastic evolution of the HPM can then be written in terms of the phasor alone as

[TABLE]

where $dS$ is additive noise corresponding to the spontaneous emissions. In the grandcanonical limit, $\gamma$ and $R$ can be treated as constants, such that for the expectation value,

[TABLE]

From the quantum regression theorem Scully and Zubairy (1997), it is then clear that

[TABLE]

where

[TABLE]

By another line of reasoning, the phase evolution is dominated by large ‘phase jumps’ when the photon number in the cavity vanishes as described in Schmitt et al. (2016); Schmitt (2015); *thesisJulianEng. From estimating the probability of having zero photons in the cavity, one obtains the ‘phase jump rate’ $\Gamma_{\text{PJ}}^{0}=\frac{B_{12}M_{\downarrow}}{\overline{n}^{\zeta}}$. In the grandcanonical limit, $\zeta=1$, so that this is consistent with (40) up to a scaling factor of order one. The picture of phase jumps would further predict that the phase evolution is suppressed if the probability for the zero-photon state vanishes ($\zeta\rightarrow\infty$), as occurs outside of the grandcanonical limit, towards the canonical regime $M_{\text{eff}}\ll\overline{n}^{2}$, where $M_{\text{eff}}$ is defined as in Eq. (11). By additional numerical simulations shown in Fig. 7, we see that this is not the case. At most, the first-order coherence time, rescaled with the molecule number, only scales by a factor two. This is not entirely unsurprising as the remaining value corresponds to the standard Schawlow-Townes dephasing Scully and Zubairy (1997) occurring in the limit where the number distribution is sharply peaked, as can already be seen from (24). We are thus led to write generically

[TABLE]

where $\xi=2$ in the grandcanonical limit and $\xi=4$ in the canonical limit. Comparing with the decay of second-order coherence (13), we find that

[TABLE]

where we have used (17) for the third equality. However, on the inset of fig. 7, we see that on sufficiently short timescales, there is always an initial decay with $\xi=2$. There is a clear intuition here. In general, $\gamma-R=\overline{\gamma}-\overline{R}+(B_{12}+B_{21})\delta n$, where $\delta n=n-\overline{n}$. When, at short timescales, $\delta n$ remains approximately constant, a similar reasoning to the grandcanonical regime yields

[TABLE]

At later times, higher-order effects set in and restrict the dynamics. It is instructive to compare these results with the ones obtained by Whittaker and Eastham Whittaker and Eastham (2009); Love et al. (2008) in a polariton context. There, a Schawlow-Townes decay is predicted to be of the form

[TABLE]

For $\Gamma_{2}\tau\ll 1$, this reduces to an exponential decay at rate $\frac{\eta^{2}}{2}\Gamma_{2}$, whereas for $\Gamma_{2}\tau\gg 1$, (46) decays exponentially at rate $\frac{\eta^{2}}{4}\Gamma_{2}$. This can be understood because at short times the number fluctuations contribute to the decay of $g^{(1)}$, but after a time $1/\Gamma_{2}$ only Schawlow-Townes phase diffusion remains. There remain a few differences between the physics of (46) and the photon condensate. First, (46) is derived in Whittaker and Eastham (2009); Love et al. (2008) under the explicit assumption that the whole ensemble has a Gaussian number distribution peaked around $\overline{n}$ 111This is different from our ansatz where we consider only individual trajectories as having a Gaussian number distribution.. This assumption is physical for the photon condensate only in the canonical regime, where the probability of having zero photons is negligible. This means that the predicted long time value of $\xi$ equals 4, in agreement with our prediction for the canonical regime.

Finally, the timescale of the transition between the two decay rates is for photon condensates not determined by $\Gamma_{2}$, but only by the second term in (13), which is responsible for the time-dependence of $R$ and $\gamma$. In the grandcanonical limit, this means that a slowing of the decay would only take place at infinitely long times unlike the prediction of (46), as is seen on Fig. 6.

V.2 Effect of Kerr-interactions

In presence of finite photon-photon interactions, the decay of first-order coherence is altered, as we see on the right panel of Fig. 6. The profile is rather Gaussian than exponential. This is characteristic for the so-called Henry mechanism Henry (1982). Whereas Schawlow-Townes decay is attributed to direct fluctuations of the phase, Henry decay results from phase diffusion as a consequence of number fluctuations causing a change of the interaction energy. For the Henry effect, Ref. Whittaker and Eastham (2009); *Whittakerpolaritons2 predicts an additional decay

[TABLE]

which reduces to a Gaussian decay with characteristic time $(\sqrt{2}/U\overline{\operatorname{\expectationvalue{\operatorname{\hat{n}}}}}\eta)$ at short timescales. As we see on figure 6, there is good agreement with our numerical results on short times, although deviations occur at later times that we can attribute to the non-Gaussian character of the number distribution.

The fact that the prediction of Gaussian decay for short timescales remains valid in the grandcanonical limit can be understood because in a frame rotating at the bare cavity frequency, the expectation value of a phasefactor of a state with $n$ photons at time $t$ is

[TABLE]

where $k_{m}$ is the $m$th cumulant of the distribution of $n$. Here, we implicitly assumed that $n$ remains approximately constant on short times.

VI Conclusions

We have studied the temporal coherence of a single-mode photon condensate both with and without weak Kerr interactions, and have also accounted for the effects of driving and dissipation. We have shown thermodynamically and dynamically how interactions reduce the number fluctuations and calculated the corresponding enhancement of the decay of second-order correlations, which we numerically verified by stochastic rate equations. The driven-dissipative nature of a realistic photon condensate causes $g^{(2)}(\tau)<1$ at large $\tau$, because of the coupling of the photon number to the slow reservoir dynamics. We have reviewed the heuristic phasor model, and shown analytically and numerically its adequacy in describing experimentally relevant quantities while being numerically efficient. For comparison, a quantum mechanical model of the condensate mode based on variational quantum trajectories was introduced, assuming only classical coupling with the dye-molecules. Exhibiting similarities with both a laser far from equilibrium Henry (1982) and an isolated atomic BEC Sinatra et al. (2009), we have observed both Schawlow-Townes (exponential) and Henry (Gaussian) contributions to the decay of phase correlations in the photon condensate, depending on the interaction strength. As the effect of ‘phase jumps’ is similar to the standard Shawlow-Townes effect, we found no qualitative differences in the coherence between canonical and grandcanonical regimes. We have shown how first- and second-order coherence times are related by the number fluctuations. An interesting open question is to what extent the picture above changes in presence of thermo-optical effects. To study the latter, longer evolution times are necessary, but as we have shown an approach based on the numerically efficient heuristic phasor model is likely to be sufficient. Another possible extension is the study of coupled condensate modes, as appear in lattices of cavities. Also within a single cavity, there can be interesting physics arising from mode competition for some parameter values Hesten et al. (2018), which we have omitted here.

Finally, we have assumed here that the absorption and emission processes are Markovian and coupling between the photons and molecules is weak.

It remains an open question to what extent our conclusions remain valid outside of these approximations de Vega and Alonso (2017). Because our conclusions in strong thermal equilibrium and weak losses are similar to a threshold laser far from equilibrium, both regarding the validity of the HPM and the shape of the autocorrelation functions, we expect these conclusions to remain valid for all usual parameter values for photon condensates.

Acknowledgements.

We acknowledge stimulating discussions with J. Schmitt, M. Weitz, F. Öztürk and H. Kurkjian. We furthermore thank J. Schmitt and M. Weitz for their feedback on our manuscript. This work was financially supported by the FWO Odysseus program. Part of the computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government – department EWI.

Appendix A The variational evolution

Following reference Verstraelen and Wouters (2018) (which has however a different sign convention for the parameter $\Delta$), the trajectory evolution (IV.1) can be recast in terms of (normalized) expectation values $\expectationvalue{\operatorname{\hat{O}}}$ to be

[TABLE]

where again $\hat{\delta_{O}}=\hat{O}-\expectationvalue{O}$. By introducing a Dirac phase through $\operatorname{\hat{a}}=:e^{i\hat{\theta}}\sqrt{\operatorname{\hat{n}}}$ satisfying $[\operatorname{\hat{n}},\operatorname{\hat{\operatorname{\theta}}}]=i$ and using Wick’s theorem, we obtain for the Gaussian correlation functions $\operatorname{\expectationvalue{\operatorname{\hat{n}}}},\operatorname{\expectationvalue*{\operatorname{\hat{\operatorname{\theta}}}}},\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{n}}\operatorname{\hat{\delta}_{n}}}}=\operatorname{\expectationvalue*{\operatorname{\hat{n}}\operatorname{\hat{n}}}}-\operatorname{\expectationvalue{\operatorname{\hat{n}}}}^{2},\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{\operatorname{\theta}}}\operatorname{\hat{\delta}_{\operatorname{\theta}}}}}=\operatorname{\expectationvalue*{\operatorname{\hat{\operatorname{\theta}}}\operatorname{\hat{\operatorname{\theta}}}}}-\operatorname{\expectationvalue*{\operatorname{\hat{\operatorname{\theta}}}}}^{2},\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{n}}\operatorname{\hat{\delta}_{\operatorname{\theta}}}}_{\text{sym}}}=\operatorname{\expectationvalue*{\operatorname{\hat{n}}\operatorname{\hat{\operatorname{\theta}}}}}/2+\expectationvalue*{\operatorname{\hat{\operatorname{\theta}}}\operatorname{\hat{n}}}/2-\operatorname{\expectationvalue{\operatorname{\hat{n}}}}\operatorname{\expectationvalue*{\operatorname{\hat{\operatorname{\theta}}}}}$ the evolution

[TABLE]

Here, the coefficients $C$ and $D$ are defined as in Verstraelen and Wouters (2018) and

[TABLE]

By considering the measurements to be perfect, we can do a restriction towards pure states, for which relation

[TABLE]

is satisfied Verstraelen and Wouters (2018). This constraint allows to compute $\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{\operatorname{\theta}}}\operatorname{\hat{\delta}_{\operatorname{\theta}}}}}$ (or another variance) directly from the others.

Appendix B Transitions between exact and variational regimes

B.1 From variational to exact

When $\operatorname{\expectationvalue{\operatorname{\hat{n}}}}$ decreases from the variational regime below the threshold $n_{\text{trans},\searrow}$, the $N\Theta$-Gaussian state must be explicitly expressed in Fock-base to continue numerically exact evolution.

A generic Gaussian density operator Gardiner and Zoller (2004) can, by definition, be written as

[TABLE]

for a quadratic $\hat{H}_{\text{eff}}$. A Gaussian state that is pure is obtained by taking the limit $\beta\rightarrow\infty$, which is equivalent with taking the lowest eigenvector (‘ground state’) of $\hat{H}_{\text{eff}}$.

In order to construct $\hat{H}_{\text{eff}}$, we need explicit matrix representations of $\operatorname{\hat{n}}$ and $\operatorname{\hat{\operatorname{\theta}}}$. For the particle number operator $\operatorname{\hat{n}}=\operatorname{\hat{a}^{\dagger}}\operatorname{\hat{a}}$ this is straightforward. Regarding $\operatorname{\hat{\operatorname{\theta}}}$ we encounter the fact that phase is no true observable, with the consequence that, strictly speaking, no hermitian phase operator exists Nieto (1993). However, as long as the Fock space is truncated (which is the case here) at level $N_{\text{max}}$ , a meaningful phase operator can be obtained through the Pegg-Barnett formalism to be

[TABLE]

Nieto (1993). One free parameter, $\theta_{0}$, remains, corresponding to the phase-cut. That is, because phase is a periodic variable, there must be a cut where (going counterclockwise) the phase sharply changes with $-2\pi$. In order to avoid secondary effects of this cut, we will use $\theta_{0}=\operatorname{\expectationvalue*{\operatorname{\hat{\operatorname{\theta}}}}}-\pi$ such that the $N\Theta$-Gaussian state is as far from the phase-cut as possible.

Using the operators $\operatorname{\hat{n}}$ and $\operatorname{\hat{\operatorname{\theta}}}:=\operatorname{\hat{\operatorname{\theta}}}_{PB}(\operatorname{\expectationvalue*{\operatorname{\hat{\operatorname{\theta}}}}}-\pi)$, an effective Hamiltonian

[TABLE]

is constructed. $\ket{\psi}$ is obtained as the eigenvector of $\hat{H}_{\text{eff}}$ corresponding to the lowest eigenvalue.

B.2 From exact to variational

When, in the exact regime, $\operatorname{\expectationvalue{\operatorname{\hat{n}}}}=\bra{\psi}\operatorname{\hat{n}}\ket{\psi}$ increases above $n_{\text{trans},\nearrow}$, the Gaussian moments must be computed. As in the previous case, this is entirely straightforward regarding $\operatorname{\expectationvalue{\operatorname{\hat{n}}}}$ and $\operatorname{\expectationvalue*{\operatorname{\hat{\delta}_{n}}\operatorname{\hat{\delta}_{n}}}}$. For the phase, we want to use an operator $\operatorname{\hat{\operatorname{\theta}}}_{PB}(\theta_{0})$ again, although it is this time not a priori clear which value of $\theta_{0}$ to use. We will therefore do an initial guess $\theta_{0}^{(0)}$ and construct a corresponding $\operatorname{\hat{\operatorname{\theta}}}_{PB}^{(0)}=\operatorname{\hat{\operatorname{\theta}}}_{PB}(\theta_{0}^{(0)})$. We then iteratively use $\theta_{0}^{i+1}=\bra{\psi}\operatorname{\hat{\operatorname{\theta}}}_{PB}^{(i)}\ket{\psi}-\pi$ and repeat the procedure self-consistently until convergence is reached. The resulting $\operatorname{\hat{\operatorname{\theta}}}_{PB}(\theta_{0}^{(f)})$ can then be used to proceed in calculating the Gaussian moments. We have verified that the above procedure results in the expected phase $\arg{\left[\alpha\right]}$ when applied to an arbitrary coherent state $\ket{\alpha}$.

After this transition, we use purity relation (57) as a numerical check and as verification for the validity of the Gaussian ansatz.

The fact that in this work two different definitions of a phase operator (Dirac and Pegg-Barnett) are used causes no problems because the particle number at the transitions is sufficiently large and the trajectory states are sufficiently well-behaved.