Statistics of radiation due to non-degenerate Josephson parametric down conversion
Lisa Arndt Fabian Hassler




lisa.arndt@rwth-aachen.de
JARA Institute for Quantum Information, RWTH Aachen University, 52056 Aachen, Germany


JARA Institute for Quantum Information, RWTH Aachen University, 52056 Aachen, Germany






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. 




Statistics of radiation due to non-degenerate Josephson parametric down conversion
    Fabian Hassler
    July 9, 2019
==================================================================================





§ INTRODUCTION

Parametric down-conversion<cit.> is an important building block of many optical quantum information applications<cit.>. It describes the process that converts a pump photon of frequency Ω into a pair of photons with frequencies Ω_a and Ω_b which add up to the frequency of the pump photon, Ω_a+Ω_b=Ω. This operation can be implemented by nonlinear oscillators with resonance frequencies Ω_a and Ω_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 pairs<cit.> as well as squeezed states of light.<cit.>

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 experiments<cit.> 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 oscillators<cit.>. 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.<cit.> Motivated by these results, we want to investigate the radiation statistics of the non-degenerate setup illustrated in Fig. <ref>. The setup is composed of a Josephson junction, dc-biased with voltage V, that is in series with two microwave resonators with resonance frequencies Ω_a and Ω_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 oscillator<cit.>. 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. <ref>. In Sec. <ref>, we discuss its Keldysh action. Moreover, we introduce the counting fields which will later provide access to the radiation statistics. In Sec. <ref>, 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. <ref> and discuss the FCS in two different limits in Sec. <ref>. In Sec. <ref>, we discuss the interpretation in terms of bursts and calculate the probability of large deviations from the equilibrium. We conclude in Sec. <ref>.




§ 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. <ref>. Each resonator k ∈{a,b} is characterized by its resonance frequency Ω_k and an impedance Z_k^0 at low frequency. The environment is modeled as resistors at temperature T with a coupling rate Γ_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 Γ_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 Γ_k=Γ_k^d+Γ_k^0. Since we are primarily interested in correlated events in both resonators, we will normalize time and frequency to the natural scale Γ=(Γ_a+Γ_b)/2. 

Close to resonance, the impedance Z_k of resonator k is given by


    Z_k(ω)=Z_k^0Ω_k/γ_k-i(ω-Ω_k),


valid for |ω-Ω_k|≪Ω_k. Here, we introduced the abbreviation γ_k=Γ_k/(Γ_a+Γ_b).

We are interested in the non-degenerate case, where both resonators have different resonance frequencies, Ω_a≠Ω_b. The superconducting phase difference across the junction depends on the superconducting phases ϕ_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 Ω_k≫γ_k. Additionally, we will assume that the impedance far from resonance Z_k^0 is small at the quantum scale, G_QZ_k^0≪1, with G_Q=e^2/πħ the conductance quantum. Then, the contribution to the phase fluctuations from low frequencies (δϕ_k)^2≃ G_QZ_k^0≪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).<cit.> The broadening can be neglected, if it is much smaller than the width of both resonators ħΓγ_k≫ k_B T G_Q(Z_a^0+Z_b^0), which restricts the allowed temperature range of our calculation.



§ KELDYSH ACTION


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


    𝒵=∫𝒟[ϕ_a^+]𝒟[ϕ_a^-]𝒟[ϕ_b^+]𝒟[ϕ_b^-] e^iS[ϕ_a^+,ϕ_a^-,ϕ_b^+,ϕ_b^-].


