Input-output description of microwave radiation in the dynamical Coulomb blockade

We study microwave radiation emitted by a small voltage-biased Josephson junction connected to a superconducting transmission line. An input-output formalism for the radiation field is established, using a perturbation expansion in the junction's critical current. Using output field operators solved up to the second order, we estimate the spectral density and the second-order coherence of the emitted field. For typical transmission line impedances and at frequencies below the main emission peak at the Josephson frequency, radiation occurs predominantly due to two-photon emission. This emission is characterized by a high degree of photon bunching if detected symmetrically around half of the Josephson frequency. Strong phase fluctuations in the transmission line make related nonclassical phase-dependent amplitude correlations short lived, and there is no steady-state two-mode squeezing. However, the radiation is shown to violate the classical Cauchy-Schwarz inequality of intensity cross-correlations, demonstrating the nonclassicality of the photon pair production in this region.


Introduction
In a resistive environment, charge tunneling across a voltage-biased Josephson junction (JJ) triggers simultaneous microwave emission (1,2) , that carries away all or some of the gained electrostatic energy. For voltages below the superconducting gap, eV < ∆, the radiation is purely due to Cooper-pair tunneling since charge tunneling via excitation of quasiparticles is energetically forbidden. This mechanism leads to spectroscopically sharp features, which can be used for probing of environmental energy levels (3,4,5) , as photon absorption by the environmental modes influences the simultaneously measurable dc-current. A novel idea is to use the voltage-biased small JJ as a source of nonclassical microwave radiation, i.e. to convert the applied dc-voltage into correlated microwave photons (6,7) . This has recently stimulated theoretical studies of the emitted microwave field (7,8,9,10) .
In this article, we investigate microwave radiation produced in a dc-voltage-biased superconducting transmission line that is terminated by a small JJ. We establish an input-output formalism for the field operators in the transmission line. In this formalism, the electric current across the Josephson junction acts as a nonlinear and time-dependent boundary condition for the microwave field (11,12) . We solve this perturbatively as a power series in the junction's critical current (7) , and give explicit expressions for the output field operators up to the second order in the critical current. Assuming thermal equilibrium of the input field we recover the limit of incoherent Cooper-pair tunneling (1,2) , where microwave emission is due to an incoherent sequence of Cooper-pair tunneling events. The photon flux has been studied recently experimentally in this regime, (6) and it was shown that the radiation at voltages below the Josephson frequency, for typical tranmission-line impedances, is due to simultaneous two-photon emission. Using the formalism established in this article, we also study the nonclassicality of the photon pair production occurring in this region.
Using field operators up to second order in the junction's critical current, we derive analytical expressions for the first and and second-order (photon) coherences for typical transmission lines. We reproduce results for the photon-flux density and simultaneous electric current, previously derived using the P (E)-theory (2,6) . Further, the emission characteristics below the Josephson frequency is shown to be highly bunched, if detected symmetrically around half the Josephson frequency, eV /h. We then further study the nonclassicality of the photon pair production occuring below the Josephson frequency (7) . Strong phase fluctuations in the transmission line make phase-dependent nonclassical amplitude correlations short lived, and lead to a vanishing two-mode squeezing, in the steady-state. We thus proceed to prove the nonclassicality of the radiation in a different way, considering the classical Cauchy-Schwarz inequality of intensity cross-correlations, a nonclassicality test that is not affected by dephasing. Using the developed methods, we derive an equivalent inequality but expressed in terms of P (E)-functions. This is used to show that the emitted photons below the Josephson frequency violate the inequality, demonstrating the nonclassicality of the photon pair production in this region. The article has the following structure. In Section 2 we introduce the model we use to describe the radiation in the transmission line, and the nonlinear and timedependent boundary condition created by the Josephson junction. In section 3, we derive the solution by establishing a perturbation expansion in the junction's critical current. In section 4, we derive results for the microwave spectral density in the used leading-order approximation, whose validity is also addressed in this section. Higherorder coherences and the nonclassicality of the output radiation are investigated in Section 5. The technical details of the calculations are given in Appendices A-D.

The system and the model
Our system consists of a dc-voltage-biased transmission line terminated by a small Josephson junction, figure 1(a). We consider explicitly two types of environments: (i) a semi-infinite transmission line (i.e. Z 0 = Z 1 ) and (ii) a semi-infinite transmission line with a λ/4 cavity (i.e. Z 1 > Z 0 ). Case (i) allows for analytical solutions, while case (ii) enhances the output radiation at the cavity resonances, which is important in experiments.

