# Positional Disorder (Random Gaussian Phase Shifts) in the Fully Frustrated Josephson Junction Array (2D XY Model)

\[

## Abstract

We consider the effect of positional disorder on a Josephson junction array with an applied magnetic field of $`f=1/2`$ flux quantum per unit cell. This is equivalent to the problem of random Gaussian phase shifts in the fully frustrated 2D XY model. Using simple analytical arguments and numerical simulations, we present evidence that the ground state vortex lattice of the pure model becomes disordered, in the thermodynamic limit, by any finite amount of positional disorder.

\]

The stability of vortex lattices to random disorder is a topic of considerable recent interest, motivated by studies of the high temperature superconductors. In two dimensions (2D), periodic arrays of Josephson junctions form a well controlled system for investigating similar issues of vortex fluctuations and disorder. Here we consider the effect of “positional” disorder on the vortex lattice of the fully frustrated Josephson array, with $`f=1/2`$ flux quantum of applied magnetic field per unit cell.

Positional disorder  was first discussed with respect to the Kosterlitz-Thouless (KT) transition for the $`f=0`$ model in zero magnetic field. Early arguments  predicting a reentrant normal phase at low temperatures have been revised by recent works  which argue that there is a finite critical disorder strength $`\sigma _c\sqrt{\pi /8}`$; for $`\sigma <\sigma _c`$ an ordered state persists for $`0TT_c(\sigma )`$. For the pure $`f=1/2`$ case on a square grid , the ordered state has two broken symmetries: the $`U(1)`$ symmetry (“KT-like” order) associated with superconducting phase coherence, and the $`Z(2)`$ symmetry (“Ising-like” order) associated with the “checkerboard” vortex lattice, in which a vortex sits on every other site. Previous works  have considered the effect of positional disorder on this $`f=1/2`$ model; all have concluded that both Ising-like and KT-like order persist for at least small disorder strengths $`\sigma `$. In this work, however, we present new arguments that suggest that, for $`f=1/2`$, the critical disorder is $`\sigma _c=0`$.

The Hamiltonian for the Josephson array is given by the “frustrated” 2D XY model ,

$$[\theta _i]=\underset{i\mu }{}U(\theta _i\theta _{i+\widehat{\mu }}A_{i\mu }),$$
(1)   
where $`i`$ are the sites of a periodic square grid with basis vectors $`\widehat{\mu }=\widehat{x}`$, $`\widehat{y}`$, the sum is over all nearest neighbor (n.n.) bonds $`i,i+\widehat{\mu }`$ and $`\theta _i\theta _{i+\widehat{\mu }}A_{i\mu }`$ is the gauge invariant phase difference across the bond, with $`A_{i\mu }=(\varphi _0/2\pi )_i^{i+\widehat{\mu }}𝐀𝑑\mathrm{}`$ the integral of the vector potential.

Positional disorder arises from random geometric distortions of the bonds of the grid, resulting in, $`A_{i\mu }=A_{i\mu }^{(0)}+\delta A_{i\mu }`$; $`A_{i\mu }^{(0)}`$ is the value in the absence of disorder, and $`\delta A_{i\mu }`$ is the random deviation. We take the $`\delta A_{i\mu }`$ to be independent Gaussian random variables with

$$[\delta A_{i\mu }]=0,\mathrm{and}[\delta A_{i\mu }\delta A_{j\nu }]=\sigma ^2\delta _{ij}\delta \mu \nu .$$
(2)   
$`[\mathrm{}]`$ denotes an average over the quenched disorder. The positionally disordered array is thus also referred to as the XY model with random Gaussian phase shifts.

When $`U(\varphi )`$ is the Villain function , the Hamiltonian (1) is equivalent to a dual “Coulomb gas” of interacting vortices ,

