Statistics of radiation due to non-degenerate Josephson parametric down conversion

TL;DR

This paper provides a theoretical analysis of the photon statistics and correlations in non-degenerate Josephson parametric down-conversion, revealing how asymmetries and thermal effects influence radiation properties.

Contribution

It introduces a comprehensive theoretical framework for the full counting statistics and coherence properties of non-degenerate Josephson parametric oscillators, including effects of asymmetry and thermal noise.

Findings

Photon pair emission is affected by mode linewidth asymmetry.

Strong violation of Cauchy-Schwarz inequality occurs at specific driving strengths.

Thermal effects modify the radiation statistics and coherence properties.

Abstract

In a process called parametric down-conversion, a dc-biased Josephson junction coupled to two microwave resonators emits photon pairs when the Josephson frequency matches the sum of the two resonance frequencies. Recent experiments have shown that such a setup permits analyzing the correlation of the radiation. Motivated by these results, we study theoretically the full counting statistics of a non-degenerate parametric oscillator below the threshold of instability. Furthermore, we analyze the second-order coherences of the radiation and discuss thermal effects on the radiation statistics. We provide results for the driving strength at which the Cauchy-Schwarz inequality is most strongly violated. Additionally, we study the impact of asymmetry in the linewidth of the modes---a distinctive property of the non-degenerate resonance effect. In particular, we find that the radiation from the mode with the larger linewidth preferably takes the total detuning, while the other mode emits photons at its resonance frequency.

Equations

$$Z_k(\omega )=\frac{Z_k^0{\Omega }_k}{{\gamma }_k-i(\omega -{\Omega }_k)},$$

$$\mathcal{Z}=\int \mathcal{D}[{\varphi }_a^{+}]\mathcal{D}[{\varphi }_a^{-}]\mathcal{D}[{\varphi }_b^{+}]\mathcal{D}[{\varphi }_b^{-}]e^{iS[{\varphi }_a^{+},{\varphi }_a^{-},{\varphi }_b^{+},{\varphi }_b^{-}]}.$$

$$\displaystyle S_{J}=\frac{E_{J}}{\hbar\Gamma}\int\!dt\Big{\{}\!\cos$$

$$\displaystyle-\cos\!\big{[}\phi_{V}(t)+\phi_{a}^{-}(t)+\phi_{b}^{-}(t)\big{]}\Big{\}},$$

$$S_k$$

$$M_k(\omega )$$

$$+[2n_T(\omega )+1]\mathrm{Re}\left[Z_k(\omega {)}^{-1}\right](1-{\alpha }_k)\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$$

$$\displaystyle\quad+\![2n_{0}(\omega)+1]\mathrm{Re}\!\left[Z_{k}(\omega)^{-1}\right]\alpha_{k}\!\begin{bmatrix}1\!&-1\\ -1\!&1\end{bmatrix}\!\Bigg{\}},$$

$${\varphi }_k^{\pm }(t)=\mathrm{Re}[e^{-i({\Omega }_k+\Delta /2)t}{\psi }_k^{\pm }(t)].$$

$$S_J=-\frac{E_J}{2\hslash \Gamma }\int dt\{s_J[{\psi }_a^{+}(t),{\psi }_b^{+}(t)]-s_J[{\psi }_a^{-}(t),{\psi }_b^{-}(t)]\},$$

$$s_J({\psi }_a,{\psi }_b)=\frac{J_1(|{\psi }_a|)J_1(|{\psi }_b|)}{|{\psi }_a||{\psi }_b|}({\psi }_a{\psi }_b+{\psi }_a^{\ast }{\psi }_b^{\ast }),$$

$$S_k=\frac{i}{16\pi G_Q}\int$$

$$Z_k^0\overline{M}{}_k(\nu )=-i$$

$$+{\gamma }_k\left[\begin{array}{cc}2n_k+1& -2n_k\\ -2(n_k+1)& 2n_k+1\end{array}\right],$$

$$\frac{d{\psi }_k}{dt}=$$

$$\displaystyle\,\,-\!\psi_{k}\frac{J_{1}(|\psi_{\bar{k}}|)J_{2}(|\psi_{k}|)}{2|\psi_{\bar{k}}||\psi_{k}|^{2}}(\psi_{a}\psi_{b}\!+\!\psi_{a}^{\ast}\psi_{b}^{\ast})\Bigg{]}.$$

$$A=\frac{\sqrt{8{ϵ}_{\ast }}\sqrt{ϵ-{ϵ}_{\ast }}{\sigma }_a}{\sqrt{{ϵ}_{\ast }^2({\sigma }_a^2+1)+2{\Delta }^2({\gamma }_b+{\sigma }_a^2{\gamma }_a)}},B=\frac{A}{{\sigma }_a},$$

$$e^{i({\xi }_a+{\xi }_b)}=\frac{\Delta -i}{{ϵ}_{\ast }}.$$

$$S_J=-\frac{E_J}{8\hslash \Gamma }\int dt[{\psi }_a^{+}(t){\psi }_b^{+}(t)-{\psi }_a^{-}(t){\psi }_b^{-}(t)+c.c.],$$