Heisenberg equations of motion and quantization of the EM field
Following (12) , we model the system using a discretized circuit model, depicted in figure 1(b). The total Lagrangian of the system can be decomposed as The cavity (0 < x < d) and the free space (x > d) Lagrangians are respectively (11) Input-output description of microwave radiation in the dynamical Coulomb blockade 4 Here Φ m (t) is the magnetic flux of node m, see figure 1(b). The Josephson junction is described by the term Here C J is the junction's capacitance, E J is the Josephson coupling energy, Φ 0 = h/2e is the flux quantum, and the dc-voltage bias results in the term Φ V = V t.
We consider the Heisenberg equations of motion in the continuum limit δx → 0. Inside each region, free space (i = 0) or cavity (i = 1), one obtains a Klein-Gordon Here Φ(x, t) is the position dependent magnetic flux. We write down a solution for the cavity region as (0 < x < d) Here Z 1 = L 1 /C 1 is the corresponding characteristic impedance and k c ω = ω √ C 1 L 1 the wave number. Similarly, we write the free-space solution as (x > d) The in-field creation operators of photons, a † (ω), and the annihilation operators, a(ω), satisfy the commutation relation (13) As a consistency check of our theory, we will show that (8) is satisfied also for the out field operators.

Boundary conditions for the EM field
The boundary conditions appear at the Josephson junction (x = 0) and at the possible discrete change of the transmission-line parameters (x = d). Generally, we will have three boundary conditions to solve, and three unknown fields [a c in (ω), a c out (ω), and a f out (ω)]. The requirements of a continuous voltage distribution and current conservation across x = d imply the linear conditions, These can be solved by Fourier transformation.
Input-output description of microwave radiation in the dynamical Coulomb blockade 5 The main challenge is to solve the boundary condition at the junction, where current conservation gives a nonlinear and time-dependent condition Here, I c = (2π/Φ 0 )E J is the junction's critical current and ω J = 2eV /h is the Josephson frequency. Classically, for Z 0 = Z 1 and T = 0 (no input), this is equivalent to RCJmodel of a Josephson junction (14) .

Perturbative input-output approach
We want to derive a solution for the outgoing free-space field, a f out (ω), as a function of the incoming field in the same region, a f in (ω). The challenge is that the boundary condition at the junction is highly nonlinear, containing all moments of the field operators a c in/out (ω) and a c in/out (ω) † . A linearization of this condition, performed in Appendix A, captures some of the essential physics, but is limited to the second order in the photonic processes and does not, for example, describe correctly the effect of lowfrequency phase fluctuations. In the case of a weakly damped resonator (Z 1 ≫ Z 0 ), a single cavity mode can also be picked out and be described as a damped oscillator (9,10) , but with a limited description of the important low-frequency modes. Here, we use a different approach and derive a solution for the continuous-mode output field operators as a power series in the junction's critical current I c .

Unperturbed solution
The starting point is the solution for I c = 0, i.e. when Cooper-pair tunneling is neglected. By a Fourier transformation (Appendix B) we obtain the linear dependence . Similarly, we can solve for the cavity out-field a c out (ω) as a function of the free-space input a f in (ω), where A(ω) gives the response of the cavity to external drive and possesses information of its resonance frequencies. With the help of this, the operator for the phase difference at the junction, φ 0 (t) ≡ 2πΦ(0, t)/Φ 0 , can be written Here,Ā(ω) = A(ω)/C * (ω) and we also note the useful relation R(ω) =Ā(ω)/Ā * (ω). The corresponding phase fluctuations are equivalent to that of the "tunneling" impedance Z t (ω) (2) , defined as where R Q = h/4e 2 is the (superconducting) resistance quantum. In the ensemble average we assume thermal equilibrium for the incoming free-space modes, used throughout this article. For the open line (Z 0 = Z 1 ) case, the impedance (15) describes a capacitor C J and a resistor Z 0 in parallel, whereas for the cavity case (Z 1 ≫ Z 0 ) it describes a capacitively shunted λ/4-resonator, with resonances approximately at ω k = (2k + 1)ω 0 , where ω 0 = π/2d √ C 1 L 1 and k ∈ [0, 1, 2, . . .].

First and second order solutions
To seek for a solution, that includes Cooper-pair tunneling (I c = 0), we multiply the right-hand side of boundary condition (11) by a formal dimensionless parameter ξ, and correspondingly write the solution for the annihilation operators of the outgoing field in the open space as a f out (ω) = ∞ n=0 ξ n a n (ω). The zeroth-order solution, a 0 (ω), corresponds to I c = 0 and was obtained in (12). We make a similar expansion for the fields inside the cavity and for the phase difference across the Josephson junction. The input field in the free space, a f in (ω), is independent of ξ, as the output in this region does not reflect back. The task is to find the other fields as a function of the known input, a f in (ω), for small I c = 0. We solve the boundary condition at the junction order by order in ξ. By a straightforward calculation we find the leading-order solution (Appendix B) We observe that the operator is a (sinusoidal) function of the zeroth-order phasedifference operator (14). In the second order for ξ (and I c ) we obtain Here, the operator z(t) (∝ I c ) is a solution to the equation φ 1 (t) = [φ 0 (t), z(t)], where φ 1 is the phase-difference operator in the first order, obtained via the leading-order solution (16), see Appendix B. Important here is that the operators φ i (t) do not commute with each other. Thus, we have obtained a solution for (11) to second order in ξ as a function of the operator z(t), which still needs to be solved. In the case of semi-infinite transmission line (Z 0 = Z 1 ), we find a simple explicit form of z(t), Here, the junction RC-time gives the high frequency cut-off ω c = 1/Z 0 C J . The solution has an apparent divergence at t = t ′ , which cancels for symmetry reasons in all measurable quantities discussed in this article. We observe that the operator z(t) is also a trigonometric function of the zeroth-order phase-difference operator. Also is a scalar function. The trigonometric functions can be decomposed into exponential operators e ±i[φ 0 (t)−ω J t] , that correspond to charge transfers of 2e across the JJ in the two possible directions. (2) Thus, we see that the solutions (16)(17) include all possible tunneling processes up to second order.
We can now study the consistency of our solution, by checking if the output radiation field satisfies the commutation relation (8). This property is vital as it, for example, secures causality in the theory (13) . We obtain up to second order, using the input field commutation relation Since |R(ω)| = 1, the first term on the right-hand side produces the desired δ-function. A straightforward calculation (Appendix B) shows that the rest of the terms, order by order in ξ, sum to zero. This confirms that the commutation relation is indeed valid, up to the order our solution allows us to check this.