$$[n_i]=\frac{1}{2}\underset{ij}{}(n_if\delta f_i)G_{ij}(n_jf\delta f_j).$$
(3)   
The sum is over all pairs of dual sites $`i`$, $`j`$, $`n_i`$ is the integer vorticity on site $`i`$, and the interaction $`G_{ij}`$ is the Green’s function for the 2D discrete Laplacian operator, $`\mathrm{\Delta }_{ik}G_{kj}=2\pi \delta _{ij}`$, where $`\mathrm{\Delta }_{ij}\delta _{i,j+\widehat{x}}+\delta _{i,j\widehat{x}}+\delta _{i,j+\widehat{y}}+\delta _{i,j\widehat{y}}4\delta _{ij}`$. For large separations, $`G_{ij}\mathrm{ln}|𝐫_i𝐫_j|`$. The $`f_if+\delta f_i`$ are $`(1/2\pi )`$ times the circulation of the $`A_{i\mu }`$ around dual site $`i`$; $`f`$ is the average applied flux, while $`\delta f_i`$ is the deviation due to the random $`\delta A_{i\mu }`$,

$$\delta f_i=\frac{1}{2\pi }\left[\delta A_{i,x}+\delta A_{i+\widehat{x},y}\delta A_{i+\widehat{y},x}\delta A_{i,y}\right].$$
(4)   
Geometrically distorting a bond increases the flux through the cell on one side of the bond, while reducing the flux through the cell on the opposite side by the same amount. The $`\delta f_i`$ are thus anticorrelated among n.n. sites. Positional disorder is thus the same as random dipole pairs of quenched charges $`\pm \delta f_i`$ . From Eqs. (2) and (4) we get,

$$[\delta f_i]=0,\mathrm{and}[\delta f_i\delta f_j]=\frac{\sigma ^2}{4\pi ^2}\mathrm{\Delta }_{ij}.$$
(5)   
The Hamiltonian (3) can be rewritten as interacting charges in a one body random potential ,

$$[q_i]=\frac{1}{2}\underset{ij}{}q_iG_{ij}q_j\underset{i}{}q_iV_i,$$
(6)   
where $`q_in_if`$, and the random potential is $`V_i=_jG_{ij}\delta f_j`$. For $`f=1/2`$, $`q_i=\pm 1/2`$. From Eq. (5),

$`[V_i]=0,`$ $`\mathrm{and}[V_iV_j]={\displaystyle \underset{k,l}{}}G_{ik}[\delta f_k\delta f_l]G_{lj}`$ (7)   
$`=`$ $`{\displaystyle \frac{\sigma ^2}{4\pi ^2}}{\displaystyle \underset{k,l}{}}G_{ik}\mathrm{\Delta }_{kl}G_{lj}={\displaystyle \frac{\sigma ^2}{2\pi }}G_{ij},`$ (8)   
The $`V_i`$ thus have logarithmic long range correlations.

We now use an Imry-Ma  type argument to estimate the stability of the doubly degenerate checkerboard ground state to the formation of a square domain of side $`L`$. The energy of such an excitation consists of a domain wall term, $`E_d`$, which is present for the pure case, and a pinning term, $`E_p`$, due to the interaction with the random $`V_i`$. $`E_d(L)`$ has the form ,

$$E_daL+c\mathrm{ln}L+d.$$
(9)   
The first term is the interfacial tension of the domain wall; the second term comes from net charge that builds up at the corners of the domain . Calculating $`E_d(L)`$ numerically for a pure system, we find an excellent fit to Eq. (9), with $`a=0.28`$, $`c=0.15`$, and $`d=0.058`$.

By Eq. (8), the average pinning energy of the domain $`𝒟`$, $`[E_p]=2_{i𝒟}q_i[V_i]=0`$, but the variance is,

$$[E_p^2]=4\underset{i,j𝒟}{}q_i[V_iV_j]q_j=\frac{4\sigma ^2}{2\pi }\underset{i,j𝒟}{}q_iG_{ij}q_j=\frac{4\sigma ^2}{\pi }E_0,$$
(10)   
where $`E_0=(\pi /32)L^2`$ is the ground state energy of the checkerboard domain . The root mean square pinning energy is thus,

