Enhanced Dynamic Casimir Effect in Temporally and Spatially Modulated Josephson Transmission Line

Real photon pairs can be created in a dynamic cavity with an oscillating boundary or temporally modulated refractive index of the constituent medium. This effect is called dynamic Casimir effect (DCE), which represents one of the most amazing predictions of quantum field theory. The DCE has been experimentally observed in Josephson metamaterials embedded in a microwave cavity. However, the efficiency of the observed DCE is extremely weak, entailing a complex external signal enhancement process to detect the signal. Here, it is shown that the DCE can be drastically enhanced in a dynamic 1D cavity consisting of a superconducting quantum interference device (SQUID)‐based Josephson transmission line with both temporal and spatial modulation on the effective inductance profile through flux‐biasing. Such a system can resonantly generate photons at driving frequencies equal to even or odd integer times of that of the fundamental cavity mode governed by the symmetry of the spatial modulation. Interesting spectral and scaling behaviors for photons excited at the band edge are further observed. The discovery introduces a new degree of freedom—spatial modulation—to enhance the efficiency of DCE.


Introduction
Quantum field theory predicts that vacuum is not purely empty but has fluctuations, which causes the so-called Casimir force when the presence of metals or dielectrics alters the vacuum expectation value of the energy of the second quantized electromagnetic field. These vacuum fluctuations can even be converted into real particles by certain dynamical external perturbations, and this process is called the dynamic Casimir effect, [1][2][3][4][5][6][7][8][9] which describes the creation of photon pairs from a quantum vacuum when the boundary conditions or the material DOI: 10.1002/lpor.201900164 properties are dynamically changed. Due to a time-dependent geometric/material configuration, the quantum field experiences parametric amplification processes over time, which is accompanied by the generation of real quanta. Among various studied DCE systems, the most frequently discussed configurations are 1D homogeneous cavities with vibrating boundaries [1,2] and cavities with fixed boundaries but dynamically modulated optical properties. [3,4] Usually, a harmonically driven configuration [5][6][7][8] is employed to enhance DCE through resonant excitation at certain frequencies. In this process, DCE predicts a direct energy conversion from mechanical kinetic energy to electromagnetic energy.
However, to observe a detectable signal of DCE, the moving speed of the geometrical configuration must be a substantial fraction of the speed of light, which becomes unrealistic for a macroscopic mechanical system. [10] Various solutions have been theoretically proposed to achieve this rapid modulation of the cavity resonances, for example, by using laser or magnetic field to modulate the reflection property of the semiconductor [11,12] or superconducting [13] mirror, or through non-adiabatic temporal modulation of the light-matter coupling strength [14][15][16] in cavity QED systems, or by using modulation of nonlinear materials, [17] cold matter material, [18] or superconductor circuits [19] to vary the refractive index of the filler inside the cavity. Recently, by using SQUID-based circuits, [19][20][21][22][23][24][25][26] the effective length or refraction property of a cavity was tuned by modulating the effective inductance through a time-varying magnetic flux at very high frequencies (>10 GHz), leading to experimental observation of squeezing signals of DCE in the microwave regime, which verified the theoretical prediction.
Thus far, most previous theories and experiments [19,20] were based on dynamic cavities with homogeneous media, which are weakly excited even at resonance. Furthermore, DCE in such a homogeneous cavity can only either involve a single photonic mode by modulating the refraction index, [19] or introduce a fixed type of interaction between different photonic modes by modulating the boundaries. [20] Herein, we introduce a new spatial modulation degree of freedom by inhomogeneously varying the effective inductance in a cavity containing a 1D SQUID-based Josephson transmission line (JTL). [27][28][29][30][31][32][33][34] We discuss the squeezing (degenerate parametric process for producing single-mode squeezing) and acceleration (nondegenerate parametric process for producing two-mode squeezing) effect caused by the variation of resonant frequencies and eigenmodes, respectively. [35,36] We numerically calculate the intensity of photon flux due to these two processes. A significant DCE occurs only for resonant excitation at even or odd integer harmonics of the fundamental resonant frequency, which depends on the spatial symmetries of the dynamic perturbations on the JTL cavity. Besides, the photon creation rate displays interesting scaling behaviors with the number of unit cells. Our results provide a powerful degree of freedom to enhance the rate of the DCE photon generation process, and open new avenues for investigating DCE in metamaterials, photonic crystals, and related systems.

