A Josephson radiation comb generator

We propose the implementation of a Josephson Radiation Comb Generator (JRCG) based on a dc superconducting quantum interference device (SQUID) driven by an external magnetic field. When the magnetic flux crosses a diffraction node of the critical current interference pattern, the superconducting phase undergoes a jump of π and a voltage pulse is generated at the extremes of the SQUID. Under periodic drive this allows one to generate a sequence of sharp, evenly spaced voltage pulses. In the frequency domain, this corresponds to a comb-like structure similar to the one exploited in optics and metrology. With this device it is possible to generate up to several hundreds of harmonics of the driving frequency. For example, a chain of 50 identical high-critical-temperature SQUIDs driven at 1 GHz can deliver up to a 0.5 nW at 200 GHz. The availability of a fully solid-state radiation comb generator such as the JRCG, easily integrable on chip, may pave the way to a number of technological applications, from metrology to sub-millimeter wave generation.


I. SOLUTION OF THE RCSJ EQUATION IN A DC SQUID
We consider a dc SQUID composed by two Josephson junctions and subject to a magnetic flux Φ. The total current though the SQUID is I J = I c1 sin ϕ 1 + I c2 sin ϕ 2 , where I ci and ϕ i are the critical current and the phase across the i-th junction, respectively. Because of the flux quantization constraint, it follows that (ϕ 1 − ϕ 2 )/2 = πΦ/Φ 0 . Introducing the phase across the SQUID ϕ = (ϕ 1 + ϕ 2 )/2, we get where φ = πΦ/Φ 0 , I + = I c1 + I c2 , r = (I c1 − I c2 )/(I c1 + I c2 ) and Φ 0 2 × 10 −15 Wb is the flux quantum.
Starting from the RCSJ model [1], we can write an equation of motion for the phase ϕ as where C is the capacitance, R is the total shunting resistance of the SQUID, I B is the external biasing current and f (ϕ, t) = I J [ϕ; φ(t)]/I + . We rescale the above equation in terms of the driving frequency ν: where δ = I B /I + , c = 2πRCν and A. Analytical solution for I B = 0 Let us consider the case of a symmetric dc SQUID (r = 0), overdamped junctions (c ≈ 0) and zero current bias (δ = 0). Then (3) reduces to where F (τ ) = cos[φ(τ )]. We notice that if the initial condition is ϕ(0) = ϕ 0 = kπ, the above equation has trivial dynamics ϕ(t) = 0. This means that, even if small, we cannot neglect the influence of theφ term. However, if ϕ = kπ, we can effectively neglect the capacitive contribution and solve analytically the differential equation (3) to obtain We suppose that |α| 1 (and α > 0). If t 0 dτ F (τ ) assumes positive and negative values, the argument of the arctangent increases or decreases exponentially depending on its sign. For simplicity, we consider t 0 dτ F (τ ) < 0 and the case of small ϕ 0 . For ϕ 0 > 0, ϕ(t) exponentially reaches π, and for ϕ 0 < 0, the evolution is similar but ϕ(t) varies between ϕ 0 and −π. Therefore, in both cases we have an exponential π jump of the phase but its direction is determined by the initial phase ϕ 0 . The rate of the exponential jump is determined by I + R: the larger I + R, the sharper the voltage pulse.
Recalling the phase particle analogy discussed in the main text, with c = 0 if the phase is initially in the minimum ϕ 0 = 2kπ it will remains in the same point even when it becomes an unstable maximum. A small shift of the initial condition induces an exponential dynamics since at t = 1/(4ν) the phase is close to (but not on) a potential maximum. The direction of this shift (and, in the case discussed, the sign of ϕ 0 ) determines the direction of the fall and the jump of the phase.
The change in time of ϕ in Eq. (6) is associated to a voltage (Notice that we have gone back the real time t = τ /(2πν) unit.) For a small oscillation around the first diffraction node it is possible to have an analytical expression for V (t). In this case, cos(πΦ/Φ 0 ) ≈ (π/2) cos(2πνt). After a straightforward calculation, we approximate the function for small times, i.e., t 1/ν to obtain The above expression has an exponential behavior (increasing and decreasing) when the system crosses a diffraction node.
The mean value of the voltage square, which is related to the delivered power, is given where T is the averaging time. Taking into account Eq. (4) and for long average time T 2/(πα ν), we obtain