$$[E_p]_{rms}=bL,b=\frac{\sigma }{2\sqrt{2}}0.35\sigma .$$
(11)   
For domains whose energy is lowered by the interaction with $`V_i`$, the typical excitation energy is $`E=E_d[E_p]_{rms}`$. Eqs. (9) and (11) imply that when $`b>a`$, i.e. when $`\sigma >\sigma _c0.8`$, $`E(L)`$ has a maximum at $`L=\xi (c/2\sqrt{2})/(\sigma \sigma _c)`$. Domains of size $`L>\xi `$ will lower their energy by increasing in size, and so disorder the system. Thus, one naively expects that when $`\sigma <\sigma _c`$ the system preserves its Ising-like order, but when $`\sigma >\sigma _c`$ the system is disordered into domains of typical size $`\xi `$.

However the leading size dependencies of Eqs. (9) and (11), $`E_d[E_p]_{rms}L`$, are exactly the same as found in the 2D n.n. random field Ising model (RFIM). For the RFIM it is known  that 2D is the lower critical dimension, that the randomness causes domains walls at $`T=0`$ always to roughen and so acquire an effective negative line tension, and that the critical disorder is $`\sigma _c=0`$, i.e. any amount of disorder, no matter how weak, destroys the Ising-like order of the pure case. By analogy, we suggest that the positionally disordered $`f=1/2`$ 2D XY model similarly has $`\sigma _c=0`$. Our conclusion, that $`[E_p]_{rms}L`$ as in the 2D RFIM, follows from a subtle cancellation between the long range interactions between charges $`q_i`$, and the long range correlations of the random potential $`V_i`$.

To check this prediction, we carry out Monte Carlo (MC) simulations of the Hamiltonian (3) with periodic boundary conditions on $`L\times L`$ square grids. Our MC procedure is as follows . One MC excitation attempt consists of the insertion of a neutral $`n=\pm 1`$ vortex pair on n.n. or next n.n. sites, which is accepted or rejected using the usual Metropolis algorithm. $`L^2`$ such attempts we call one MC pass. At each temperature we typically used $`4000`$ MC passes to equilibrate the system, followed by $`128,000`$ MC passes to compute averages. Every $`100`$ passes we attempt a global excitation reversing the sign of all the charges, $`q_iq_i`$. For each disorder realization we cooled down two distinct “replicas”, starting with different random charge configurations and using different random number sequences. In only about $`3\%`$ of the cases did the two replicas fail to give reasonable agreement.

To test for Ising-like order we define an order parameter analogously to an Ising antiferromagnet,

$$M=\frac{1}{L^2}\underset{i}{}q_i(1)^{x_i+y_i}.$$
(12)   
We first consider $`\sigma =0.3`$, smaller than both the naive estimate of $`\sigma _c=0.8`$, and the $`\sigma _c=\sqrt{\pi /8}0.63`$ of the $`f=0`$ model. Fig. 1 plots $`[M^2]`$ vs. $`T`$, averaged over $`200`$ disorder realizations, for sizes $`L=10`$, $`14`$ and $`20`$. All curves start to increase from zero near $`T0.13`$, which is $`T_c(\sigma =0)`$ of the pure model. However $`[M^2]`$ at low $`T`$ decreases steadily with increasing $`L`$. The reason for this becomes clearer if we consider the histogram of values of $`M^2`$ that occur as we sample the different realizations of disorder. We show such histograms in Figs. 1b-d, for the lowest temperature $`T=0.02`$. As $`L`$ increases, the statistical weight shifts from predominantly ordered systems ($`M^2=1/4`$), to predominantly disordered systems ($`M^2=0`$). Assuming that this trend continues, we expect that as $`L\mathrm{}`$, $`[M^2]0`$.

To measure the “random field correlation length” $`\xi `$, we consider the vortex correlation function,