$${\gamma }_k\left[\begin{array}{cc}2n_k+1& -2n_k\\ -2(n_k+1)[1+f_k(e^{i{\chi }_k}-1)]& 2n_k+1\end{array}\right].$$

$$P(N_k)=\int \frac{d{\chi }_k}{2\pi }\mathcal{Z}({\chi }_k)e^{-i{\chi }_k{N}_k}.$$

$$S=\frac{i}{16\pi G_Q}\int$$

$$A(\nu ,{\chi }_a,{\chi }_b)=$$

$$J=\frac{i2\pi E_J{G}_Q}{\hslash \Gamma }\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right].$$

$$g_{kl}^{(2)}(\tau )=\frac{⟨:I_l(\tau )I_k(0):⟩}{⟨I_k⟩⟨I_l⟩},k,l\in \{a,b\}$$

$${\bar{I}}_k=\frac{f_k{\gamma }_k(n_k+1)}{8\pi G_Q{Z}_k^0}⟨{\psi }_k^{-\ast }{\psi }_k^{+}⟩,$$

$${\bar{I}}_k=2f_k{\gamma }_k\frac{1+n_k}{{ϵ}_{\ast }^2-{ϵ}^2}\left\{n_k{ϵ}_{\ast }^2+{ϵ}^2\left[(1+n_{\bar{k}}){\gamma }_{\bar{k}}-n_k{\gamma }_k\right]\right\}.$$

$${\bar{I}}_k=$$

$$p(\nu )=$$

$$r$$

$$\gamma ={\gamma }_a-{\gamma }_b,|\gamma |<1.$$

$$⟨:I_l(\tau )I_k(0):⟩$$

$$\propto ⟨{\psi }_k^{-}(0{)}^{\ast }{\psi }_l^{-}(\tau {)}^{\ast }{\psi }_l^{+}(\tau ){\psi }_k^{+}(0)⟩.$$

$${\tau }_{ϵ}=\frac{1+{\gamma }^2{\Delta }^2}{{ϵ}_{\ast }({ϵ}_{\ast }-ϵ)(1-{\gamma }^2)}.$$

$$g_{aa}^{(2)}(0)=g_{bb}^{(2)}(0)=2,$$

$$g_{ab}^{(2)}(0)\stackrel{\text{ϵ≪ϵˉ}}{=}1+\frac{{ϵ}^2}{{ϵ}_{\ast }^2}\frac{(1+n_a+n_b{)}^2{\gamma }_a{\gamma }_b}{n_a{n}_b},$$

$$g_{ab}^{(2)}(0)\stackrel{\text{ϵ≈ϵ∗}}{=}2+\frac{2({ϵ}_{\ast }-ϵ)}{{ϵ}_{\ast }(1+n_a+n_b)}\left(1-n_a\frac{{\gamma }_a}{{\gamma }_b}-n_b\frac{{\gamma }_b}{{\gamma }_a}\right).$$

Full text

I Introduction

Parametric down-conversionguckenheimer is an important building block of many optical quantum information applicationsGisin:02 ; Scarani:05 . It describes the process that converts a pump photon of frequency $\Omega$ into a pair of photons with frequencies $\Omega_{a}$ and $\Omega_{b}$ which add up to the frequency of the pump photon, $\Omega_{a}+\Omega_{b}=\Omega$. This operation can be implemented by nonlinear oscillators with resonance frequencies $\Omega_{a}$ and $\Omega_{b}$ that are ac-driven below an instability threshold set by the nonlinear parameters of the oscillator. Parametric down-conversion has been employed to generate entangled photon pairsOu:92 as well as squeezed states of light.Wu:86

In most optical experiments, only a small fraction of the emitted pairs is detected due to a limited collection efficiency. Therefore, the unconditional emission events appear uncorrelated. However, recent experimentsAstafiev:07 ; Groot:10 ; hofheinz:11 ; Hoi:11 ; Liu:14 ; Chen:14 have advanced quantum non-linear optics into the microwave frequency range. In particular, the energy pumped into an electromagnetic environment by a dc-biased Josephson junction is of special interest, as all of the energy associated with a tunneling Cooper pair is confined to transmission lines and electrical oscillatorsDrummond:90 ; Padurariu:12 ; Leppakangas:13 ; Armour:13 ; Gramich:13 ; Souquet:14 ; Jin:15 ; Armour:15 ; Kubala:15 ; Leppakangas:16 ; Wustmann:17 . This constitutes a large technical advantage compared to optical experiments as the concentration of the radiation leads to almost perfect detection efficiencies.

Recently, accurate measurements of the radiation emitted by a dc-biased Josephson junction in series with two high-quality microwave resonators of different frequencies have demonstrated that such a setup permits the study of correlations.Westig:17 Motivated by these results, we want to investigate the radiation statistics of the non-degenerate setup illustrated in Fig. 1. The setup is composed of a Josephson junction, dc-biased with voltage $V$, that is in series with two microwave resonators with resonance frequencies $\Omega_{a}$ and $\Omega_{b}$. The Josephson frequency can be tuned to the sum of the two resonance frequencies, leading to parametric resonance.

