# Theoretical study of dissociative recombination and vibrational excitation
of the BF${}_{2}^{+}$ ion by an electron impact

V. Kokoouline1<EMAIL_ADDRESS>M. Ayouz2 J. Zs Mezei3,4,5
<EMAIL_ADDRESS>K. Hassouni4 I. F. Schneider5,6 1Department of Physics,
University of Central Florida, Orlando, FL 32816, USA  2LGPM, CentraleSupelec,
Université Paris-Saclay, 91190 Gif-sur-Yvette, France  3HUN-REN Institute for
Nuclear Research (ATOMKI), H-4001 Debrecen, Hungary  4LSPM CNRS-UPR3407,
Université Paris 13, 93430 Villetaneuse, France  5LOMC-UMR6294, CNRS,
Université Le Havre Normandie, 76600 Le Havre, France  6LAC-UMR9188, CNRS
Université Paris-Saclay, F-91405 Orsay, France

###### Abstract

Cross sections for dissociative recombination and electron-impact vibrational
excitation of the BF${}_{2}^{+}$ molecular ion are computed using a
theoretical approach that combines the normal modes approximation for the
vibrational states of the target ion and the UK R-matrix code employed to
evaluate electron-ion scattering matrices for fixed geometries of the ion.
Thermally-averaged rate coefficients are obtained from the cross sections for
temperatures $10-3000$ K.

###### pacs:

33.80. -b, 42.50. Hz

††preprint: APS/123-QED

Non-equilibrium plasmas produced with electrical discharges in BF3 containing
feed gas present a continuously increasing interest for a large number of
applications. In particular, BF3 is very often the boron carrier when plasmas
are used for material processing Sennikov et al. (2017); Torigoe et al.
(2016); Gonzatti et al. (2010). Basically, BF3 plasmas are used either for (i)
the synthesis of ultra-hard boron compounds, e.g. boron carbides Sennikov et
al. (2017), (ii) the deposition of boron nitride, an advanced material with a
large number of functionalities Torigoe et al. (2016), and (iii) p-type doping
by boron in the semi-conductor and photovoltaic industries Gonzatti et al.
(2010); Duchaine et al. (2012). As far as doping applications are concerned,
plasma immersion ion implantation (PIII) processes are probably one of the
most promising in terms of cost and technical performance Duchaine et al.
(2012); Young et al. (2016). These processes make use of very low pressure,
very high density magnetized plasmas generated in a BF3 containing feed gas.
Depending on the level of the power coupled to the plasma, several cations -
BF${}_{3}^{+}$, BF${}_{2}^{+}$, BF+, B+ \- and anions - e.g. F- \- may be
produced M. Farber and Srivastava (1984); Yong-Ki and Irikura (2000). The
positive ions are then extracted from the source to an implantation chamber
where the processed silicon substrate is submitted to very high negative
voltage pulses. These pulses result in a large acceleration of the positive
ions that are implanted in the substrate, which results in the doping of this
later. The implantation depth depends on the nature and the energy
distribution of the ions impinging the substrate, while the doping level
depends on the ion flux and, consequently, on the plasma density. The plasma
density and the relative predominance of the different ions are determined by
the ionization kinetics in the source region. The plasma sources used in PIII
processes are usually magnetized Stewart and Lieberman (1991). The ambipolar
diffusion in the radial direction that is perpendicular to the magnetic field
is strongly reduced with respect to the parallel diffusion. Under typical
plasmas conditions used in PIII processes, i.e., $B\approx 100-500$ Gauss and
$p\approx 1$ Pa, the electron cyclotron frequency is the GHz range, while the
electron-neutral momentum transfer collision frequency is in the MHz range,
according to Eq. (5.4.5) of Lieberman and Lichtenberg (2005) the diffusion
coefficient in the radial direction is reduced by approximately $6$ orders of
magnitude with respect to the parallel diffusion coefficient. The
characteristic time estimated for perpendicular diffusion for typical PIII
reactor geometries is at the order of a hundred second, while the dissociative
recombination characteristic time is around $0.1$ s. Following a similar
procedure, one can easily show that mutual neutralization is also likely to
dominate diffusion losses provided the negative ion density is of the same
order of magnitude as electron density, i.e. moderately electronegative
plasmas. It appears therefore that under PIII discharge conditions, positive
ions losses at the reactor wall are dominated by their dissociative
recombination with electrons and by their neutralization by collisions with
negative ions. The investigation of dissociative recombination of molecular
ions is therefore of major interest for these processes. This is especially
the case of BF${}_{2}^{+}$, which is often the major ion in BF3 containing
plasma in discharge conditions corresponding to PIII process Young et al.
(2016). This study is a continuation of a previous work performed in the
theoretical framework of the multichannel quantum defect theory on the
dissociative recombination and competitive processes of BF+ Mezei et al.
(2016).

The article is organized as follows. After the Introduction, Section II
describes the theoretical approach used in the present calculation. In Section
II, the obtained cross sections and the corresponding rate coefficients are
displayed and discussed. Section III concludes the study.

##  I Theoretical approach

###  I.1 Dissociative recombination and vibrational excitation cross section
formulas

Since the basic formalism used in our model is presented in detail in Refs.
Douguet et al. (2011); Fonseca dos Santos et al. (2014), we restrict ourselves
here to underline its major ideas. The studied molecular cation,
BF${}^{+}_{2}$, is linear in its equilibrium geometry.

Table 1: Vibrational frequencies (in cm-1) obtained in this study and compared with previous data available in literature. Mode | symmetric stretching | bending | asymmetric stretching |  |
---|---|---|---|---|---
Symmetry | $\Sigma_{g}^{+}$ |  $\Pi_{u}^{+}$ and $\Pi_{u}^{-}$ | $\Sigma_{u}^{+}$ |  |
| $\omega_{1}$ | $\omega_{2}$ | $\omega_{3}$ | Method | Reference
| 1062.8 | 469.2 | 2146.1 | CI/ cc-pVTZ | This work
| 1023 | 443 | 2088 | MP2 |  Ref. Pyykko and Zhao (1990)
| 1030 | 450 |  | CI |  Ref. Perić and Peyerimhoff (1993)
|  |  |  2026.1$\pm$ 0.2 | Exp. |  Ref. Jacox and Thompson (1995)

The theoretical model starts with the following assumptions (see for example
Fonseca dos Santos et al. (2014)): (i) the rotation of the molecule is
neglected, (ii) the cross section is averaged over the autoionizing
resonances, (iii) the autoinization lifetime is assumed to be much longer than
the predissociation lifetime , and (iv) the harmonic approximation is used to
describe the vibrational state of the core ion. Using (i)-(iv) and applying
the frame transformation, the DR cross section is given by Eq.(13) of Fonseca
dos Santos et al. (2014), in which the scattering matrix elements were
expanded to first order in the normal coordinates. The cross section for
vibrational excitation (VE) of the mode $i$ writes

$\displaystyle\sigma^{VE}_{i}(E_{el})=\frac{\pi\hbar^{2}}{4mE_{el}}g_{i}\sum_{ll^{\prime}\lambda\lambda^{\prime}}\left|\frac{\partial
S_{l\lambda,l^{\prime}\lambda^{\prime}}}{\partial{q}_{i}}\right|^{2}\theta(E_{el}-\hbar\omega_{i})\,.$
(1)

Here $q_{i}$, $\hbar\omega_{i}$ and $g_{i}$ $(i=1-3)$ are respectively the
dimensionless coordinate, the energy and the degeneracy of the mode $i$. The
degeneracy is two for the bending mode 2 and one for the symmetric and
asymmetric stretching modes. $S_{l\lambda,l^{\prime}\lambda^{\prime}}$ is an
element of the fixed-nuclei scattering matrix for electron–BF${}_{2}^{+}$
collisions with the initial channel ($\lambda l$) and the exit channel
($\lambda^{\prime}l^{\prime}$), $l$ being the electron angular momentum and
$\lambda$ its projections on the molecular axis. And finally, $m$ is the
reduced mass of the electron-ion system, $E_{el}$ the incident energy of the
electron, $\theta$ the Heaviside step function. The present theoretical
approach can describe the (de-)excitation process changing only one quanta in
each normal mode of the target ion. (De-)excitation cross sections for the
changing two or more quanta in a mode, are much smaller (the propensity rule)
and are neglected in this study. In the present case, the initial state of the
ion is the ground vibrational level, so the electron can only be captured into
the first excited vibrational state of each normal mode.

