Inhomogeneous Josephson junction chains for superinductance optimization

We report a theoretical study of the low-frequency impedance of a Josephson junction chain whose parameters vary in space. Our goal is to find the optimal spatial profile which maximizes the total inductance of the chain without shrinking the low-frequency window where the chain behaves as an inductor. If the spatial modulation is introduced by varying the junction areas, we find that the best result is obtained for a spatially homogeneous chain, reported earlier in the literature. An improvement over the homogeneous result can be obtained by representing the junctions by SQUIDs with different loop areas, so the inductances can be varied by applying a magnetic field. Still, we find that this improvement becomes less important for longer chains.


Introduction
Quantum engineering in superconducting nanocircuits is a rapidly developing field, due to progress in sample fabrication techniques which has been occurring in the past decade [1] Complex circuits with many elements can be routinely fabricated on a chip nowadays.Due to superconductivity, electromagnetic signals propagate in such circuits with extremely low losses, and the circuit properties can be tuned by applying an external magnetic field.Highly inductive elements are often needed in such nanocircuits, to realize a large non-dissipative impedance.Applications of large inductances include protection of fluxonium qubits from the charge noise [2], tunable microwave impedance matching [3], or a potential implementation of the electrical current standard in quantum metrology based on Bloch oscillations [4,5,6].
Because any geometrical inductor (a coil being the standard textbook example) also necessarily possesses a parasitic self-capacitance which starts to dominate at high frequencies, its non-dissipative impedance is limited by the vacuum impedance, ∼ µ 0 / 0 = 4αR Q , where α ≈ 1/137 is the fine structure constant, and R Q ≈ 13 kΩ is the resistance quantum [7].Indeed, the inductance of a geometrical inductor is due to the magnetic field produced by the current, which acts on the current itself.The relativistic nature of this effect is the intrinsic reason for its weakness.This limitation can be overcome by using superconducting materials whose inductance is due to the kinetic energy of the Cooper pair condensate [8], and thus is of non-relativistic origin.The term "superinductance" Correspondence to: denis.basko@lpmmc.cnrs.fr is often used to denote such superconductivity-based inductance.
Several structures, based on Josephson junctions (JJs), have been reported to work as superinductors [9,10].In the first one, a large inductance was obtained by putting N Josephson junctions in series, which gave the total inductance N L (L is the inductance of a single junction).In Ref. [10], magnetic-field-induced frustration was used to increase the inductance, which then exhibited a strong nonlinearity.Here, we focus on the linear case, and analyze structures analogous to that of Ref. [9].
A simple strategy to increase the total inductance of a JJ chain would then be to make L and/or N as large as possible.However, in either case one faces some limitations.In the first case, the JJ inductance L is inversely proportional to the Josephson energy of the junction, E J = ( /2e) 2 (1/L).To work as an inductor, the junction must be in the superconducting regime, E J E C , where the charging energy E C = (2e) 2 /(2C) is determined by the junction capacitance.This condition sets a lower limit on E J , or, equivalently, an upper limit L < L max , or a lower limit on the junction area A, as both E J , C ∝ A.
Limitations on the junction number N arise from the dependence of the chain response on the frequency ω.The phase slip rate, although exponentially suppressed for E J E C , grows with N , giving rise to a finite dc resistance, which spoils the purely inductive response of the chain at low frequencies.From the high-frequency side, the effective bandwidth of the inductive response is restricted by electromagnetic modes supported by the chain, ω ω 1 (the lowest mode frequency).Crucially, besides the capacitance C of the junction between neighboring supercon- LC acoustic-like region flat region ducting islands, each island has a small capacitance C g to the ground.This capacitance gives rise to screening of the Coulomb interaction between the islands on a length scale λ = C/C g and produces an acoustic-like region of the mode dispersion ω(q) = (LC) −1/2 2 (q)/[ 2 (q) + λ −2 ] of spatially homogeneous chains, where (q) = 2 sin(q/2), and q is the wavenumber, 0 q π (Fig. 1).The first mode corresponds to q = π/(N + 1), so for large N πλ, the frequency of the lowest mode ω 1 ∝ 1/N , and the inductive response bandwidth shrinks with increasing N .This was the main limitation for the device studied in Ref. [9], where a special effort was made to decrease the parasitic ground capacitance C g .
The above argumentation works for spatially homogeneous chains, whose total inductance is determined by just two parameters, the single-junction inductance L and their number N , if L is assumed to be the same for all junctions.This, however, need not be the case, since an arbitrary spatial profile of junction sizes along the chain can be produced during the sample fabrication.A spatial modulation of junction parameters modifies the normal modes of the chain, and can manifest itself in various situations.For example, Josephson energy renormalization by coupling to the normal modes was shown to be affected by a modulation of the chain parameters [11].Effect of the normal mode structure on dephasing of the fluxonium qubit was discussed in Ref. [12].For the present problem, one can try to optimize the total inductance and the operation bandwidth of the chain using many more degrees of freedom than just L and N , because the parameters of each of the N junctions can be treated as optimization variables.To study, whether one can take advantage of this large number of variables and improve the homogeneous chain result of Ref. [9] by carefully choosing the spatial profile of the junction parameters, is the purpose of the present work.
In this paper, we consider two ways to introduce a spatial inhomogeneity into the structure.One is to vary the area A n of each junction n (assuming the island area to be already optimized to minimize the ground capacitance

S n
A n A n Fig. 2. A schematic representation of a SQUID (top view): two superconducting islands, top and bottom, are connected by two junctions forming a loop.At zero magnetic field, the SQUID inductance Ln(0) is determined by the total area of the junctions An, shown by hatching.When a magnetic field B is applied, the inductance Ln(B) = Ln(0)/| cos(πBSn/Φ0)| is determined by the magnetic flux BSn through the SQUID loop area Sn, represented by the white circular region in the center.
as was done in Refs.[2,9,13]).This leads to a simultaneous variation of the junction inductances L n and capacitances C n , such that their product L n C n = const.Optimizing over all areas {A n }, we find that the best result is still achieved for a homogeneous configuration.
The second way to introduce a spatial variation of the junction parameters is to represent each junction by a SQUID (superconducting quantum interference device).When subject to a magnetic field B, a SQUID behaves like an effective Josephson junction with a field-dependent Josephson energy E J (B) = E J (0)| cos(πBS/Φ 0 )|, where Φ 0 = 2π /(2e) is the superconducting flux quantum, and S is the SQUID loop area which determines the magnetic flux BS through the SQUID (Fig. 2).Then, if all SQUIDs have different areas S n , the inductance of each junction of the chain, L n (B) = L n (0)/| cos(πBS n /Φ 0 )|, varies in space, and this variation is independent of the variation of the capacitance C n (the latter is controlled by the junction area A n , independent of the loop area S n ).In this case, we show that one can indeed improve over the homogeneous result, by placing SQUIDs with larger loop area (higher inductance) near the ends of the chain.Still, the obtained improvement over the homogeneous result turns out to decrease with the increasing chain length.
The paper is organized as follows.In the next section we specify the model and formally pose the optimization problem.In Sec. 3 we analyze the case when only the junction areas A n vary in space.In Sec. 4 we study variation of the SQUID loop areas S n .In Sec. 5 we give our conclusions.