In particular, we want to study the second-order coherences as well as the full counting statistics (FCS) of the resulting radiation. Note that the FCS of the non-degenerate optical parametric oscillator below the instability threshold has been studied previously in Ref. Vyas:92, for the specific case of zero detuning between the Josephson frequency and the sum of the two resonance frequencies as well as zero asymmetry in the linewidth of the modes. We revisit the FCS in the limit of the long measurement times employing a different framework based on the Keldysh path integral formalism previously applied in the context of the degenerate parametric oscillatorPadurariu:12 . We extend the earlier results to the case of an arbitrary mismatch between driving frequency and resonance frequencies as well as an arbitrary asymmetry in the linewidth of the modes. We discuss the qualitative differences between the degenerate and the non-degenerate case caused by the interplay of asymmetry and detuning. In particular, we find that the radiation from the mode with the larger linewidth preferably takes the total detuning, while the other mode emits photons at its resonance frequency. Furthermore, we analyze the second-order coherences of the radiation and investigate how thermal effects impact the radiation statistics.

The paper is organized as follows. We describe the setup in Sec. II. In Sec. III, we discuss its Keldysh action. Moreover, we introduce the counting fields which will later provide access to the radiation statistics. In Sec. IV, we analyze the second-order coherences and discuss the frequency spectrum of the radiation as well as the Fano factors. We include thermal effects in our discussion and calculate the driving strength for which the Cauchy-Schwarz inequality is most strongly violated. We calculate the FCS in Sec. V and discuss the FCS in two different limits in Sec. VI. In Sec. VII, we discuss the interpretation in terms of bursts and calculate the probability of large deviations from the equilibrium. We conclude in Sec. VIII.

II Setup

The studied setup is composed of a Josephson junction, biased by a dc-voltage source, that is in series with two microwave resonators. The detailed setup is illustrated in Fig. 1. Each resonator $k\in\{a,b\}$ is characterized by its resonance frequency $\Omega_{k}$ and an impedance $Z_{k}^{0}$ at low frequency. The environment is modeled as resistors at temperature $T$ with a coupling rate $\Gamma_{k}^{0}$. Additionally, each resonator is connected to a detector. The detectors are assumed to be well isolated such that they effectively operate at $T=0$. They count the radiation emitted by the resonators at a detector rate $\Gamma_{k}^{d}$. The photons dissipated by the resistors are not detected causing a finite detection efficiency $f_{k}$. In total, the emission of photons in each resonator occurs on the frequency scale $\Gamma_{k}=\Gamma_{k}^{d}+\Gamma_{k}^{0}$. Since we are primarily interested in correlated events in both resonators, we will normalize time and frequency to the natural scale $\Gamma=(\Gamma_{a}+\Gamma_{b})/2$.

Close to resonance, the impedance $Z_{k}$ of resonator $k$ is given by

[TABLE]

valid for $|\omega-\Omega_{k}|\ll\Omega_{k}$. Here, we introduced the abbreviation $\gamma_{k}=\Gamma_{k}/(\Gamma_{a}+\Gamma_{b})$.

We are interested in the non-degenerate case, where both resonators have different resonance frequencies, $\Omega_{a}\neq\Omega_{b}$. The superconducting phase difference across the junction depends on the superconducting phases $\phi_{k}$ across each resonator as well as the applied dc voltage $V$. For the following calculations, we require that both resonators have a high quality factor $\Omega_{k}\gg\gamma_{k}$. Additionally, we will assume that the impedance far from resonance $Z_{k}^{0}$ is small at the quantum scale, $G_{Q}Z_{k}^{0}\ll 1$, with $G_{Q}=e^{2}/\pi\hbar$ the conductance quantum. Then, the contribution to the phase fluctuations from low frequencies $(\delta\phi_{k})^{2}\simeq G_{Q}Z_{k}^{0}\ll 1$ can be safely neglected. In this context, we also want to neglect the broadening of the Josephson emission line by low-frequency phase noise. This broadening can be estimated as $k_{B}TG_{Q}(Z_{a}^{0}+Z_{b}^{0})$.likharev ; Arndt:18 The broadening can be neglected, if it is much smaller than the width of both resonators $\hbar\Gamma\gamma_{k}\gg k_{B}TG_{Q}(Z_{a}^{0}+Z_{b}^{0})$, which restricts the allowed temperature range of our calculation.

III Keldysh action

We describe the quantum mechanical system using the Keldysh path integral formalism. The partition function $\mathcal{Z}$ is given by the path integral over the Keldysh contour weighted by the factor $e^{iS}$, where $S$ corresponds to the quantum action

[TABLE]

Here, $\phi_{k}^{\pm}(t)$ refers to the superconducting phase difference across the resonators along the forward (backward) part of the Keldysh contour. We divide the action into the contribution from the junction and the resonators, $S=S_{J}+S_{a}+S_{b}$. The action of the Josephson junction is given by an integral of the phase-dependent Josephson energy $U(\phi)=-E_{J}\cos\phi$ and reads

[TABLE]

with $\phi_{V}(t)=2eVt/\hbar\Gamma$ the phase difference across the voltage source.

In frequency space, with $\phi_{k}(\omega)=\int dt\,e^{i\omega t}\phi_{k}(t)$, the action of resonator $k$ with impedance $Z_{k}(\omega)$ is given by