Emission characteristics I: Photon-flux density
Thus, having derived explicit expressions for the outgoing field-operators to second order in the junction's critical current, we go on to study properties of the output radiation. We first investigate general relations for the amplitude correlations, and after this consider their explicit forms. The truncation of the power series to the leading order can be made for small transparency JJs, i.e. for small I c . The exact definition of "small" is then addressed in Section 4.4. In later parts of the article, Sections 5 and 6, we discuss results of similar calculations but done for higher-order correlations.

General properties for the amplitude correlations
By a direct calculation, we obtain for the amplitude correlations related to the photonflux and power-spectral densities, We use here the notation a out (ω) ≡ a f out (ω). The function f (ω) is identified as the the photon-flux density (15) . This diagonal form is a result of the finite phase-coherence time, present already in the zeroth-order phase-difference (14). The phase difference performs a quantum Brownian motion in time, (16) and it follows that expectation values of type and for a single value ω J = 140 µeV/h (right panels). The Josephson radiation is seen as a diagonal resonance ω = ω J in the photon flux density (left panels). The flux is asymmetric around the diagonal, as for ω > ω J the emission is suppressed by the temperature, but for ω < ω J multi-photon production results in extra emission. When the pair production of photons dominates, the photon flux becomes symmetric around half of the Jospehson frequency ω J , as seen for the single picked value ω J = 140 µeV/h (right panels). In (b) this is approximately the sum of the two resonance frequencies in the cavity, and emission to these modes is enhanced.
e iφ 0 (t) e iφ 0 (t ′ ) are zero. This also implies that the amplitude correlations related to possible squeezing are zero, In general, due to the random phase fluctuations there is on average no phase coherence in the output radiation. Further, only even powers of the critical current are present in the power series of the photon-flux density, or of any higher-order correlator considered in this article, where functions F n (ω) are independent of I c . This follows again from the phase fluctuations, namely because Π m e ±iφ 0 (tm) = 0, for odd integers m.
We further divide the leading-order result for the output photon-flux-density, (19), as Here the term f t (ω) is the photon-flux density created by inelastic Cooper-pair tunneling, and the part f th (ω) describes reflection of the incoming thermal photons, being finite only for T = 0. In the following we consider the explicit forms of f t (ω) and f th (ω).

Emission from inelastic Cooper-pair tunneling f t (ω)
Radiation due to inelastic Cooper-pair tunneling is obtained by inserting the leadingorder solution for both operators a ( †) out in (19). A straightforward calculation gives (Appendix C), Here, the function P (E) is the probability to exchange energy E with the electromagnetic environment, in this case with the transmission line, defined as where the phase correlator function, , is a measure of phase fluctuations in the zeroth order. Equation (23) was obtained first in reference (6) by applying the theory of inelastic Cooper-pair tunneling (2) , i.e. P (E)-theory, and keeping track of the simultaneously emitted photons.
We will now analyze the obtained photon-flux density in more detail and compare it with the classical solution, which consists of continuous radiation at the Josephson frequency ω J , broadened by low-frequency phase fluctuations. The classical power spectral density, defined as S(ω) =hωf (ω), has the approximative form (14) S cl (ω) =h where we assume a small Γ = 4πk B T Z 0 /R Q ≡ 2hD compared tohω J . The same result is obtained also from (23) by inserting phase correlations of classical (thermal) phase diffusion (17) , J(t) = −D|t|. Especially, for T = 0 one has J(t) = 0 and P (E) = δ(E). Then all the radiation is emitted at ω = ω J with the total power I 2 c Re[Z t (ω J )]/2. In an exact classical solution also higher harmonics and a change in the dc-voltage across the junction exist, but the main picture remains.
In the quantum-mechanical treatment of the EM fields, two qualitative differences appear when T → 0: (i) the linewidth remains finite due to shot noise in the charge transport (14,18,19,20,21,22,23,24) and (ii) radiation below ω J has a finite tail due to multi-photon emission (6,7) . Whereas property (i) is not captured by the leading-order perturbation theory done here (except for the derivation of the zero-frequency shot noise, see section 4.4), property (ii) is seen as an asymmetric broadening of the P (E)-function. At zero temperature and for Z 0 = Z 1 , the P (E)-function has an approximative form (E > 0) (1,2) , where γ E is the Euler constant, E C J = e 2 /2C J is the junction charging energy, and ρ = Z 0 /R Q is the dimensionless resistance of the transmission line.
For E < 0 one has P (E) = 0, i.e. no energy can be extracted from the environment. This is obtained by using (1,16,25) . This zero temperature long-time behaviour is a good approximation also at finite temperatures for frequencies ω < ω J − k B T /h. For a typical low-Ohmic transmission line, ρ ≪ 1, the resulting power density is peaked at the Josephson frequency ω J with the magnitude This is symmetric around half the Josephson frequency, ω J /2, indicating that the radiation results from photons created in pairs (7) whose frequencies ω a and ω b add up to the Josephson frequency ω J . This result can also be derived also by straight linearization of the boundary condition, which includes maximally two photon emission processes, as demonstrated in Appendix A. Similar results hold also for the cavity configuration, Z 1 > Z 0 . Especially, if the Josephson frequency matches with the sum of the frequency of two modes, strong pair production to these modes is observed. Numerical results for the photon-flux density for the free-space and cavity configurations are presented in figure 2.