Formal setting of the optimization problem
We consider a chain of N + 1 superconducting islands.Each island is connected to its nearest neighbors by Josephson junctions, so the chain has N junctions (Fig. 3).We assume N 1.When the junctions are in the superconducting regime, E J E C , the oscillations of the superconducting phase ϕ n on each island are small.Then, the Josephson current through the nth junction from island n to (n + 1) can be written as is the junction critical current.Thus, the voltage drop across the junction can be determined by using the Josephson relation: This expression shows that the junction behaves as a linear inductor, and the Josephson kinetic inductance is given by L n = /(2eI c n ).In addition, we assume that the dissipation is very small and can be neglected.Then, an isolated chain is equivalent to the electric circuit shown in Fig. 3(b), where C n is the capacitance formed by the neighboring superconducting islands, and C g n is the capacitance of each island to ground.We define the complex impedance Z(ω) of the chain at frequency ω as the ratio of the voltage V ω e −iωt on an external ac voltage source, connected to the islands n = 1 and n = N + 1, to the current I ω e −iωt through this source (Fig. 3).
To determine the normal mode frequencies of this circuit, one can apply the Kirchhoff's law at each of the N +1 nodes of this circuit.This gives the following system of linear equations for the voltages V n : where the junction admittance is defined as System ( 2) can be written in the matrix form, Υ (ω)V = I, in terms of the column vectors ) and I T = (I ω , 0, . . ., 0, −I ω ), as well as the corresponding matrix Υ (ω).Then the chain impedance Z(ω) can be expressed in terms of the matrix elements of the inverse . At low frequencies, the admittances are dominated by the inductive part, so the impedance is given by Z(ω → 0) = −iωL tot , where L tot is the total inductance of the chain, The approximation Z(ω) ≈ −iωL tot is valid as long as ω ω 1 , where ω 1 is the lowest normal mode frequency, for which det Υ (ω) = 0.
As discussed in the introduction, ideally one would like to increase both L tot and ω 1 , but these two requirements are in conflict.Thus, one can try to maximize L tot at fixed ω 1 , or maximize ω 1 while keeping L tot fixed.We prefer the second option, as the constraint expressed by Eq. ( 4) is much easier to resolve than the constraint ω 1 = const.Thus, our optimization problem is formulated as follows: find the spatial profile of L n , C n , C g n which maximizes ω 1 while keeping L tot fixed.To complete the formulation of the problem, we have to specify the independent variables over which the optimization is performed.
The shape and size of the superconducting islands and of the junctions between them can be well controlled in the fabrication process.It is easy to notice that while the parameters L n , C n are mostly determined by the junction areas, the parasitic ground capacitances C g n are mostly determined by the island sizes.Thus, the first obvious step is to minimize the island sizes as much as possible while keeping constant the junction areas, as any part of the island area which does not participate in the junctions, does not contribute to the inductance, but decreases ω 1 .This optimization was performed in Refs.[2,9,13,14].Then, the ground capacitance of the nth island becomes a function of the areas of the junctions in which it particpates, n − 1 and n.This function was calculated numerically in Ref. [9], and the resulting dependence resembles a weak power law or a logarithm.We will assume that this first optimization step has been performed.Then, our first setting corresponds to independent variation of all junction areas, which are allowed to vary in a certain range.In the fabrication process, quite a wide range of sizes can be achieved, and the restriction on the areas rather comes from physical considerations.One restriction is that for too small areas, the condition E J E C is violated, and then the classical description of small phase oscillations is no longer valid.Indeed, the amplitude of a quantum phase slip, ∝ e −(2/π)R Q /Zn [15] is exponentially suppressed only for small junction impedances, Z n ≡ L n /C n R Q (we remind that R Q ≈ 13 kΩ denotes the resistance quantum), so too large impedances are not allowed.The junction impedance is inversely proportional to its area, so the area cannot be made too small.On the other hand, if the junction area is too large, the junction can no longer be treated as a zero-dimensional object, because the frequency of its own electromagnetic modes becomes too low.
Let us choose the smallest allowed junction area as the unit of area.Then the largest allowed junction area A max 1 is an independent dimensionless parameter of the problem.The junction inductance and capacitance at the smallest area, L max and C min , can be chosen as the units of inductance and capacitance, respectively.Thus, we have N dimensionless variables A n , allowed to vary in the range They determine the inductance and the capacitance of each junction as and Eq. ( 4) thus imposes a constraint on the set {A n }.
In this case, the plasma frequency of each junction is un- Finally, for the ground capacitances we use a simple form where g(x) is some function, growing sublinearly with x (a power law or a logarithm).All qualitative arguments given below are not sensitive to the specific dependence g(x); in the numerical calculations, we set g(x) = √ x, as mentioned in Ref. [14].To define Eq. (5c) at the ends, we set A 0 ≡ A 1 , A N +1 ≡ A N .Thus, the first optimization problem is fully defined as maximization of ω 1 determined from Eqs. (2), whose coefficients are expressed by Eqs.(5b) and (5c) in terms of the dimensionless areas A n .The optmization variables are the areas A n in the allowed range (5a) and subject to constraint (4), as well as the number of the junctions N itself.Note that constraint (4) and inequalities (5a) restrict the number of junctions N to the interval The second way of producing a spatial variation of the JJ chain parameters is to replace each junction by a SQUID.Each SQUID is characterized by its loop area S n , independent of the junction area A n (Fig. 2).By applying a magnetic field B, one can change the SQUID inductance as where the zero-field inductance L n (0) is determined by the junction area A n .This way of tuning the properties of the JJ by magnetic field is routinely used in experiments (see, e. g., Ref. [16]).Here, it is crucial for us that the spatial variation of inductance is independent of that of capacitance, which was not the case in the previous model, since in Eq. (5b) the product L n C n remained fixed.Thus, instead of the optimization problem defined by Eqs.(5a)-(5c) via variables A 1 , . . ., A N , we consider another problem defined via variables F 1 , . . ., F N : 1 All junction areas are assumed to be the same, A n = 1, so the plasma frequency of each SQUID is modulated as