$$S(𝐤)=\frac{1}{L^2}\underset{i,j}{}e^{i𝐤(𝐫_i𝐫_j)}n_in_j.$$
(13)   
For the pure case, $`S(𝐤)`$ in the ordered phase has singular Bragg peaks at $`𝐊=\pm \pi \widehat{x}\pm \pi \widehat{y}`$. If the vortex lattice is disordered, these peaks will broaden, and their finite width provides a measure of $`\xi `$. Writing $`𝐤=𝐊+\delta 𝐤`$, and assuming a Lorentzian shape for the disorder averaged peak, $`[S(𝐤)]1/(\delta k^2+\xi ^2)`$, we determine $`\xi `$ by fitting to this form for $`\delta k=0`$, and $`\delta k=2\pi /L`$ . In Fig. 2 we show $`\xi `$ vs. $`\sigma `$ at our lowest $`T=0.02`$, for several system sizes $`L`$. Only for our smallest value $`\sigma =0.25`$ does a finite size effect remain. In this case, however, $`\xi `$ decreases as $`L`$ increases. This is in contrast to the increase of $`\xi `$ with $`L`$ that one would expect if one were approaching a second order transition. This behavior is consistent with that seen in Figs. 1b-d, where as $`L`$ increases, a greater fraction of the disorder realizations result in disordered states.

We next fit our results for $`\xi (\sigma )`$ to several possible scaling expressions: (i) $`\xi e^{C/(\sigma \sigma _c)^2}`$, (ii) $`\xi e^{C/(\sigma \sigma _c)}`$, and (iii) $`\xi |\sigma \sigma _c|^p`$. The first has been suggested by Binder  for the 2D RFIM. While in Binder’s expression $`\sigma _c=0`$, here we leave it as an arbitrary parameter to be determined from the fit. The second has been suggested for the positionally disordered $`f=0`$ model , in which $`\sigma _c>0`$. The third is the familiar power law form. Using data for only the largest $`L`$ for each $`\sigma `$, the results of these fits are shown in Fig. 2. The value of $`\sigma _c`$ and the $`\chi ^2`$ of the fit for each case is (i) $`\sigma _c=0.0046\pm 0.050`$, $`\chi ^2=67`$;(ii) $`\sigma _c=0.0134\pm 0.055`$, $`\chi ^2=67`$; (iii) $`\sigma _c=0.0013\pm 0.098`$, $`p=2.86\pm 0.84`$, $`\chi ^2=7.6`$. The power law (iii) gives a significantly better fit than (i) or (ii), however all give $`\sigma _c=0`$ within the estimated error. Given the rather limited range of the data, the above fits should be treated with caution. However they do indicate that the data contains no suggestion of a diverging $`\xi `$ at a finite $`\sigma `$. Coupled with our Imry-Ma argument, we thus find a consistent picture suggesting that $`\sigma _c=0`$ for the $`f=1/2`$ 2D XY model.

Returning to the case $`\sigma =0.3`$, where Ising-like order has been lost, we now consider whether the system may still have a finite temperature “spin glass” transition to a disordered but frozen vortex state. To test for this we measure the self and cross overlaps , $`Q_{\mathrm{self}}`$ and $`Q_{\mathrm{cross}}`$,

$`Q_{\mathrm{self}\alpha }`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}{\displaystyle \underset{i}{}}n_i^{(\alpha )}(t)n_i^{(\alpha )}(t+\tau )`$ (14)   
$`Q_{\mathrm{cross}}`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}{\displaystyle \underset{i}{}}n_i^{(\alpha )}(t)n_i^{(\beta )}(t).`$ (15)   
$`\alpha `$ and $`\beta `$ index the two independent replicas. For $`\tau `$ sufficiently large we expect $`Q_{\mathrm{self}\mathrm{\hspace{0.17em}1}}=Q_{\mathrm{self}\mathrm{\hspace{0.17em}2}}=Q_{\mathrm{cross}}`$, if the system is well equilibrated. Averaging Eq. (15) over several values of $`\tau 2000`$ to improve our statistics, we plot $`[Q_{\mathrm{self}\mathrm{\hspace{0.17em}1}}]`$, $`[Q_{\mathrm{self}\mathrm{\hspace{0.17em}2}}]`$ and $`[Q_{\mathrm{cross}}]`$ vs. $`T`$ in Fig. 3a. We see that our system is fairly well equilibrated down to the lowest $`T`$ we study. To test for a spin glass transition, we measure the overlap susceptibility,