Elastic and inelastic reflection of thermal photons, f th (ω)
In addition to the radiation created by inelastic Cooper-pair tunneling, the leadingorder result (22) has a term proportional to the Bose factor, which we further divide as . The part f 0 (ω) describes the zeroth order (elastic) reflection of photons at the junction, The inelastic term f in (ω) comes from correlators between the zeroth and second-order operators. For the free-space configuration, We interpret this as inelastic reflection of thermal photons, exchanging energy with a Cooper-pair tunneling either direction. The term proportional to P (±hω J −hω) contributes as photon emission to the frequency ω, and the term ∝ P (±hω J +hω) as photon absorption from this frequency. Such processes do not contribute to the net current, and are a small correction to f t (ω) for the situations considered in this article.

Convergence
So far we have found out that it is the phase fluctuations across the JJ that describe Cooper-pair tunneling and simultaneous photon emission in the leading order. To study the convergence of the perturbation expansion, we then investigate the spectrum of the phase fluctuations at the junction. In particular, we compare the magnitude of the zeroth-order contribution with the magnitude of the leading-order contribution. For a rapidly converging perturbation expansion, the latter should be much smaller than the first. This should hold for all frequencies, since the right-hand side of the boundary condition at the junction (11), mixes all combinations of frequency terms summing up to ω J . This leads to the comparison Here