B. Solution in presence of a current bias
We consider now the case of a small biasing current I B in a symmetric SQUID (r = 0).
By "small biasing current", we mean that its effect must be negligible with respect to the driven dynamics, i.e., δ α, but must dominate the capacitor dynamics, i.e., δ c.
When ϕ ≈ kπ, the Josephson current contribution αF (τ ) sin ϕ is small and the dynamics is determined only by I B . Equation (3) reduces to Away from ϕ ≈ kπ, the Josephson current contribution dominates and, therefore, the equation for motion approximatively reads Since the behavior of ϕ(t) under these two different dynamics is drastically different (linear versus exponential change), we can tune the system parameters in order to effectively separate the two regimes governed by Eqs. (10) and (11) and generate the sequence of exponential phase jumps.
In other terms, the presence of a current bias has two effects. The first is to transport the system away from the region ϕ ≈ kπ in which the drive contribution is small and the dynamics is dominated by the capacitive term. The second is to breaks the symmetry of the system inducing the jumps always in the same direction.
The effect of a current bias can be interpreted in terms of the tilted-washboard potential discussed in the main text.

C. Asymmetric SQUID case
From Eq. (1), the energy potential reads (we neglect the current bias and set δ = 0) [1] To find the position of the maxima and the minima we derive E J (t) with respect to ϕ and equal it to zero. The corresponding equation reads tan ϕ = r tan φ. For r = 0, we see that the values of ϕ satisfying the above equation do not depend on time. On the contrary, since φ = πΦ/Φ 0 , they depend on time through Φ for any r = 0 as discussed in the main text.
For a periodic drive, the particle is found alternatively on the left and on the right of the maximum and it rolls in alternate directions producing the alternate pattern of the voltage pulses.

II. FINITE-SIZE EFFECTS
We consider a long chain of SQUIDs located at x = x k , (k = 1, 2, ..., N ) and coupled to a transmission line of length L. We assume that the first SQUID is at point x 0 = 0 and that x k = ka, where a is the distance between the SQUIDs. Every SQUID emits radiation with the spectrum v(ω). The time dependent voltage drop across it is where ω 0 is the driving frequency. Then the voltage at the end of the transmission line, i.e., at x = L, can be estimated as a sum of the signals coming from each SQUID: where c is the speed of electromagnetic waves propagating along the transmission line. For simplicity, we assumed that the speed is the same both in the transmission line and in the chain of SQUIDs. Performing summation over k in Eq. (13) and considering that N 1, one finds Thus, the intensity of n-th harmonics at The power of the n-th harmonics scales with the number of junction as N 2 if N anω/(2c) ≤ 1.
For the parameters used in the simulation, considering a wavelength of λ = 6 × 10 −4 m corresponding to 200 GHz in Silicon, and a SQUID distance of a = 1 µm, we have that for N = 50 the array can be considered as a lumped element.

III. POWER DELIVERED TO THE LOAD RESISTOR
As discussed in the main text, current conservation through the SQUID implies that dynamics of each SQUID is independent from the rest of the chain. As a consequence, the voltage at the extremes of the chain scales as the number N of SQUIDs. Accordingly, the intrinsic power generated by the device (that is, the power delivered to an infinite load) scales as N 2 . In this section we discuss the effect of a finite load on the JRCG. Let the SQUIDs be identical and V (t) be the voltage drop across each SQUID. Then the voltage at the load is V tot (t) = N V (t) and the current flowing through it is is the (real) impedance of the load. This current flows back to the SQUID array and adds up, with opposite sign, to the bias current I B . As a result, (2) must be rewritten as Recalling that V i (t) = ( /2e)φ, we obtain with The effective change in the shunt resistance modifies the SQUID dynamics. In particular, it also reduces the power P = N 2 V 2 /R L that can be delivered to the load. From Eq. (9), we see that for a single SQUID V 2 ∝ R eff . Using (18), we find that P scales as N 2 for

IV. EFFECT OF THE THERMAL NOISE
To estimate the effect of the thermal noise we use the Langevin equation where ξ(t) is the white noise with correlation function We have numerically solved the associated stochastic differential equation in case of a YBCO SQUID with 1 GHz drive. We have considered a symmetric YBCO SQUID with negligible capacitance. The noise source has been taken at temperature of 4.2 K. With these parameters, the dynamics of the SQUID shown in Fig. 1 is essentially identical to that without noise source.
The thermal noise results in a small broadening of the resonances at high frequency (see Fig. 2). The signal to noise ratio is still about 10 3 and it can be further increased by decreasing the working temperature of the device.
The effect of noise is maximum when the energy barrier is shallow and undesired transitions are most likely to occur. Therefore, we expect an increased noise influence for slow frequency drive since the system remains in a shallow barrier potential for a longer time.
However, even in this situation these noise effects can be reduced by increasing I B in order to restore the privileged direction of the dynamics.