Temporally and Spatially Modulated SQUID-Based JTL Cavity
The 1D SQUID-based JTL cavity is shown in Figure 1, which consists of an array of DC-SQUIDs embedded in a coplanar transmission line. Each DC-SQUID consists of two Josephson junctions connected in parallel to form a loop. By assuming that all junctions have approximately the same fabrication parameters, a weakly excited DC-SQUID junction with working frequency well below its plasma frequency (ω p /2π = 36 GHz in ref. [21] see the Material section) can behave as a linear inductor with a local Josephson inductance: L J ( ext ) = ( 0 /2π) 2 · [2E J cos(π ext / 0 )] −1 , where E J is the Josephson energy, 0 = h/2e is the magnetic flux quantum, and ext is the external magnetic flux. [20][21][22]29,30] This SQUID array can couple to the external transmission line through capacitors. The capacitor enforces the current mode that flows through the SQUID and can be considered as an imperfect open boundary (Neumann boundary condition) with small output current. In our design, a DC-signal is used to generate a static bias-flux, while another AC-signal is used to dynamically tune the external flux and effective inductance of the SQUID circuit. We design a special on-chip AC flux line-an S-shaped wiggling line consisting of connected u-shape/n-shape semi-square unit of size d with opposite orientations. It is placed underneath the superconducting circuit layer and separated by a dielectric layer. The different winding directions between neighboring semi-squares, marked as "A" and "B," generate opposite directions of the perturbed magnetic flux. Therefore, the external magnetic flux creates both temporal and spatial modulation: ext (x, t) = dc + ac (x, t) = dc + δ · h(x) sin(ω d t − γ d ), with ω d the modulation frequency, γ d the phase of modulation field, and dc /δ the strength of DC/AC flux. The distribution function h(x), in the ideal case, is a periodic rectangular function, with the periodic length of each unit cell D and the length of each segment equivalent to the size d of the semi-square AC line. This function h(x) is equal to 1 for "A" segment and −1 for "B" segment, which is fully determined by the AC flux line setup. As shown later, it is the symmetry of the perturbation distribution that affects the excitation, and the periodic rectangular waveform h(x) provides a good approximation for the external magnetic flux. Furthermore, when the amplitude of the AC-signal flux is much smaller than the effective DC-signal flux, and both the two fluxes are much smaller than the magnetic flux quantum: ac mod( dc , 0 ) 0 (mod is short for modulo operation), we can realize the temporally and spatially modulated effective inductance.
with the time-dependent term δ L (t) = δ L 0 · sin(ω d t − γ d ) and the amplitude term δ L 0 = π δ / 0 · tan(π dc / 0 ) 1. Herein, we will focus on two setups-spatially symmetric ABBA and antisymmetric AABB external flux, respectively. The flux distribution (equivalent to the inductance distribution) consists of N c repeated unit cells with a total length of 4N c ·d. The parameters of a typical SQUID are detailed in the Material section. Following the standard techniques of field quantization of transmission line, [37,38] we quantize the field in such a SQUIDsbased JTL and obtain the Hamiltonian as H = n ω n a + n a n + 1 2 where a n and a + n are, respectively, the canonical annihilation and creation operators for the nth instantaneous eigenmode u n (x, t) in a cavity with resonant frequency ω n (t), and M nm (t) = n|ṁ describes the coupling strength between any two eigenmodes (see more detailed analysis in the Supporting Information). From the Hamiltonian and the definition of the instantaneous annihilation/creation operator defined in Equation (S1.11, Supporting Information), the evolution equation for the nth canonical annihilation operator can be obtained from d dt a n = ∂ ∂t a n + 1 i [ a n , H], which is where ξ n = 1 2 · ∂ log ωn ∂t and μ ± nm = 1 2 · ( ωn ωm ± ωm ωn ) · M nm . On the right hand side, the first term is the conventional harmonic term, while the other two terms represent the two DCE contributions responsible for the creation of photon pairs: the squeezing term ξ n arising from the variation of the resonant frequency, and the acceleration term μ ± nm caused mainly by the variation of the eigenfunction. [39][40][41] The concrete forms of these two terms will fully determine the behavior of DCE, which, as shown later, depend on the spatial symmetry of the AC flux line configuration.
It is noted that Equations (2) and (3) only describe DCE with perfect Neumann boundary condition. However, for a system implemented with realistic capacitors, the imperfect boundaries should be modeled by the standard input-output formalism discussed later. Since the Q-factor of a SQUID-based microwave cavity can easily reach 10 5 , [27] the intensity of radiation calculated from this approximation process can be sufficiently accurate.