and each variable F n represents the ratio
where Φ max is some maximal magnetic flux allowed to pierce the SQUID loops in order for the device to remain in the superconducting regime E J E C .Clearly, {F n } are independent variables, because {S n } are independent, and additional freedom is introduced by the magnetic field.Just like before, the only constraint on F n is Eq. ( 4), and it restricts the chain length N to the interval The two optimization problems, defined by Eqs.(5a)-(5c) and by Eqs.(8a)-(8b), will be studied in the next two sections, respectively.

Junction area modulations
Before we proceed with optimization for inhomogeneous JJ chains, it is useful to see what can be achieved in the homogeneous case, for future reference.For the problem (5a)-(5c), with all A n = A, we have only two variables, A and N .Constraint (4) fixes It is convenient to denote the first mode frequency for this homogeneous chain by Ω N .It is given by [9] for N 1.This is a decreasing function of x ≡ N/N 0 for any g(x) growing slower than linearly with x.Thus, ω 1 is maximized by taking N = N 0 , all L n = L max .We denote the corresponding value of ω 1 by Ω N0 .
To improve this result using an inhomogeneous chain, one should take some N > N 0 [a smaller one would be incompatible with the constraint (4)], and hope that the gain in ω 1 from the inhomogeneiety would overcome the loss due to the length increase.A qualitative idea of the best spatial profile A n can be obtained from the perturbation theory for system (2), developed in Ref. [17].Let us use the homogeneous chain of length N with all A n = A = N/N 0 and the first mode frequency Ω N as the zero approximation.If we now modify each junction area by a small amount ∆A n , the first-order frequency shift is given by [17] The dependence of α n on n is quite simple (sin 2 + const), and α n is the largest for n = (N + 1)/2, in the middle of the chain.The value at the maximum α (N +1)/2 > 0 as long as [2Ag (A)/g(A)] sin 2 [π/(N + 1)] < 1, which is the case for any sublinear g(x) and N > 4. Thus, the center of the chain contributes the most to the increase of ω 1 .
Let us take N = N 0 + 1.Then, the largest increase of the areas near the center, allowed by constraint (4), is obtained by keeping N 0 −1 junctions with A n = 1, and two more junctions with A n = 2, to be put in the center.(Note that it is impossible to keep N 0 junctions with A n = 1, as the constraint would require the remaining one to have A n = ∞).As the area change for the central junctions is not small, the perturbative Eq. (12a) is not sufficient to describe this situation.Still, ω 1 for this structure can be found analytically.The result of this straightforward but bulky calculation, given in Appendix A:, is that the resulting frequency is always smaller than Ω N0 .
The full optimization of all junction areas {A n }, subject to constraint (4), can be performed numerically.For any N > N 0 , we maximize ω 1 as a function of all the areas, calculated numerically from the eigenvalue equation det Υ (ω) = 0.The resulting maximum ω 1 is plotted versus N in Fig. 4 for several values of C min /C g min and A max .The analytical result of Appendix A: shows that the curve starts to bend down at N = N 0 + 1, and the numerics shows that the same trend is followed for all N .Thus the optimal ω 1 at N > N 0 is always below the best value for the homogeneous chain, Ω N0 .In Fig. 5 we show the optimal spatial profile {A n }, corresponding to one of the points in Fig. 4. Indeed, the best ω 1 for a fixed N is obtained by placing the largest junctions in the middle of the chain.Still, the resulting gain in ω 1 is smaller than the loss due to the increase of the chain length from N 0 to N .