The situation is similar if the electron energy is not sufficient to excite
the ion and then to leave it. In such a situation, the present model suggests
that the excitation probability of ion by the electron is described by the
same physics, but instead of leaving the vibrationally excited ion, the
electron is captured in a Rydberg resonances attached to that vibrational
state, excited by the electron. If the electron is captured by the ion, the
system will most likely dissociate, rather then autoionize. Correspondingly,
the cross section for dissociative recombination (DR) is then obtained Fonseca
dos Santos et al. (2014) as

$\displaystyle\sigma^{DR}(E_{el})=\frac{\pi\hbar^{2}}{4mE_{el}}\sum_{i=1}^{3}g_{i}\sum_{ll^{\prime}\lambda\lambda^{\prime}}\left|\frac{\partial
S_{l\lambda,l^{\prime}\lambda^{\prime}}}{\partial{q}_{i}}\right|^{2}\theta(\hbar\omega_{i}-E_{el})\,.$
(2)

Here $i$ runs over all three modes: two stretching modes $\nu_{1}$ and
$\nu_{3}$ (symmetric and asymmetric) with respective frequencies $\omega_{1}$
and $\omega_{3}$ and corresponding coordinates $q_{1}$ and $q_{3}$, and a
doubly degenerate transverse mode $\nu_{2}$ with a lower frequency
$\omega_{2}$ and coordinates $(q_{2x},q_{2y})$.

To calculate the derivative of the scattering matrix $\partial
S_{l\lambda,l^{\prime}\lambda^{\prime}}/\partial q_{i}$ with respect to the
normal coordinate $q_{i}$, the scattering matrix is evaluated for two values
of $q_{i}$ keeping the other normal coordinates $q_{i^{\prime}}$ fixed at
$q_{i^{\prime}}=0$.

The elements $S_{l\lambda,l^{\prime}\lambda^{\prime}}(\vec{q})$ of the
scattering matrix $\hat{S}(\vec{q})$ at a given geometry $\vec{q}$ specified
by four normal coordinates
$\vec{q}=\\{q_{1},q_{2},q_{3},q_{4}\\}=\\{q_{1},q_{2x},q_{3},q_{2y}\\}$ are
computed from the reactance matrix $\hat{K}$, obtained numerically as
discussed below:

$\displaystyle\hat{S}=\frac{\hat{1}+i\hat{K}}{\hat{1}-i\hat{K}}\,,$ (3)

where $\hat{1}$ is the identity matrix.

###  I.2 The properties of the BF${}_{2}^{+}$ ion and the scattering
calculations

The main electronic ground state configuration ${}^{1}\Sigma_{g}^{+}$ of the
ion in its natural point group symmetry $D_{\infty h}$ is
$\left(1\sigma_{u}^{+}\right)^{2}\left(1\sigma_{g}^{+}\right)^{2}\left(2\sigma_{g}^{+}\right)^{2}\left(2\sigma_{u}^{+}\right)^{2}\left(3\sigma_{g}^{+}\right)^{2}\left(4\sigma_{g}^{+}\right)^{2}\left(3\sigma_{u}^{+}\right)^{2}\left(1\pi_{u}^{-}\right)^{2}$
$\left(1\pi_{u}^{+}\right)^{2}\left(1\pi_{g}^{-}\right)^{2}\left(1\pi_{g}^{+}\right)^{2}$.
The normal coordinates and the related frequencies are obtained using the cc-
pVTZ basis set centered on each atom and including $s$, $p$ and $d$ orbitals.
Performing configuration interaction (CI) calculation in the $D_{2h}$ symmetry
group, using the MOLPRO suite of codes Werner et al. (2012), we found
equilibrium geometry of the ion with the B-F distance of $1.2215$ Å. Table 1
gives the vibrational frequencies obtained from CI calculations using cc-pVTZ
basis set, and compares with the data available in literature.

Figure 1:  Eigenphase sums as a function of the electron scattering energy
$E_{el}$ for $q_{i}=0.01$ (dimensionless) for the symmetric stretching (left
panel), bending (right upper panel), and asymmetric stretching (right bottom
panel) modes. The curves of different color correspond to different symmetries
of the $e^{-}$–BF${}_{2}^{+}$ system.

The electron scattering calculations were performed using the UK R-Matrix code
Tennyson (2010); Carr et al. (2012) with the Quantemol-N interface Tennyson et
al. (2007). The calculations were performed in the abelian subgroup $D_{2h}$
and the target BF${}_{2}^{+}$ ion was assumed to be in its ground electronic
state. In our CI model, we freeze 10 electrons in the core
$1a_{g},2a_{g},3a_{g},1b_{2u},2b_{2u}$, while the remaining 12 electrons are
kept free in the active space of
$4a_{g},5a_{g},1b_{3u},2b_{3u},3b_{2u},4b_{2u},1b_{1g},1b_{1u},2b_{1u}$,
$1b_{3g}$ molecular orbitals. A total of eight electronic excited target
states are represented by 1844 configuration state functions (CSFs) for the
ground state. All the generated states up to 16 eV were retained in the final
close-coupling calculation. We used an R-matrix sphere of radius 10 bohrs and
continuum Gaussian-type orbitals with partial waves up to $l<4$. In the
following, this calculation will be referred to CAS1.

Initially, several basis sets - including DZP and cc-pVTZ - were tested to
investigate the stability of the target properties such as polarizability and
ground state energy and, finally, we chose the cc-pVTZ basis set in order to
perform the scattering calculations. Since the first electronically excited
state ${}^{3}\Sigma_{g}^{-}$ is approximately 13 eV above the dissociation
limit for the ground state, this latter state is essentially isolated and the
non-adiabatic effects are expected to be small. Therefore, for low electron
energy collisions, i.e. bellow 10 eV, only the ground electronic state is open
for ionization in e–BF${}_{2}^{+}$ collisions, and the dimension of the
geometry-fixed scattering matrix remains the same at low collision energies.

At low collision energies the fixed-nuclei scattering matrix depends only
weakly on energy. A sharper energy-dependence is observed at certain
relatively-high energies, corresponding to positions of Rydberg states
attached to the excited electronic states of the ion. A quantity convenient
for the analysis of the energy dependence of the scattering matrix is the
eigenphase sum. Figure 1 displays the eigenphase sum for the three different
geometries corresponding to a small displacement from equilibrium along each
normal mode of the BF${}_{2}^{+}$ ion. Bending and asymmetric stretching mode
calculations were performed in the $C_{2v}$ abelian subgroup, while for the
symmetric mode the group $D_{2h}$ was used in the calculations. The variation
of the eigenphase sums is smooth for energies below 6 eV. Above this value, a
sharp energy dependence at certain energies is observed due to the presence of
electronic Rydberg resonances attached to closed ionization limits.

##  II Cross sections and rate coefficients

Figure 2:  Vibrational excitation of BF${}_{2}^{+}$: probabilities
corresponding to the normal vibrational modes of the target ion.

For convenience, we introduce the quantities

$P_{i}=\frac{g_{i}}{2}\sum_{ll^{\prime}\lambda\lambda^{\prime}}\left|\frac{\partial
S_{l\lambda,l^{\prime}\lambda^{\prime}}}{\partial{q}_{i}}\right|^{2}\,,$ (4)

which could be interpreted as the probability for excitation of the
vibrational mode $i$, are listed in table 2. Figure 2 shows the weakly
dependence of that quantity on energy and, therefore, could be used as a
constant in calculations of the thermally-averaged rate coefficients.

Table 2: Parameters of Eqs. (8) and (9) calculated at $E=0.1$ eV collision energy. Mode $i$ | $P_{i}$
---|---
symmetric stretching | 0.053
bending | 0.20
asymmetric stretching | 0.05

Using $P_{i}$, the cross sections of Eqs. (1) and (2) are written as

$\displaystyle\sigma^{VE}_{i}(E_{el})=\frac{\pi\hbar^{2}}{2mE_{el}}P_{i}\
\theta(E_{el}-\hbar\omega_{i})\,,$ (5)
$\displaystyle\sigma^{DR}(E_{el})=\frac{\pi\hbar^{2}}{2mE_{el}}\sum_{i=1}^{3}P_{i}\
\theta(\hbar\omega_{i}-E_{el})\,.$ (6)

