Magnetic monopoles and superinsulation in Josephson junction arrays

Electric-magnetic duality or S-duality, extending the symmetry of Maxwell's equations by including the symmetry between Noether electric charges and topological magnetic monopoles, is one of the most fundamental concepts of modern physics. In two-dimensional systems harboring Cooper pairs, S-duality manifests in the emergence of superinsulation, a state dual to superconductivity, which exhibits an infinite resistance at finite temperatures. The mechanism behind this infinite resistance is the linear charge confinement by a magnetic monopole plasma. This plasma constricts electric field lines connecting the charge-anti-charge pairs into electric strings, in analogy to quarks within hadrons. Yet the origin of the monopole plasma remains an open question. Here we consider a two-dimensional Josephson junction array (JJA) and reveal that the magnetic monopole plasma arises as quantum instantons, thus establishing the underlying mechanism of superinsulation as two-dimensional quantum tunneling events. We calculate the string tension and the dimension of an electric pion determining the minimal size of a system capable of hosting superinsulation. Our findings pave the way for study of fundamental S-duality in desktop experiments on JJA and superconducting films.

The superinsulating state, dual to superconductivity 1-7 is a remarkable manifestation of S-duality 8 in condensed matter physics. Superinsulators exhibit infinite resistance at finite temperatures, mirroring the infinite conductance of superconductors. The mechanism preventing charge transport is the linear charge confinement 7 of both Cooper pairs and normal excitations by a magnetic monopole plasma. This plasma constricts electric field lines connecting the charge-anticharge pairs into electric strings, in analogy to quarks within hadrons 9 .
Maxwell equations in vacuum are symmetric under the duality transformations interchanging the electric and magnetic fields E → B and B → −E (we use hereafter natural units c=1, =1, ε 0 =1). This duality holds in the presence of field sources, provided magnetic monopoles 8 are included along with electric charges. While despite intensive searches 10 , no elementary particles with a net magnetic charge have ever been observed, monopoles emerge and are detected as topological excitations in strongly correlated systems, see, for example 11,12 . Notably, these monopoles emerge as classical particles that freeze out upon cooling down the system. A drastically different class of phenomena arises if monopoles form a monopole plasma as a result of multiple instanton quantum tunneling events. In this case, a monopole plasma offers an ideal screening mechanism for electric fields, and the system harboring the monopole plasma makes a perfect dielectric with zero static dielectric constant, ε = 0, as long as the electric field does not exceed some threshold value 13 . Now, as the perfect diamagnetism is associated with an infinite conductance, i.e. super-conductivity, the perfect dielectricity should correspond to dual superconductors possessing an infinite resistance, i.e. superinsulators 1,5 .
Dual superconductivity was introduced in the 1970s by 't Hooft as a Gedankenexperiment for quark confinement 9 . The idea was exactly that perfect dielectricity, in analogy to the Meissner effect in superconductors, would squeeze electric fields into thin flux tubes with quarks at their ends. When quarks are pulled apart, it is energetically more favourable to pull out of the vacuum additional quark-antiquark pairs and to form several short strings instead of a long one. As a consequence, the colour charge can never be observed at distances above the fundamental length scale, 1/Λ QCD and quarks are confined. Only colour-neutral hadron jets can be observed in collider events.
Superinsulation as an emergent condensed matter state was first proposed in 1 on the basis of electric-magnetic duality and independently reinvented in 5 on the basis of duality between two different symmetry realizations of the uncertainty principle. Experiments reporting superinsulation detected it in films experiencing superconductor-insulator transition 3,4 . Both considerations 1,5 involved the symmetric interchange of charges and vortices in 2D systems, Finally, the topological gauge theory of superinsulation put forth in 7 revealed that the relevant fundamental duality is the one relating charges and magnetic monopoles rather than vortices. Accordingly, superinsulation is the result of the proliferation of the monopole plasma and represents the Abelian realization of dual superconductivity 14 in condensed matter. The experimental implications, including the Berezinskii-Kosterlitz-Thouless (BKT) criticality of the deconfinement transition and the electric field-induced breakdown of confinement, were observed in NbTiN films 13,15 . Yet, the monopoles of 7 has emerged in the framework of a long-distance effective field theory of thin films. Here we complete the description of superinsulation and consider a Josephson junction array (JJA), which, in particular, represents a "microscopic" model for a superconducting film 16 , and develop an exact magnetic monopole theory of superinsulation in JJA.