Symmetry-Dependent DCE Contributions
With perfect Neumann boundary condition I = 0, at an arbitrary time t, we can solve the instantaneous eigenmodes u n (x, t) of JTL cavity through the transfer matrix method, in which both the eigenfrequency and eigenfunction vary harmonically in time.  Figure 2e,f show the band structures and the typical field distributions for these two systems: ABBA and AABB with N c repeated unit cells. In the case without the time-dependent perturbation δ L = 0, the system functions just the same as a homogeneous vacuum cavity with standing wave eigenmodes and its band structure is located at a series of equally separated points (K n,0 , ω n,0 ) = (n s π /N c D, nπ v dc /N c D), in which n s denotes folded index of the eigenmode index n to the 1st Brillouin zone, and v dc denotes the speed of wave with zero AC signal. Clearly, there are N c eigenmodes in each band which equals to the number of repeated unit cells. With a small perturbation δ L , the band structure exhibits only a noticeable deformation close to the band edge denoted by the K-point in Figure 2b,g. A significant energy shift occurs only at the Brillouin zone boundary because of the breaking of the twofold degeneracy at K-point due to the perturbation. But for the two systems being investigated, completely different behaviors are observed near the K-point: the system with symmetric ABBA perturbation has an obvious frequency shift, while the one with antisymmetric AABB perturbation exhibits a very slight frequency shift but the presence of an edge eigenmode, as shown in Figure 2b,g. This difference will result in very distinct behaviors of the DCE photon creation process in these two systems.
For the given harmonic perturbation described in Equation (1), coefficients of the two terms in Equation (3) that lead to DCE can be expressed as which are both modulated at frequency ω d . With δ L << 1, the region that we are interested in, those coefficients can be analytically solved through the standard time-independent perturbation method. [42][43][44] Specifically, we start from the exact and analytical standing-wave eigenmodes in the configuration with zero AC flux and then treat the variation of the effective inductance as a perturbation term. We can then analytically solve the corrections of each eigenmode, and obtain two amplitude functions of Equation (4).
which perfectly matches the numerical results (see details in the Supporting Information). In Equation (5), h nm is a timeindependent symmetric matrix with h nm = h mn , which is defined as The form of matrix h nm is fully determined by the distribution function h(x), or by a more decisive property: the symmetry of h(x), which, as we will show, leads to the interesting selection rules in the DCE process.
To understand the influences of symmetry in different JTL configurations, we analyze the concrete forms of the two amplitude functions w n and C nm . They are related to the Fourier components , which strongly depends on the symmetry of the distribution function h(x). In particular, w n of the squeezing term ξ n is In the ABBA configuration, an even form of h(x) renders h n a nonzero value only when n = n K ࣕ (2s+1)·N c (s is an arbitrary integer), which corresponds to the K-eigenmodes at the edge of the Brillouin zone. This is consistent with the observed frequency shift at those K-points in the ABBA system. On the contrary, in the AABB configuration, an odd h(x) function forces h n (h nn and w n ) zero for all modes, which is consistent with the nearly δ L -independent eigenfrequency in the AABB system. The more accurate numerical results of w n presented in Figure 2c,g for the two configurations with five repeated unit cells verify our perturbative analysis.
The amplitude function C nm of the acceleration contribution µ nm ± is an anti-symmetric real matrix with C nm = -C mn = Re(C nm ). For n ± m = 2l (n and m are both even or odd), C nm has nonzero value only when the lth Fourier component h l is nonzero. Since h l is nonzero when l = n K in the ABBA system, only the eigenmode pairs satisfying n ± m = 2n K can couple to each other. These eigenmode pairs have the same wave vector and nonzero C nm coupling elements. On the other hand, for n ± m = 2l + 1, C nm couples two eigenmodes with different parities. Therefore, C nm is always zero in a symmetric system, whereas there is no constraint on its value in an anti-symmetric system. Based on the discussion above, we can obtain all the acceleration matrix elements for the two systems, as shown in Figure 2d,h. They are both sparse matrices with nonzero elements only when the index pair (n, m) has the same (different) parity in the symmetric (anti-symmetric) configuration.
Worthy of mentioning, the position where the Fourier component h l becomes zero is solely determined by the symmetry of the perturbation. Even in the case of N c >> 1, where the ABBA and AABB configurations appear to be almost identical except for the boundary position, the significant symmetry differences between the two systems still result in totally different behaviors. A more complicate distribution function h(x) due to more detailed features in the magnetic flux can only slightly alter the value of h nm , but does not affect the resonant excitation conditions of DCE contributions. Therefore, the spatial symmetry of the AC flux line is a decisive factor that determines the distinct spectra and resonance conditions of the photon creation process in these two systems.