SQUID loop area modulations
As in the previous section, we start by a straighforward study of the homogeneous case.Constraint (4) fixes F = N/N 0 , so for the chain with 11) we have where λ = C g min /C min is the screening length, which does not depend on N and {F n } for problem (8b)-(8a).Expression ( 13) reaches maximum at N = πλ, so we have to consider three cases for the position of this value with respect to the interval (10).
(i) In the case πλ > F max N 0 , the frequency is maximized by taking the longest possible chain, N = F max N 0 .This case corresponds to the regime when for all allowed N the first mode is on the flat part of the mode dispersion curve (Fig. 1).This means that we have demanded a value of L tot which is too small; a larger inductance can be obtained by simply increasing the length at almost no cost in ω 1 .So, this case has no practical relevance.
(ii) When N 0 πλ F max N 0 , the frequency is maximized at N = πλ.This corresponds to the first mode frequency roughly at the boundary between the flat part of the mode dispersion curve and its acoustic part.
(iii) In the case πλ < N 0 , the frequency is maximized by taking the shortest possible chain.This regime corresponds to demanding such a large inductance L tot that the first mode necessarily belongs to the acoustic part of the dispersion curve.This is the regime where the competition between L tot and ω 1 is the most severe; it is in this regime that a gain in ω 1 by introducing a spatial variation of F n would be the most interesting for practical purposes.
The perturbation theory in small modulations ∆F n with respect to a homogeneous chain with N junctions gives a result, similar to Eq. (12a): which again tells us that inductance modulations in the center of the chain contribute the most to the increase in ω 1 .As in the previous section, we now consider a chain of length N = N 0 + 1 with inductances of two junctions in the center smaller by a factor F = 2.The explicit calculation of given in Appendix A: shows that this chain has ω 1 > Ω N0 , and thus one can indeed improve over the homogeneous result.However, for long chains, N 0 πλ, the gain is quite small: Is it possible to gain more in ω 1 by choosing a chain length N significantly exceeding N 0 ?As a trial spatial profile, let us consider a long chain with a central region of length N − 2N 1 1 where the inductances are smaller by a factor F than in the surrounding (although this piecewise profile does not coincide with the true optimal one, found numerically below, it allows for a simple analytical solution): Constraint ( 4) then fixes For N − 2N 1 1, we can study the problem in the continuum limit, replacing the junction number n by a continuous variable x.In addition, let us focus on the most interesting case of long chains N 0 πλ, then one can approximate the mode dispersion by the acoustic one, ω(q) ≈ q/ √ LC g .Then, Eqs. ( 2) are transformed into the Helmholtz equation with von Neumann boundary conditions at the ends of the chain, (18) For the piecewise function L(x), given by Eq. ( 16), and for a given frequency ω, the wavenumbers in the outer regions and in the central region are given by q = ω L max C g min and by q/ √ F, respectively.Thus, taking advantage of the symmetry of L(x) with respect to x → N − x, we seek V (x) in the form (the first mode is odd) The requirement of continuity of V and (1/L)(∂V /∂x) at x = N 1 , N − N 1 yields the following equation for q: For all F > 1, upon increasing N 1 from 0 to N 0 /2 (that is, upon decreasing N from N 0 F to N 0 ), the solution monotonically rises from q = π/(N 0 √ F) to q = π/N 0 (Fig. 6), the highest frequency being achieved in the shortest homogeneous chain.This means that in the limit N 0 πλ the gain in ω 1 is so small that it is not captured by the acoustic approximation.
To check these considerations numerically, we perform the full optimization of all {F n }, subject to constraint (4).As in the previous section, for any N > N 0 , we maximize ω 1 as a function of all the areas, calculated numerically from the eigenvalue equation det Υ (ω) = 0.The resulting maximum ω 1 is plotted versus N in Fig. 7  values of λ and F max .The optimal spatial profile of the inductance is shown in Fig. 8; as in the previous section, it corresponds to putting the small-inductance junctions in the middle of the chain, and the large-inductance ones near the ends.From the analytical arguments above, we do not expect the first mode frequency for the optimal inhomogeneous chain of optimal length to be much larger than for the shortest homogeneous chain.This is checked numerically in Fig. 9(a), where we plot the two frequencies as a function of N 0 (we remind that at fixed F max , N 0 parametrizes the desired total inductance).For long chains, the improvement due to spatial modulation is indeed negligible.The optimal length of the modulated chain is close to N 0 at large N 0 (up to a constant offset), as shown in Fig. 9(b).