[TABLE]

with $n_{T}(\omega)=[\exp(\hbar\omega\Gamma/k_{B}T)-1]^{-1}$ the Bose-Einstein occupation at temperature $T$. Additionally, we introduced the factor $\alpha_{k}=\Gamma_{k}^{d}/\Gamma_{k}$, which originates from the current division between the detector and the resistor.

To observe the parametric resonance effect, we tune the Josephson frequency $\Omega_{J}$ close to the sum of the resonator frequencies by setting the dc-bias voltage to $V=\hbar\Gamma\Omega_{J}/2e$, with $\Omega_{J}=\Omega_{a}+\Omega_{b}+\Delta$. The detuning $\Delta$ is assumed to be small, $\Delta\ll\Omega_{k}$. Under this assumptioncommut , we can perform a rotating-wave approximation by introducing the slow, complex variables $\psi_{k}$

[TABLE]

Note, that the choice of how we distribute the detuning $\Delta$ between both rotations has no physical consequence and corresponds simply to a total frequency shift in our calculations.

By inserting Eq. (6) in Eq. (III) and neglecting all fast-oscillating terms, we obtain a new action for the slow variables $\psi_{k}(t)$

[TABLE]

with the Bessel functions $J_{m}(x)$. Similarly, the action of resonator $k$ can be expressed as

[TABLE]

with $n_{k}=n_{T}(\Omega_{k})(1-\alpha_{k})$ the average number of photons in the mode in thermal equilibrium.

At $T=0$, the ’classical’ equations of motion for the slow, complex variables $\psi_{k}$ correspond to the saddle-point solution of the total action $S$ for $n_{a}=n_{b}=0$ and read

[TABLE]

Here, $\bar{k}$ corresponds to the other resonator. Additionally, we introduced the dimensionless driving strength $\epsilon=2\pi E_{J}G_{Q}\sqrt{Z_{a}^{0}Z_{b}^{0}}/\hbar\Gamma\sqrt{\gamma_{a}\gamma_{b}}$ as well as the abbreviation $\sigma_{k}=(Z_{k}^{0}\gamma_{\bar{k}}/Z_{\bar{k}}^{0}\gamma_{k})^{1/2}$. By linearizing Eq. (III), we find that the classical equations have only nontrivial, nonzero solutions $\psi_{k}\neq 0$ if the Josephson energy exceeds the threshold $\epsilon\geq\epsilon_{\ast}=\sqrt{1+\Delta^{2}}$.Wustmann:17 Above this threshold, self-sustained parametric oscillations emerge, leading to coherent radiation from the two oscillators. The two radiation frequencies are shifted from the resonance frequencies due to the detuning $\Delta$. Just above the threshold, the radiation from resonator $k$ is detuned by $\gamma_{k}\Delta$. Therefore, the mode with the larger linewidth is influenced by the detuning more strongly. We will come back to this property later, when we discuss the frequency distribution of photons below the threshold.

To better understand the behavior of the system above the threshold, we want to calculate the equilibrium positions that correspond to the self-sustained oscillations. Using analogous steps to Ref. Armour:15, , we introduce amplitude-phase coordinates for both oscillators, $\psi_{a}=Ae^{i\xi_{a}}$ and $\psi_{b}=Be^{i\xi_{b}}$, and expand Eq. (III) to leading order in $\epsilon-\epsilon_{*}$. In addition to the trivial solution $A=B=0$, we obtain the fixed point positions

[TABLE]

The fixed point value for the relative phase $\xi_{a}-\xi_{b}$ is arbitrary, demonstrating that the phase-locking between the Josephson oscillations and the resonator modes only locks the total phase. This property is qualitatively different to the degenerate parametric oscillatorWustmann:13 , where the degeneracy of the phase is only two-fold and the phase can take two values differing by $\pi$.

While classical, coherent radiation is no longer possible below the threshold, quantum fluctuations enable the emission of photon pairs. Due to the low impedance environment, phase fluctuations typically remain small, $(\delta\psi_{k})^{2}\simeq Z_{k}^{0}G_{Q}\ll 1$. Upon approaching the threshold, the fluctuations increase. However, for sufficiently small impedances, the crossover region where the fluctuations become of the order of $1$ remains narrow and can be estimated as $|\epsilon_{*}-\epsilon|\simeq Z_{k}^{0}G_{Q}\gamma_{\bar{k}}\epsilon_{*}\ll\epsilon_{*}$. Therefore, we can expand the Josephson action in Eq. (7) up to leading order in $\psi_{k}$

[TABLE]

valid up to a narrow interval below the threshold.

Next, we want to introduce the counting fields $\chi_{k}(t)$ which will later provide access to the radiation statistics. The counting fields count only the photons emitted at the detectors, which are effectively at zero temperature. A common choice is a piecewise-constant counting field, with $\chi_{k}(t)=\chi_{k}$ during the time interval when the respective detector counts the emitted photons and $\chi_{k}(t)=0$ at all other times. We can include the counting fields in our calculation by modifying the second term in Eq. (9) intochi

[TABLE]