I. RESULTS
We start with the notion that, contrary to charges, vortices are topological excitations, characterized by a topological quantum number. The configuration space of the theory of vortices decomposes into so-called superselection sectors, characterized by the integer total vortex number, which are connected via instantons, nonperturbative configurations representing quantum tunneling events between topological vacua 17 . As a consequence, charges are conserved but vortices are not and can "appear" and"disappear" via quantum tunneling events forming the instantons. In two spatial di-mensions (2D), these instantons are nothing but magnetic monopoles 18 . The instantons are known to make a noticeable impact on the low-temperature physics of one-dimensional (1D) system. In particular, the global O(2) model, representing the physics of 1D superconducting quantum wires with screened Coulomb interactions, admits instantons representing quantum phase slips 18 . These quantum phase slips cause a superconductor-tometal quantum transition 19,20 at zero temperature, an insulating phase possibly emerging in finite systems coupled to the environment 21 . Remarkably, in 2D 7,13 , the monopole instantons manifest a much more profound and striking action, governing not only metallic but superinsulating behavior.
We consider a square Josephson junction array (JJA) with the spacing comprising superconducting islands with the nearest-neighbor Josephson coupling of the strength E J . Each island has a self-capacitance C 0 and mutual capacitances C to its nearest neighbors. The corresponding charging energies are E C 0 = e 2 /2C 0 and E C = e 2 /2C. The degrees of freedom of the array are the integer multiples of the fundamental charge unit 2e of the Cooper pair on each island, q x ∈ Z, and the quantummechanically conjugated phases ϕ x ∈ [0, 2π]. The partition function for such a JJA 16 is given by (see Methods) where S is the Euclidean action and the sum runs over the 3D Euclidean lattice with spacing 0 in the "time" direction, which, as we will show below, represents the (inverse) tunneling frequency. Here, ∆ i and∆ i are forward and backward finite differences, ∆ ≡∆ i ∆ i is the corresponding 2D finite difference Laplacian and ∆ 0 and ∆ 0 are forward and backward finite time differences (see Methods). The integer charges q x interact via the twodimensional Yukawa potential with the mass C 0 /C/ . In the experimentally accessible nearest-neighbors capacitance limit C C 0 , this implies a two-dimensional Coulomb law at distances smaller than the electrostatic screening length Λ = (C/C 0 ). Then the charging energy E C and the Josephson coupling E J are the two relevant energy scales which can be further traded for one energy parameter ω P = √ 8E C E J , the Josephson plasma frequency, and one numerical parameter g = π 2 E J /2E C , the dimensionless conductance. In the following we will consider the physics of JJA at energies much below the plasma frequency, which takes the role of the natural ultraviolet (UV) cutoff in the theory.
In the limit C 0 = 0, which we will henceforth consider, the partition function of the JJA can be mapped exactly 1 onto the partition function of a topological Chern-Simons gauge theory 23 (see Methods), A magnetic monopole instanton depicted by the red ball is assigned to the center of the 3D cube in the 3D Euclidean lattice with the spatial spacing , comprising a single JJA plaquette and the elemental unit 0 along the quantized Euclidean time. An elemental fluxon carrying the phase 2π of a single vortex located in the JJA plaquette splits into four parts each carried away through the vertical plaquettes, and the original JJA plaquette one unit time later does not contain the vortex anymore. The monopole tunneling event interpolates between two states differing by one unit of the topological quantum number. The quantity 0 2 can be considered as the volume of the monopole. The monopole itself is anisotropic, with no flux coming out in the third direction, because of the deep non-relativistic limit of the effective compact QED action in JJA.
where K µν is the lattice Chern-Simons operator 1 (see Methods). Here, a µ and b µ are fictitious gauge fields representing conserved charge and vortex fluctuations by their dual field strengths, j µ = K µν b ν and φ µ = K µν a ν , respectively. The first term in the action is the topological mixed Chern-Simons term 23 between these two types of dual fluctuations. The integers Q i are the electric topological excitations of the system, the integers M i are the magnetic topological excitations. Together with the vortex number M 0 , the latter form a 3-current M µ which is conserved due to the gauge invariance in the b µ gauge sector,∆ µ M µ = 0. Due to this constraint, only the two integers M i are the independent degrees of freedom. From the point of view of the original Minkowski space-time, the 3-current M µ describes events in which one vortex disappears from the array, the flux being "car-ried away" by the spatial vortex currents M i . From the Euclidean space-time point of view, however, M µ are the components of a 3D magnetic field. A configuration such as the one in Fig. 1 represents thus a unit magnetic monopole, the JJA vortex on the lower plaquette playing the role of the Dirac string 8 . The integer monopole charge is m = ∆ i M i . The asymmetry of a monopole, whose flux flows out only in the spatial directions, but not over a whole 3D lattice cube, is due to the deep nonrelativistic limit of the JJA gauge theory. One sees in Fig. 1 how the flux of the JJA vortex is divided up into four parts and is carried away by the M i in the spatial directions. As a consequence, on the upper plaquette of the cube, representing the same JJA plaquette one quantum of time later, there is no more any vortex. Thus, the magnetic monopole m expresses the tunneling of the system between two different topological vacua.
We now discuss the implications of the monopole plasma proliferation. This occurs in the phase where the electric topological excitations Q i are suppressed because of their large energy, one sets Q i = 0. To establish the nature of the monopole plasma phase, we derive the electromagnetic response of the system by coupling the charge current j µ = K µν b ν to the real physical electromagnetic potential A µ where S is the Chern-Simons gauge theory action intro-duced in Eq. (2). Integrating out the matter fields, and taking the limit 0 ω P 1, one arrives at the effective action for electromagnetic fields in the monopole plasma phase as where F i are the spatial components of the dual electromagnetic vector strength F µ =K µν A ν and where the magnetic topological excitations, encoding the monopoles, are now additional dynamical variables which have to be summed over in the partition function. This is a deep non-relativistic version of Polyakov's compact QED action 14,18 in which only electric fields survive. Its form shows that the action is periodic under shifts F i → F i + 2πM i , with integer M i and that the gauge fields are thus indeed compact, i.e. angular variables defined on the interval [−π, +π].
One is used to the fact that electromagnetic fields mediate Coulomb forces between static charges, a 1/|x| potential in 3D, or a log|x| potential in 2D. Monopoles in compact electromagnetism drastically change this, as we now show. We consider two external probe charges of strength ±q ext and compute how their interaction potential is changed by the monopoles. To do so we consider the expectation value of the Wilson loop operator W (C), where C is a closed loop in 3D Euclidean space-time (a factor is absorbed into the gauge field A µ to make it dimensionless), When the loop C is restricted to the plane formed by the Euclidean time and one of the space coordinates, W (C) measures the potential between two external probe charges ±q ext . A perimeter law indicates a shortrange potential, while an area-law is tantamount to a linear interaction between the probe charges 17,18 . For couplings g/ 0 ω P large enough, the action peaks near F i = M i , allowing for the saddle-point approximation to compute the Wilson loop. Using the lattice Stoke's theorem, one rewrites Eq. (5) as For q ext = 1, i.e. Cooper pair probes, the result is (see Methods) where A is the area of the surface S enclosed by the loop C. This area law indicates a linear potential between test Cooper pairs, with the string tension where z is the instanton fugacity, G(0) is the value of the 2D lattice Coulomb potential at coinciding points and we have reinstated physical units. The string binds together charges, prevents charge transport on arrays of a sufficient size and is the origin of the infinite resistance characterizing superinsulation. For single electron probes, q ext = 1/2 in our units, the string tension is Single electrons are thus also confined, hence the absence of charge transport mediated by thermally excited normal quasiparticles in the superinsulating state, that has remained a tantalizing puzzle ever since the experimental discovery of the superinsulation 5 . The string tension contains two factors. The first, ω P / depends solely on the "classical" array parameters, the lattice spacing and plasma frequency. The second factor depends on the quantum characteristics of the array, the dimensionless conductance and the ratio between the long tunneling time 0 and the short phase oscillation period 1/ω P . To estimate the typical string size string = c /σ we will take the following typical values for experimental JJA, = 100 nm and ω P = 10 GHz. Then, the first contribution to the string size amounts to 0 string / = c/ ω P ≈ 550. This number is reduced by the second factor in (9). However, even for 0 ω P = O(1000), we still have string / ≈ 150, at the border of the total size of typical, presently manufactured JJA. In other words, In field theory parlance these JJA are too far from the infrared confining fixed point 22 and due to asymptotic freedom (for a review see 17 ) only the screened Coulomb forces within an electric "meson" can be observed. In order to detect superinsulation on JJA one must thus design an array with sufficiently high plasma frequency, sufficiently small interdistance between centres of superconducting islands and with the linear size sufficiently large to fit an entire string, presumably in the thousands of lattice spacings. A promising platform for highly controllable JJA, capable to host superinsulation, is offered by proximity arrays that can be driven through the SIT by either a gate voltage 24 or a magnetic field 25 . Note that in the system adopted in 25 , the dual twin to a Cooper pair Mott insulator, the vortex Mott insulator, has already been observed at the superconducting side of the SIT.
Finally, a comment about the role of disorder is due. In Ref. 26 the origin of superinsulation was related solely to disorder localizing charges. This mix up emerged since, in 26 , superinsulation was confused with manybody-localization, which was introduced in the seminal paper 27 . Our results conclusively show that the origin of superinsulation lies in the proliferation of quantum tunneling events (magnetic monopole instantons), which can be viewed as the 2D generalization of 1D quantum slips in wires, with no role of disorder involved.