Here, ϕ_k^±(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(ϕ)=-E_Jcosϕ and reads


    S_J=E_J/ħΓ∫dt{cos   [ϕ_V(t)+ϕ_a^+(t)+ϕ_b^+(t)]
       -cos[ϕ_V(t)+ϕ_a^-(t)+ϕ_b^-(t)]},


with ϕ_V(t)=2eVt/ħΓ the phase difference across the voltage source.

In frequency space, with ϕ_k(ω)=∫ dt  e^iω tϕ_k(t), the action of resonator k with impedance Z_k(ω) is given by


    S_k   =i/8π G_Q∫dω/2π∑_α,β=±ϕ_k^α(ω)^∗M_k^αβ(ω) ϕ_k^β(ω),
    
    M_k(ω)   =ω {Re[Z_k(ω)^-1][  0  1; -1  0 ]+iIm[Z_k(ω)^-1][  1  0;  0 -1 ]
         +[2n_T(ω)+1]
    Re[Z_k(ω)^-1] (1-α_k) [  1 -1; -1  1 ]
         +[2n_0(ω)+1]Re[Z_k(ω)^-1] α_k [  1 -1; -1  1 ]} ,


with n_T(ω)=[exp (ħωΓ/k_B T)-1]^-1 the Bose-Einstein occupation at temperature T. Additionally, we introduced the factor α_k=Γ_k^d/Γ_k, which originates from the current division between the detector and the resistor. 

To observe the parametric resonance effect, we tune the Josephson frequency Ω_J close to the sum of the resonator frequencies by setting the dc-bias voltage to V=ħΓΩ_J/2e, with Ω_J=Ω_a+Ω_b+Δ. The detuning Δ is assumed to be small, Δ≪Ω_k. Under this assumption<cit.>, we can perform a rotating-wave approximation by introducing the slow, complex variables ψ_k


    ϕ_k^±(t)=Re[e^-i(Ω_k+Δ/2) tψ_k^±(t)].


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

By inserting Eq. (<ref>) in Eq. (<ref>) and neglecting all fast-oscillating terms, we obtain a new action for the slow variables ψ_k(t)


    S_J=-E_J/2ħΓ∫dt { s_J[ψ_a^+(t),ψ_b^+(t)]-s_J[ψ_a^-(t),ψ_b^-(t)]},
       s_J(ψ_a,ψ_b)=J_1(|ψ_a|)J_1(|ψ_b|)/|ψ_a||ψ_b|(ψ_aψ_b+ψ_a^∗ψ_b^∗),


with the Bessel functions J_m(x). Similarly, the action of resonator k can be expressed as 


    S_k=i/16π G_Q∫   dν/2π∑_α,β=±ψ_k^α(ν)^∗M_k^αβ(ν) ψ_k^β(ν),
    
    Z_k^0M_k(ν)=-i   (ν+Δ/2) [  1  0;  0 -1 ]
           +γ_k
     [    2n_k+1     -2n_k; -2(n_k+1)    2n_k+1 ] ,


with n_k=n_T(Ω_k)(1-α_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 ψ_k correspond to the saddle-point solution of the total action S for n_a=n_b=0 and read


    dψ_k/dt=   ( iΔ/2
    -γ_k)ψ_k-4iϵγ_kσ_k[ψ_k̅^∗J_1(|ψ_a|)J_1(|ψ_b|)/|ψ_a||ψ_b|
         -ψ_kJ_1(|ψ_k̅|)J_2(|ψ_k|)/2|ψ_k̅||ψ_k|^2(ψ_aψ_b+ψ_a^∗ψ_b^∗)].


Here, k̅ corresponds to the other resonator. Additionally, we introduced the dimensionless driving strength ϵ=2π E_J G_Q √(Z_a^0Z_b^0)/ħΓ√(γ_aγ_b) as well as the abbreviation σ_k=(Z_k^0γ_k̅/Z_k̅^0γ_k)^1/2. By linearizing Eq. (<ref>), we find that the classical equations have only nontrivial, nonzero solutions ψ_k≠ 0 if the Josephson energy exceeds the threshold ϵ≥ϵ_∗=√(1+Δ^2).<cit.> 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 Δ. Just above the threshold, the radiation from resonator k is detuned by γ_kΔ. 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, ψ_a=Ae^iξ_a and ψ_b=Be^iξ_b, and expand Eq. (<ref>) to leading order in ϵ-ϵ_*. In addition to the trivial solution A=B=0, we obtain the fixed point positions


    A=√(8ϵ_*)√(ϵ-ϵ_*)σ_a/√(ϵ_*^2(σ_a^2+1)+2Δ^2(γ_b+σ_a^2γ_a)),   B=A/σ_a,
       e^i(ξ_a+ξ_b)=Δ-i/ϵ_*.


The fixed point value for the relative phase ξ_a-ξ_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 oscillator<cit.>, where the degeneracy of the phase is only two-fold and the phase can take two values differing by π.

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, (δψ_k)^2≃ Z_k^0G_Q≪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 |ϵ_*-ϵ|≃ Z_k^0G_Qγ_k̅ϵ_*≪ϵ_*. Therefore, we can expand the Josephson action in Eq. (<ref>) up to leading order in ψ_k 


    S_J=-E_J/8ħΓ∫dt [ψ_a^+(t)ψ_b^+(t)-ψ_a^-(t)ψ_b^-(t)+c.c.],


valid up to a narrow interval below the threshold.

Next, we want to introduce the counting fields χ_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 χ_k(t)=χ_k during the time interval when the respective detector counts the emitted photons and χ_k(t)=0 at all other times. We can include the counting fields in our calculation by modifying the second term in Eq. (<ref>) into<cit.>


    γ_k
     [                     2n_k+1                      -2n_k; -2(n_k+1)[1+f_k(e^iχ_k-1)]                     2n_k+1 ].


Here, we have introduced the counting efficiency f_k=α_k/(n_k+1) which is generally finite due to additional dissipation at the resistors at temperature T. With the modification in Eq. (<ref>), 𝒵(χ_k) represents the characteristic function of the probability distribution of detecting N_k photons within the chosen time interval


    P(N_k)=∫dχ_k/2π𝒵(χ_k)e^-iχ_k N_k.


The total action of the linearized system can be written in compact form by introducing the vector Ψ(ν)=[ψ_a^+(ν),ψ_a^-(ν),ψ_b^+(-ν)^∗,ψ_b^-(-ν)^∗]^T. In this basis, the action reads


    S=i/16π G_Q∫   dν/2π∑_α,βΨ^α(ν)^∗A^αβ(ν,χ_a,χ_b) Ψ^β(ν),
    
    A(ν,χ_a,χ_b)=   [    M_a(ν,χ_a)             J;             J M_b^T(-ν,χ_b) ].


Here, M_k(ν,χ_k) is given by Eq. (<ref>) with the modification in Eq. (<ref>) and 


    J=i2π E_J G_Q/ħΓ[  1  0;  0 -1 ].







§ SECOND-ORDER COHERENCES

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


    g^(2)_kl(τ)=⟨:I_l(τ)I_k(0):⟩/⟨ I_k⟩⟨ I_l⟩,   k,l∈{ a,b}


of the radiation. Here, we introduced the photon current I_k(τ) 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.<cit.> 

Let us choose the time-dependent, piecewise constant counting fields as follows. During the time interval (0,Δ t), we count the photons emitted by resonator k and therefore set χ_k(t)=χ_k. Within the time interval (τ,τ+Δ t), we count the photons emitted by resonator l and set χ_l(t)=χ_l. At all other times, the counting fields are set to zero. We assume that the measurement time Δ t is small, such that the average number of photons emitted is small, Δ t ≪1. Additionally, we assume that the time distance between the measurement intervals is large, τ≫Δ t. 

We can calculate the number of photons Δ⟨ N_k⟩ detected in the small time interval Δ t with the characteristic function via Δ⟨ N_k⟩=d 𝒵/d (iχ_k)|_χ_k=0. We obtain


    I̅_k= f_k γ_k(n_k+1)/8π G_Q Z_k^0⟨ψ_k^-∗ψ_k^+⟩ ,


with I̅_k=Δ⟨ N_k⟩/Δ 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. <ref>. We insert the correlator in Eq. (<ref>) and obtain


    I̅_k=2 f_kγ_k1+n_k/ϵ_*^2-ϵ^2{ n_kϵ_*^2
    +ϵ^2[(1+n_k̅)γ_k̅-n_kγ_k]}.


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


    I̅_k=     f_k∫dν/2π ρ(ν),    ρ(ν)=(1-γ^2)^2ϵ^2/4p(ν),
    
    p(ν)=     ν ^4 +1+r/2ν
      ^2 - νγΔ + (1-r)^2+4 γ ^2 Δ
     ^2/16.


Here, we have introduced the parameter


    r    = (1-γ^2) ϵ^2 + γ ^2 -Δ ^2,


where r ranges from γ^2-Δ^2 at zero drive to 1-γ^2Δ^2 at the threshold. Additionally, we have introduced the asymmetry 


    γ = γ_a-γ_b,     |γ|
     < 1.




Figure <ref> displays the emission spectrum as a function of frequency. At low driving strength, ϵ≪ϵ_*, the resonators are most likely to emit photons either at their resonance frequency Ω_k or at the detuned frequency Ω_k+Δ. 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, ϵ≈ϵ_*, the detuning is split between both resonators such that the radiation from resonator k is detuned by γ_kΔ. This behavior coincides with the classical results discussed in the previous section. 


 

Analogously to Eq. (<ref>), we can calculate the second-order correlator


    ⟨:I_l(τ)I_k(0):⟩   =1/(Δ t)^2d^2 𝒵/d (iχ_k)d (iχ_l)|_χ_k,χ_l=0
       ∝⟨ψ_k^- (0)^∗ψ_l^-(τ)^∗ψ_l^+(τ)ψ_k^+(0)⟩ .


By applying Wick's theorem, we can calculate the correlator using the results from App. <ref>. Close to threshold, ϵ≈ϵ_*, the coherences g_aa, g_bb, and g_ab all demonstrate the same asymptotic behavior. We obtain g_kl^(2)(τ)=1+exp (-|τ|/τ_ϵ) with


    τ_ϵ=1+γ^2Δ^2/ϵ_*(ϵ_*-ϵ)(1-γ^2).


The expression suggests the emergence of a new long time scale τ_ϵ in the vicinity of the threshold which diverges upon approaching the threshold. This behavior is similar to that of the degenerate oscillator <cit.>. 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 τ=0 are also of interest. Figure <ref> displays the second-order coherences for τ=0 as a function of driving strength. Analytically, we obtain


    g_aa^(2)(0)=g_bb^(2)(0)=2,
       g_ab^(2)(0)ϵ≪ϵ̅=1+ϵ^2/ϵ_*^2(1+n_a+n_b)^2γ_aγ_b/n_a n_b,
       g_ab^(2)(0)ϵ≈ϵ_*=2+2(ϵ_*-ϵ)/ϵ_*(1+n_a+n_b)(1-n_aγ_a/γ_b-n_bγ_b/γ_a).


Here, ϵ̅ is given by


    ϵ̅=ϵ_*{n_a n_b/[(1+n_a)γ_a-n_bγ_b][(1+n_b)γ_b-n_aγ_a]}^1/4


and corresponds to the drive at which g_ab^(2)(0) reaches its maximum value of


    g_ab^(2)(0)ϵ=ϵ̅=1+γ_aγ_b/(√(γ_a n_a)+√(γ_b n_b) )^2,


valid for small temperature such that n_a,n_b≪1. Figure <ref> 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


    √(g_aa^(2)(0)g_bb^(2)(0))≥ g_ab^(2)(0)


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 ϵ≥ϵ̅^2/ϵ_* to generate nonclassical radiation. Furthermore, Eq. (<ref>) indicated that for n_aγ_a/γ_b+n_bγ_b/γ_a≥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 via<cit.>


    F_kk=1+I̅_k∫ dτ[g_kk^(2)(τ)-1]


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


    F_kkϵ≪ϵ̅=    1+2 f_kn_k+2ϵ^2f_kγ_k̅/ϵ_*^4[n_kϵ_*^2
            +(1+2n_k+n_k̅)(4γ_k+ϵ_*^2γ_k̅-ϵ_*^2γ_k)],
    
    F_kkϵ≈ϵ_*=    f_k(1+γ^2Δ^2)/2(ϵ_*-ϵ)^2(1+2n_k+n_k̅),


valid for low temperature, such that n_a,n_b≪1. In the limit ϵ→ 0, we obtain the Fano factor expected from thermal radiation<cit.>. Upon approaching the threshold, the Fano factor diverges quadratically. Furthermore, we know from Eq. (<ref>) 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 τ_ϵ in the vicinity of the threshold, see also Sec. <ref> and Ref. Padurariu:12.

Additionally, we can calculate the Fano factor F_ab, given by<cit.>


    F_ab=√(I̅_aI̅_b)∫ dτ[g_ab^(2)(τ)-1],


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


    F_abϵ≪ϵ̅=    √(f_af_bγ_aγ_b/n_a n_b)ϵ^2/ϵ_*^2,
    
    F_abϵ≈ϵ_*=    √(f_af_b)(1+γ^2Δ^2)/2(ϵ_*-ϵ)^2(1+3/2n_a+3/2n_b),


valid for n_a,n_b≪1. Figure <ref> 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, ϵ≪ϵ̅, 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≪1, even a small drive can strongly increase F_ab. In the limit T→ 0, the number of correlated photons jumps from 0 at zero drive to √(f_af_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 k̅. However, the jump is very sensitive to thermal effects and easily washed out as displayed in Fig. <ref>.




§ 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,τ). Analogously to the previous section, we set χ_k(t)=χ_k during the measurement interval and χ_k(t)=0 at all other times. Additionally, we introduce discrete frequencies spaced equally apart with 2π/τ. Since the action separates for each discrete frequency, evaluating the path integrals amounts to calculating the determinant of the matrix A(ν,χ_a,χ_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


    𝒵(χ_a,χ_b)=e^λ (χ_a,χ_b).


Here, the cumulant-generating function is given by


    λ(χ_a,χ_b)   =- τ∫dν/2πln[ A(ν,χ_a,χ_b)/ A(ν,χ_a=χ_b=0)]
       =- τ∫dν/2πln[1 - ρ(ν)s(χ_a,χ_b)],


with


    s(χ_a,χ_b) =     (e^iχ_a -1)f_a(1-f_b)+(e^iχ_b -1)f_b(1-f_a) 
       + (e^i(χ_a+χ_b) -1)f_af_b


and the frequency spectrum of the radiation ρ(ν) defined in Eq. (<ref>). The structure of Eq. (<ref>) 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 k̅. However, the detection of a single, seemingly uncorrelated photon is possible due to a finite detection efficiency. 

Following the steps presented in App. <ref>, the integral in Eq. (<ref>) can be evaluated to obtain


    λ(χ_a,χ_b) =  τ(1-√(m)),


where m is the solution of


    (m-1) [ m (r-m) + γ^2 Δ^2 ] = (1-γ^2)^2 ϵ^2s m


that coincides with 1 for s=s(χ_a,χ_b)=0. This is one of the main results of this paper. 
The full analytic expression of m is given in App. <ref>. 

We can use Eq. (<ref>) to obtain the derivatives


    dm/ds|_s=0   = - (1-γ^2) ϵ^2/ϵ_*^2-ϵ ^2
    d^2m/ds^2|_s=0   = -
      2(1+γ ^2 Δ ^2) (1-γ^2) ϵ^4/ (ϵ_* ^2-ϵ ^2)^3.


These derivatives are useful to obtain explicit expressions for the various cumulants. The average intensity I̅_k =⟨ N_k⟩/τ is given by


    ⟨ N_k⟩/τ=1/τd λ/d (iχ_k)|_χ_a,χ_b=0= f_k ϵ^2 (1-γ^2)/2(ϵ_*^2-ϵ^2),


which coincides with our previous results in Eq. (<ref>) at T=0. The intensity noise S_k=⟨⟨ N_k^2⟩⟩/τ is given by 


    ⟨⟨ N_k^2⟩⟩/τ   =1/τd^2 λ/d (iχ_k)^2|_χ_a,χ_b=0= ⟨ N_k⟩/τ
                +  f_k^2ϵ^4(5-r+3γ^2Δ^2)(1-γ^2)/4(ϵ_*^2-ϵ^2)^3.


Analogously, the cross-resonator noise S_ab=⟨⟨ N_a N_b⟩⟩/τ is given by 


    ⟨⟨ N_a N_b⟩⟩/τ   =1/τd^2 λ/d (iχ_a)d (iχ_b)|_χ_a,χ_b=0= √(f_af_b⟨ N_a⟩⟨ N_b⟩)/τ
              +  f_a f_bϵ^4(5-r+3γ^2Δ^2)(1-γ^2)/4(ϵ_*^2-ϵ^2)^3.


The Fano factors previously discussed in Sec. <ref> can also be obtained directly from these results with F_kk=S_k/I̅_k and F_ab=S_ab/(I̅_aI̅_b)^1/2. 



§ LIMITS

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


 §.§ No detuning or no asymmetry

In the cases where either Δ =0 (no detuning) or γ=0 (no asymmetry), Eq. (<ref>) simplifies and we find the simple explicit expression


    m = 1/2[1 + r + √((1-r)^2
     -4 (1-γ^2)^2 ϵ^2 s)].


The resulting expression for the cumulant-generating function coincides with the results for the optical, non-degenerate oscillator in Ref. Vyas:92 for τ→∞. 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 τ→∞. 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. <ref>, when we discussed how the detuning is distributed between both resonators.



 §.§ Close to threshold

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


    d^nm/ds^n|_s=0= - (2n)!/2 (2n-1) n!(1-γ^2)(1+γ^2Δ^2)^n-1ϵ_*^2n/ (ϵ_*
     ^2-ϵ ^2)^2n -1.


The corresponding function m is given by


    m=1   + (1-γ^2)/2(1+γ^2 Δ^2)(ϵ^2 - ϵ_*^2 
             +√((ϵ_*^2-ϵ^2)^2 -4 (1+γ^2 Δ^2)ϵ_*^2 s)).


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

    (m-1) [ r+2γ^2 Δ^2 -m (1+ γ^2 Δ^2) ] = (1-γ^2)^2 ϵ_*^2
     s.


This expansion stems from the fact that, at the threshold, the cumulant-generating function becomes non-analytic for s=0 (m=1), see App. <ref>. 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.



§ BURSTS AND LARGE DEVIATIONS


In this section, we will assume a perfect counting efficiency, f_a=f_b=1. In this limit, Eq. (<ref>) is reduced to s(χ)=e^iχ -1 with χ=χ_a+χ_b. This shows that photons are always detected in pairs with I̅_a=I̅_b=I̅. The form of λ (χ) in Eq. (<ref>) allows for an expansion in e^iχ. We obtain


    λ (χ)=∑_j=1^∞λ_j( e^ijχ-1),


where the coefficients are given by


    λ_j=τ/j∫dν/2π[ ρ(ν)/1+ρ(ν)]^j.


In the vicinity of the threshold, we can calculate the asymptotic behavior of the coefficients analytically.
In the limit j≫1, we can approximate Eq. (<ref>) by a Gaussian integral and obtain


    λ_j≫ 1≃τ(1-γ^2)ϵ_*/4j^3/2√(π(1+γ^2Δ^2))e^-jμ_0,


where μ_0 corresponds to the minimum of ln[1+1/ρ(ν)] as a function of frequency ν. In the vicinity of the threshold μ_0 is given by


    μ_0=(ϵ_*-ϵ)^2/1+γ^2Δ^2=(1-γ^2)^2ϵ_*^2/16I̅^2(1+γ^2Δ^2).


A more general expression for μ_0 is given in App. <ref>. 

The form of the cumulant-generating function in Eq. (<ref>) 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 oscillator<cit.>, 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 I̅. Using Eq. (<ref>), we can calculate the corresponding cumulant generating function λ(χ) by approximating the sum in Eq. (<ref>) by an integral. Then, we evaluate the probability in Eq. (<ref>) at N=Iτ using the saddle-point approximation. To exponential accuracy, we obtain


    P(I)   =∫dχ/2πexp [-iχ Iτ+λ(χ)]
       ∝exp[-τ/4τ_ϵ(I/I̅+I̅/I-2)],


with τ_ϵ previously defined in Eq. (<ref>). The form of this expression is identical to the degenerate oscillator case<cit.> suggesting a universality of both systems that is connected to the diverging time scale τ_ϵ. Figure <ref> displays the probability of large deviations at different driving strengths. Asymmetry increases the probability of large deviations.





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



§ DERIVATION OF CORRELATORS

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


    G_αβ(τ)=⟨Ψ^α(τ)Ψ^β(0)^∗⟩=16π G_Q(A^αβ)^-1(τ),


where the last step follows from Gaussian integration of Eq. (<ref>) with (A^αβ)^-1(τ) defined via inverse Fourier transform


    (A^αβ)^-1(τ)=∫dν/2πe^-iντ(A^αβ)^-1(ν).


Here, (A^αβ)^-1(ν) is the inverse of the action matrix defined in Eq. (<ref>) with χ_a=χ_b=0. The determinant of the action matrix is given by p(ν) in Eq. (<ref>). In general, p(ν) has four complex roots ν_1, ν_1^∗, ν_2, ν_2^∗ with


    ν_1=1/2(i-√(-r-2iγΔ)),
       ν_2=1/2(i+√(-r-2iγΔ)).



In order to perform the inverse Fourier transform, we have to separate the advanced (τ>0) and retarded (τ<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


    G^A(τ)=2π G_Q[ G_1^A e^-iν_1^*τ+G_2^Ae^-iν_2^*τ],


with 


    G_1^A=B(ν_1^*)/ir+2γΔ+√(-r+2iγΔ)(1-iγΔ),
       G_2^A=B(ν_2^*)/ir+2γΔ-√(-r+2iγΔ)(1-iγΔ).


Analogously, we obtain for the retarded propagator G^R 


    G^R(τ)=2π G_Q[ G_1^R e^-iν_1τ+G_2^R e^-iν_2τ],


with 


    G_1^R=B(ν_1)/ir-2γΔ-√(-r-2iγΔ)(1+iγΔ),
       G_2^R=B(ν_2)/ir-2γΔ+√(-r-2iγΔ)(1+iγΔ).


Here, the 4 × 4 matrix B(ν) is given by



    B(ν)=    n_a[      f_21^a(ν,γ)      f_14(-ν,-γ); -f_14(-ν^*,-γ)^*       f_12^b(-γ) ]⊗[ 1 1; 1 1 ]+n_b[      f_12^a(γ)      f_14(ν,γ); -f_14(ν^*,γ)^*  f_21^b(-ν,-γ) ]⊗[ 1 1; 1 1 ]
       +[        f_11^a(ν,γ)          f_12^a(γ)          f_13(ν,γ)          f_14(ν,γ);        f_21^a(ν,γ)   -f_11^a(ν^*,γ)^*        f_14(-ν,-γ)          f_24(ν,γ);     -f_24(ν^*,γ)^*     -f_14(ν^*,γ)^*      f_11^b(-ν,-γ)      f_21^b(-ν,-γ);   -f_14(-ν^*,-γ)^*     -f_13(ν^*,γ)^*         f_12^b(-γ) -f_11^b(-ν^*,-γ)^* ],


where we use the abbreviations

    g(ν,γ)=i-iγ+Δ-2ν,
       h(ν,γ)=(i-Δ)^2-ϵ^2+γ^2(1+ϵ^2)+4iγν-4ν^2,
       f_11^k(ν,γ)=Z_k^0 g(ν,γ)h(ν,γ),
       f_12^k(γ)=-2iZ_k^0(1-γ)^2(1+γ) ϵ^2,
       f_21^k(ν,γ)=-2iZ_k^0(1+γ)g(ν,γ)g(ν^*,γ)^*,
       f_13(ν,γ)=√(Z_a^0Z_b^0)√(1-γ^2) h(ν,γ)ϵ,
       f_14(ν,γ)=-2i√(Z_a^0Z_b^0)√(1-γ^2)(1-γ) g(-ν^*,-γ)^*ϵ,
       f_24(ν,γ)=-√(Z_a^0Z_b^0)√(1-γ^2) [h(ν^*,γ)^*+4(1-γ^2)]ϵ.




§ FCS INTEGRAL

We show how one can calculate integrals of the type


    I(ζ) = -∫dω/2πln(1-ζ/g(ω))


with g(ω) a polynomial with real coefficients and g(ω) ∝ω^2n, n∈ℕ for ω→∞ such that the integral converges. In the following, we consider ζ∈ℝ. However, the results can be analytically continued to complex
ζ.

The first step is to integrate by parts with the result


    I(ζ) = ∫dω/2πζω g'(ω)/g(ω)
     [g(ω ) - ζ].


We denote with a_1, ⋯, a_n, a_1^*, ⋯, a_n^* the (simple) roots of
the polynomial g(ω) that are chosen such that Im a_j
>0. Similar, b_1(ζ), ⋯, b_n(ζ), b_1^*(ζ), ⋯,
b_n^*(ζ) denote the roots of the polynomial g(ω) -ζ which coincide
with a_j for ζ =0.

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


    I(ζ) = i∑_j=1^n b_j(ζ) - i ∑_j=1^n a_j  .


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


    g(ω) = ω^4 + p ω^2 + q ω +r(ζ),


the negative square of sum of the roots in the upper half plane, m(ζ) =
-[b_1(ζ)+b_2(ζ)]^2, is a root of the resolvent cubic


    m^3 - 2 p m^2 +[p^2 -4 r(ζ)] m +q^2 .


With this property, Eq. (<ref>) reduces to


    I(ζ) = √(m(ζ=0)) -√(m(ζ)) .




§ FULL SOLUTION OF THE CUMULANT-GENERATING FUNCTION

The cumulant-generating function is given by Eq. (<ref>) with m the solution of Eq. (<ref>) 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



    m=1/3{ 1+r+2√(1-r+r^2-3s̃+3γ^2Δ^2)cos[1/3arccos((1+r)(2-5r+2r^2-9s̃)+9γ^2Δ^2(r-2)/2(1-r+r^2-3s̃+3γ^2Δ^2)^3/2)] },



with s̃=(1-γ^2)^2 ϵ^2 s. However, in many situations, e.g. the calculation of the cumulants, the implicit solution of m in Eq. (<ref>) is considerably easier to use.



§ NON-ANALYTICITY OF CUMULANT-GENERATING FUNCTION


For the extreme statistics in Sec. <ref>, we need to find the minimum value μ_0>0 of ln[1+1/ρ(ν)] as a function of frequency ν. We assume a perfect counting efficiency, f_a=f_b=1. Then, Eq. (<ref>) is reduced to s(χ)=e^iχ -1 with χ=χ_a+χ_b and the cumulant-generating function in Eq. (<ref>) is given by


    λ(χ)=- τ∫dν/2πln[1 - ρ(ν)(e^iχ -1)].


This cumulant-generating function has a singularity at χ=-iμ_0. Therefore, we can use our result for λ(χ) in Eq. (<ref>) to learn more about μ_0=ln (s_0+1). From Eq. (<ref>), it is clear that m can only coincide with 0 for γΔ =0. Therefore, the square root in Eq. (<ref>) 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. (<ref>) and obtain a cubic equation for s̃=(1-γ^2)^2 ϵ^2 s


    (r-1+γ^2Δ^2)(r^2+4γ^2Δ^2)
    =4s̃^3-s̃^2(1-10r+r^2+12γ^2Δ^2)
       -2s̃[r^3-4r^2+r-γ^2Δ^2(r^2-r+10+6γ^2Δ^2)].


The relevant point of non-analyticity of the solution m is given by the solution of Eq. (<ref>) that coincides with 0 at the threshold where r=1-γ^2Δ^2. In the limits previously discussed, we can express the result analytically.


 §.§ No detuning or no asymmetry

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


    μ_0= ln[1-r/(1-γ^2)^2 ϵ^2]    if    r < -1 
    ln[1+(r-1)^2/4(1-γ^2)^2 ϵ^2]    if    r > -1,


which qualitatively corresponds to the solution for the degenerate oscillator<cit.>. For γΔ≠ 0, the discontinuity of the second-order derivative is washed out. This is a qualitative difference between the non-degenerate and the degenerate oscillator.



 §.§ Close to threshold

In the vicinity of the threshold, we obtain


    μ_0=(ϵ_*-ϵ)^2/1+γ^2Δ^2,


which corresponds to the μ_0 obtained in Eq. (<ref>).

10


guckenheimer
J. Guckenheimer and P. Holmes,
 Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector
  Fields
 (Springer-Verlag, New York, 1983).

Gisin:02
N. Gisin, G. Ribordy, W. Tittel, and H. Zbinden,
 http://dx.doi.org/10.1103/RevModPhys.74.145
 Rev. Mod. Phys. 74, 145 (2002).

Scarani:05
V. Scarani, S. Iblisdir, N. Gisin, and A. Acín,
 http://dx.doi.org/10.1103/RevModPhys.77.1225
 Rev. Mod. Phys. 77, 1225 (2005).

Ou:92
Z. Y. Ou, S. F. Pereira, H. J. Kimble, and K. C. Peng,
 http://dx.doi.org/10.1103/PhysRevLett.68.3663
 Phys. Rev. Lett. 68, 3663 (1992).

Wu:86
L.-A. Wu, H. J. Kimble, J. L. Hall, and H. Wu,
 http://dx.doi.org/10.1103/PhysRevLett.57.2520
 Phys. Rev. Lett. 57, 2520 (1986).

Astafiev:07
O. Astafiev, K. Inomata, A. O. Niskanen, T. Yamamoto, Y. A. Pashkin, Y.
  Nakamura, and J. S. Tsai,
 Nature (London) 449, 588 (2007).

Groot:10
P. C. de Groot, J. Lisenfeld, R. N. Schouten, S. Ashhab, A. Lupaşcu, C. J.
  P. M. Harmans, and J. E. Mooij,
 Nature Physics 6, 763 EP  (2010).

hofheinz:11
M. Hofheinz, F. Portier, Q. Baudouin, P. Joyez, D. Vion, P. Bertet, P. Roche,
  and D. Esteve,
 http://dx.doi.org/10.1103/PhysRevLett.106.217005
 Phys. Rev. Lett. 106, 217005 (2011).

Hoi:11
I.-C. Hoi, C. M. Wilson, G. Johansson, T. Palomaki, B. Peropadre, and P.
  Delsing,
 http://dx.doi.org/10.1103/PhysRevLett.107.073601
 Phys. Rev. Lett. 107, 073601 (2011).

Liu:14
Y.-Y. Liu, K. D. Petersson, J. Stehlik, J. M. Taylor, and J. R. Petta,
 http://dx.doi.org/10.1103/PhysRevLett.113.036801
 Phys. Rev. Lett. 113, 036801 (2014).

Chen:14
F. Chen, J. Li, A. D. Armour, E. Brahimi, J. Stettenheim, A. J. Sirois, R. W.
  Simmonds, M. P. Blencowe, and A. J. Rimberg,
 http://dx.doi.org/10.1103/PhysRevB.90.020506
 Phys. Rev. B 90, 020506 (2014).

Drummond:90
P. D. Drummond and M. D. Reid,
 http://dx.doi.org/10.1103/PhysRevA.41.3930
 Phys. Rev. A 41, 3930 (1990).

Padurariu:12
C. Padurariu, F. Hassler, and Y. V. Nazarov,
 http://dx.doi.org/10.1103/PhysRevB.86.054514
 Phys. Rev. B 86, 054514 (2012).

Leppakangas:13
J. Leppäkangas, G. Johansson, M. Marthaler, and M. Fogelström,
 http://dx.doi.org/10.1103/PhysRevLett.110.267004
 Phys. Rev. Lett. 110, 267004 (2013).

Armour:13
A. D. Armour, M. P. Blencowe, E. Brahimi, and A. J. Rimberg,
 http://dx.doi.org/10.1103/PhysRevLett.111.247001
 Phys. Rev. Lett. 111, 247001 (2013).

Gramich:13
V. Gramich, B. Kubala, S. Rohrer, and J. Ankerhold,
 http://dx.doi.org/10.1103/PhysRevLett.111.247002
 Phys. Rev. Lett. 111, 247002 (2013).

Souquet:14
J. R. Souquet, M. J. Woolley, J. Gabelli, P. Simon, and A. A. Clerk,
 Nature Communications 5, 5562 (2014).

Jin:15
J. Jin, M. Marthaler, and G. Schön,
 http://dx.doi.org/10.1103/PhysRevB.91.085421
 Phys. Rev. B 91, 085421 (2015).

Armour:15
A. D. Armour, B. Kubala, and J. Ankerhold,
 http://dx.doi.org/10.1103/PhysRevB.91.184508
 Phys. Rev. B 91, 184508 (2015).

Kubala:15
B. Kubala, V. Gramich, and J. Ankerhold,
 http://dx.doi.org/10.1088/0031-8949/2015/t165/014029
 Physica Scripta T165, 014029 (2015).

Leppakangas:16
J. Leppäkangas, M. Fogelström, M. Marthaler, and G. Johansson,
 http://dx.doi.org/10.1103/PhysRevB.93.014506
 Phys. Rev. B 93, 014506 (2016).

Wustmann:17
W. Wustmann and V. Shumeiko,
 http://dx.doi.org/10.1103/PhysRevApplied.8.024018
 Phys. Rev. Applied 8, 024018 (2017).

Westig:17
M. Westig, B. Kubala, O. Parlavecchio, Y. Mukharsky, C. Altimiras, P. Joyez, D.
  Vion, P. Roche, D. Esteve, M. Hofheinz, M. Trif, P. Simon, J. Ankerhold, and
  F. Portier,
 http://dx.doi.org/10.1103/PhysRevLett.119.137001
 Phys. Rev. Lett. 119, 137001 (2017).

Vyas:92
R. Vyas,
 http://dx.doi.org/10.1103/PhysRevA.46.395
 Phys. Rev. A 46, 395 (1992).

likharev
K. K. Likharev,
 Dynamics of Josephson Junctions and Circuits
 (Gordon and Breach Science Publishers, 1986).

Arndt:18
L. Arndt, A. Roy, and F. Hassler,
 http://dx.doi.org/10.1103/PhysRevB.98.014525
 Phys. Rev. B 98, 014525 (2018).

commut
To avoid mutual excitations between both resonators, the ratio
  Ω_a/Ω_b should be hard to approximate by the ratio of two small
  integer numbers. This is provided, e.g., by the golden ratio.

Wustmann:13
W. Wustmann and V. Shumeiko,
 http://dx.doi.org/10.1103/PhysRevB.87.184501
 Phys. Rev. B 87, 184501 (2013).

chi
Note that this modification is done in time space with the piecewise-constant
  counting fields χ_k(t).

beenakker:01
C. W. J. Beenakker and H. Schomerus,
 http://dx.doi.org/10.1103/PhysRevLett.86.700
 Phys. Rev. Lett. 86, 700 (2001).

hassler:15
F. Hassler and D. Otten,
 http://dx.doi.org/10.1103/PhysRevB.92.195417
 Phys. Rev. B 92, 195417 (2015).

Brange:19
F. Brange, P. Menczel, and C. Flindt,
 http://dx.doi.org/10.1103/PhysRevB.99.085418
 Phys. Rev. B 99, 085418 (2019).

fano
Note that the summand 1 is missing in the definition of F_ab compared to
  the definition of F_kk since there is no 'self'-correlation.

Vyas:89
R. Vyas and S. Singh,
 http://dx.doi.org/10.1103/PhysRevA.40.5147
 Phys. Rev. A 40, 5147 (1989).