Here, we have introduced the counting efficiency $f_{k}=\alpha_{k}/(n_{k}+1)$ which is generally finite due to additional dissipation at the resistors at temperature $T$. With the modification in Eq. (14), $\mathcal{Z}(\chi_{k})$ represents the characteristic function of the probability distribution of detecting $N_{k}$ photons within the chosen time interval

[TABLE]

The total action of the linearized system can be written in compact form by introducing the vector $\Psi(\nu)=[\psi_{a}^{+}(\nu),\psi_{a}^{-}(\nu),\psi_{b}^{+}(-\nu)^{\ast},\psi_{b}^{-}(-\nu)^{\ast}]^{T}$. In this basis, the action reads

[TABLE]

Here, $\mkern 4.0mu\overline{\mkern-4.0muM\mkern-4.0mu}\mkern 4.0mu_{k}(\nu,\chi_{k})$ is given by Eq. (9) with the modification in Eq. (14) and

[TABLE]

IV Second-order coherences

In this section, we want to calculate the second-order coherences

[TABLE]

of the radiation. Here, we introduced the photon current $I_{k}(\tau)$ in the detector from resonator $k$. In our calculation, the Keldysh path integral formalism guarantees the contour-ordering of the operators due to the inherent time-ordering along the Keldysh contour.beenakker:01 ; hassler:15

Let us choose the time-dependent, piecewise constant counting fields as follows. During the time interval $(0,\Delta t)$, we count the photons emitted by resonator $k$ and therefore set $\chi_{k}(t)=\chi_{k}$. Within the time interval $(\tau,\tau+\Delta t)$, we count the photons emitted by resonator $l$ and set $\chi_{l}(t)=\chi_{l}$. At all other times, the counting fields are set to zero. We assume that the measurement time $\Delta t$ is small, such that the average number of photons emitted is small, $\Delta t\ll 1$. Additionally, we assume that the time distance between the measurement intervals is large, $\tau\gg\Delta t$.

We can calculate the number of photons $\Delta\langle N_{k}\rangle$ detected in the small time interval $\Delta t$ with the characteristic function via $\Delta\langle N_{k}\rangle=d\mathcal{Z}/d(i\chi_{k})|_{\chi_{k}=0}$. We obtain

[TABLE]

with $\bar{I}_{k}=\Delta\langle N_{k}\rangle/\Delta t$. Since the action is Gaussian, the calculation of the correlators is straightforward. To avoid cluttering the main article, the detailed derivation can be found in App. A. We insert the correlator in Eq. (20) and obtain

[TABLE]

It is instructive to express this intensity as an integral of the frequency spectrum of the radiation $\rho(\nu)$. At $T=0$, we obtain

[TABLE]

Here, we have introduced the parameter

[TABLE]

where $r$ ranges from $\gamma^{2}-\Delta^{2}$ at zero drive to $1-\gamma^{2}\Delta^{2}$ at the threshold. Additionally, we have introduced the asymmetry

[TABLE]

Figure 2 displays the emission spectrum as a function of frequency. At low driving strength, $\epsilon\ll\epsilon_{*}$, the resonators are most likely to emit photons either at their resonance frequency $\Omega_{k}$ or at the detuned frequency $\Omega_{k}+\Delta$. The detuning is not split between the resonators. Which resonator preferably takes the total detuning is determined by asymmetry. At zero asymmetry, both resonators are equally likely to emit photons at the detuned frequency. However, in the presence of asymmetry, the radiation from the mode with the larger linewidth is more likely to contain the total detuning. This is a distinct feature of the non-degenerate oscillator.

Close to the threshold, $\epsilon\approx\epsilon_{*}$, the detuning is split between both resonators such that the radiation from resonator $k$ is detuned by $\gamma_{k}\Delta$. This behavior coincides with the classical results discussed in the previous section.

Analogously to Eq. (20), we can calculate the second-order correlator

[TABLE]

By applying Wick’s theorem, we can calculate the correlator using the results from App. A. Close to threshold, $\epsilon\approx\epsilon_{*}$, the coherences $g_{aa}$, $g_{bb}$, and $g_{ab}$ all demonstrate the same asymptotic behavior. We obtain $g_{kl}^{(2)}(\tau)=1+\exp(-|\tau|/\tau_{\epsilon})$ with

[TABLE]

The expression suggests the emergence of a new long time scale $\tau_{\epsilon}$ in the vicinity of the threshold which diverges upon approaching the threshold. This behavior is similar to that of the degenerate oscillator Padurariu:12 . However, for the non-degenerate oscillator asymmetry can have a strong influence on the time scale.

In addition to the exponential behavior, the second-order coherences for $\tau=0$ are also of interest. Figure 3 displays the second-order coherences for $\tau=0$ as a function of driving strength. Analytically, we obtain

[TABLE]

Here, $\bar{\epsilon}$ is given by

[TABLE]

and corresponds to the drive at which $g_{ab}^{(2)}(0)$ reaches its maximum value of

[TABLE]

valid for small temperature such that $n_{a},n_{b}\ll 1$. Figure 4 displays the maximum value of $g_{ab}^{(2)}(0)$ as a function of asymmetry. If the occupation numbers of both resonators in thermal equilibrium coincide, asymmetry decreases coherence. However, since the resonance frequencies of the resonators differ, the occupation numbers in thermal equilibrium differ, too. In this case, asymmetry can increase the coherence. To achieve maximal coherence, the resonator with the larger resonance frequency and, thus, the lower photon occupation in thermal equilibrium should feature the larger linewidth.