Conclusions
In this work, we explored the possibility to optimize the frequency range where a JJ chain can work as a superinductor, by a careful choice of the spatial profile of the junction parameters.In the case when junction areas are varied, the best result is still obtained for a spatially homogeneous chain, as in Ref. [9].Another way to introduce a spatial variation is to represent the junctions by SQUIDs whose loop areas are different.Then, by applying a magnetic field, one can vary the junction inductance independently from its capacitance.We show that this strategy can indeed give an improvement with respect to the homogeneous case, if the most inductive junctions are placed near the ends of the chain, and the least inductive ones in the middle.Still, we find that this improvement becomes less important for longer chains.The qualitative difference between the cases of junction area and SQUID loop area modulations stems from the fact that in the first case, the plasma frequency of each junction, 1/ √ L n C n , remains fixed.In the second case, the junction inductance can be decreased independently from the capacitance, which leads to an increase of the local plasma frequency, and to a certain degree increases the overall frequency scale.The N0 dependence of the chain length Nopt, at which the optimal value of ω1 is obtained.On both panels the blue circles and the red squares correspond to Fmax = 2 and 10, respectively, and we took λ = 20.
We are grateful to W. Guichard and F. Hekking for stimulating discussions and critical reading of the manuscript.We also acknowledge support from the European Research Council (Grant No. 306731).with yet unknown q which will be determined by matching the solutions in the middle of the chain.Note that as q is related to the frequency by the dispersion relation (A.3), which is a monotonically increasing function, it is sufficient to check whether the value of q, obtained by matching the solutions, is larger or smaller than the one corresponding to the shortest homogeneous chain, q 0 = π/(N 0 + 1) = π/(2N 1 ).