II. DISCUSSION
We demonstrated that magnetic monopoles appearing in JJA are a deep non-relativistic version of the monopoles introduced by Polyakov in the framework of compact QED 14,18 and derived how these instantons dominate the JJA dynamics at low energies, far below the JJA plasma frequency. The tension of the string binding charges into neutral "mesons" is expressed through JJA parameters, the distance between superconducting islands, , the plasma frequency, ω P , and the dimensionless conductance g. We found that both Cooper pairs and normal excitations, are confined by monopoles, thereby resolving the enigma of the absence of current due to single-charge excitations in superinsulators. One of the experimental implications of our results is that the typical JJA used so far are far too coarse and small to accommodate an entire electric string. In field theory parlance they are too far from the infrared confining fixed point 22 and due to asymptotic freedom 17 only the screened Coulomb forces within an electric "meson" can be observed. The large size of the electric mesons reflects the fact that the electromagnetic interaction is much weaker than the strong force.This explains the paradox-ical enigma why superinsulators where experimentally seen in films but not yet in the paradigmatic JJA system for which they were first derived 5,7 and indicates the direction for further experimental research. Devising large-size JJA and proximity arrays will open an opportunity of observing superinsulation in highly controllable and tuneable systems and of exploring the fundamental properties of S-duality via the desktop experiments.
The Hamiltonian for a planar JJA of spacing , with nearest neighbours Josephson couplings E J , ground capacitances C 0 and nearest neighbours capacitances C is given by 16,33 H = where boldface characters denote the sites of the twodimensional array, < xy > indicates nearest neighbours, V x is the electric potential of the island at x and ϕ x the phase of its order parameter (we shall use natural units c = 1, = 1, ε 0 = 1). The Hamiltonian (11) can be rewritten as where ∆ i and∆ i are forward and backward finite differences and ∆ ≡∆ i ∆ i is the two-dimensional finite difference Laplacian. The phases ϕ x are quantummechanically conjugated to the charges E x on the islands: these are quantized in integer multiples of 2e (Cooper pairs), E x = 2eq x , q x ∈ Z, where e is the electron charge. The Hamiltonian (12) can be expressed in terms of charges and phases by noting that the electric potentials V x are determined by the charges E x via a discrete version of Poisson's equation: Using this in (12) we get where E C ≡ e 2 /2C. The integer charges q x interact via a two-dimensional Yukawa potential of mass C 0 /C/ . In the nearest-neighbours capacitance limit C C 0 , which is accessible experimentally, this becomes essentially a two-dimensional Coulomb law.
The partition function of the JJA admits a phase-space path-integral representation 33 where β = 1/T is the inverse temperature. In (15) time has to be considered also as discrete, as generally appropriate when degrees of freedom can change only in integer steps. We introduce thus a discrete time step 0 whose inverse represents the ultraviolet (UV) energy cutoff in the model. The interval 0 represents the minimal time interval on which the dynamics is still governed by the horizontal Hamiltonian (14). For frequencies above 1/ 0 new modes can be excited. We thus substitute the time integrals and space sums over a lattice with nodes x by a sum over space-time lattice nodes x, with x 0 = t denoting the discrete time direction. Denoting by ∆ 0 the (forward) finite time differences we obtain Eq. (1) of the main text.