Then, the Cauchy-Schwarz inequality

[TABLE]

is most strongly violated. The violation of the inequality is linked to nonclassical correlations between the radiation from the two resonators and demonstrates the quantum character of the radiation. At $T=0$, the Cauchy-Schwarz inequality is violated everywhere below the threshold. Due to thermal processes at finite temperature, the inequality is no longer violated for all driving strength. Instead, the driving strength needs to exceed $\epsilon\geq\bar{\epsilon}^{2}/\epsilon_{*}$ to generate nonclassical radiation. Furthermore, Eq. (30) indicated that for $n_{a}\gamma_{a}/\gamma_{b}+n_{b}\gamma_{b}/\gamma_{a}\geq 1$, the Cauchy-Schwarz inequality is always fulfilled for any driving strength.

The second-order coherences can also be used to calculate the Fano factor $F_{kk}$. The Fano factor is connected to the second-order coherence viahassler:15

[TABLE]

and corresponds to the number of correlated photons from either resonator. We obtain

[TABLE]

valid for low temperature, such that $n_{a},n_{b}\ll 1$. In the limit $\epsilon\rightarrow 0$, we obtain the Fano factor expected from thermal radiationBrange:19 . Upon approaching the threshold, the Fano factor diverges quadratically. Furthermore, we know from Eq. (20) that the number of photons present in either resonator diverges only linearly upon approaching the threshold. Thus, the number of correlated photons exceeds by far the number of photons present in either resonator. This behavior is connected to the long times scale $\tau_{\epsilon}$ in the vicinity of the threshold, see also Sec. VII and Ref. Padurariu:12, .

Additionally, we can calculate the Fano factor $F_{ab}$, given byfano

[TABLE]

which corresponds to the number of correlated photons across resonators. We obtain

[TABLE]

valid for $n_{a},n_{b}\ll 1$. Figure 5 displays both Fano factors as a function of driving strength.

Close to threshold, the behavior of $F_{ab}$ is similar to that of $F_{kk}$. In fact, for $f_{a}=f_{b}$, we find $F_{ab}=(F_{aa}+F_{bb})/2$. In contrast, the behavior of $F_{ab}$ at small drive, $\epsilon\ll\bar{\epsilon}$, differs from the behavior of $F_{kk}$ and depends strongly on thermal effects. As expected, at zero drive, there is no correlation between resonator $a$ and $b$. However, at small temperatures such that $n_{a},n_{b}\ll 1$, even a small drive can strongly increase $F_{ab}$. In the limit $T\rightarrow 0$, the number of correlated photons jumps from [math] at zero drive to $\sqrt{f_{a}f_{b}}$ at finite drive. Thus, at $T=0$ and finite drive, a photon emitted from resonator $k$ is always correlated with at least one photon emitted from resonator $\bar{k}$. However, the jump is very sensitive to thermal effects and easily washed out as displayed in Fig. 5.

V Full counting statistics

This section focuses on the calculation of the full counting statistics (FCS). Since we already discussed the influence of thermal radiation in the previous section, we will focus on the limit $T=0$. We count the photons emitted from either resonator in the measurement interval $(0,\tau)$. Analogously to the previous section, we set $\chi_{k}(t)=\chi_{k}$ during the measurement interval and $\chi_{k}(t)=0$ at all other times. Additionally, we introduce discrete frequencies spaced equally apart with $2\pi/\tau$. Since the action separates for each discrete frequency, evaluating the path integrals amounts to calculating the determinant of the matrix $A(\nu,\chi_{a},\chi_{b})$ at each frequency. We turn the resulting product of determinants into the exponent of a sum. In the long time limit, we can transform the sum back to a frequency integral and obtain

[TABLE]

Here, the cumulant-generating function is given by

[TABLE]

with

[TABLE]

and the frequency spectrum of the radiation $\rho(\nu)$ defined in Eq. (22). The structure of Eq. (42) suggests that photons are always emitted in pairs. This is consistent with the results from the previous section. A photon emitted from resonator $k$ is always correlated with at least one photon emitted from resonator $\bar{k}$. However, the detection of a single, seemingly uncorrelated photon is possible due to a finite detection efficiency.

Following the steps presented in App. B, the integral in Eq. (41) can be evaluated to obtain

[TABLE]

where $m$ is the solution of

[TABLE]

that coincides with $1$ for $s=s(\chi_{a},\chi_{b})=0$. This is one of the main results of this paper. The full analytic expression of $m$ is given in App. C.

We can use Eq. (44) to obtain the derivatives

[TABLE]

These derivatives are useful to obtain explicit expressions for the various cumulants. The average intensity $\bar{I}_{k}=\langle N_{k}\rangle/\tau$ is given by

[TABLE]

which coincides with our previous results in Eq. (21) at $T=0$. The intensity noise $S_{k}=\langle\langle N_{k}^{2}\rangle\rangle/\tau$ is given by

[TABLE]

Analogously, the cross-resonator noise $S_{ab}=\langle\langle N_{a}N_{b}\rangle\rangle/\tau$ is given by