Input-Output Model for the Radiation Spectrum
To compute the output photon flux spectrum, we use the framework of the standard input-output theory, by allowing for an external field a in to couple into the cavity. The evolution equation, cf. Equation (3), is now modified as in refs. [19][20][21]25,26]: where κ n denotes the coupling strength. Based on the radiation property of a Fabry-Pérot cavity, [45] we can assume that these coupling strengths κ n are the same for all eigenmodes with κ n = κ = v dc π/(2D · Q K ), where Q K is the Q-factor of the first edge eigenmode (the K-mode shown in Figure 2b,f). To obtain a time-independent form for the Hamiltonian, further progress is needed. Indeed, both the squeezing/acceleration terms and the creation/annihilation operators contain fast oscillatory phase. By averaging these time-dependent terms and ignoring those fast oscillating ones, a set of closed equations can be obtained to describe the DCE process in the JTL cavity, which solves the output canonical operators in terms of the input canonical operators (see details in the Supporting Information). Furthermore, for the radiation on a frequency lower than the modulation frequency ω d = g·π v dc /N c D = gω 1,0 , by ignoring the small scattering correction μ + nmâ m , [21] we can obtain I 2×2Nmax is the direct product of a 2-by-2 identity matrix and a 1-by-N max vector of ones. N max is the truncated eigenmodes index during the numerical calculation process, and the summation of this index indicates that the output photon flux consists of contributions from all the modes. In this expression, the N max -by-N max diagonal matrix χ with χ n (ω) = [κ − i(ω − ω n,0 )] −1 describes the cavity susceptibility. The other N max -by-N max matrix F =F II −F I describes the effective DCE contribution, where the diagonal matrix part F I with F nn I (ω d ) = 1 4 δ L 0 ω d · w n δ(2n − g ) and the symmetric off-diagonal matrix part F II with F nm I I (ω d ) = 1 4 δ L 0 ω d · ( √ n/m − √ m/n)C nm · δ(n + m − g ) indicate the effective squeezing and acceleration contributions, respectively. The δ-function in F matrix from the averaging process describes the conservation of energy in the DCE process. Therefore, for frequency from zero to the modulation frequency ω d , with a general solution formulaâ out (ω) = S 11 (ω) · a in (ω) + S 21 (ω) · [â in (ω d − ω)] + , the spectrum of output radiation n out (ω) =< a + out (ω) a out (ω) > can be numerically solved as with averaged thermal input photon density n in (ω, [21,25] The first two terms are of classical origin and the last term is a purely quantum mechanical effect which originates from the vacuum quantum fluctuations, or in other words, from the DCE. We note that the vacuum fluctuation dominates the output flux with n in (ω, T ) ≈ 0 at low temperature.

Radiation Spectrum at a Different Modulation Frequency
Consider the DCE process at zero temperature. For a general modulation frequency withg = 2s in the ABBA configuration or g = 2s + 1 in the AABB configuration (s is an arbitrary integer number), both DCE contributions F nn I (ω d ) and F nm I I (ω d ) are zero. In this case, the equations have a trivial solution:â out (ω) = [1 − n 2κχ n (ω)] ·â in (ω). The zero intensity S 21 (ω) indicates that the radiated photons are generated only from classical thermal fluctuations with resonance enhancement, and there are no contributions from quantum vacuum fluctuations.
For cases with nonzero F matrix element, the JTL cavity produces nontrivial radiation spectrum with peaks occurring near some typical eigenfrequencies whose indexes are strongly connected to the index pairs (n, m) with nonzero F nm . The condition when there is a nonzero contribution from the squeezing term F nn I is that g must be an even integer with g = 2n, which leads to resonant excitation of photon pairs at the nth eigenmodes and produces peaks in the radiation spectrum near the frequency ω n,0 . Accompanied with the preceding analysis of the amplitude function w n based on the parity of perturbation h(x), this condition implies that that photons can be created by the squeezing effect only in the symmetric ABBA system with g = 2n K for those K-points eigenmodes, while there are no excitations arising from the squeezing effect in the anti-symmetric AABB system.
A similar argument can be applied to the resonant excitation due to the acceleration term F nm II . Photons in the two modes indexed n and m can be created/annihilated simultaneously for g = m + n. In the ABBA system, since only those modes satisfying m ± n = 2n k have nonzero C nm matrix elements, there exist two possibilities for the resonant excitation with an even g parameter with g ࣙ 2N c . In the first case with g ࣔ 2n K , only specific mode pairs indexed m = (g/2) + n K and n = (g/2) − n K can be excited. While in the second case with g = 2n k , all mode pairs satisfying g = m + n can be excited simultaneously, as shown in the zero-temperature radiation photon spectrum n out (ω) at different modulation frequency ω d (g parameter) in the ABBA system in Figure 3a. Meanwhile, in the AABB system, since only those pairs with m ± n = 2s +1 have nonzero C nm matrix elements, resonant excitation can occur for any odd integer value g = 2s + 1, in which condition many (n, m) pairs with different parity can be excited simultaneously, as shown in Figure 3b can also be scattered from the nth to the mth eigenmode or vice versa for g = |m -n|, which will induce output peaks at frequency larger than the modulation frequency ω d and can be negligible at low temperature and for weakly excited cases. In addition to the previous spectrum analysis, we now discuss the intensity of the output photon flux N out = 1/2π · ω d 0 dω · n out (ω) in different cases. In the symmetric ABBA system, only those with even g parameters can resonantly create photon pairs. Due to the contributions from both the squeezing term F I n K and the acceleration terms F I I g ±n,n , the excitations with g = 2n K show very different features from the rest, which exhibits a broadband radiation spectrum in Figure 3a and the peaks in N out in Figure 3c. On the other hand, for even g with gࣔ2n K , the contribution solely comes from the acceleration term F I I g /2+n K ,g /2−n K , which corresponds to a sparse radiation spectrum and a lower value in N tot . Meanwhile, in the AABB system, only those odd g parameters can resonantly create photon pairs, as shown in Figure 3b,c. Some stronger peaks exist near the K-points, which is caused by the relatively stronger eigenmode variation near the Brillouin zone boundary. As a comparison, we present in Figure 3c the result of a traditional homogeneous cavity marked as "AAAA," which has zero acceleration contribution and can resonantly radiate photon pairs near the nth eigenmode only with even g parameters g = 2n and n out (ω) = 1 16 δ L 0 2 ω d 2 h g /2,g /2 2 κ 2 · |χ g /2 (ω)χ g /2 * (ω d − ω)| 2 (with F â κ). The spectrum has a peak 1 16 δ L 0 2 ω d 2 h g /2,g /2 2 κ −2 at ω = ω n = ω d /2 which is proportional to the square of the Q-factor and has a full width at half maximum 2κ. The output photon flux N out in this AAAA system does not have nontrivial peaks, but is only dependent on the square of the modulation frequency, due to the linear relationship between F I and ω d . All these evolutions are consistent with the selective excitation conditions discussed in Section 4.2.

Scaling Law for the Radiation Spectrum
Furthermore, for the two types of resonances in the symmetric ABBA system, we find interesting scaling laws based on their different resonance sources. Here we discuss the effect of a repeated number of unit cell N c on the output photon flux N out . As we know, the decay rate of a homogeneous Fabry-Pérot cavity should be inversely proportional to the duration for photons traveling in a round trip between two boundaries. Thus, we assume Q K -factor being proportional to the length of our JTL cavity, that is, Q K ∝ 1/κ ∝ N c . [45] For different resonant excitation cases with g equal/unequal to 2n K in the symmetric ABBA system, we find totally different scaling behaviors between the output photon flux N out and the number of unit cells N c , as shown in Figure 4b.
For JTL cavities with different N c but the same boundaries, increasing N c introduces more internal eigenmodes in a given frequency region, because the number of eigenmodes in each band  a) The zero-temperature harmonic radiation photon spectra n out (ω) in the ABBA JTL configuration, with same boundary condition but a different repeated number of unit cells N c . The normalization modulation frequency index g s = g/N c is 5 and 6, respectively, b) shows the total photon flux N tot in the ABBA JTL cavity with different g s parameter, and in the AAAA system with g s = 6, with same boundary condition but different N c . In all of these calculations, we choose Q K = 2.5N c , T = 0K and solve the spectrum by using the analytical method.
is the same as the number of unit cell N c whereas the eigenfrequency of band edge mode is almost independent of N c , as shown in Figure 2. However, for those eigenmode pairs (n·N c , m·N c ) with the same eigenfrequency in different N c setups, due to the repeated property of the distribution function h(x), it is easy to verify that they have the same DCE coupling matrix F n·Nc ,m·Nc (ω d ) and resonant conditions. For gࣔ2n K in the ABBA system, only some discrete internal eigenmode pairs are selectively excited, which will not be enhanced with an increasing eigenmode density. Therefore, the only effect of increasing N c on this resonant condition is a varied coupling strength κ ∝ 1/N c , which induces a narrower resonance with a line width 2κ and a higher radiation spectrum peak that grows quadratically with the Q-factor, as shown in Figure 4a, just like that in an AAAA system. As a result, the output photon flux N out is linearly scaled with N c , as shown in Figure 4b with normalized modulation frequency index g s = g/N c = 4, 5, 8 in ABBA system and that in the AAAA system. However, for g = 2n K in ABBA system, a broadband eigenmode with all possible eigenmode pairs satisfying ω n ± ω m = ω d can be excited. Therefore, an increasing N c brings an additional enhancement factor due to the increasing density of cavity modes. The total DCE photons flux N tot is quadratically related to N c , as shown in Figure 4b with g s = 2, 6, 10. This explains different scaling behaviors between the total DCE photons flux N tot and the number of unit cells N c . The DCE photon flux N tot exhibits a linear growth for gࣔ2n K and a quadratic growth for g = 2n K with an increasing number of unit cells N c , which will be easier to be detected by increasing the unit cells number N c in a real JTL cavity. It is worth mentioning that a homogenous cavity with a vibrating boundary can also show additional growth by increasing the cavity length while keeping a constant vibration strengtḣ l /l . However, this is at the cost of increasing its vibration velocity, which is practically impossible due to the restriction of the speed of light. In our setup, however, we effectively allocate the variations into each unit, which can strongly enhance the DCE through directly increasing the sample length and is more feasible for experiments. Indeed, similar to the cavity with a vibrating mirror, the overall DCE efficiency in our system can also be significantly enhanced by increasing the perturbation strength δ L0 with N out ∝ δ L0 2 . This perturbation strength δ L0 can reach as high as 25% [27] given a sufficiently strong AC signal. Such an electrical modulation scheme is far more superior to the mechanical modulation approach. Overall, the intensity of photon flux in our setup can be detected under the existing experimental conditions, with much better feasibility than existing approaches. [19,20]