This leads us to the general condition
This is a relation for the smallness of the junction's critical current I c . For frequencies below the cut-off frequency ω c , , and in the following we will replace Z f (ω) by Z t (ω). We now examine the condition (31) explicitly at zero frequency, at the Josephson frequency, and then finally for frequencies between these special frequencies.
Let us investigate the zero frequency limit by multiplying each side of (31) byhω. We get Using the leading-order result for the simultaneous electric current This compares Johnson-Nyquist current noise (left-hand side) with the shot-noise coming from the charge transport (right-hand side), the noise considered also in Refs. (14,18,19,20,21,22,23,24) . For an Ohmic impedance (Z 0 = Z 1 ) this implies The second condition is calculated at the Josephson frequency ω = ω J . The contribution from ω J would act back to the low-frequency spectrum in the next perturbation round and would affect, for example, to a possible shift in the average voltage across the junction. Analysis at the Josephson frequency gives us To derive this, we have used the approximation P (0) ≈ (1/Γπ). We assume now, that the Josephson frequency does not match with a resonance frequency in the cavity (but it still can match a sum of two). We have then Re[Z t (ω J )] ≈ Re[Z t (0)], and we get the condition For an Ohmic impedance this is usually slightly more strict as (34), as typically ρ ∼ 1/20 > 1/2π. It is also independent of Z 0 . For the Ohmic case this can be then converted to a demand that thermal dephasing has to be faster than inelastic Cooper-pair tunneling, since However, if the Josephson frequency is exactly at the resonance, we get Here we have used the result for the Q-factors of the resonance modes Q n = (2n For the analysis at the middle frequencies For a resonant environment this sets a limit between the critical current and the sharpness (Q-factor) of the mode. One gets then the condition I c Z 1 Q 0 ≪ V . For an Ohmic impedance we get the condition I c Z 0 ≪ V , a known convergence condition for the higher orders of P (E)-theory obtained in reference (24) .

Emission characteristics II: Second-order coherence
To study statistics of the emitted photons in more detail, we investigate the secondorder coherence G (2) (τ ) for the output radiation, i.e. the probability to detect a pair of photons with time interval τ . The possibility for multi-photon emission implies bunching of the outgoing photons, meaning an increased probability of detecting photon pairs simultaneously. In this section, we consider our results for G (2) (τ ), obtained by including the leading contributions up to the fourth order in the critical current I c .

Photon coherences
We start with the first-order coherence, G (1) (τ ), for a continuous-mode field defined as (26) Similarly as before, we estimate this up to second order in I c . We can relate this to the photon-flux density, equation (22), and obtain Input-output description of microwave radiation in the dynamical Coulomb blockade 13 , 10 −2 (green), and 10 −1 (blue) c = (πE J / ω J ) 2 Figure 3. The normalized second-order coherence for a JJ connected directly to the free space, as estimated from (39) and (41) using the P (E)-function (26). We use here ρ = 10 −3 , 10 −2 , and 10 −1 , from the bottom to top, and we have chosen E C =hω J . We see that the bunching, g t (0), is close to the estimate (43) when ρ ≪ 1. The second-order coherence decays due to finite bandwidth, and has analogous form as the intensity pattern in single-slit diffraction.
In the following, we are interested in the contribution due to inelastic Cooper-pair tunneling, f t (ω), Here, we have neglected the vanishing contribution due to backward Cooper-pair tunneling, ∝ P (−hω J −hω). The second-order coherence, gives information of correlations between the emitted photons. This is defined for a continuous-mode field as (26) The leading-order contribution for (40) comes again from the second order in I c , which describes the effect of single-Cooper-pair tunneling. To obtain analytical results we calculate G t (τ ) for the JJ connected directly to the free space, Z 0 = Z 1 , at very low temperatures (for a more general expression see Appendix D). After a straightforward calculation we get Here, we have neglected terms proportional to the Bose factor, i.e. ∝ f 0 (ω). In the following we use this result to study photon bunching in the output radiation.
Input-output description of microwave radiation in the dynamical Coulomb blockade 14

Bunching
An important quantity describing photon emission is the relation between the first and second-order coherences, This basically compares probabilities for single and two photon detection. If g (2) (0) < 1, the field is called antibunched, and if g (2) (0) > 1 the field is bunched. For a Poissonian process the result is g (2) (0) = 1, while for thermal radiation g (2) (0) = 2.
Arbitrarily high bunching is possible also classically whereas antibunching is a pure sign of nonclassicality (26) .
With the analytical values of the first and the second-order coherences we can immediately get an estimate for the bunching in the free-space configuration (Z 0 = Z 1 ). We consider a typical transmission line (ρ ≪ 1, ω J ≪ ω c ) and solution (26), and obtain (Appendix D) This can be made arbitrary large by decreasing the critical current, i.e., the emitted power. This property is typical for a pair production of photons. Notably is that result (43) is independent of ρ, even though the power is proportional to ρ. In figure 3, we visualize the time dependence of g t (τ ). As Cooper-pair tunneling is also accompanied by an emission of low-energy photons, describing a simultaneous change in the voltage across the junction, it is more clear for the interpretation of the results not to include frequencies in the neighborhood of ω = 0 or ω = ω J . We consider then a small interval ∆ω of frequencies around half the Josephson frequency ω J /2, i.e. ω J /2 − ∆ω/2 < ω < ω J /2 + ∆ω/2, which in an experiment corresponds to a filtering of the output radiation (27) . One obtains for the corresponding second-order coherence (Appendix D) The related first-order coherence, within the same approximation, is G (1) (0) ≈ ρI 2 c Z 2 0 ∆ω/ω J . Therefore, we obtain the bunching, if measured in a small frequency interval around ω J /2, This is proportional to 1/ρ 2 , and for the considered case of small ρ, is much larger than result (43) for the complete output field. If we consider detection in a small frequency interval completely above half the Josephson frequency ω J /2, we obtain the leading-order result g (2) t (0) ≈ 0 (Appendix D), and exactly zero at the zero temperature. This is because the photon production at  2 , multiplied by the photon-flux densities at the two frequencies, i.e. N × f t (ω a )f t (ω b ). The colour scale is normalised to the maximal value, which is 2 × 10 5 higher for the cavity configuration (right) compared to the JJ connected directly to the free space (left). Negative values are a sign of nonclassicality and the parameters correspond to figure 2. The violation occurs around the condition of photon pair production ω a + ω b = ω J = 140 µeV/h. For the cavity configuration the observed nonclassicality is enhanced when the two frequencies match the two lowest modes of the cavity, maximising the photon pair production.
these frequencies occurs through multi-photon processes, and emission of two (or more) photons above ω J /2 is not possible from a single-Cooper-pair process. However, the result g (2) t (0) = 0 does not imply that the field is antibunched, since contributions from higher orders is neglected. The next-order contribution for G (2) t (0) comes from the fourth order, which has a special meaning as this is also the leading order of |G (1) t (0)| 2 . This order is also the first one to describe photon emission from two Cooper-pair tunnelings. Analytical results can be obtained for the Ohmic environment in the considered case ρ ≪ 1. The most important contribution becomes from a term describing two photon emission processes due to two (correlated) Cooper-pair tunnelings, a † 1 a † 1 a 1 a 1 . For small ρ and approximation J(t) = −D|t| − iπρSgn(t) (17) , we get through a lengthy analytical calculation a contribution g (2) t (0) = 2−B, whereB ∈ (1, 2). At the Josephson frequency and for a bandwidth larger than thermal dephasing D, we getB ≈ 1. Therefore the photon characteristics nearby this frequency is close to a Poissonian. Well below ω J the calculation givesB << 1.

Nonclassicality
The electromagnetic field is nonclassical if it cannot be described by the classical theory of electromagnetism. One example is a quadrature squeezed state of single mode, for which the width of the Wigner quasiprobability distribution in one of the quadratures is smaller than the width of a coherent state, i.e. the quantum description of a classical coherent signal. (28) The quadrature squeezing is measured through amplitude auto-and cross-correlations, which are phase-sensitive quantities. In our system, the Josephson junction is driven by a dc voltage, which suffers from both thermal and transport noise. As we will see, this leads to a rather short phase coherence time and no steadystate squeezing. There exists however a number of other relations, that are satisfied by a classical field, but can be violated by a general quantum mechanical field. (29) These are useful in our system, if they are immune to dephasing. An example of such a nonclassicality test is the Cauchy-Schwarz inequality for intensity auto-and cross-correlations, that is known to be violated maximally for a field created through parametric down conversion (28) .
In this section, we first address the question of quadrature squeezing in the output radiation, and then go on to derive a Cauchy-Schwarz inequality for photonflux correlations in the leading-order approximation, which we find to be an optimal way of detecting nonclassicality in the considered system.

Quadrature squeezing and dephasing
The pair production of photons implies quadrature squeezing, (28) which is characterized by correlators of type a out (ω)a out (ω ′ ) . The result (20), however, means that such nonclassical correlations do not exist on average, due to dephasing of the phase difference across the junction. This can be qualitatively visualized as a diffusion of the angle of quadrature squeezing. The situation is analogous to a parametric down conversion with an nonideal drive (26) .
To investigate how long it takes for the squeezing angle to be randomized, if one would know its value (or distribution) at time t = 0, we consider the phase-coherence function, In the long-time limit and for finite temperatures its behaviour is defined by the zerofrequency impedance Z 0 , via J(t) ∼ −D|t|, where D = 2πk B T ρ/h. Therefore, we identify D as the dephasing rate of quadrature squeezing. For typical values for the low-frequency impedance Z 0 = 50 Ω and T = 20 mK, one has 1/D ≈ 8 ns. Such dephasing times are therefore a very relevant property of a voltage-driven system, and a challenge for a measurement of phase-dependent system properties.

Cauchy-Schwarz inequality for intensity cross-correlations
A nonclassicality test that is not affected by phase fluctuations must be of higher order in the operators a ( †) out . The logical thing to do is to add two more operators to the ensemble average that characterizes squeezing, a out a out . Basically we have two possibilities to consider, the second-order coherence, of type a † out a † out a out a out , or the intensity correlator, of type a † out a out a † out a out . The second-order coherence was considered in Section 5, and was found to reveal a high degree of bunching, as a result of photon pair production. However, only antibunching would be a proof of nonclassicality. Therefore, we will now investigate a nonclassicality condition based on intensity cross-correlations. In the countable-mode case, a suitable Cauchy-Schwarz inequality is of the form (29) In the following, we apply this condition to the considered continous-mode case.
In the case of continuum of modes, we practically estimate G (2) (0) over a small frequency range ∆ω around ω 1 or ω 2 , and similarly for the corresponding cross correlator. Through a straightforward calculation we obtain for the considered auto-correlator (Appendix D) To keep the notation short we mark here a i ≡ a(ω i ). G (2) (0) at ω a is calculated up to the second-order in I c and similarly for the contribution at ω b . On the other hand, the intensity cross-correlations between the two frequencies have the form (when |ω a − ω b | > ∆ω) With the results (47-48) we get then the Cauchy-Schwarz inequality (in the limit ∆ω → 0 and calculated up to second-order in I c ) This result is valid for both the free-space and the cavity configuration. The inequality (49) is defined only via the P (E)-function (24). The left-hand side of (49) has a maximum when ω a + ω b = ω J , i.e. when the argument goes to zero. If at the same time |ω a − ω b | ≫ k B T , the right-hand side is close to zero, as one of the P (E)-functions has a large negative argument compared to the temperature. In this case the inequality becomes violated, which we visualize in figure 4. The violation is due to nonclassical photon pair production. Generally at ω a = ω b the nonclassicality cannot be tested with this inequality since the two sides are equal by definition. The use of a resonant environment (Z 1 ≫ Z 0 ) does not change the violation of the inequality (49) qualitatively. However, it significantly increase the photon emission rate at certain frequencies, which facilitates experimental detection (7) . As the contribution from thermal radiation is neglected here, the tested frequencies ω a(b) should be well above k B T /h. Also, in an experiment a detection over a finite bandwidth is used, whose effect should be carefully analyzed.

Conclusions and outlook
In conclusion, we have derived a continuous-mode solution for microwave radiation in a transmission line with a dc-voltage bias and which is terminated by a small Josephson junction. This is done by a perturbative treatment of the boundary condition that describes Cooper-pair tunneling across the Josephson junction. We showed that the method reproduces the previously derived expression for the created photon-flux density, obtained by applying the P (E)-theory. We extended this first by determining the corresponding second-order coherence. We found that the emitted microwave field has a high degree of bunching due to photon pair production at frequencies below the Josephson frequency. We then addressed the question of nonclassicality of the emitted radiation in this region and showed that the photon pair production violates the classical Cauchy-Schwarz inequality for intensity cross-correlations.
The established method opens a possibility for further detailed study of the radiation characteristics in this system. For example, calculations in higher order access the question of the effect of correlations between consecutive tunneling Cooper pairs. For the considered case of low-Ohmic environment, we obtained bunching in the output radiation. On the other hand, when the transmission-line impedance is increased beyond the resistance quantum, antibunching of the Cooper-pair tunneling is expected, due to Coulomb blockade. (2) In this regime, the output photons should also be antibunched. Also, summation to all orders can be feasible, if known summation methods for this type of perturbation expansions work. Overall, this system is very rich in physics, covering the limit of dynamical Coulomb blockade at low impedances, to Coulomb blockade in the high-impedance limit. The question of the detailed form and properties of the related output radiation makes this system very interesting for future works. This is motivated also by the technical development towards simultaneous measurements of both microwaves and electrical currents.
Appendix A: Linearization of the boundary condition at the junction A straightforward way to solve for the out-field is to linearize the boundary condition (11) and Fourier transform the problem. The silent assumption is the small fluctuations of the phase difference, φ(t). This is actually not usually the case, since the (zeroth order) phase difference performs a quantum Brownian motion in time (2,16) . However, the linearization turns out to give correct results for frequencies ω ∼ ω J /2 (ρ ≪ 1), where such fluctuations stay small. This is consistent with pair production of photons in this frequency range.
In the following we consider linearization in the case of the free-space configuration, Z 0 = Z 1 , and take the limit C J → 0. To do this properly, we rewrite the right-hand side of (11) using the identity Expanding the right-hand side of this up to linear order in φ(t), and Fourier Input-output description of microwave radiation in the dynamical Coulomb blockade 19 transforming, we get for this Here, for simplicity, we have introduced negative frequencies as a(−ω) ≡ a † (ω). The first two terms represent radiation at the Josephson frequency ω J , while the other terms describe mixing of this with additional photonic process where, for example, ω J is splitted into two frequencies.
We continue by solving the corresponding boundary condition, This can be done by writing the equation into the matrix form Here the frequency vector A is constructed as (11) where |δω| < ω J . This form is possible since the boundary condition mixes only frequencies differing by ω J . We have also introduced a cut-off Ω = Nω J . The matrices M have diagonals 1 and first nondiagonals (in the n:th row) d ± [(−N − 1 + n)ω J + δω], where the plus sign corresponds to the term M n,n+1 and we have and d in = −d out . Equation (.1) has to be solved generally numerically. For small d approximative analytical solution can be sought with the ansats We then find a solution in the lowest order for I c S ± (ω) = −2d ± (ω).

(.3)
We can now estimate the photon flux density below the Josephson frequency. Using (.3) one gets for ω < ω J and at T = 0 .

Input-output description of microwave radiation in the dynamical Coulomb blockade 20
Appendix B: Perturbative input-output approach

(.9)
Here the output (input) cavity field in the n:th order is labelled as b n (c n ). We rewrite the boundary condition at the junction as b n (x = 0) = c n (x = 0)e iθ(ω) + (.10) Here the formal operation {·} n picks out n:th order contribution. It follows then a n (ω) = i Z 1 hωπ In the leading order, we can only include zeroth-order phase difference in the Taylor expansion of operators e ±iφ 0 +ξφ 1 +ξ 2 ... , and we immedately obtain (16). Generally, one can solve for the phase difference at the junction in the n:th order, √ ω B(ω)a n (ω)e −iωt + H.c., (.12) This is self consistent, since the same order result for a n (ω) depends only on the previousorder phase-differences.

Input-output description of microwave radiation in the dynamical Coulomb blockade 21
The phase difference in the leading order has a central role when constructing the general solution. We get for this We calculate the explicit form for Z 0 = Z 1 , We aim to expand the term ξI c sin[φ(t) − ω J t] to second order in ξ. We have formally φ(t) = φ 0 (t) + ξφ 1 (t) + . . ., and we need to properly expand the functions to first order in ξ. We know that up to first order in ξ we can include only operators φ 0 (t) and ξφ 1 (t) in the Taylor expansion. Therefore we can put φ(t) = φ 0 (t) + ξφ 1 (t).
We now define an operator z(t) through the relation φ 1 (t) = [φ 0 (t), z(t)]. Through a direct calculation of the Taylor expansion one gets for the first-order contribution To solve z we first evaluate commutator We get for the special case Z 0 = Z 1 , We then investigate the commutator Here we have used that [φ 0 (t), e ±iφ 0 (t ′ ) ] = ±iC(t − t ′ )e ±iφ 0 (t ′ ) . Comparing this with solution (.12) we deduce .

(.20)
For the open-space configuration (Z 0 = Z 1 ) we get theñ which leads to (17). This expansion (with methods shown in Appendix C) can be also used to rederive the P (E)-theory net current across the JJ, used in section 4.4.

Derivation of the term f t (ω)
Using the leading order solutions for both operators a ( †) out in the expression a † out (ω)a out (ω ′ ) , we get a contribution for the photon flux Here we have used the fact that expectation values of the form e iφ 0 (t) e iφ 0 (t ′ ) are zero due to random phase fluctuations. Also contributions such as a 0 (ω)a 1 (ω ′ ) vanish due to the same reason. We use now the following property of bosonic operators (2) , where . We then perform a change of variables, x = t − t ′ , y = (t + t ′ )/2, do integrations over x and y, and obtain Derivation of the term f in (ω) We derive now the inelastic reflection of thermal photons, f in (ω). To do this we take use of the zeroth-order solution (.22) and the second-order solution (17). We get The next step is to calculate the ensemble average. By applying the Wick's theorem we get We have also e ±iφ 0 (t) e ±iφ 0 (t ′ ) = 0. These relations lead to We can perform integration over t ′′ by using The term inside the last parentheses can be put into the form Using these relation, we obtain For the last two time integrations we do a change of variables x = t − t ′ and y = (t+t ′ )/2. The resulting y-dependence is in a factor exp[iy(ω ′′ −ω)]. The integration over y leads to the factor 2πδ(ω − ω ′′ ). Performing integration over ω ′′ and division by 2π (to obtain f in ), one gets Adding this with the contribution a † 0 (ω)a 2 (ω ′′ ) = a † 2 (ω ′′ )a 0 (ω) * , and using ω = ω ′′ , one sees that only two times the real part of (.32) survives. Thus, We know that J(−x) = J * (x), and therefore Im[e J(−x) ] = −Im[e J(x) ]. Therefore the part ∝ Sgn(x) cancels out due to symmetry reasons. Because only e J(x) is a complex number, the result does not change if we take the imaginary part over the whole expression [without the part ∝ Sgn(x)], instead of only over e J(x) . This leads to result (29).
Appendix D: Estimating higher-order coherences up to second order in I c We want to calculate expressions such out . We neglect the contribution if using once the second order solution a † 2 , as this is proportional to f 0 (ω). We will also neglect backward directed Cooper-pair tunneling, i.e. we take the approximation a 1 (ω) ∝ dte i(ω−ω J )t e iφ 0 (t) and neglect terms of type ∝ dte i(ω+ω J )t e −iφ 0 (t) . Such tunneling against the voltage is well suppressed by the temperature.
We take use of the following result for the ensemble average of bosonic operators φ, Here we use the notation φ ′ = φ 0 (t ′ ) and similarly for others. The result can be derived by expressing the exponential functions as powers series and applying the Wick's theorem. Important here is that the order of the operators φ i stays the same when contracted into the pairs. This is then also valid for permutations of the initial operators. The result (.35) is also immune to exchanging the signs in the exponents.

Intensity cross-correlations
We consider first the correlator Once we obtain expression for this, the other orderings of the operators a ( †) out can be deduced by using general relations for the field operators, equation (8). We need to calculate the sum of all the orderings, However, it turns out that only the first of these four terms is important, as the other terms are proportional to f 0 (ω), and can be neglected. This can be understood by rewriting the ensemble average with the help of a formal density matrixρ of the unperturbed system, Tr a(ω ′ )a † (ω ′′ )a(ω ′′′ )ρa † (ω) : tunneling with radiation (∝ e ±iφ 0 (t) ) has to be inserted around the density matrix, otherwise the result is zero at To calculate (.36), the difficulty is to perform integration over all times of the term In result (.35) the first term inside the brackets is the easiest to calculate, as the timedependent terms are only functions of t − t ′′′ or t ′ − t ′′ . We do twice similar change of variables as when calculating f t (ω) (Appendix C) and obtain the first contribution for the term (.37), Here we mark p(t ′ − t ′′ ) ≡ φ(t ′ )φ(t ′′ ) and it is determined by Eq. (15). We have neglected a contribution proportional to f 0 (ω). Only the last term, proportional to φφ ′′ φ ′ φ ′′′ , gives another finite contribution at T = 0, We note that in this case the operators of same type are paired [a(ω ′ ) ↔ a(ω ′′′ ), a † (ω) ↔ a † (ω ′′ )]. To obtain expression for (.36) we sum up these two results and multiply them withÃ ω ′ ω ′′ ωω ′′′ Ic 8eπ 2 , wherẽ . (.40) We obtain for the intensity cross correlations up to second order in I c , Second-order coherence G (2) We consider now the expectation value a † (ω)a † (ω ′ )a(ω ′′ )a(ω ′′′ ) . By using the identity [a out (ω), a † out (ω ′ )] = δ(ω − ω ′ ) and doing the exchange ω ′ ↔ ω ′′ in term (.41), we get the result a † (ω)a † (ω ′ )a(ω ′′ )a(ω ′′′ ) = This can also be derived through a direct calculation, as was done for the intensity cross correlator.

Cauchy-Schwarz inequality
The Cauchy-Schwarz inequality compares intensity correlations with the second-order coherence. We calculate these in a small frequency interval ∆ω around the frequencies ω a and ω b . This assumes a filtering of the measured output signal into these frequencies (intensity-correlations).