Figure 3 illustrates the DR cross section $\sigma^{DR}(E_{el})$ computed based
on Eq. (5). At very low scattering energies, i.e. below 0.02 eV, the DR cross
section is a smooth function inversely proportional to the incident energy of
the electron, as predicted by the Wigner threshold law, whereas at higher
energies, it exhibits a drop at each vibrational threshold.

Figure 3:  Cross section for the dissociative recombination of BF${}_{2}^{+}$.
The lowest vibrational threshold of each normal mode are indicated by arrows.

Due to the simple analytical form of the cross sections, the corresponding
rate coefficients are easily evaluated from the general expression

$\displaystyle\alpha(T)=\frac{8\pi}{\left(2\pi
k_{b}T\right)^{3/2}}\int_{0}^{\infty}\sigma(E_{el})\exp\left(-\frac{E_{el}}{k_{b}T}\right)E_{el}dE_{el}\,,$
(7)

giving

$\displaystyle\alpha^{VE}_{i}(T)=\sqrt{\frac{2\pi}{k_{b}T}}\frac{\hbar^{2}}{m^{3/2}}\
P_{i}\ \exp\left(-\frac{\hbar\omega_{i}}{k_{b}T}\right)\,,$ (8)
$\displaystyle\alpha^{DR}(T)=\sqrt{\frac{2\pi}{k_{b}T}}\frac{\hbar^{2}}{m^{3/2}}\
\sum_{i=1}^{3}P_{i}\
\left[1-\exp\left(-\frac{\hbar\omega_{i}}{k_{b}T}\right)\right]\,.$ (9)

where $k_{b}$ is the Boltzmann coefficient and $T$ is the temperature. The
thermally averaged rate coefficients for dissociative recombination and
vibrational excitation are shown in Figure 4.

In order to access the uncertainty of the present theoretical model, we have
performed a complete calculation of the DR and VE rate coefficients using
different basis sets and orbital space in the electron-scattering
calculations. The calculations were performed for three sets of parameters (1)
the CAS1 with the cc-pVTZ basis set, mentioned above; (2) a calculation
(referred here as cc-pVTZ CAS2) similar to (1) but with a smaller orbital
space, where 8 electrons are kept free in the active space; and (3) a
calculation (referred here as cc-pVQZ CAS1) similar to (1) but with the larger
basis cc-pVQZ. The results are in shown in Fig. 4. The difference between the
rate coefficients produced in the three calculations for the DR process and
the VE of the asymmetric stretching and bending modes is about $4\%$. The
uncertainty of the rate coefficient for the VE of the symmetric stretching
mode is larger, varying in the interval $10-40\%$ for different temperatures.
Notice that the overall probability for the symmetric stretching excitation is
much smaller than the probabilities for other modes and DR.

Figure 4:  Dissociative recombination and vibrational excitation of
BF${}_{2}^{+}$: rate coefficients - Eqs. (8) and (9). To give an idea about
the uncertainty of the present results, we also plotted results of two
calculations with the cc-pVTZ CAS2 (dashed line) and the cc-pVQZ CAS1 (dotted
line) sets of parameters of the model. For the vibrational excitation of the
symmetric stretching mode, the cc-pVTZ CAS2 and cc-pVQZ CAS1 curves are
indistinguishable and slightly below the cc-pVTZ CAS1 curve. For the remaining
three processes, the three calculations produce the curves almost
indistinguishable in the figure. Figure 5:  Relative importance of the
dissociative recombination and vibrational excitation of BF${}_{2}^{+}$, with
respect to those of BF+.

##  III Conclusions and discussions

In this study, cross sections and rate coefficients for dissociative
recombination and vibrational excitation of BF${}_{2}^{+}$ by electron-impact
were obtained using a theoretical approach that combines the normal modes
approximation for the vibrational states of the target ion, the vibrational
frame transformation, and the UK R-matrix code. The thermally averaged rate
coefficients have a simple analytical form.

The obtained thermally-averaged rate coefficients are relevant for the kinetic
modeling of molecule based cold non-equilibrium plasmas, in the context of
complete lack of other theoretical or experimental data on these processes for
this cation, and are ready to be used in the modeling of the fluorine/boron
plasma for etching or implantation processes. Indeed, we are presently able to
make important statements on the relative importance of BF${}_{2}^{+}$ with
respect to BF+ on the population and excitation balance. In particular, as
shown in Figure 5, the dissociative recombination of BF${}_{2}^{+}$ strongly
dominates that of BF+ below 7000 K, and the vibrational excitation displays
the same feature below 5000 K.

The rotational structure of the target ion and of the neutral molecule was
neglected in the present approach, which implies that the obtained cross
sections and rate coefficients should be viewed as averaged over initial
rotational states and summed over final rotational states of the corresponding
initial and final vibrational levels (for vibrational excitation) or
dissociative states (for dissociative recombination). Purely rotational
transitions, i.e. without changing the vibrational state, might be useful to
model very cold environment, below 40 K, which is not the case for the
presently investigated BF3 plasma. BF${}_{2}^{+}$ has no permanent dipole, so
the rotational transitions are likely to have very small cross sections.

## Acknowledgments

This work was supported by the National Science Foundation, Grant No
PHY-18-06915. Ioan F. Schneider thanks for generous financial support from La
Région Haute-Normandie, via the GRR Electronique, Energie et Matériaux and the
project BIOENGINE, from the Fédération de Recherche ”Energie, Propulsion,
Environnement”, and from the LabEx EMC3 and FEDER via the projects PicoLIBS
(ANR-10-LABEX-09-01), EMoPlaF and CO2-VIRIDIS. Ioan F. Schneider and János
Zsolt Mezei acknowledge support from the CNRS via the GdR THEMS, IAEA (Vienna)
via the Coordinated Research Project ”Light Element Atom, Molecule and Radical
Behaviour in the Divertor and Edge Plasma Regions”, the Programme National
”Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP
co-funded by CEA and CNES, and Fédération de Recherche ”Fusion par Confinement
Magnétique” (CNRS and CEA). János Zsolt Mezei and Khalid Hassouni acknowledge
support from USPC via ENUMPP and Labex SEAM. Ioan F. Schneider thanks
Laboratoire Aimé Cotton for outstanding hospitality.

## Data availability

Upon a reasonable request, the data supporting this article will be provided
by the corresponding author.

## References

  * Sennikov et al. (2017) P. G. Sennikov, R. A. Kornev, and A. I. Shishkin, Plasma Chem. Plasma Process. 37, 997 (2017).
  * Torigoe et al. (2016) M. Torigoe, K. Teii, and S. Matsumoto, IEEE Trans. Plasma Sci. 44, 3219 (2016).
  * Gonzatti et al. (2010) F. Gonzatti, F. Milesi, V. Delaye, J. Duchaine, F. Torregrosa, H. Etienne, and K. Yckache, AIP Conf. Proc. 1321, 27 (2010).
  * Duchaine et al. (2012) J. Duchaine, F. Milesi, R. Coquand, S. Barraud, S. Reboh, F. Gonzatti, F. Mazen, and F. Torregrosa, AIP Conf. Proc. 1496, 71 (2012).
  * Young et al. (2016) D. L. Young, W. Nemetti, V. La Salvia, M. R. Page, S. Theingi, J. Aguiar, B. G. Lee, and P. Stardins, Sol. Energy Mater. Sol. Cells 158, 68 (2016).
  * M. Farber and Srivastava (1984) M. M. Farber and R. D. Srivastava, J. Chem. Phys 81, 241 (1984).
  * Yong-Ki and Irikura (2000) K. Yong-Ki and K. K. Irikura, AIP Conf. Proc. 543, 220 (2000).
  * Stewart and Lieberman (1991) R. A. Stewart and M. A. Lieberman, J. Applied Physics 70, 3481 (1991).
  * Lieberman and Lichtenberg (2005) M. A. Lieberman and A. J. Lichtenberg, _Principles of plasma discharges and materials processing_ (John Wiley & Sons, Inc., Hoboken, New Jersey, 2005).
  * Mezei et al. (2016) J. Z. Mezei, F. Colboc, N. N Pop, S. Ilie, K. Chakrabarti, S. Niyonzima, M. Lepers, A. Bultel, O. Dulieu, O. Motapon, et al., Plasma Sources Sci. Technol. 25, 055022 (2016).
  * Douguet et al. (2011) N. Douguet, A. Orel, I. Mikhailov, I. F. Schneider, C. H. Greene, and V. Kokoouline, J. Phys.: Conf. Series 300, 012015 (2011).
  * Fonseca dos Santos et al. (2014) S. Fonseca dos Santos, N. Douguet, V. Kokoouline, and A. E. Orel, J. Chem. Phys. 140, 164308 (2014).
  * Werner et al. (2012) H.-J. Werner, P. J. Knowles, R. Lindh, F. R. Manby, M. Schütz, et al., WIREs Comput Mol Sci 2, 242 (2012).
  * Pyykko and Zhao (1990) P. Pyykko and Y. Zhao, J. Chem. Phys. 94, 7753 (1990).
  * Perić and Peyerimhoff (1993) M. Perić and S. Peyerimhoff, Mol. Phys. 78, 877 (1993).
  * Jacox and Thompson (1995) M. E. Jacox and W. E. Thompson, J. Chem. Phys. 102, 4747 (1995).
  * Tennyson (2010) J. Tennyson, Phys. Rep. 491, 29 (2010).
  * Carr et al. (2012) J. Carr, P. Galiatsatos, J. Gorfinkiel, A. Harvey, M. Lysaght, D. Madden, Z. Mašín, M. Plummer, J. Tennyson, and H. Varambhia, Euro. Phys. J. D 66, 58 (2012).
  * Tennyson et al. (2007) J. Tennyson, D. B. Brown, J. J. Munro, I. Rozum, H. N. Varambhia, and N. Vinci, J. Phys. Conf. Series 86, 012001 (2007).