[TABLE]

The Fano factors previously discussed in Sec. IV can also be obtained directly from these results with $F_{kk}=S_{k}/\bar{I}_{k}$ and $F_{ab}=S_{ab}/(\bar{I}_{a}\bar{I}_{b})^{1/2}$.

VI Limits

In this section, we briefly discuss the FCS in two different limits.

VI.1 No detuning or no asymmetry

In the cases where either $\Delta=0$ (no detuning) or $\gamma=0$ (no asymmetry), Eq. (44) simplifies and we find the simple explicit expression

[TABLE]

The resulting expression for the cumulant-generating function coincides with the results for the optical, non-degenerate oscillator in Ref. Vyas:92, for $\tau\to\infty$. Additionally, the result is of similar form as the FCS of the degenerate oscillator in Ref. Padurariu:12, as well as the optical result in Ref. Vyas:89, for $\tau\to\infty$. In this limit, the non-degenerate oscillator displays the same qualitative characteristics as the degenerate oscillator. This indicates that asymmetry can only have a qualitative effect on the statistics of the radiation, if the detuning is finite. We already discovered such an effect in Sec. IV, when we discussed how the detuning is distributed between both resonators.

VI.2 Close to threshold

In the vicinity of the threshold, we can obtain a general expression for the derivatives of $m$ using Eq. (44)

[TABLE]

The corresponding function $m$ is given by

[TABLE]

The result can in fact also be obtained directly by expanding Eq. (44) around $m=1$ with the result

[TABLE]

This expansion stems from the fact that, at the threshold, the cumulant-generating function becomes non-analytic for $s=0$ ($m=1$), see App. D. The non-analytic point of the cumulant-generating function is also of interest in the next section when we discuss the probability of big deviations.

VII Bursts and large deviations

In this section, we will assume a perfect counting efficiency, $f_{a}=f_{b}=1$. In this limit, Eq. (42) is reduced to $s(\chi)=e^{i\chi}-1$ with $\chi=\chi_{a}+\chi_{b}$. This shows that photons are always detected in pairs with $\bar{I}_{a}=\bar{I}_{b}=\bar{I}$. The form of $\lambda(\chi)$ in Eq. (41) allows for an expansion in $e^{i\chi}$. We obtain

[TABLE]

where the coefficients are given by

[TABLE]

In the vicinity of the threshold, we can calculate the asymptotic behavior of the coefficients analytically. In the limit $j\gg 1$, we can approximate Eq. (55) by a Gaussian integral and obtain

[TABLE]

where $\mu_{0}$ corresponds to the minimum of $\ln[1+1/\rho(\nu)]$ as a function of frequency $\nu$. In the vicinity of the threshold $\mu_{0}$ is given by

[TABLE]

A more general expression for $\mu_{0}$ is given in App. D.

The form of the cumulant-generating function in Eq. (55) suggests that photons are emitted in uncorrelated bursts of $j$ photon pairs. Since this behavior has been previously discussed in the context of the degenerate parametric oscillatorPadurariu:12 , we will not go into more details here.

Instead, we want calculate the probability of large deviations of the photon current $I$ from its average value $\bar{I}$. Using Eq. (56), we can calculate the corresponding cumulant generating function $\lambda(\chi)$ by approximating the sum in Eq. (54) by an integral. Then, we evaluate the probability in Eq. (15) at $N=I\tau$ using the saddle-point approximation. To exponential accuracy, we obtain

[TABLE]

with $\tau_{\epsilon}$ previously defined in Eq. (27). The form of this expression is identical to the degenerate oscillator casePadurariu:12 suggesting a universality of both systems that is connected to the diverging time scale $\tau_{\epsilon}$. Figure 6 displays the probability of large deviations at different driving strengths. Asymmetry increases the probability of large deviations.

VIII Conclusion

In view of recent advances in experimental setups that enable the measurement of correlations, we have studied the statistics of the radiation due to non-degenerate Josephson parametric resonance. We discussed the impact of detuning and asymmetry on the frequency spectrum of the radiation. We found that at low driving strength the detuning is not split between both resonators. Instead, the photons from the resonator with the larger linewidth contain the total detuning. With increasing driving strength, the detuning is split between the resonators resembling the classical behavior above the instability threshold.

Next, we investigated how thermal effects and asymmetry impact the second-order coherence as well as the Fano factor of the emitted radiation. In particular, we were able to obtain the driving strength for which the Cauchy-Schwarz inequality is most strongly violated. Furthermore, we found that the cross-coherence can be increased, if the resonator with the larger resonance frequency has the larger linewidth. Both properties add up to an optimal parameter regime for demonstrating non-classical correlations.

In addition, we calculated the full counting statistics of the radiation. We discovered that in the limit of no asymmetry or no detuning the full counting statistics of the non-degenerate oscillator can be mapped back to the degenerate oscillator. Our results are valid below the threshold where phase fluctuations are small. In the future, it would be interesting to investigate the radiation statistics in the narrow crossover region at the instability threshold where phase fluctuations increase.

Appendix A Derivation of correlators

For the calculations presented in the main text, we need the propagators

[TABLE]