Fig. 1 .
Fig. 1.Dispersion curve of a JJ chain with C/C g = 25 (solid curve) and modes of a chain with N = 100 junctions (filled circles).The mode frequency ω is measured in the units of the junction plasma frequency 1/ √ LC.

Fig. 3 .
Fig. 3. (a) A schematic view of the Josephson junction chain and its impedance definition.(b) Linear circuit, equivalent to the chain shown in (a), described by Eqs.(2).

Fig. 4 .
Fig. 4. The first mode frequency ω1 (in units of the plasma frequency ωp ≡ 1/ √ LmaxCmin) obtained by full numerical optimization of all junction areas {An}, subject to constraint (4).We take N0 = 25 for all curves, while λ 2 ≡ Cmin/C g min = 400 and 16 for panels (a) and (b), respectively.Two values of Amax = 3 and 10 were chosen, shown by the blue and red symbols (lower and upper curves), respectively, on each panel.The solid curve shows ΩN , the first mode frequency for the homogeneous chain with Ln = LmaxN0/N , Cn = CminN/N0, C g n = C g min N/N0, and the dashed horizontal line shows the best homogeneous result ΩN 0 .

Fig. 7 .
Fig. 7.The first mode frequency ω1, in the units of the plasma frequency ωp ≡ 1/ √ LmaxCmin, obtained by full numerical optimization of all {Fn}, subject to constraint (4), shown by symbols for Fmax = 2 and 10 (blue circles and red squares, respectively).The solid curve shows the first mode frequency ΩN for the homogeneous chain with Ln = LmaxN0/N , Cn = Cmin, C g n = C g min .We take N0 = 25 for all curves, while λ = 20 and 4 for panels (a) and (b), respectively.The dashed horizontal line shows the best homogeneous result.

Fig. 9 .
Fig. 9. (a) The N0 dependence of the first mode frequency for the optimal inhomogeneous chain having the optimal length (symbols) and for the shortest homogeneous chain with N = N0 (solid curve).The frequencies are measured in the units of the plasma frequency ωp ≡ 1/ √ LmaxCmin.(b)The N0 dependence of the chain length Nopt, at which the optimal value of ω1 is obtained.On both panels the blue circles and the red squares correspond to Fmax = 2 and 10, respectively, and we took λ = 20.

Fig. 10 .
Fig. 10.A schematic representation of a JJ chain with two modified central junctions.