11institutetext:  Laboratoire de Physique et Modélisation des Milieux
Condensés, Université de Grenoble-Alpes and CNRS, 25 Avenue des Martyrs, 38042
Grenoble, France 22institutetext: Institut Néel, Université de Grenoble-Alpes
and CNRS, 25 Avenue des Martyrs, 38042 Grenoble, France

# Inhomogeneous Josephson junction chains: a superconducting meta-material for
superinductance optimization

D. V. Nguyen 1122 D. M. Basko<EMAIL_ADDRESS>1 1

###### Abstract

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.

##  1 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 Jung2014  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 Manucharyan2009 ,
tunable microwave impedance matching Altimiras2013 , or a potential
implementation of the electrical current standard in quantum metrology based
on Bloch oscillations Likharev1985 ; Mooij2006 ; Guichard2010 .

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, $\sim\sqrt{\mu_{0}/\epsilon_{0}}=4\alpha R_{Q}$, where
$\alpha\approx{1}/137$ is the fine structure constant, and $R_{Q}\approx
13\>\mathrm{k}\Omega$ is the resistance quantum Feynman . 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 Tinkham , and thus is of non-
relativistic origin. The term “superinductance” is often used to denote such
superconductivity-based inductance.

Several structures, based on Josephson junctions (JJs), have been reported to
work as superinductors Masluk2012 ; Bell2012 . In the first one, a large
inductance was obtained by putting $N$ Josephson junctions in series, which
gave the total inductance $NL$ ($L$ is the inductance of a single junction).
In Ref. Bell2012 , 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. Masluk2012 .

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}=(\hbar/2e)^{2}(1/L)$. To work as an inductor, the junction must be in
the superconducting regime, $E_{J}\gg{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_{\mathrm{max}}$, or a lower limit on the junction area $\mathcal{A}$, as
both $E_{J},C\propto\mathcal{A}$.

Limitations on the junction number $N$ arise from the dependence of the chain
response on the frequency $\omega$. The phase slip rate, although
exponentially suppressed for $E_{J}\gg{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, $\omega\ll\omega_{1}$ (the lowest mode frequency).
Crucially, besides the capacitance $C$ of the junction between neighboring
superconducting islands, each island has a small capacitance $C^{\mathrm{g}}$
to the ground. This capacitance gives rise to screening of the Coulomb
interaction between the islands on a length scale
$\lambda=\sqrt{C/C^{\mathrm{g}}}$ and produces an acoustic-like region of the
mode dispersion
$\omega(q)=(LC)^{-1/2}\sqrt{\epsilon^{2}(q)/[\epsilon^{2}(q)+\lambda^{-2}]}$
of spatially homogeneous chains, where $\epsilon(q)=2\sin(q/2)$, and $q$ is
the wavenumber, $0\leqslant{q}\leqslant\pi$ (Fig. 1). The first mode
corresponds to $q=\pi/(N+1)$, so for large $N\gg\pi\lambda$, the frequency of
the lowest mode $\omega_{1}\propto{1}/N$, and the inductive response bandwidth
shrinks with increasing $N$. This was the main limitation for the device
studied in Ref. Masluk2012 , where a special effort was made to decrease the
parasitic ground capacitance $C^{\mathrm{g}}$.

Figure 1: Dispersion curve of a JJ chain with $C/C^{\mathrm{g}}=25$ (solid
curve) and modes of a chain with $N=100$ junctions (filled circles). The mode
frequency $\omega$ is measured in the units of the junction plasma frequency
$1/\sqrt{LC}$.

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 Taguchi2015 .
Effect of the normal mode structure on dephasing of the fluxonium qubit was
discussed in Ref. Viola2015 . 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. Masluk2012  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 $\mathcal{A}_{n}$ of each junction $n$
(assuming the island area to be already optimized to minimize the ground
capacitance as was done in Refs. Manucharyan2009 ; Masluk2012 ;
Manucharyan2012 ). This leads to a simultaneous variation of the junction
inductances $L_{n}$ and capacitances $C_{n}$, such that their product
$L_{n}C_{n}=\mathrm{const}$. Optimizing over all areas
$\\{\mathcal{A}_{n}\\}$, we find that the best result is still achieved for a
homogeneous configuration.

Figure 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 $L_{n}(0)$ is
determined by the total area of the junctions $\mathcal{A}_{n}$, shown by
hatching. When a magnetic field $B$ is applied, the inductance
$L_{n}(B)=L_{n}(0)/|\cos(\pi{B}S_{n}/\Phi_{0})|$ is determined by the magnetic
flux $BS_{n}$ through the SQUID loop area $S_{n}$, represented by the white
circular region in the center.

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(\pi{B}S/\Phi_{0})|$, where $\Phi_{0}=2\pi\hbar/(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(\pi{B}S_{n}/\Phi_{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 $\mathcal{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 $\mathcal{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.

##  2 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\gg{1}$. When the junctions are in the
superconducting regime, $E_{J}\gg{E}_{C}$, the oscillations of the
superconducting phase $\varphi_{n}$ on each island are small. Then, the
Josephson current through the $n$th junction from island $n$ to $(n+1)$ can be
written as $I_{n}=I_{n}^{c}\sin(\varphi_{n+1}-\varphi_{n})\approx
I_{n}^{c}(\varphi_{n+1}-\varphi_{n})$ Tinkham . Here,
$I_{n}^{c}=\hbar/(2eL_{n})=(2e/\hbar)E_{J,n}$ is the junction critical
current. Thus, the voltage drop across the junction can be determined by using
the Josephson relation:

$V_{n+1}-V_{n}=\frac{\hbar}{2e}\left(\frac{d\varphi_{n+1}}{dt}-\frac{d\varphi_{n}}{dt}\right)=\frac{\hbar}{2eI_{n}^{c}}\frac{dI_{n}}{dt}.$
(1)

This expression shows that the junction behaves as a linear inductor, and the
Josephson kinetic inductance is given by $L_{n}=\hbar/(2eI_{n}^{c})$. 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^{\mathrm{g}}_{n}$ is the capacitance of each
island to ground. We define the complex impedance $Z(\omega)$ of the chain at
frequency $\omega$ as the ratio of the voltage $V_{\omega}{e}^{-i\omega{t}}$
on an external ac voltage source, connected to the islands $n=1$ and $n=N+1$,
to the current $I_{\omega}{e}^{-i\omega{t}}$ through this source (Fig. 3).

Figure 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).

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}$:

$\displaystyle
Y_{1}({V_{1}}-{V_{2}})-i\omega{C}^{\mathrm{g}}_{1}V_{1}=I_{\omega},$ (2a)
$\displaystyle
Y_{n}(V_{n}-V_{n+1})+Y_{n-1}(V_{n}-V_{n-1})-i\omega{C}^{\mathrm{g}}_{n}V_{n}=0$
$\displaystyle\hskip 142.26378pt(n=2,\dots,N),$ (2b) $\displaystyle
Y_{N}(V_{N+1}-V_{N})-i\omega{C}^{\mathrm{g}}_{N+1}V_{N+1}=-I_{\omega},$ (2c)

where the junction admittance is defined as

$Y_{n}(\omega)=-i\omega{C_{n}}-\frac{1}{i\omega{L_{n}}}.$ (3)

System (2) can be written in the matrix form,
$\hat{\Upsilon}(\omega)\mathbf{V}=\mathbf{I}$, in terms of the column vectors
$\mathbf{V}^{T}=(V_{1},\ldots,V_{N+1})$ and
$\mathbf{I}^{T}=(I_{\omega},0,\ldots,0,-I_{\omega})$, as well as the
corresponding matrix $\hat{\Upsilon}(\omega)$. Then the chain impedance
$Z(\omega)$ can be expressed in terms of the matrix elements of the inverse
$\hat{\Upsilon}^{-1}(\omega)$ as
$Z=(\hat{\Upsilon}^{-1})_{11}-(\hat{\Upsilon}^{-1})_{1,N+1}-(\hat{\Upsilon}^{-1})_{N+1,1}+(\hat{\Upsilon}^{-1})_{N+1,N+1}$.
At low frequencies, the admittances are dominated by the inductive part, so
the impedance is given by $Z(\omega\to{0})=-i\omega{L}_{\mathrm{tot}}$, where
$L_{\mathrm{tot}}$ is the total inductance of the chain,

$L_{\mathrm{tot}}=\sum_{n=1}^{N}L_{n}.$ (4)

The approximation $Z(\omega)\approx-i\omega{L}_{\mathrm{tot}}$ is valid as
long as $\omega\ll\omega_{1}$, where $\omega_{1}$ is the lowest normal mode
frequency, for which $\det\hat{\Upsilon}(\omega)=0$.

As discussed in the introduction, ideally one would like to increase both
$L_{\mathrm{tot}}$ and $\omega_{1}$, but these two requirements are in
conflict. Thus, one can try to maximize $L_{\mathrm{tot}}$ at fixed
$\omega_{1}$, or maximize $\omega_{1}$ while keeping $L_{\mathrm{tot}}$ fixed.
We prefer the second option, as the constraint expressed by Eq. (4) is much
easier to resolve than the constraint $\omega_{1}=\mathrm{const}$. Thus, our
optimization problem is formulated as follows: find the spatial profile of
$L_{n},C_{n},C_{n}^{\mathrm{g}}$ which maximizes $\omega_{1}$ while keeping
$L_{\mathrm{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_{n}^{\mathrm{g}}$ 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 $\omega_{1}$.
This optimization was performed in Refs. Manucharyan2009 ; Masluk2012 ;
Manucharyan2012 ; ManucharyanThesis . Then, the ground capacitance of the
$n$th 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.
Masluk2012 , 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}\gg{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, $\propto{e}^{-(2/\pi)R_{Q}/Z_{n}}$
Matveev2002  is exponentially suppressed only for small junction impedances,
$Z_{n}\equiv\sqrt{L_{n}/C_{n}}\ll{R}_{Q}$ (we remind that $R_{Q}\approx
13\>\mathrm{k}\Omega$ 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 $\mathcal{A}_{\mathrm{max}}\gg{1}$ is an
independent dimensionless parameter of the problem. The junction inductance
and capacitance at the smallest area, $L_{\mathrm{max}}$ and
$C_{\mathrm{min}}$, can be chosen as the units of inductance and capacitance,
respectively. Thus, we have $N$ dimensionless variables $\mathcal{A}_{n}$,
allowed to vary in the range

$1\leqslant\mathcal{A}_{n}\leqslant\mathcal{A}_{\mathrm{max}}.$ (5a) They
determine the inductance and the capacitance of each junction as
$L_{n}=\frac{L_{\mathrm{max}}}{\mathcal{A}_{n}},\quad
C_{n}=C_{\mathrm{min}}\mathcal{A}_{n},$ (5b) and Eq. (4) thus imposes a
constraint on the set $\\{\mathcal{A}_{n}\\}$. In this case, the plasma
frequency of each junction is unchanged,
$1/\sqrt{{L_{n}}{C_{n}}}=1/\sqrt{{L_{\max}}{C_{\min}}}\equiv{\omega_{p}}$.
Finally, for the ground capacitances we use a simple form
$C^{\mathrm{g}}_{n}=C^{\mathrm{g}}_{\mathrm{min}}\,g(\mathcal{A}_{n-1}/2+\mathcal{A}_{n}/2),$
(5c)

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)=\sqrt{x}$, as mentioned in Ref. ManucharyanThesis . To define Eq. (5c)
at the ends, we set $\mathcal{A}_{0}\equiv\mathcal{A}_{1}$,
$\mathcal{A}_{N+1}\equiv\mathcal{A}_{N}$. Thus, the first optimization problem
is fully defined as maximization of $\omega_{1}$ determined from Eqs. (2),
whose coefficients are expressed by Eqs. (5b) and (5c) in terms of the
dimensionless areas $\mathcal{A}_{n}$. The optmization variables are the areas
$\mathcal{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

$N_{0}\equiv\frac{L_{\mathrm{tot}}}{L_{\mathrm{max}}}\leqslant{N}\leqslant{N}_{0}\mathcal{A}_{\mathrm{max}}.$
(6)

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 $\mathcal{A}_{n}$ (Fig. 2). By
applying a magnetic field $B$, one can change the SQUID inductance as

$L_{n}(B)=\frac{L_{n}(0)}{|\cos(\pi{B}S_{n}/\Phi_{0})|}$ (7)

where the zero-field inductance $L_{n}(0)$ is determined by the junction area
$\mathcal{A}_{n}$. This way of tuning the properties of the JJ by magnetic
field is routinely used in experiments (see, e. g., Ref. Weissl2015 ). 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
$\mathcal{A}_{1},\ldots,\mathcal{A}_{N}$, we consider another problem defined
via variables $\mathcal{F}_{1},\ldots,\mathcal{F}_{N}$:

$\displaystyle 1\leqslant\mathcal{F}_{n}\leqslant\mathcal{F}_{\mathrm{max}},$
(8a) $\displaystyle L_{n}=\frac{L_{\mathrm{max}}}{\mathcal{F}_{n}},\quad
C_{n}=C_{\mathrm{min}},\quad
C_{n}^{\mathrm{g}}=C^{\mathrm{g}}_{\mathrm{min}},$ (8b)

All junction areas are assumed to be the same, $\mathcal{A}_{n}=1$, so the
plasma frequency of each SQUID is modulated as
$1/\sqrt{L_{n}C_{n}}={\omega_{p}}\sqrt{\mathcal{F}_{n}}$, and each variable
$\mathcal{F}_{n}$ represents the ratio

$\mathcal{F}_{n}=\frac{|\cos(\pi{B}S_{n}/\Phi_{0})|}{\cos(\pi\Phi_{\mathrm{max}}/\Phi_{0})},\quad\frac{1}{\mathcal{F}_{\mathrm{max}}}\equiv\cos\frac{\pi\Phi_{\mathrm{max}}}{\Phi_{0}},$
(9)

where $\Phi_{\mathrm{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}\gg{E}_{C}$. Clearly, $\\{\mathcal{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
$\mathcal{F}_{n}$ is Eq. (4), and it restricts the chain length $N$ to the
interval

$N_{0}\equiv\frac{L_{\mathrm{tot}}}{L_{\mathrm{max}}}\leqslant{N}\leqslant\mathcal{F}_{\mathrm{max}}N_{0}.$
(10)

The two optimization problems, defined by Eqs. (5a)–(5c) and by Eqs.
(8a)–(8b), will be studied in the next two sections, respectively.

##  3 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 $\mathcal{A}_{n}=\mathcal{A}$, we have only
two variables, $\mathcal{A}$ and $N$. Constraint (4) fixes
$\mathcal{A}=N/N_{0}$, $L_{n}=L_{\mathrm{max}}N_{0}/N$,
$C_{n}=C_{\mathrm{min}}N/N_{0}$,
$C^{\mathrm{g}}_{n}=C^{\mathrm{g}}_{\mathrm{min}}\,g(N/N_{0})$. It is
convenient to denote the first mode frequency for this homogeneous chain by
$\Omega_{N}$. It is given by Masluk2012

$\displaystyle\Omega_{N}^{2}={}$
$\displaystyle{}\frac{1}{LC}\,\frac{1-\cos[\pi/(N+1)]}{1-\cos[\pi/(N+1)]+C^{\mathrm{g}}/(2C)}\approx$
$\displaystyle\approx{}$
$\displaystyle{}\frac{(L_{\mathrm{max}}C_{\mathrm{min}})^{-1}}{1+(C^{\mathrm{g}}_{\mathrm{min}}/C_{\mathrm{min}})(N_{0}/\pi)^{2}[x/g(x)]^{2}}$
(11)

for $N\gg{1}$. This is a decreasing function of $x\equiv{N}/N_{0}$ for any
$g(x)$ growing slower than linearly with $x$. Thus, $\omega_{1}$ is maximized
by taking $N=N_{0}$, all $L_{n}=L_{\mathrm{max}}$. We denote the corresponding
value of $\omega_{1}$ by $\Omega_{N_{0}}$.

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 $\omega_{1}$ from the inhomogeneiety would overcome the
loss due to the length increase. A qualitative idea of the best spatial
profile $\mathcal{A}_{n}$ can be obtained from the perturbation theory for
system (2), developed in Ref. Basko2013 . Let us use the homogeneous chain of
length $N$ with all $\mathcal{A}_{n}=\mathcal{A}=N/N_{0}$ and the first mode
frequency $\Omega_{N}$ as the zero approximation. If we now modify each
junction area by a small amount $\Delta\mathcal{A}_{n}$, the first-order
frequency shift is given by Basko2013

$\displaystyle\frac{\Delta\omega_{1}}{\Omega_{N}}=\frac{1}{N+1}\left(1-\Omega_{N}^{2}L_{\mathrm{max}}C_{\mathrm{min}}\right)\sum\limits_{n=1}^{N}\alpha_{n}\,\frac{\Delta\mathcal{A}_{n}}{\mathcal{A}},$
(12a)
$\displaystyle\alpha_{n}=\sin^{2}\frac{\pi{n}}{N+1}+\frac{2{\mathcal{A}}g^{\prime}(\mathcal{A})}{g(\mathcal{A})}\times{}$
$\displaystyle\qquad{}\times\left(\cos\frac{\pi}{N+1}\sin^{2}\frac{\pi{n}}{N+1}-\cos^{2}\frac{\pi/2}{N+1}\right).$
(12b)

The dependence of $\alpha_{n}$ on $n$ is quite simple
($\sin^{2}+\>\mathrm{const}$), and $\alpha_{n}$ is the largest for
$n=(N+1)/2$, in the middle of the chain. The value at the maximum
$\alpha_{(N+1)/2}>0$ as long as
$[2\mathcal{A}g^{\prime}(\mathcal{A})/g(\mathcal{A})]\sin^{2}[\pi/(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 $\omega_{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 $\mathcal{A}_{n}=1$, and two more junctions with $\mathcal{A}_{n}=2$, to
be put in the center. (Note that it is impossible to keep $N_{0}$ junctions
with $\mathcal{A}_{n}=1$, as the constraint would require the remaining one to
have $\mathcal{A}_{n}=\infty$). As the area change for the central junctions
is not small, the perturbative Eq. (12a) is not sufficient to describe this
situation. Still, $\omega_{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 $\Omega_{N_{0}}$.

The full optimization of all junction areas $\\{\mathcal{A}_{n}\\}$, subject
to constraint (4), can be performed numerically. For any $N>N_{0}$, we
maximize $\omega_{1}$ as a function of all the areas, calculated numerically
from the eigenvalue equation $\det\hat{\Upsilon}(\omega)=0$. The resulting
maximum $\omega_{1}$ is plotted versus $N$ in Fig. 4 for several values of
$C_{\mathrm{min}}/C^{\mathrm{g}}_{\mathrm{min}}$ and
$\mathcal{A}_{\mathrm{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 $\omega_{1}$ at $N>N_{0}$
is always below the best value for the homogeneous chain, $\Omega_{N_{0}}$. In
Fig. 5 we show the optimal spatial profile $\\{\mathcal{A}_{n}\\}$,
corresponding to one of the points in Fig. 4. Indeed, the best $\omega_{1}$
for a fixed $N$ is obtained by placing the largest junctions in the middle of
the chain. Still, the resulting gain in $\omega_{1}$ is smaller than the loss
due to the increase of the chain length from $N_{0}$ to $N$.

Figure 4: The first mode frequency $\omega_{1}$ (in units of the plasma
frequency
$\omega_{\mathrm{p}}\equiv{1}/\sqrt{L_{\mathrm{max}}C_{\mathrm{min}}}$)
obtained by full numerical optimization of all junction areas
$\\{\mathcal{A}_{n}\\}$, subject to constraint (4). We take $N_{0}=25$ for all
curves, while $\lambda^{2}\equiv
C_{\mathrm{min}}/C^{\mathrm{g}}_{\mathrm{min}}=400$ and 16 for panels (a) and
(b), respectively. Two values of $\mathcal{A}_{\mathrm{max}}=3$ and 10 were
chosen, shown by the blue and red symbols (lower and upper curves),
respectively, on each panel. The solid curve shows $\Omega_{N}$, the first
mode frequency for the homogeneous chain with $L_{n}=L_{\mathrm{max}}N_{0}/N$,
$C_{n}=C_{\mathrm{min}}N/N_{0}$,
$C^{\mathrm{g}}_{n}=C^{\mathrm{g}}_{\mathrm{min}}\sqrt{N/N_{0}}$, and the
dashed horizontal line shows the best homogeneous result $\Omega_{N_{0}}$.
Figure 5: The optimal spatial profile $\\{\mathcal{A}_{n}\\}$, giving the
largest $\omega_{1}$ for $N_{0}=25$, $N=50$, $\mathcal{A}_{\mathrm{max}}=10$,
$C_{\mathrm{min}}/C^{\mathrm{g}}_{\mathrm{min}}=400$.

##  4 SQUID loop area modulations

As in the previous section, we start by a straighforward study of the
homogeneous case. Constraint (4) fixes $\mathcal{F}=N/N_{0}$, so for the chain
with $L_{n}=L_{\mathrm{max}}N_{0}/N$, $C_{n}=C_{\mathrm{min}}$,
$C_{n}^{\mathrm{g}}=C_{\mathrm{min}}^{\mathrm{g}}$, assuming $N\gg{1}$,
instead of Eq. (3) we have

$\Omega_{N}^{2}\approx\frac{\pi\lambda}{L_{\mathrm{tot}}C_{\mathrm{min}}}\,\frac{N/(\pi\lambda)}{1+N^{2}/(\pi\lambda)^{2}},$
(13)

where $\lambda=\sqrt{C_{\mathrm{min}}^{\mathrm{g}}/C_{\mathrm{min}}}$ is the
screening length, which does not depend on $N$ and $\\{\mathcal{F}_{n}\\}$ for
problem (8b)–(8a). Expression (13) reaches maximum at $N=\pi\lambda$, so we
have to consider three cases for the position of this value with respect to
the interval (10).

(i) In the case $\pi\lambda>\mathcal{F}_{\mathrm{max}}N_{0}$, the frequency is
maximized by taking the longest possible chain,
$N=\mathcal{F}_{\mathrm{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_{\mathrm{tot}}$
which is too small; a larger inductance can be obtained by simply increasing
the length at almost no cost in $\omega_{1}$. So, this case has no practical
relevance.

(ii) When $N_{0}\leqslant\pi\lambda\leqslant\mathcal{F}_{\mathrm{max}}N_{0}$,
the frequency is maximized at $N=\pi\lambda$. 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 $\pi\lambda<N_{0}$, the frequency is maximized by taking the
shortest possible chain. This regime corresponds to demanding such a large
inductance $L_{\mathrm{tot}}$ that the first mode necessarily belongs to the
acoustic part of the dispersion curve. This is the regime where the
competition between $L_{\mathrm{tot}}$ and $\omega_{1}$ is the most severe; it
is in this regime that a gain in $\omega_{1}$ by introducing a spatial
variation of $\mathcal{F}_{n}$ would be the most interesting for practical
purposes.

The perturbation theory in small modulations $\Delta\mathcal{F}_{n}$ with
respect to a homogeneous chain with $N$ junctions gives a result, similar to
Eq. (12a):

$\frac{\Delta\omega_{1}}{\Omega_{N}}=\frac{1}{N+1}\sum\limits_{n=1}^{N}\frac{\Delta\mathcal{F}_{n}}{\mathcal{F}}\sin^{2}\frac{\pi{n}}{N+1},$
(14)

which again tells us that inductance modulations in the center of the chain
contribute the most to the increase in $\omega_{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 $\mathcal{F}=2$. The explicit
calculation of given in Appendix A: shows that this chain has
$\omega_{1}>\Omega_{N_{0}}$, and thus one can indeed improve over the
homogeneous result. However, for long chains, $N_{0}\gg\pi\lambda$, the gain
is quite small:

$\omega_{1}-\Omega_{N_{0}}\approx\frac{1}{2N_{0}\sqrt{LC}}\left(\frac{\pi\lambda}{N_{0}}\right)^{3}.$
(15)

Is it possible to gain more in $\omega_{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}\gg{1}$ where the
inductances are smaller by a factor $\mathcal{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):

$L_{n}=\left\\{\begin{array}[]{ll}L_{\mathrm{max}},&1\leqslant{n}\leqslant{N}_{1},\\\
L_{\mathrm{max}}/\mathcal{F},&N_{1}<n<N-N_{1},\\\
L_{\mathrm{max}},&N-N_{1}\leqslant{n}\leqslant{N}.\end{array}\right.$ (16)

Constraint (4) then fixes

$N_{1}=\frac{\mathcal{F}N_{0}-N}{2(\mathcal{F}-1)}.$ (17)

For $N-2N_{1}\gg{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}\gg\pi\lambda$,
then one can approximate the mode dispersion by the acoustic one,
$\omega(q)\approx{q}/\sqrt{LC^{\mathrm{g}}}$. Then, Eqs. (2) are transformed
into the Helmholtz equation with von Neumann boundary conditions at the ends
of the chain,

$\left(\frac{\partial}{\partial{x}}\frac{1}{L(x)}\frac{\partial}{\partial{x}}+\omega^{2}C^{\mathrm{g}}\right)V(x)=0,\quad\left.\frac{\partial{V}}{\partial{x}}\right|_{x=0,N}=0.$
(18)

For the piecewise function $L(x)$, given by Eq. (16), and for a given
frequency $\omega$, the wavenumbers in the outer regions and in the central
region are given by
$q=\omega\sqrt{L_{\mathrm{max}}C^{\mathrm{g}}_{\mathrm{min}}}$ and by
$q/\sqrt{\mathcal{F}}$, respectively. Thus, taking advantage of the symmetry
of $L(x)$ with respect to $x\to{N}-x$, we seek $V(x)$ in the form (the first
mode is odd)

$V(x)=\left\\{\begin{array}[]{ll}A\cos{q}x,&0<x<N_{1},\\\
A^{\prime}\sin[q(N/2-x)/\sqrt{\mathcal{F}}],&N_{1}<x<N-N_{1},\\\ -A\cos(qN-
qx),&N-N_{1}<x<N.\end{array}\right.$ (19)

The requirement of continuity of $V$ and $(1/L)(\partial{V}/\partial{x})$ at
$x=N_{1},N-N_{1}$ yields the following equation for $q$:

$\tan\left[\sqrt{\mathcal{F}}\,q\left(\frac{N_{0}}{2}-N_{1}\right)\right]=\sqrt{\mathcal{F}}\cot{q}N_{1}.$
(20)

For all $\mathcal{F}>1$, upon increasing $N_{1}$ from 0 to $N_{0}/2$ (that is,
upon decreasing $N$ from $N_{0}\mathcal{F}$ to $N_{0}$), the solution
monotonically rises from $q=\pi/(N_{0}\sqrt{\mathcal{F}})$ to $q=\pi/N_{0}$
(Fig. 6), the highest frequency being achieved in the shortest homogeneous
chain. This means that in the limit $N_{0}\gg\pi\lambda$ the gain in
$\omega_{1}$ is so small that it is not captured by the acoustic
approximation.

Figure 6: Solution of Eq. (20) as a function of $N_{1}$ for different values
of $\mathcal{F}=2,4,8,16$ (from the upper to the lower curve, respectively).

To check these considerations numerically, we perform the full optimization of
all $\\{\mathcal{F}_{n}\\}$, subject to constraint (4). As in the previous
section, for any $N>N_{0}$, we maximize $\omega_{1}$ as a function of all the
areas, calculated numerically from the eigenvalue equation
$\det\hat{\Upsilon}(\omega)=0$. The resulting maximum $\omega_{1}$ is plotted
versus $N$ in Fig. 7 for several values of $\lambda$ and
$\mathcal{F}_{\mathrm{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 $\mathcal{F}_{\mathrm{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).

Figure 7: The first mode frequency $\omega_{1}$, in the units of the plasma
frequency
$\omega_{\mathrm{p}}\equiv{1}/\sqrt{L_{\mathrm{max}}C_{\mathrm{min}}}$,
obtained by full numerical optimization of all $\\{\mathcal{F}_{n}\\}$,
subject to constraint (4), shown by symbols for $\mathcal{F}_{\mathrm{max}}=2$
and 10 (blue circles and red squares, respectively). The solid curve shows the
first mode frequency $\Omega_{N}$ for the homogeneous chain with
$L_{n}=L_{\mathrm{max}}N_{0}/N$, $C_{n}=C_{\mathrm{min}}$,
$C_{n}^{\mathrm{g}}=C_{\mathrm{min}}^{\mathrm{g}}$. We take $N_{0}=25$ for all
curves, while $\lambda=20$ and 4 for panels (a) and (b), respectively. The
dashed horizontal line shows the best homogeneous result.  Figure 8: The
optimal spatial profile of the inductance, $L_{n}/L_{\mathrm{max}}$, giving
the largest $\omega_{1}$ for $N=46$, $\mathcal{F}_{\mathrm{max}}=2$ (blue
circles) and for $N=111$, $\mathcal{F}_{\mathrm{max}}=10$ (red squares). Other
parameters are $N_{0}=25$, $\lambda=20$.

Figure 9: (a) The $N_{0}$ dependence of the first mode frequency for the
optimal inhomogeneous chain having the optimal length (symbols) and for the
shortest homogeneous chain with $N=N_{0}$ (solid curve). The frequencies are
measured in the units of the plasma frequency
$\omega_{\mathrm{p}}\equiv{1}/\sqrt{L_{\mathrm{max}}C_{\mathrm{min}}}$. (b)
The $N_{0}$ dependence of the chain length $N_{\mathrm{opt}}$, at which the
optimal value of $\omega_{1}$ is obtained. On both panels the blue circles and
the red squares correspond to $\mathcal{F}_{\mathrm{max}}=2$ and 10,
respectively, and we took $\lambda=20$.

##  5 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. Masluk2012 . 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/\sqrt{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.

###### Acknowledgements.

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

##  Appendix A: Chain with two central junctions modified

Let us start by deriving the dispersion relation for a homogeneous chain of
$N$ junctions with parameters $L_{1},\ldots,L_{N}=L$, $C_{1},\ldots,C_{N}=C$,
$C^{\mathrm{g}}_{1},\ldots,C^{\mathrm{g}}_{N+1}=C^{\mathrm{g}}$, and
$\sqrt{C/C^{\mathrm{g}}}\equiv\lambda$. Eq. (2) becomes

$\frac{1-\omega^{2}LC}{\omega^{2}LC^{\mathrm{g}}}\left(2{V_{n}}-{V_{n+1}}-{V_{n-1}}\right)-V_{n}=0.$
(A.1)

A plane wave, $V_{n}=A_{\pm}{e}^{\pm{i}qn}$, with any $A_{\pm}$ and $q$
satisfies this equation, provided that

$\frac{1-\omega^{2}LC}{\omega^{2}LC^{\mathrm{g}}}\,2(1-\cos{q})-1=0,$ (A.2)

which gives the usual dispersion relation Masluk2012 ,

$\omega(q)=\frac{1}{{\sqrt{LC}}}\sqrt{\frac{{1-\cos q}}{1-\cos
q+1/(2\lambda^{2})}}.$ (A.3)

For a given $\omega$, we seek the solution in the form
$A_{+}e^{iqn}+A_{-}e^{-iqn}$, and substitute it into Eqs. (2a), (2c) at the
ends of the chain, which play the role of the boundary conditions. These give,
respectively,

$\displaystyle A_{+}e^{iq}(1-e^{-iq})+A_{-}e^{-iq}(1-e^{iq})=0,$ (A.4a)
$\displaystyle A_{+}e^{iq(N+1)}(1-e^{iq})+A_{-}e^{-iq(N+1)}(1-e^{-iq})=0.$
(A.4b)

The first of these equations requires the solution to have the form

$V_{n}=A\cos[q(n-1/2)],$ (A.5a) while the second one imposes the form
$V_{n}=B\cos[q(n-N-3/2)],$ (A.5b)

with some $A$ and $B$. Matching these expressions in the bulk of the chain, we
obtain two possibilities, corresponding to even and odd modes with respect to
reflection $n\to{N}+2-n$:

$\displaystyle A=B,\quad q(n-1/2)=q(n-N-3/2)+2\pi{k},$ (A.6a) $\displaystyle
A=-B,\quad q(n-1/2)=q(n-N-3/2)+2\pi{k}+\pi,$ (A.6b)

where $k$ is an integer. Thus, the even modes have $q=2k\pi/(N+1)$, and the
odd ones $q=(2k+1)\pi/(N+1)$. Note that the first mode is odd.

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

Now, let us consider a chain with two central junctions modified. We assume
$N=N_{0}+1$ to be even, $N=2N_{1}$, then we can again take advantage of the
reflection symmetry, $n\to{N}+2-n$. In the homogeneous part of the chain we
assume $L_{n}=L$, $C_{n}=C$ and $C^{\mathrm{g}}_{n}=C^{\mathrm{g}}$, while in
the central region we set

$\displaystyle L_{N_{1}-}=L_{N_{1}+1}=\frac{L}{\xi},$ (A.7a) $\displaystyle
C_{N_{1}}=C_{N_{1}+1}=\zeta{C},$ (A.7b) $\displaystyle
C^{\mathrm{g}}_{N_{1}}=C^{\mathrm{g}}_{N_{1}+2}=\eta{C}^{\mathrm{g}},\quad
C^{\mathrm{g}}_{N_{1}+1}=\eta^{\prime}{C}^{\mathrm{g}}.$ (A.7c)

For the junction areas’ variation, considered in Sec. 3, we have to set
$\xi=\zeta=2$, $\eta=g(3/2)$, $\eta^{\prime}=g(2)$. For the loop areas’
variation (Sec. 4), we have $\xi=2$, $\zeta=\eta=\eta^{\prime}=1$. The
reflection symmetry is preserved, so the modes can still be classified as even
or odd, and by continuity we know that the first mode is odd. Thus, similarly
to Eqs. (A.5a), (A.5b), we can seek $V_{n}$ in the form

$V_{n}=\left\\{\begin{array}[]{ll}A\cos[q(n-1/2)],&n\leqslant{N}_{1},\\\
0,&n=N_{1}+1,\\\ -A\cos[q(n-N-3/2)],&n\geqslant{N}_{1}+2,\end{array}\right.$
(A.8)

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}=\pi/(N_{0}+1)=\pi/(2N_{1})$.

$V_{n}$ in the form (A.8) automatically satisfy the Kirchhoff laws for the
nodes $n=1,\ldots,N_{1}-1,N_{1}+1,N+3,\ldots,2N_{1}+1$. The Kirchhoff laws for
the remaining $n=N_{1},N_{1}+2$ are identical, so we have one independent
equation which determines $q$:

$\begin{split}&\frac{1-\omega^{2}LC}{\omega^{2}LC^{\mathrm{g}}}\left[\cos(qN_{1}-q/2)-\cos(qN_{1}-3q/2)\right]+{}\\\
&{}+\left(\frac{\xi-\zeta\omega^{2}LC}{\omega^{2}LC^{\mathrm{g}}}-\eta\right)\cos(qN_{1}-q/2)=0,\end{split}$
(A.9)

where the frequency $\omega$ is related to $q$ by Eq. (A.3). Note that
$\eta^{\prime}$ dropped out of the equation as $V_{N_{1}+1}=0$. Eq.
(LABEL:cosqN1=) can be identically rewritten as

$\displaystyle S(q)\cos(qN_{1}-q/2)-\cos(qN_{1}-3/2)=0,$ $\displaystyle
S(q)\equiv
1+\xi+\left[\lambda^{2}(\xi-\zeta)-\eta\right]4\sin^{2}\frac{q}{2},$

or, equivalently, as

$\displaystyle\cot{q}N_{1}{}{}={}$
$\displaystyle{}\frac{2-\xi-[1-\eta+\lambda^{2}(\xi-\zeta)]4\sin^{2}(q/2)}{\xi+[1-\eta+\lambda^{2}(\xi-\zeta)]4\sin^{2}(q/2)}\tan\frac{q}{2}.$
(A.10)

The left-hand side of this equation passes through zero precisely at
$q=\pi/(2N_{1})=q_{0}\ll{1}$, with a large negative slope. Thus, to find out
whether the solution $q=q_{*}$ is larger or smaller than $q_{0}$, we just need
to check the sign of the right-hand side at $q=q_{0}$.

When only junction areas $\mathcal{A}_{N_{1}},\mathcal{A}_{N_{1}+1}$, are
varied, that is, $\xi=\zeta=2$ and $\eta=g(3/2)>1$, the large factor
$\lambda^{2}$ drops out, so the right-hand side of Eq. (A.10) is necessarily
positive, and thus $q_{*}<q_{0}$. For the variation of SQUID loop areas only,
we have $\xi=2$, $\zeta=\eta=1$, which leads to $q_{*}>q_{0}$. Note, however,
that the difference $q_{*}-q_{0}$ is quite small:

$q_{*}-q_{0}\approx\frac{\pi}{2N_{0}^{2}}\left(\frac{\pi\lambda}{N_{0}}\right)^{2}.$
(A.11)

## References

  * (1) P. Jung, A. V. Ustinov, S. M. Anlage, Supercond. Sci. Technol. 27, 073001 (2014)
  * (2) V. E. Manucharyan, J. Koch, L. I. Glazman, M. H. Devoret, Science 326, 113 (2009)
  * (3) C. Altimiras, O. Parlavecchio, P. Joyez, D. Vion, P. Roche, D. Esteve, F. Portier, Appl. Phys. Lett. 103, 212601 (2013)
  * (4) K. K. Likharev, A. B. Zorin, J. Low Temp. Phys. 59, 347 (1985)
  * (5) J. E. Mooij, Y. V. Nazarov, Nat. Phys. 2, 169 (2006)
  * (6) W. Guichard, F. W. J. Hekking, Phys. Rev. B 81, 064508 (2010)
  * (7) R. P. Feynman, R. B. Leighton, M. Sands, The Feynman Lectures on Physics, Vol. 2 (Addison-Wesley, Reading, MA, 1964)
  * (8) M. Tinkham, Introduction to Superconductivity (McGraw-Hill Book Co., New York, 2004)
  * (9) N. A. Masluk, I. M. Pop, A. Kamal, Z. K. Minev, M. H. Devoret, Phys. Rev. Lett. 109, 137002 (2012)
  * (10) M. T. Bell, I. A. Sadovskyy, L. B. Ioffe, A. Y. Kitaev, and M. E. Gershenson, Phys. Rev. Lett. 109, 137003 (2012)
  * (11) M. Taguchi, D. M. Basko, F. W. J. Hekking, Phys. Rev. B 92, 024507 (2015)
  * (12) G. Viola, G. Catelani, Phys. Rev. B 92, 224511 (2015)
  * (13) V. E. Manucharyan, N. A. Masluk, A. Kamal, J. Koch, L. I. Glazman, M. H. Devoret, Phys. Rev. B 85, 024521 (2012).
  * (14) V. E. Manucharyan, Ph. D. Thesis (Yale University, 2012).
  * (15) K. A. Matveev, A. I. Larkin, L. I. Glazman, Phys. Rev. Lett. 89, 096802 (2002).
  * (16) T. Weißl, G. Rastelli, I. Matei, I. M. Pop, O. Buisson, F. W. J. Hekking, W. Guichard, Phys. Rev. B 91, 014507 (2015)
  * (17) D. M. Basko, F. W. J. Hekking, Phys. Rev. B 88, 094507 (2013)