$$\chi _Q=L^2\left\{[Q_{\mathrm{cross}}^2][Q_{\mathrm{cross}}^2]\right\},$$
(16)   
which we plot vs. $`T`$ in Fig. 3b for various system sizes. The peak in $`\chi _Q`$ near $`T0.06`$ shows no noticeable increase as $`L`$ increases, thus suggesting that there is no finite temperature spin glass transition.

If the vortices are not frozen, but are free to diffuse, one expects that superconducting phase coherence is also destroyed. To explicitly test this we measure the helicity modulus. The Hamiltonian (3) can viewed as representing the XY model with “fluctuating twist” boundary conditions . Using the method of Ref. , we determine the dependence of the total free energy $`F`$ of the corresponding XY model, as a function of the twist $`(\mathrm{\Delta }_x,\mathrm{\Delta }_y)`$ which is applied in a “fixed twist” boundary condition. We then determine the $`(\mathrm{\Delta }_{x0},\mathrm{\Delta }_{y0})`$ that minimizes $`F`$; the helicity modulus tensor is then the curvature of $`F`$ at the minimizing twist, $`\mathrm{{\rm Y}}_{\mu \nu }=^2F/\mathrm{\Delta }_\mu \mathrm{\Delta }_\nu `$. In Fig. 4a we plot $`\mathrm{{\rm Y}}_1`$, the largest of the two eigenvalues of $`\mathrm{{\rm Y}}_{\mu \nu }`$, vs. $`T`$, for $`\sigma =0.3`$ and sizes $`L=10,14,20`$. At all $`T`$, $`\mathrm{{\rm Y}}_1`$ continues to decrease as $`L`$ increases, giving no suggestion of a finite temperature transition. In Figs. 4b-d we plot histograms of the minimizing twist $`\mathrm{\Delta }_0`$ for the three sizes $`L`$. Note, in choosing our random phase shifts $`\delta A_{i\mu }`$, we impose the constraint $`_i\delta A_{i\mu }=0`$ in order to remove one trivial source of $`\mathrm{\Delta }_00`$. We see that the width of the distributions of $`\mathrm{\Delta }_0`$ steadily increases with increasing $`L`$, suggesting  that the strength of the random disorder is renormalizing to greater values on larger length scales.

To conclude, our results suggest that Ising-like order is destroyed for any finite amount of positional disorder. Further, we found in one specific case that when the Ising-like order vanished, no spin glass order or phase coherence existed either. We speculate that this remains true as well for any finite disorder strength. Although $`\sigma _c=0`$, the finite $`\xi (\sigma )`$ nevertheless can become extremely large for small values of $`\sigma `$. When $`\xi `$ exceeds the size of the experimental or numerical sample, the system will indeed look ordered. We believe this explains previous numerical work on this problem which reported the persistence of Ising-like order at small $`\sigma `$. In the most recent of these works, Cataudella  reports at $`\sigma 0.113`$ a finite $`T_c`$ to an Ising-like ordered state. The correlation length exponent that he finds is $`\nu 1.7`$, clearly different from that of the pure model. Using our scaling form (iii) we can estimate that at this value of $`\sigma `$, $`\xi 120`$, much larger than Cataudella’s largest system size of $`L=36`$. His results may thus be reflecting a cross over region at $`L<\xi `$, rather than a true transition.

We thank Prof. Y. Shapir for many valuable discussions. This work has been supported by DOE grant DE-FG02-89ER14017.