Effects of Finite Temperature and Loss
Finally, we discuss the effects of imperfections on the radiation spectrum, including finite temperature and intrinsic loss. To quantify the effect of the thermal excitation, we calculate the ratio of the radiation signal intensity of the classical to the  quantum DCE origin based on Equation (9) and plot the critical temperature point with a 1% classical signal threshold in the ABBA system with different modulation frequencies. As shown in Figure 5a, the radiations for almost all resonant excitation points are dominated by quantum DCE when the temperature is below 50 mK, which can be implemented with 3He-4He dilution refrigeration technology. Low temperature is also essential for the operation of the SQUIDs. Next, we consider the influence of loss dissipation in the realistic circuit. A dissipation factor κ α is introduced to the inputoutput theory to describe a relatively low intrinsic loss in SQUID circuits. This dissipation factor will only modify the expression of matrix χ (ω) in Equation (8) by χ n (ω) = [κ + κ a,n − i(ω − ω n0 )] −1 . Similar to the homogeneous Fabry-Pérot cavity, κ α,n is set to be proportional to the mode index n with κ α,n = κ α,K · n/N c , [45] with κ α,K defined as the dissipation factor of the first K-eigenmode. This dissipation rapidly reduces the intensity of the output photon flux, as shown in Figure 5b for a typical resonant excitation in ABBA system with g s = 6 at 10mK. Therefore, a low-loss SQUID circuit is necessary for the DCE signal to be detected.

Conclusion
To summarize, we have studied DCE in a 1D inhomogeneous JTL cavity with either symmetric or antisymmetric spatial modulation of the inductance. We have shown that the symmetry plays a very important role in selecting the DCE driving frequencies, and the corresponding radiation photon spectrum exhibits an intriguing scaling law. Our result may inspire further experimental effort for achieving enhanced DCE in inhomogeneous systems.

Material
Herein to solve the DCE efficiency we use the representative SQUID parameters [21] as: size a 0 = 450 nm, capacitances C J1 = C J2 = 45 fF, inductance L J ( dc ) = 21.7 nH and plasma frequency ω p /2π = 36 GHz, and the perturbation amplitude for the effective inductance is set to δ L0 = 0.001. Each AC segment covers ten SQUID units with size d = 10a 0 = 4.5 µm and the whole length of a typical design with repeated unit cells N c = 20 is N c D = 0.36 mm with size of each unit cell D = 4d = 18 µm.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.