where the last step follows from Gaussian integration of Eq. (16) with $(A^{\alpha\beta})^{-1}(\tau)$ defined via inverse Fourier transform

[TABLE]

Here, $(A^{\alpha\beta})^{-1}(\nu)$ is the inverse of the action matrix defined in Eq. (17) with $\chi_{a}=\chi_{b}=0$. The determinant of the action matrix is given by $p(\nu)$ in Eq. (23). In general, $p(\nu)$ has four complex roots $\nu_{1},\,\nu_{1}^{\ast},\,\nu_{2},\,\nu_{2}^{\ast}$ with

[TABLE]

In order to perform the inverse Fourier transform, we have to separate the advanced ($\tau>0$) and retarded ($\tau<0$) parts of the propagator. Then, we can perform the integral by closing the contour in the lower (upper) halve of the complex plane. For the advanced propagator $G^{A}$ we obtain

[TABLE]

with

[TABLE]

Analogously, we obtain for the retarded propagator $G^{R}$

[TABLE]

with

[TABLE]

Here, the $4\,\times\,4$ matrix $B(\nu)$ is given by

[TABLE]

where we use the abbreviations

[TABLE]

Appendix B FCS integral

We show how one can calculate integrals of the type

[TABLE]

with $g(\omega)$ a polynomial with real coefficients and $g(\omega)\propto\omega^{2n}$, $n\in\mathbb{N}$ for $\omega\to\infty$ such that the integral converges. In the following, we consider $\zeta\in\mathbb{R}$. However, the results can be analytically continued to complex $\zeta$.

The first step is to integrate by parts with the result

[TABLE]

We denote with $a_{1},\cdots,a_{n},a_{1}^{*},\cdots,a_{n}^{*}$ the (simple) roots of the polynomial $g(\omega)$ that are chosen such that $\operatorname{Im}a_{j}>0$. Similar, $b_{1}(\zeta),\cdots,b_{n}(\zeta),b_{1}^{*}(\zeta),\cdots,b_{n}^{*}(\zeta)$ denote the roots of the polynomial $g(\omega)-\zeta$ which coincide with $a_{j}$ for $\zeta=0$.

The integral (72) can then be evaluated with the residue theorem. In particular, we obtain

[TABLE]

For $n=2$, we can use the following property of a depressed quartic polynomial to calculate the sum: given a polynomial of the form

[TABLE]

the negative square of sum of the roots in the upper half plane, $m(\zeta)=-[b_{1}(\zeta)+b_{2}(\zeta)]^{2}$, is a root of the resolvent cubic

[TABLE]

With this property, Eq. (73) reduces to

[TABLE]

Appendix C Full solution of the cumulant-generating function

The cumulant-generating function is given by Eq. (43) with $m$ the solution of Eq. (44) that coincides with $1$ for $s=0$. We use the trigonometric method to calculate an analytic expression for the relevant solution of $m$ and obtain

[TABLE]

with $\tilde{s}=(1-\gamma^{2})^{2}\epsilon^{2}s$. However, in many situations, e.g. the calculation of the cumulants, the implicit solution of $m$ in Eq. (44) is considerably easier to use.

Appendix D Non-analyticity of cumulant-generating function

For the extreme statistics in Sec. VII, we need to find the minimum value $\mu_{0}>0$ of $\ln[1+1/\rho(\nu)]$ as a function of frequency $\nu$. We assume a perfect counting efficiency, $f_{a}=f_{b}=1$. Then, Eq. (42) is reduced to $s(\chi)=e^{i\chi}-1$ with $\chi=\chi_{a}+\chi_{b}$ and the cumulant-generating function in Eq. (41) is given by

[TABLE]

This cumulant-generating function has a singularity at $\chi=-i\mu_{0}$. Therefore, we can use our result for $\lambda(\chi)$ in Eq. (43) to learn more about $\mu_{0}=\ln(s_{0}+1)$. From Eq. (44), it is clear that $m$ can only coincide with [math] for $\gamma\Delta=0$. Therefore, the square root in Eq. (43) can only cause the non-analyticity in this special case. However, as the solution for $s_{0}$ needs to be continuous, it is sufficient to find the point of non-analyticity of the relevant solution $m$. We take the discriminant of Eq. (44) and obtain a cubic equation for $\tilde{s}=(1-\gamma^{2})^{2}\epsilon^{2}s$

[TABLE]

The relevant point of non-analyticity of the solution $m$ is given by the solution of Eq. (D) that coincides with [math] at the threshold where $r=1-\gamma^{2}\Delta^{2}$. In the limits previously discussed, we can express the result analytically.

D.1 No detuning or no asymmetry

For $\gamma\Delta=0$, the solution has a non-analytic point (second-order derivative is discontinuous) at $r=-1$. The solution is given by

[TABLE]

which qualitatively corresponds to the solution for the degenerate oscillatorPadurariu:12 . For $\gamma\Delta\neq 0$, the discontinuity of the second-order derivative is washed out. This is a qualitative difference between the non-degenerate and the degenerate oscillator.

D.2 Close to threshold

In the vicinity of the threshold, we obtain

[TABLE]

which corresponds to the $\mu_{0}$ obtained in Eq. (57).