B. Lattice Chern-Simons operator
Formulating a lattice version of the Chern-Simons operator µαν ∂ α requires some care, if gauge invariance has to be properly implemented 1 . We introduce first the forward and backward finite difference and shift operators on the 3D Euclidean lattice with sites denoted by {x} and directions indicated by Greek letters and lattice spacing, , whereμ denotes a unit vector in direction µ and d = in the spatial directions, d = 0 in the Euclidean time direction. Summation by parts on the lattice interchanges both the two finite differences (with a minus sign) and the two shift operators. Gauge transformations are defined by using the forward finite differences. In terms of these operators one can then define two lattice Chern-Simons terms where no summation is implied over the equal indices µ and ν. Summation by parts on the lattice interchanges also these two operators (without any minus sign). Gauge invariance is then guaranteed by the relations Note that the product of the two Chern-Simons terms gives the lattice Maxwell operator where ∆ =∆ µ ∆ µ is the 3D Laplace operator.

C. Deriving the Chern-Simons gauge theory
We first use the Villain representation (for a review see 34 ) to express the cosine interaction in (15) in terms of a set of integer link variables a i . By introducing real charge currents j i and assuming that the size of the system is much smaller than Λ, so that we can safely set C 0 → 0 from now on, we arrive at where we have dropped the underscripts referring to the lattice positions of the variables and where the summation over equal Greek 3D lattice direction indices is implied. We have also introduced the notation j 0 for the integer charges.
To proceed further, we stress that Eq. (20), derived from what is viewed as the standard JJA Hamiltonian, misses a crucial piece. This missing contribution is a kinetic term proportional to the vortex mass. The omission of this term is common practice when considering overdamped junctions. However, this omission is a priori not justified in arrays, in which collective effects may lead to a renormalization of the vortex mass, no matter how large its bare value may be. It is known since the very early days of JJA that integrating over charge fluctuations leads anyway to a vortex kinetic term 33 . It is a general principle in field theory that whatever is induced by fluctuations must be included at bare level. It is also known that dissipation is substantially reduced when the Coulomb interaction becomes long-range 35 , exactly the regime we are interested in. Finally, ballistic motion of vortices has indeed been observed experimentally 36 . This calls for the necessity of adding a vortex kinetic term to the action. Such a vortex kinetic term would involve the time derivative of a i , the integer variable conjugated to the charge currents and would represent phase slips corresponding to one vortex moving from one plaquette to a neighbouring one. When the coefficient of this kinetic term (∂ 0 a i ) 2 takes the value π 2 /4 0 E C we can introduce a real Lagrange multiplier a 0 and a fictitious electric field φ i = K iµ a µ and write the vortex kinetic term and the charge Coulomb interaction compactly as In this representation, the Coulomb interaction between the charges follows from the Gauss law constraint associated with the Lagrange multiplier a 0 . This procedure works only at a particular value of the vortex mass that makes the model self-dual under the interchange of charge and vortex variables while simultaneously interchanging π 2 E J ↔ 2E C , or alternatively, g ↔ 1/g. As a consequence it was called the self-dual approximation in 1 . We shall adopt this approximation here too. At this point we note that the charge current j µ is conserved and, hence, it can be represented as the field strength associated to a second fictitious gauge field b µ as where b 0 is a real variable, while b i are integers. We then use Poisson's formula, turning a sum over integers {n µ } into an integral over real variables, to make all components of the gauge fields a µ and b µ real, at the price of introducing integer link variables Q i and M i , Finally we note that the quantities (1/2π)K µν ∆ ν ϕ are the circulations of the array phases around the plaquettes orthogonal to the direction µ in 3D Euclidean space-time and are thus quantized as 2π integers. We can thus absorb the quantities K i0 ∆ 0 ϕ +K ij ∆ j ϕ in a redefinition of the integers M i and defineK 0i ∆ i ϕ = 2πM 0 . The original integral over the phases ϕ can then be traded for a sum over the vortex numbers M 0 , which is the gauge theory partition function in the main text.

