Radiation statistics of a degenerate parametric oscillator at threshold

We thank the referee for reviewing our work and for the overall positive feedback. In the following we will address the points of concern raised by the referee.

The referee writes:

The authors start by presenting a presumably general classical Langevin equation for parametric oscillators (1), which is specified in Eq (2), where the function f is defined. The authors comment that an integer parameter m in this function takes value 2 for a “generic oscillator” and value 4 for the Duffing oscillator. However, all the results in the paper are derived from Eq. (3), which is not a corresponding FPE but is an effective equation for slow dynamics derived from the quantum Lindblad equation (11).

Our response:

The referee has misunderstood this part of our manuscript by assuming that Eq. (3) is derived from Eq. (11). In fact Eq. (3) directly follows from Eq. (2) by the equivalence of a Fokker-Planck equation to its Langevin description. We have added a new appendix A to make this important point more clear and provide the necessary steps of the derivation. For the argument of the universality of the Liouvillian, it is crucial that (3) [as well as (2)] simply follow from symmetry arguments and the fact that the nonlinearity is weak.

The derivation is not presented in the main text but is given in the Supplement. Here the authors explicitly write the Hamiltonian for two specific cases, voltage-biased Josephson junction and Duffing oscillator. One would assume, logically, that the former corresponds to a “generic oscillator” in (2) with m=2, and the latter to Duffing oscillator with m=4. However, this is not the case, both fall in the class m=2. A possible reason for this could be a different symmetry of the driving term in (11) compared to the common one due to the parametric driving of the oscillator frequency.

We are afraid but we cannot follow the referee at this point. We never state that the Duffing oscillator corresponds to m=2. In fact it is stated at several points (the first being the first paragraph of the introduction) that the Duffing oscillator is not generic (with m=4) due to the additional symmetry. This result is not new but rather well-known in the literature, see e.g. Ref. [4].

Furthermore, the presence of dissipative constant in the nonlinear terms in (11) and (2) indicates that the source of this nonlinearity is a multiphoton coupling to the environment rather than, commonly, a nonlinear energy of the oscillator. Whether this is taken care of in Eq. (11) remains unclear since the quantity $L$ in the jump operator is never defined.

The fact that $\gamma$ appears in our expression of $f(x)$ is simply due to the fact that we would like to measure the nonlinearity with a dimensionless parameter $\alpha$. Close to threshold, we have $\epsilon \approx \gamma$ as the relevant bare frequency scale. We have added a new footnote 3 in the revised manuscript in order to explain this fact. Note that in Eq. (11) the (generic) jump operator $L$ is specified to $a$ and $a^\dagger$ as the terms describing the emission and absorption of photons.

The authors should explain in the very beginning what physical oscillators are in mind.

We have added two sentences at the end of the third paragraph of the introduction listing physical examples, where our results can be realized. The systems cover a broad range of physical implementations.

An exact solution for a parametrically driven quantum nonlinear (Duffing) oscillator was derived in Opt. Commun. 127, 230, 1996 (not quoted). There the cumulants at the threshold were computed, and it would be useful to make a comparison with the results obtained in this paper.

We would like to thank the referee for bringing the reference to our attention. However, we would also like to point out that the results of that paper, as far as we understand, only concern the stationary state inside the cavity. However, our results are valid for the photon-current leaving the cavity. For the latter, the dynamics of the system is important, not only the stationary state. In particular, the Fano factor of the photon-counting depends on the second-order coherence $g^{(2)}(\tau)$ for all times, whereas the paper mentioned by the referee only discusses $g^{(2)}(0)$. Because of this, it is not possible to compare the results to ours. Still, we have included the reference at the appropriate point in the introduction.

It remained unclear to me whether Eq. (3) describes classical or quantum fluctuations. In the latter case, the results should refer to intrinsic fluctuations of the oscillator since the definition of the output photon flux adopted in the paper does not include the quantum correlations.

Equation (3) describes both quantum and classical fluctuations. These correspond to $\bar n$ (classical) and $\frac12$ in Eq. (2). The fluctuations are included in a quasiclassical fashion that is possible for an oscillator with a weak nonlinearity (see Ref. [31] for details). The referee is correct that one has to be careful about proper normal ordering of the output photon flux in order to exclude vacuum photons in the counting. Note however that close to threshold the typical number of photons $N_0$ is large such that normal ordering is only an irrelevant correction. We have added a sentence to footnote 6 to explain this important point.