D. Summing over the monopole instanton plasma
We start from the instanton plasma representation of the Wilson loop expectation value, where ∆ 2 is the 2D Laplacian. Following Polyakov we introduce a scalar field χ and we rewrite this as where the angle η = 2π∆ i S i /(−∆ 2 ) represents a dipole sheet on the Wilson surface S and the monopole fugacity z is determined by the self-interaction as with G(0) being the inverse of the 2D Laplacian at coinciding arguments. In (27) we also adopted the dilute instanton approximation, valid at low g, in which one takes into account only single monopoles m = ±1. The sum can now be explicitly performed, with the result, x ∆iχ∆iχ+ 8πg 0 ω P z(1−cos(χ+qextη)) , (28) By shifting the field χ by −q ext and introducing µ 2 = 4πgz/ 0 ω P we can rewrite this as W (C) = 1 Z χˆD χ e − 0 ω P 2πg x 1 2 ∆i(χ−qextη)∆i(χ−qextη)+µ 2 (1−cos(χ)) .
For small g this integral is dominated by the classical solution to the equation of motion where ∆ 2 is the 2D spatial Laplacian. Far from the boundaries of the Wilson surface S, where S i = 0, this reduces to a one-dimensional equation. Taking the Wilson surface S to lie in the (t, x 2 ) plane, this gives, Following 18 we solve this equation in the continuum limit, For q ext = 1 (corresponding to Cooper pairs in our case), the classical solution with the boundary conditions χ cl → 0 for |x 1 | → ∞ is Inserting this back in (29) we get formula (8) in the main text. The solution with the same boundary conditions for q ext = 1/2 (corresponding to single electrons in our case), instead is The action of this solution leads to formula (10) in the main text.