Summarizing, the paper meets the Expectation criterium 3: Opens a new pathway in an existing research direction, with clear potential for multipronged follow-up work. The paper meets all the Acceptance criteria except for (partially) item 3, see the report. The mentioned shortcomings must be explained in the revised manuscript

The revised manuscript has considerably improved the presentation of our results. We are confident that it now provides sufficient details such that the arguments and derivations can be reproduced.

We thank the referee for reviewing our work and for the positive feedback. In the following we will address the points of concern raised by the referee.

The referee writes:

1) The abstract is not sufficiently detailed. At least, there should be a short mentioning of the method used to obtain the results. Otherwise, the paper fails to fulfill point 2 of the general acceptance criteria and cannot be published.

Our response:

We have added two sentences to the abstract in order to explain our method.

2) While the mathematics of the paper is presented in a sound and rigorous way, the authors could perhaps elaborate a bit more on the physical interpretation of the results. For instance, could the authors please comment on why the Fano factor reaches its maximum for a driving stronger than the threshold? Is this an expected/unexpected result?

We would like to point out that close to the threshold both the nonlinearity as well as the quantum/classical fluctuations are important. In this fluctuation dominated regime, it is notoriously difficult to make a simple physical intuition work. The main point of our paper is a data-collapse; i.e., the numerical solution of the universal Liouvillian (Fig. 3) can be used to "predict" the Fano factor in all parametric systems close to the threshold. The microscopic details only enter via the scale factors $x_\ast, \tau_\ast$.

Nevertheless, we would like to offer to the referee the way we understand the fact that the maximum of Fig. 3 is for positive $\beta$: it is rather clear, that the Fano factor (starting at $F=2$ for low driving and approaching a coherent state with $F=1$ above threshold) has a maximum in between. A good description of the physics corresponding to the $\mathbb{Z}_2$ symmetry-breaking is a quantum particle in a double-well potential $-\beta x^2 + x^4$ where $\beta$ is tuned (see also the answer to point 4). Below the threshold ($\beta<0$), only a single minimum is present. Increasing $\beta$, the variance of the position increases (which corresponds to larger noise and Fano factor). However, the variance does not have a maximum at the threshold ($\beta=0$). Making $\beta$ slightly positive, the particle still gets more "room" to fluctuate as the confinement becomes weaker. Only at a larger value of $\beta$, the particle is confined to one of the two minima at finite $x_s$ and the variance (which is now around this minima) decreases.

3) Do the authors expect higher-order cumulants to fulfill universal ratios similar to the one that they identify for the first three cumulants [see Eq. (9)]? Why/why not?

There is an "infinite" number of other ratios that are also universal. We have given the one involving the lowest-cumulants as these are the easiest to measure. One can find other ratios from Eq. (8). The rule of the game is to come up with ratios where the factors $N_0$ and $\tau_*$ drop out. Other such ratios are for example $\frac{\langle N\rangle \langle\langle N^4\rangle\rangle}{ \langle\langle N^2 \rangle\rangle \langle\langle N^3 \rangle\rangle}= \frac{c_1 c_4}{c_2 c_3}$ or even $\frac{\langle N\rangle^2 \langle\langle N^4\rangle\rangle}{ \langle\langle N^2 \rangle\rangle^3 }= \frac{c_1^2 c_4}{c_2^3}$. Note that these other ratios involve higher-order cumulants.

4) Higher-order cumulants play an important role for determining Lee-Yang zeros and thus predicting the occurrence of phase transitions. Could the authors comment on the importance of their own results for further investigations on the behavior of Lee-Yang zeros for driven systems? What new insights about the physics at threshold can be gained from the results in this regard?

As far as we understand the theory of Lee-Yang zeros for open systems, the question is about the zeros of $Z(s) = \exp( \lambda(s) \tau/\tau_*)$. It is possible to map our Eq. (3) onto the Schrödinger equation of a (quantum) particle in a potential (see e.g. https://doi.org/10.1016/j.physa.2017.10.021). Using this mapping, we believe that we can show that $\lambda(s)$ remains finite for all $s\in\mathbb{C}$. The reason is that the potential, due to the finite nonlinearity, is confining for all $s$. Therefore, our result corresponds to the fact that $Z(s)$ is an analytical function without any zeros. The physical reason is that all the cumulants remain finite, due to the inclusion of the nonlinearity. The universal scaling of the cumulants (in the long-time limit at threshold) is then given by the scaling as a function of $\alpha$ as the cumulants diverge only for $\alpha \to 0$.