text:
# Domain Wall Renormalization Group Study of XY Model with Quenched Random Phase Shifts

## I Introduction

The $`XY`$ model with quenched random phase shifts as a model for a superconducting glass has been intensively investigated over the last decade, focusing on the so-called gauge glass model which corresponds to the case with maximal disorder. Since a transport current exerts a force on a flux lattice, it tends to move in response which causes dissipation of the current. The existence of disorder, which destroys the flux lattice structure, is essential to pin the vortices in order for a superconducting phase to exist in a high $`T_c`$ superconductor vg-fisher ; vg-blatter ; vg-natter . Although there exists no proof whether or not the gauge glass and vortex glass are in the same universality class, it is of interest as the simplest model of a disordered superconductor and is still not understood despite all the effort expended on it.
From numericalgg-hs ; gg-rtyf ; gg-g ; gg-fty ; gg-ks ; gg-mg  and experimental experiment1  studies, it is believed that the gauge glass has no ordered phase at any finite temperature in two dimensions. In three dimensions, numerical domain wall renormalization group (DWRG) studiesDWRG1 ; DWRG2  indicate that the lower critical dimension seems to be close to three. However the situation is less conclusive, since the simulations are limited to small system sizes. Finite temperature Monte Carlo studies yield a transition temperature $`T_c/JO(1)`$gg-hs ; gg-rtyf ; gg-wy  which is difficult to reconcile with DWRG studiesgg-rtyf ; gg-g ; gg-ks ; gg-mg  as these studies imply that the lower critical dimension for superconducting glass order is close to three. Experimentally there is also some evidence for a finite temperature phase transition to a superconducting glass phase experiment2 ; experiment3 .
The Hamiltonian of the XY model with random quenched disorder can be written as

$$H=\underset{<ij>}{}V(\theta _i\theta _jA_{ij})$$
(1)
where $`V(\varphi )`$ is an even, $`2\pi `$ periodic function of $`\varphi `$ with a maximum at $`\varphi =\pi `$ and minimum at $`\varphi =0`$, usually taken as $`V(\varphi _{ij})=J_{ij}\mathrm{cos}(\varphi _{ij})`$. The sum is over all nearest neighbor pairs of sites and the coupling constants, $`J_{ij}`$, are assumed uniform, $`J_{ij}=J>0`$. The random bond variables $`A_{ij}`$, which are responsible for the randomness and frustration, are taken to be independent and uniformly distributed in $`(\alpha \pi ,\alpha \pi ]`$ with $`0\alpha 1`$. For a gauge glass, $`\theta _i`$ is the phase of the superconducting order parameter at site $`i`$ of a square lattice in $`2D`$ and a simple cubic lattice in $`3D`$. The random bond variables $`A_{ij}`$ are taken to correspond to maximal disorder with $`\alpha =1`$. An external field applied to an extreme type II superconductor induces a uniform component $`A_{ij}^0=(2\pi /\mathrm{\Phi }_0)_i^j𝐀𝑑𝐥`$ where $`𝐀`$ is the vector potential of the applied field and $`\mathrm{\Phi }_0=hc/2e`$ is the quantum of flux. In this work, we take $`A_{ij}^0=0`$, corresponding to zero applied field. Unless explicitly stated, we consider an unscreened system with $`\alpha =1`$ corresponding to maximal disorder. The Hamiltonian of Eq. (1) also describes the $`XY`$ magnet with random Dzyaloshinski-Moriya interactions xy-magnet  and also a Josephson junction array with positional disorder. gk1 ; gk2  These studies xy-magnet ; gk1 ; gk2  showed that the existence of weak disorder $`(\alpha 1)`$ does not destroy an ordered phase at intermediate temperature but predict a re-entrant transition to a disordered phase at low temperature in two dimensions. However, recent analytic xyrdm-nsk ; xyrdm-kn ; xyrdm-s ; xyrdm-cf  and numerical gg-ks  studies suggest the absence of a re-entrant transition and that there exists an ordered phase for $`T<T_c(\alpha )`$ when $`\alpha <\alpha _c`$.
When the random bond variables $`A_{ij}`$ are restricted to $`0`$ or $`\pi `$ with equal probabilities, this model reduces to the $`\pm J`$ $`XY`$ spin glass, which is believed to be in a different universality class due to the additional reflection symmetry villain-vsg  which is absent in the case of uniformly or Gaussian distributed $`A_{ij}`$. An $`XY`$ spin glass may have both spin and chiral glass order associated with rotational and reflection symmetry, respectively. It has been suggested that, in $`2D`$ and $`3D`$, spin and chiral variables decouple at long distances and order independently sg-mc-kt ; sg-dw-kt ; sg-dw-k ; sg-dw-mg , and the lower critical dimensions are $`d_l4`$ for spin glass order and $`d_l<3`$ for chiral glass order. However, the decoupling scenario contradicts the analytic studies on a ladder lattice nhm , on a tube lattice tnh , and on a $`2D`$ lattice with a special choice of disorder nh . Recent numerical simulations sg-dw-mg  also suggests $`d_l`$ for a spin glass order may be close to three.
In this paper, we re-investigate the possibility of an ordered phase at small but finite temperature $`T`$ by a numerical domain wall renormalization group (DWRG)DWRG1 ; DWRG2 , or defect energy scaling. The domain wall or defect energy of the system is computed by using the Hamiltonian in the Coulomb gas (vortex) representation, which is more convenient for numerical work as it eliminates spin wave contributions to the energy. Although the conventional DWRG method can handle only the scaling of the coupling constant $`J(L)`$ at scale $`L`$, which is proportional to the domain wall or defect energy, our method enables us to study the flows of both the coupling constant and the disorder strength, $`A(L)`$, at length scale $`L`$ gg-ks . We apply this to the case of general disorder strength, $`0\alpha 1`$. The outline of the paper is as follows. In Section II we discuss the DWRG method and also our strategy. In Section III, we explicitly perform the transformation of the $`3D`$ Hamiltonian of Eq. (1) from the phase to the Coulomb gas representation. Our numerical method is explained in Section IV. Finally we discuss our numerical results in Section V and give a brief discussion of some of the effects of weak disorder, $`\alpha <1`$, and of finite screening of vortex - vortex interactions.

## II Strategy

The general idea behind a DWRG is to compute, analytically or numerically, the energy $`\mathrm{\Delta }E(L)`$ of a domain wall in a system of linear size $`L`$ and fit this to a finite size scaling form

$$\mathrm{\Delta }E(L)L^\theta $$
(2)
where $`\theta `$ is a stiffness exponent, whose sign is of fundamental importance. If $`\theta <0`$, $`\mathrm{\Delta }E(L)`$ vanishes in the thermodynamic limit. The energy of the domain wall or defect excitation vanishes which implies that, for $`T>0`$, the probability of the defect $`P_Le^{\mathrm{\Delta }E(L)/kT}1`$ as $`L\mathrm{}`$. This in turn implies that the density of such defects is finite when $`T>0`$ and there will be no resistance to an infinitesimal applied force and the system has no order. This is analogous to the vanishing of the shear modulus in a liquid, the superfluid density in a superfluid or superconductor and the spin stiffness constant in an isotropic magnet when $`T>T_c`$. On the other hand, if $`\theta >0`$, such defects will have zero probability when $`L=\mathrm{}`$ and the system will have finite stiffness and will be ordered at sufficiently small $`T>0`$.
In a uniform system without disorder, the definition of the energy of a domain wall of size $`L`$, $`\mathrm{\Delta }E(L)`$, is intuitively obvious. The first step is to find the ground state (GS) energy of a system of size $`L`$, which requires applying boundary conditions (BC) which are compatible with the GS configuration. For a ferromagnet, this is straightforward to implement as the GS configuration is known to be one with all spins parallel and periodic BC are compatible with this. To impose a spin domain wall perpendicular to the $`\widehat{𝐱}`$ direction, one simply changes the BC to antiperiodic along $`\widehat{𝐱}`$ and periodic in the other $`d1`$ directions. Then it immediately follows that

$$\mathrm{\Delta }E(L)=E_{ap}(L)E_p(L)L^{dn}$$
(3)
where $`n=1`$ for an Ising model and $`n=2`$ for a system with a continuous symmetry such as $`XY`$ and Heisenberg models.
One would like to use the same strategy for random systems described by Eq. (1), as suggested by Anderson for Ising spin glassespwa . However, it is not so clear how to proceed because, for a particular sample (realization of disorder), neither the GS configuration nor compatible BC is known so computing the defect energy $`\mathrm{\Delta }E(L)`$ is problematical. Assuming $`\mathrm{\Delta }E(L)`$ can be calculated, the stiffness exponent is defined by the scaling ansatz

$$\mathrm{\Delta }E(L)L^\theta $$
(4)
where $`\mathrm{}`$ denotes an average over realizations of disorder. To our knowledge, it is not known how to calculate analytically either the GS energy $`E_0(L)`$ or the energy $`E_D(L)`$ of the system containing a defect relative to this GS which means that one must proceed numerically. A number of conceptual and technical difficulties are apparent. The first, and most important, is the technical problem of computing the energy difference $`\mathrm{\Delta }E(L)`$ between the energies of the system subject to two different BC. We ultimately want the disorder averaged defect energy $`\mathrm{\Delta }E(L)=E_a(L)E_b(L)`$ where $`E_a(L)`$ is the lowest energy of a particular sample subject to BC denoted by $`a`$ and $`E_b(L)`$ with BC denoted by $`b`$. We need the individual energies $`E_a(L)`$ and $`E_b(L)`$ essentially exactly because the uncertainty in $`\mathrm{\Delta }E(L)`$ must be kept as small as possible. Also, to our knowledge, there is no proof that the scaling ansatz of Eq. (4) is a correct description and, even if it is, the only thing we can be sure of is $`\theta (d2)/2`$. All results are based on fitting data to the scaling form of Eq. (4) so one is attempting both to verify the scaling ansatz and to estimate a numerical value of $`\theta `$. For any conclusion to be believable, the data must have both very small errors and fit Eq. (4) extremely well. The first requirement of highly accurate data is the most important as the estimate of $`\theta `$ depends on this. Assuming that $`E_a(L)`$ and $`E_b(L)`$ can be determined exactly for each sample, then $`\mathrm{\Delta }E(L)`$ is also known exactly for each sample and the errors in $`\mathrm{\Delta }E(L)`$ are $`O(N^{1/2}L^{d1})`$ where $`N`$ is the number of samples of size $`L`$ in $`d`$ dimensions. If the energy minima $`E_{a,b}(L)`$ are not found exactly, a crude estimate of the errors in $`\mathrm{\Delta }E(L)`$ is $`O(N^{1/2}L^d)`$ but this is certainly too low as failure of the algorithm to find the true minima because of being trapped in a metastable state of energy $`E>E_0`$ will cause systematic errors of unknown magnitude. Empirically, we find that this can readily cause errors larger than $`\mathrm{\Delta }E(L)`$ which makes the data point useless. This is most likely to happen for large $`L`$ because the CPU time required grows uncontrollably, as do the errors, so the large $`L`$ data becomes unreliable. This technical difficulty limits the accessible sizes $`L`$ to small values as one must keep errors in individual data points small.
We are forced to conclude that the accessible sizes $`L`$ are limited by the necessity of finding essentially exact global energy minima of each of a number $`N`$ of samples subject to certain, yet to be defined, BC. To our knowledge, there is no algorithm applicable to the systems of interest which will find exact minima in polynomial time, such as the branch and cut algorithmrieger  for the $`2D`$ Ising spin glass or numerically exact combinatorial optimization algorithms combin  for gauge and vortex glass models in the infinite screening limit, so we have to live with the fact that our problem is NP complete and the required CPU time explodes as $`L`$ increases. We use simulated annealinganneal1 ; anneal2  to estimate the lowest energies, which seems considerably more efficient than simple quenching to $`T=0`$, but we are unable to go beyond $`L=7`$ in $`3D`$ and $`L=10`$ in $`2D`$. We wish to extract the stiffness exponent $`\theta `$ from the scaling ansatz of Eq. (4) with a single power law and this makes sense only if the errors on individual data points are very small and the fit to the assumed scaling form is extremely good. In our opinion, the only sensible strategy is to obtain very accurate estimates of $`\mathrm{\Delta }E(L)`$ for the limited sizes $`L`$ which are feasible for the computer power available.
In the phase representation of Eq. (1), the configuration space to be searched for the global energy minima $`E_{a,b}(L)`$ is rather large as the phases $`\theta _i(0,2\pi ]`$ are continuous variables. Searching this space in finite time is not feasible as most of the allowed configurations of the $`\theta _i`$ are not even local energy minima. It is well known that random $`XY`$ models with a Hamiltonian of Eq. (1) can be written in a Coulomb gas (CG) or vortex representation via a duality transformationduality1 ; duality2 ; nh  which leaves the partition function invariant. This expresses the Hamiltonian in terms of charge or vortex configurations which are already local energy minima kt ; jmk . Thus, a reformulation of the Hamiltonian of Eq. (1) as a CG performs a partial minimization. A further minimization of the CG Hamiltonian corresponds to searching the much smaller space of local minima. Reformulating the problem of Eq. (1) including the BC in CG language is clearly a worthwhile exercise as it dramatically reduces the number of configurations over which we have to minimize, despite introducing long-ranged Coulomb interactions between vortices. The transformation is carried out in Section III for the model of Eq. (1) in $`3D`$.
The final problem is to define what is meant by a domain wall and the BC needed to induce a wall in a finite system of size $`L`$ in $`d`$ dimensions. We imagine the system of Eq. (1) on a torus in $`2D`$ or a hypertorus in $`3D`$, which corresponds to imposing periodic BC in the phases $`\theta _{𝐢+L\widehat{𝐞}_\mu }=\theta _𝐢`$ where $`\widehat{𝐞}_\mu `$ is a unit vector in the direction $`\mu =x,y,\mathrm{},d`$ and $`𝐢=(i_x,\mathrm{},i_d)`$ with $`i_\mu =(1,2,\mathrm{},L)`$. The phases at corresponding sites $`(i,j)`$ on opposite faces are coupled by some interaction $`\stackrel{~}{V}(\theta _i,\theta _j,A_{ij})`$ which may be regarded as defining the BC. In principle, the GS is obtained by minimizing the energy with respect to the $`L^d`$ bulk variables $`\theta _i`$ and all forms of $`\stackrel{~}{V}`$. This program is beyond our ability and we restrict ourselves to those $`\stackrel{~}{V}`$ which induce a spin or chiral defect, which are related to the continuous and reflection symmetry respectively. To impose a spin defect, we choose $`\stackrel{~}{V}=V(\theta _i\theta _jA_{ij})`$. The plaquettes between the opposite faces are indistinguishable from the others and play no special role. We therefore keep fixed the frustrations $`f_𝐫=_{P(𝐫)}A_{ij}/2\pi `$ where the sum is over the bonds in a clockwise direction of the elementary plaquette centered at $`𝐫`$. We still have the freedom to add $`\mathrm{\Delta }_\mu `$ to every bond in the $`\mu `$ direction between opposite faces which imposes a global twist $`\mathrm{\Delta }_\mu `$ in the phase round a loop circling the hypertorus in the direction $`\mu `$. This is equivalent to a gauge transformation $`A_{ij}A_{ij}+\mathrm{\Delta }_\mu /L`$ on every bond $`ij`$ in the direction $`\mu `$. The lowest energy, $`E_0(\mathrm{\Delta }_\mu )`$, is $`2\pi `$ periodic in $`\mathrm{\Delta }_\mu `$ with a minimum at some $`\mathrm{\Delta }_\mu ^0`$ which depends on the sample. To induce a spin domain wall normal to $`𝐱`$, one changes the twists from their best twist (BT) values $`\mathrm{\Delta }_\mu ^0\mathrm{\Delta }_\mu ^0+\pi \delta _{\mu ,x}`$. The minimum energy subject to this constraint gives $`E_{sD}(L)`$, the energy of the system of size $`L`$ containing an extra spin defect. Note that $`E_{sD}E_0`$ for every sample but $`E_0`$ is not necessarily the absolute minimum as some other functional form of $`\stackrel{~}{V}`$ may give a lower energy. However, even if $`E_0`$ is not the true GS energy but is the energy of a state with some excitation from the GS, this method of inducing a spin defect ensures that any excitation in the BT configuration will also be present in the state with an extra spin domain wall so that $`\mathrm{\Delta }E_s^{BT}(L)E_{sD}(L)E_0(L)0`$ is not affected by these. It is convenient, but not necessary, to define the spin defect energy by a twist of $`\pi `$ from the BT value $`\mathrm{\Delta }_x^0`$. This choice yields the maximum defect energy $`\mathrm{\Delta }E(L)`$. Any other choice $`0<ϵ\pi `$ yields the same spin stiffness exponent $`\theta _s^{BT}`$ defined by

$$\mathrm{\Delta }E_s^{BT}(L,ϵ)=𝒜(ϵ)L^{\theta _s^{BT}}$$
(5)
The size $`ϵ`$ of the twist from the BT value $`\mathrm{\Delta }_x^0`$ affects only the amplitude $`𝒜(ϵ)`$ which is a maximum at $`ϵ=\pi `$.
A chiral domain wall is induced by imposing reflective BCsg-dw-kt  which means that corresponding sites $`(ij)`$ on opposite faces are connected by interactions $`\stackrel{~}{V}=V(\theta _i+\theta _jA_{ij})`$ which is equivalent to a reflection of the spins about some axis. In principle, one follows the procedure for a spin domain wall to obtain the chiral defect energy $`\mathrm{\Delta }E_c=E_{cD}E_0`$ where $`E_{cD}`$ is the minimum energy of the system with these modified interactions $`\stackrel{~}{V}`$ connecting opposite faces. However, there is no reason to expect $`E_{cD}>E_0`$ as the BC defining $`E_0`$ may trap a chiral defect in some samples and, in such cases, the modified interactions $`\stackrel{~}{V}`$ will cancel the chiral defect to give $`E_{cD}<E_0`$. This phenomenon has been observed previously in numerical simulations of the $`XY`$ spin glasssg-dw-kt ; sg-dw-mg . We therefore define the chiral defect energy as $`\mathrm{\Delta }E_c^{BT}(L)|E_{cD}(L)E_0(L)|`$ and the chiral stiffness exponent $`\theta _c^{BT}`$ by a finite size scaling ansatz analogous to Eq. (5)

$$\mathrm{\Delta }E_c^{BT}(L)L^{\theta _c^{BT}}$$
(6)
Note that the Hamiltonian of Eq. (1) is truly invariant under reflection $`\theta _i\theta _i`$ in the $`XY`$ spin glass case when $`A_{ij}=0,\pi `$ as $`A_{ij}=\pm \pi `$ are equivalent. The Hamiltonian with uniform distribution of $`A_{ij}`$ such as a gauge glass is not truly invariant under reflection which would also require $`A_{ij}A_{ij}`$ for the Hamiltonian to be invariant as it lacks the reflection symmetry.
However, in an $`XY`$ spin glass there are two possible types of order, spin glass order and a chiral glass order each with their own stiffness exponent of Eq. (5) and Eq. (6). Recently, an important prediction was made that $`\theta _s=\theta _c<0`$ for an $`XY`$ spin glass in dimension $`d<d_l`$ where $`d_l`$ is the lower critical dimensionnh . Although not rigorous, the arguments are very plausible and supported by analytic calculations on simple one dimensional systems in which $`\theta _s=\theta _c`$ exactlynhm ; tnh . There is a notable lack of analytic results in this field with which to test numerical simulations and to our knowledge this is the only one existing at present. We have checked our numerical method in $`d=2<d_l`$ and get agreement with the analytic prediction that $`\theta _s=\theta _c=0.37\pm 0.015`$ to within numerical uncertaintyxysg-ka . Assuming the conjecture nh  is correct, this agreement gives some confidence in our definition of domain wall energies as discussed above and in our numerical method in $`d=3`$ using the CG representation. There is no analogous equality of the stiffness exponents in a $`XY`$ spin glass for $`d=3>d_l`$ so we do not attempt to estimate $`\theta _c`$ in $`3D`$ but concentrate on the spin stiffness exponent $`\theta _s`$. Also, at present, we are unable to derive an expression for a CG Hamiltonian with reflective boundary condition in $`d=3`$. All previous work on the $`XY`$ gauge glassgg-rtyf ; gg-g ; gg-ks ; gg-mg  and on the $`XY`$ spin glasssg-dw-kt ; sg-dw-k ; sg-dw-mg  using the $`T=0`$ DWRG method have used different definitions for domain wall energies. Minimization with respect to the global twists $`\mathrm{\Delta }_\mu `$ is omitted, the lowest energy with $`\mathrm{\Delta }_\mu =0`$ is called $`E_p`$ and the lowest energy with $`\mathrm{\Delta }_\mu =\pi `$ is called $`E_{ap}`$. Neither of these BC is compatible with the GS configuration as both must induce some excitation from $`E_0`$. Nevertheless, the spin defect energy is defined by $`\mathrm{\Delta }E_s^{RT}|E_{ap}E_p|`$ and the spin stiffness exponent $`\theta _s^{RT}`$ by

$$\mathrm{\Delta }E_s^{RT}(L)L^{\theta _s^{RT}}$$
(7)
We call this a random twist (RT) measurement as, for a particular sample, the twists $`\mathrm{\Delta }_\mu =0,\pi `$ are two arbitrary random choices relative to the best twist $`\mathrm{\Delta }_\mu ^0`$, which is the twist which yields the lowest energy. In a uniform ferromagnet, $`\mathrm{\Delta }_\mu ^0=0`$ which is realized by periodic BC and $`\mathrm{\Delta }_\mu ^0+\pi `$ by antiperiodic BC. For a particular realization of randomness, $`\mathrm{\Delta }_\mu ^0`$ is the analogue of periodic BC in a uniform ferromagnet.

## III Transformation to Coulomb Gas Representation

In this section, we discuss the CG representation of the Hamiltonian of Eq. (1) including all finite size contributions. This representation parameterizes the energy in terms of the topological excitations on a torus in $`2D`$ and a hypertorus in $`3D`$ and includes global excitations which wind round the whole hypertorus. These latter excitations are very important for a finite system and are vital for finite size scaling considerations when one is limited to small system sizes $`L`$. Also, every allowed configuration of topological excitations is a local energy minimum kt ; jmk  as spin wave excitations decouple from the vortex excitations, which allows us to obtain a more accurate estimate of energy minima than using the phase representation of Eq. (1) with the limited CPU time available. The transformation of the two dimensional $`XY`$ model to the CG representation including boundary contributions has been discussed in detail in earlier worksnh ; duality1 ; duality2 . In this section, we use the method of Ney-Nifle and Hilhorstnh  to transform to the CG representation in $`3D`$.
We first replace the potential $`V(\varphi )`$ in Eq. (1) by a piecewise parabolic potential which is equivalent to a Villainvillain  potential at $`T=0`$. The partition function for a $`L\times L\times L`$ gauge glass model in $`3D`$ is

$$Z=_\pi ^{+\pi }\underset{i}{}d\theta _i\underset{\left\{n_{ij}\right\}}{}\mathrm{exp}\left[\beta J\underset{<ij>}{}\left(\theta _{ij}A_{ij}\right)^2\right]$$
(8)
where $`\theta _{ij}\theta _i\theta _j2\pi n_{ij}`$ and where $`n_{ij}=n_{ji}`$ are integers on the bond $`<ij>`$. By choosing one phase, $`\theta _0`$, as a reference, the partition function can be written as

$`Z`$ $`={\displaystyle _\pi ^\pi }𝑑\theta _0{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \underset{<ij>}{}}d\theta _{ij}\mathrm{exp}\left[\beta J\left(\theta _{ij}A_{ij}\right)^2\right]`$ (11)
$`\times {\displaystyle \underset{𝐫}{}}\delta ({\displaystyle \underset{P(𝐫_{xy}}{}})\theta _{ij}mod2\pi )\delta ({\displaystyle \underset{P(𝐫_{yz}}{}})\theta _{ij}mod2\pi )`$
$`\times {\displaystyle \underset{𝐫}{}}\delta ({\displaystyle \underset{P(𝐫_{zx}}{}})\theta _{ij}mod2\pi \left)\delta \right({\displaystyle \underset{x}{}}\theta _{ij}mod2\pi )`$
$`\times {\displaystyle \underset{𝐫}{}}\delta ({\displaystyle \underset{y}{}}\theta _{ij}mod2\pi )\delta ({\displaystyle \underset{z}{}}\theta _{ij}mod2\pi )`$
Here, $`𝐫`$ is the coordinate of the center of an elementary cube of the original lattice which corresponds to the coordinate of a dual lattice site. $`𝐫_{xy}`$ is the coordinate of the center of the elementary plaquette in the $`xy`$ plane and similarly for $`𝐫_{yz}`$ and $`𝐫_{zx}`$. Note that for a $`2D`$ system in the $`xy`$ plane, $`𝐫_{xy}`$ are the dual lattice sites $`𝐫`$. Since each cube has six faces (plaquettes), each of which is shared by two adjacent cubes, to each dual lattice site $`𝐫`$ we assign three independent plaquettes with centers at $`𝐫_{xy}`$, $`𝐫_{yz}`$, and $`𝐫_{zx}`$ as shown in Figs. (1) and (2). $`_{P(𝐫_{xy}})\theta _{ij}`$ is the circulation of $`\theta _{ij}`$ round the plaquette in the $`xy`$ plane of the cube at $`𝐫`$,

$$\underset{P(𝐫_{xy})}{}\theta _{ij}\left[\theta _{ij_x}+\theta _{j_xi^{^{}}}\theta _{j_yi^{^{}}}\theta _{ij_y}\right]$$
(12)
and $`_x\theta _{ij}`$ is the circulation of $`\theta _{ij}`$ along an arbitrary loop in the $`x`$-direction round the hypertorus,

$$\underset{x}{}\theta _{ij}\underset{i_x=1}{\overset{L}{}}\theta _{(i_x,i_y,i_z),(i_x+1,i_y,i_z)}$$
(13)
where $`i_\mu `$ with $`\mu =x,y,z`$ is the coordinate of the original lattice site and $`i_y`$ and $`i_z`$ are fixed. Other summations are defined similarly. Note that one needs to consider only one global loop on the hypertorus in each direction. Circulations round other global loops can be expressed in terms of circulations round any three chosen global loops and round elementary plaquettes. It is clear from the definition of $`\theta _{ij}`$ and periodic boundary condition imposed on the $`\theta _i`$, these circulations are integer multiples of $`2\pi `$. Since the delta functions can be rewritten as follows,

$`\delta `$ $`\left({\displaystyle \underset{P(𝐫_{xy})}{}}\theta _{ij}mod2\pi \right)={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{n_𝐫^z=\mathrm{}}{\overset{\mathrm{}}{}}}\mathrm{exp}\left[in_𝐫^z{\displaystyle \underset{P(𝐫_{xy})}{}}\theta _{ij}\right]`$ (14)
$`\delta `$ $`\left({\displaystyle \underset{x}{}}\theta _{ij}mod2\pi \right)={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{n_x=\mathrm{}}{\overset{\mathrm{}}{}}}\mathrm{exp}\left[in_x{\displaystyle \underset{x}{}}\theta _{ij}\right]`$
the partition function now becomes

$`Z`$ $`=Z_0{\displaystyle \underset{𝐧}{}}{\displaystyle \underset{\{𝐧_𝐫\}}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \underset{<ij>}{}}d\theta _{ij}\mathrm{exp}\left[\beta J(\theta _{ij}A_{ij})^2\right]`$ (15)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \underset{𝐫}{}}\left(in_𝐫^x{\displaystyle \underset{P(𝐫_{yz})}{}}\theta _{ij}+in_𝐫^y{\displaystyle \underset{P(𝐫_{zx})}{}}\theta _{ij}+in_𝐫^z{\displaystyle \underset{P(𝐫_{xy})}{}}\theta _{ij}\right)\right]`$ (16)
$`\times `$ $`\mathrm{exp}\left[in_x{\displaystyle \underset{x}{}}\theta _{ij}+in_y{\displaystyle \underset{y}{}}\theta _{ij}+in_z{\displaystyle \underset{z}{}}\theta _{ij}\right]`$ (17)
where $`n_𝐫^\mu `$ and $`n_\mu `$ with $`\mu =x,y,z`$ are integers, and $`_𝐧_\mu _{n_\mu }`$ and $`_{\{𝐧_𝐫\}}_\mu _{n_𝐫^\mu }`$ The sum $`_𝐫`$ is over the dual lattice sites $`𝐫`$ at the centers of the elementary cubes and $`Z_0`$ is an unimportant constant. We use the notation $`_{P(𝐫_{xy})}`$ to denote a sum over the bonds, in a clockwise direction, of the plaquette in the $`xy`$ plane centered at $`𝐫_{xy}`$ and $`_x`$ to denote a sum over bonds on a global loop round the hypertorus in the $`x`$ direction. To perform the integrations over $`\{\theta _{ij}\}`$, we choose the three global loops around the hypertorus as shown in Fig. (3).

To deal with these global loops mathematically, we introduce the following quantities,

$$\delta _𝐫^x\{\begin{array}{cc}1\hfill & \text{if }𝐫=(x,1,1)\text{ with }x=1,2,\mathrm{},L\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(18)
similarly for $`\delta _𝐫^y`$ and $`\delta _𝐫^z`$. For example, the cubes at $`𝐫=`$ $`(x,1,1)`$ have a part of the global loop in the $`x`$ direction.
The definition of $`A_𝐫^\mu `$ associated with the cube at $`𝐫`$ is also shown in Fig. (1). As the plaquettes, one can assign three independent $`A_𝐫^\mu `$ to each cube. After relabeling $`A_{ij}`$ by $`A_𝐫^\mu `$, performing the integrations in Eq. (17) over $`\theta _{ij}`$, the partition function becomes

$`Z`$ $`=Z_0{\displaystyle \underset{𝐧,𝐧_𝐫}{}}\mathrm{exp}\left[{\displaystyle \frac{1}{4\beta J}}{\displaystyle \underset{𝐫}{}}{\displaystyle \underset{\alpha }{}}\left\{\left(\times 𝐧_𝐫\right)_\alpha +\delta _𝐫^\alpha n_\alpha \right\}^2\right]`$ (19)
$`\times `$ $`\mathrm{exp}\left[i{\displaystyle \underset{𝐫}{}}{\displaystyle \underset{\alpha }{}}\left\{\left(\times 𝐧_𝐫\right)_\alpha +\delta _𝐫^\alpha n_\alpha \right\}A_𝐫^\alpha \right]`$ (20)
We use the following notation for the discrete derivative,

$$(\times 𝐧_𝐫)_x(n_𝐫^yn_{𝐫\widehat{𝐳}}^y)(n_𝐫^zn_{𝐫\widehat{𝐲}}^z)$$
(21)
and similarly for other components of $`\times 𝐧_𝐫`$. Note that, when $`𝐫=(1,y,z)`$, then $`𝐫\widehat{𝐱}=(L,y,z)`$ due to periodic boundary conditions, and similarly for $`𝐫\widehat{𝐲}`$ and $`𝐫\widehat{𝐳}`$. Applying the Poisson summation formula to $`n_𝐫^\mu `$,

$$\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}f(n)=\underset{q=\mathrm{}}{\overset{\mathrm{}}{}}_{\mathrm{}}^{\mathrm{}}𝑑ue^{2\pi iqu}f(u),$$
(22)
the partition function of Eq. (20) becomes

$`Z`$ $`=Z_0{\displaystyle \underset{𝐧}{}}{\displaystyle \underset{\{𝐪_𝐫\}}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \underset{𝐫}{}}d𝐮_𝐫\mathrm{exp}\left[2\pi i{\displaystyle \underset{𝐫}{}}𝐪_𝐫𝐮_𝐫\right]`$ (23)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{4\beta J}}{\displaystyle \underset{𝐫}{}}{\displaystyle \underset{\alpha =x,y,z}{}}\left\{\left(\times 𝐮_𝐫\right)_\alpha +\delta _𝐫^\alpha n_\alpha \right\}^2\right]`$
$`\times `$ $`\mathrm{exp}\left[i{\displaystyle \underset{𝐫}{}}{\displaystyle \underset{\alpha =x,y,z}{}}\left\{\left(\times 𝐮_𝐫\right)_\alpha +\delta _𝐫^\alpha n_\alpha \right\}A_𝐫^\alpha \right]`$
It is convenient to change integration variables to $`v_𝐫^\mu `$ to make Eq. (23) more symmetric

$`v_𝐫^x`$ $``$ $`u_𝐫^x{\displaystyle \frac{y}{2L}}n_z+{\displaystyle \frac{z}{2L}}n_y`$
$`v_𝐫^y`$ $``$ $`u_𝐫^y{\displaystyle \frac{z}{2L}}n_x+{\displaystyle \frac{x}{2L}}n_z`$ (24)
$`v_𝐫^z`$ $``$ $`u_𝐫^z{\displaystyle \frac{x}{2L}}n_y+{\displaystyle \frac{y}{2L}}n_x`$
where $`𝐫=(x,y,z)`$. The partition function of Eq. (23) now becomes

$`Z`$ $`=Z_0{\displaystyle \underset{𝐧}{}}{\displaystyle \underset{\{𝐪_𝐫\}}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \underset{𝐫}{}}d𝐯_𝐫\mathrm{exp}\left[2\pi i{\displaystyle \underset{𝐫}{}}𝐪_𝐫𝐯_𝐫\right]`$ (25)
$`\times `$ $`\mathrm{exp}[{\displaystyle \frac{\pi i}{L}}{\displaystyle \underset{𝐫}{}}𝐪_𝐫(𝐫\times 𝐧\}]`$
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{4\beta J}}{\displaystyle \underset{𝐫}{}}{\displaystyle \underset{\alpha =x,y,z}{}}\left\{\left(\times 𝐯_𝐫\right)_\alpha +{\displaystyle \frac{n_\alpha }{L}}\right\}^2\right]`$
$`\times `$ $`\mathrm{exp}\left[i{\displaystyle \underset{𝐫}{}}{\displaystyle \underset{\alpha =x,y,z}{}}\left\{\left(\times 𝐯_𝐫\right)_\alpha +{\displaystyle \frac{n_\alpha }{L}}\right\}A_𝐫^\alpha \right]`$
The terms linear in $`n_\mu `$ in the exponent of Eq. (25) vanish due to the periodic boundary condition $`_𝐫\left(v_𝐫^xv_{𝐫\widehat{𝐲}}^x\right)=0`$. Eq. (25) can be simplified by introducing frustration variables $`f_𝐫^\mu `$ at site $`𝐫`$

$`f_𝐫^x`$ $``$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{P(𝐫_{yz})}{}}A_{ij}`$ (26)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}\left(A_𝐫^yA_{𝐫+\widehat{𝐳}}^y+A_{𝐫+\widehat{𝐲}}^zA_𝐫^z\right)`$
from which $`f_𝐫^y`$ and $`f_𝐫^z`$ are obtained by cyclic permutation of $`xyz`$. After some algebra, the partition function of Eq. (25) becomes a form suitable for integration over $`v_𝐫^\mu `$

$`Z`$ $`=`$ $`Z_0{\displaystyle \underset{𝐧}{}}{\displaystyle \underset{\{𝐪_𝐫\}}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \underset{𝐫}{}}d𝐯_𝐫^x\mathrm{exp}\left[{\displaystyle \frac{2\pi i}{2L}}\left\{𝐪_𝐫\left(𝐫\times 𝐧\right)\right\}\right]`$ (27)
$`\times \mathrm{exp}\left[{\displaystyle \frac{1}{4\beta J}}{\displaystyle \underset{𝐫}{}}{\displaystyle \underset{\alpha =x,y,z}{}}\left\{\left(\times 𝐯_𝐫\right)_\alpha \right\}^2\right]`$
$`\times \mathrm{exp}\left[2\pi i{\displaystyle \underset{𝐫}{}}(𝐪_𝐫𝐟_𝐫)𝐯_𝐫\right]`$
To evaluate the integrals over $`v_𝐫^\mu `$ in Eq. (27) it is convenient to take the Fourier transform

$$𝐯_𝐫=L^{\frac{3}{2}}\underset{𝐤}{}e^{i𝐤𝐫}\stackrel{~}{𝐯}(𝐤)$$
(28)
where $`𝐤=(k_x,k_y,k_z)`$ and $`k_i=\frac{2\pi }{L}m_i`$ with $`m_i=0,1,\mathrm{},L1`$. We also decompose $`\stackrel{~}{𝐯}(𝐤)`$ into longitudinal, $`\stackrel{~}{𝐯}_L(𝐤)`$, and transverse, $`\stackrel{~}{𝐯}_T(𝐤)`$, components where

$`\stackrel{~}{v}_L^\alpha (𝐤)`$ $`=`$ $`{\displaystyle \underset{\beta }{}}{\displaystyle \frac{(1e^{ik_\alpha })(1e^{ik_\beta })}{\lambda _𝐤}}\stackrel{~}{v}^\beta (𝐤)`$
$``$ $`{\displaystyle \underset{\beta }{}}L_{\alpha \beta }\stackrel{~}{v}^\beta (𝐤)`$
$`\stackrel{~}{v}_T^\alpha (𝐤)`$ $``$ $`{\displaystyle \underset{\beta }{}}(\delta _{\alpha \beta }L_{\alpha \beta })\stackrel{~}{v}^\beta (𝐤)`$
$``$ $`{\displaystyle \underset{\beta }{}}T_{\alpha \beta }\stackrel{~}{v}^\beta (𝐤)`$
$`\lambda _𝐤`$ $`=`$ $`62\mathrm{cos}k_x2\mathrm{cos}k_y2\mathrm{cos}k_z`$ (29)
define the longitudinal and transverse projection operators $`L_{\alpha \beta }(𝐤)`$ and $`T_{\alpha \beta }(𝐤)`$. In Fourier space, the integration over $`𝐯_𝐫`$ in Eq. (27) is straightforward

$`{\displaystyle \underset{𝐤}{}d\stackrel{~}{𝐯}_L(𝐤)d\stackrel{~}{𝐯}_T(𝐤)\mathrm{exp}\left[\frac{1}{4\beta J}\underset{𝐤}{}\lambda _𝐤\stackrel{~}{𝐯}_T(𝐤)^2\right]}`$
$`\times \mathrm{exp}\left[2\pi i{\displaystyle \underset{𝐤}{}}\left\{\stackrel{~}{𝐩}_T(𝐤)\stackrel{~}{𝐯}_T(𝐤)+\stackrel{~}{𝐩}_L(𝐤)\stackrel{~}{𝐯}_L(𝐤)\right\}\right]`$ (30)
where $`\stackrel{~}{𝐩}(𝐤)=𝐪(𝐤)𝐟(𝐤)`$. Integration over $`\stackrel{~}{𝐯}_L(𝐤)`$ gives $`\stackrel{~}{𝐩}_L(𝐤)=0`$. In real space, this is the condition that the discrete divergence of the vorticity (charge) $`𝐪_𝐫`$ at any dual lattice site $`𝐫`$ obeys $`𝐪_𝐫=0`$ since $`𝐟_𝐫=0`$. The integration over the transverse components $`\stackrel{~}{𝐯}_T(𝐤)`$ are simple Gaussian integrals and are easily performed. Integration over $`\stackrel{~}{𝐯}_L(0)`$ and $`\stackrel{~}{𝐯}_T(0)`$ yield the neutrality condition

$$\stackrel{~}{𝐩}(0)=\underset{𝐫}{}𝐩_𝐫=0$$
(31)
The final step in this rather technical derivation of the Hamiltonian in the CG representation is to apply the Poisson summation formula of Eq. (22) to eliminate the $`n_\mu `$ in favor of global vortices or charges $`q_{\mu 1}`$. The Gaussian integrations yield the partition function

$$Z=Z_0\underset{\{𝐪_𝐫\}}{}\underset{q_\mu }{}\mathrm{exp}\left[\beta H(𝐪_𝐫,𝐟_𝐫,q_{\mu 1},f_{\mu 1})\right]$$
(32)
The Hamiltonian $`H`$ is identified as

$`H`$ $`=`$ $`(2\pi )^2J{\displaystyle \underset{𝐫,𝐫^{}}{}}(𝐪_𝐫𝐟_𝐫)(𝐪_𝐫^{}𝐟_𝐫^{})G(𝐫𝐫^{})`$ (33)
$`+`$ $`{\displaystyle \frac{J}{2L}}\left\{{\displaystyle \frac{\pi }{L}}{\displaystyle \underset{𝐫}{}}\left(zp_𝐫^yyp_𝐫^z\right)+Q_x\right\}^2`$ (34)
$`+`$ $`{\displaystyle \frac{J}{2L}}\left\{{\displaystyle \frac{\pi }{L}}{\displaystyle \underset{𝐫}{}}\left(xp_𝐫^zzp_𝐫^x\right)+Q_y\right\}^2`$ (35)
$`+`$ $`{\displaystyle \frac{J}{2L}}\left\{{\displaystyle \frac{\pi }{L}}{\displaystyle \underset{𝐫}{}}\left(yp_𝐫^xxp_𝐫^y\right)+Q_z\right\}^2`$ (36)
where $`G(𝐫)=L^3_{𝐤0}(\mathrm{exp}(i𝐤𝐫)1)/\lambda _𝐤`$ is the lattice Green’s function on a simple cubic lattice in $`3D`$ with periodic BC and the quantities $`Q_\mu `$ are

$$Q_x=\pi \underset{𝐫}{}(zp_𝐫^y\delta _{y,1}yp_𝐫^z\delta _{z,1})+2\pi L(q_{x1}f_{x1})$$
(37)
with $`Q_y`$ and $`Q_z`$ obtained from $`Q_x`$ by cyclic permutations of $`xyz`$. In Eq. (37), $`f_{x1}`$ is the circulation of $`A_{ij}`$ along the chosen global loop round the hypertorus in the $`x`$ direction

$$2\pi f_{x1}=\underset{x=1}{\overset{L}{}}A_{𝐫=(x,1,1)}^x$$
(38)
and similarly for $`f_{y1}`$ and $`f_{z1}`$. The integers $`q_{\mu 1}`$ are interpreted as circulations of the phase round the three independent global loops encircling the hypertorus.

Eqs. (36) and (37) are the main results of this section. These expressions give the energy of a system of size $`L`$ on a hypertorus. To find the minimum energy $`E_0`$ the Hamiltonian $`H=H(𝐪_𝐫,𝐟_𝐫,q_{\mu 1},f_{\mu 1})`$ is minimized with respect to the bulk integer valued vector charges $`𝐪_𝐫`$, the integer valued global winding numbers $`q_{\mu 1}`$ and the global frustrations $`f_{\mu 1}`$. In the case of the $`XY`$ spin glass the $`q_{\mu 1}f_{\mu 1}`$ of Eq. (37) are restricted to be integer or half integer while for the gauge glass this can have any real value. Note that minimizing with respect to $`q_{x1}`$ and $`f_{x1}`$ is exactly minimizing with respect to the twist $`\mathrm{\Delta }_x`$ and the best twist $`\mathrm{\Delta }_\mu ^0`$ corresponds to $`f_{\mu 1}^0`$, the value of $`f_{\mu 1}`$ at the energy minimum. Adding a global twist $`\mathrm{\Delta }_\mu `$ to the phases is exactly equivalent to changing the global frustration from its original value $`f_{\mu 1}f_{\mu 1}+\mathrm{\Delta }_\mu /2\pi `$. As discussed in Section II, a spin domain wall is induced by $`\mathrm{\Delta }_x^0\mathrm{\Delta }_x^0+\pi `$ which, in the CG representation, is $`f_{x1}^0f_{x1}^0+1/2`$. Note that the frustrations $`𝐟_𝐫`$ are given in terms of the random $`A_{ij}`$ and are kept fixed during the minimization.
To confront our numerical results with the only existing analytic predictionnh  we need expressions for both spin and chiral domain wall energies in dimension $`d=2<d_l`$. The CG representation of the $`XY`$ spin and gauge glasses including all finite size corrections have been known for some time in $`2D`$duality2 ; nh  and, for the sake of completeness, we quote the necessary results below. To study the spin stiffness, we join corresponding sites on opposite faces by $`\stackrel{~}{V}=V(\theta _i\theta _jA_{ij})`$ and the CG Hamiltonian is

$`H`$ $`=`$ $`(2\pi )^2J{\displaystyle \underset{\mathrm{𝐫𝐫}^{}}{}}(q_𝐫f_𝐫)G(𝐫𝐫^{})(q_𝐫^{}f_𝐫^{})`$ (39)
$`+`$ $`{\displaystyle \frac{J}{L^2}}(\sigma _x^2+\sigma _y^2)`$ (40)
where

$`\sigma _x`$ $`=`$ $`2\pi \left[L(q_{x1}f_{x1})+{\displaystyle \underset{𝐫}{}}(q_𝐫f_𝐫)y\right]`$ (41)
$`\sigma _y`$ $`=`$ $`2\pi \left[L(q_{y1}f_{y1}){\displaystyle \underset{𝐫}{}}(q_𝐫f_𝐫)x\right]`$ (42)
$`G(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}{\displaystyle \underset{𝐤0}{}}{\displaystyle \frac{e^{i𝐤𝐫}1}{42\mathrm{cos}k_x2\mathrm{cos}k_y}}`$ (43)
Here, as in $`3D`$, $`𝐫=(x,y)`$ represents the coordinates of the dual lattice sites and $`G(𝐫)`$ is the lattice Green’s function with periodic BC on the $`\theta _i`$. From Eq. (8), we see that the difference of a factor $`2`$ in the prefactor of Eq. (40) from other works is in the coupling constant $`J`$. Also note that the Hamiltonian of Eq. (40) describes the $`XY`$ spin glass when the frustrations $`f_𝐫`$ and the global frustrations $`f_{\mu 1}`$ are restricted to $`(0,1/2)`$. In the gauge glass, they can have any real value. As in $`3D`$, the topological charges $`q_𝐫`$ are integers as are $`q_{\mu 1}`$ because of the periodic BC in the $`\theta _i`$.
The last piece of information we need is the CG representation of the Hamiltonian of the $`XY`$ spin glass in $`2D`$ with a chiral domain wall imposed. This is to be found in the paper of Ney-Nifle and Hilhorstnh . A single chiral domain wall is induced by joining opposite faces by interactions $`\stackrel{~}{V}=V(\theta _i+\theta _jA_{ij})`$ which is equivalent to imposing reflective BCsg-dw-kt ; nh . In turn, this is equivalent to doubling the size in (say) the $`x`$ direction to a $`2L\times L`$ lattice in which the extra half is a charge conjugated image of the other. This system has two chiral domain walls with Hamiltoniannh

$$H_R=2\pi ^2J\underset{𝐫,𝐫^{}}{}(q_𝐫f_𝐫)\stackrel{~}{G}(𝐫𝐫^{})(q_𝐫^{}f_𝐫^{})$$
(44)
where $`\stackrel{~}{G}(𝐫)`$ is the Green’s function for a $`2L\times L`$ square lattice with periodic BC and also with $`q_{𝐫+L\widehat{𝐱}}=q_𝐫`$ and $`f_{𝐫+L\widehat{𝐱}}=f_𝐫`$. Note that the sign reversal of the frustrations $`f_𝐫`$ is not necessary for the spin glass because $`f_𝐫=\pm 1/2`$ are equivalent.
The form of the energy of Eq. (44) is used in the simulations to estimate $`\theta _c`$ the chiral stiffness exponent as it is intuitively more transparent and more convenient than the corresponding expression with a single chiral domain wallnh . Unfortunately, we have been unable to derive the analogous expression in $`3D`$ to Eq. (44) so we have no independent estimate of $`\theta _c`$ sg-dw-mg  for the $`3D`$ $`XY`$ spin glass.

## IV Numerical Method

### IV.1 Minimization Algorithm

In Section II we argued that it is necessary to find the energies $`E_0(L)`$ and $`E_D(L)`$ essentially exactly for every sample to control the errors in the small domain wall energy $`\mathrm{\Delta }E(L)`$. To our knowledge, for the systems of interest no algorithm exists which can locate the global energy minima in polynomial time. We are left with two methods: (i) repeated simple quenches from an arbitrary initial configuration to $`T=0`$ followed by a downward slide to the nearest local minimum and (ii) simulated annealinganneal1  which is considerably more efficient anneal2 . By this we mean that, for the same CPU time, simulated annealing finds a lower energy than simple quenching.
We start with the system in some randomly chosen configuration and quench to some $`T_0`$, determined by trial and error, for each size $`L`$. Then we do a Monte Carlo (MC) sweep through the system by inserting charges $`q=1`$ and $`q^{}=1`$ on an arbitrary pair of nearest neighbor sites of the dual lattice in $`2D`$ and accepting or rejecting the move according to conventional MC rules. This has the effect of inserting new charges, annihilating charges or moving a charge by one lattice spacing while maintaining charge neutrality $`q_𝐫=0`$. In $`3D`$, the elementary excitation is a loop of charge round an elementary square with vertices at dual lattice sites. This maintains $`𝐪_𝐫=0`$. The closed loop of charge can lie in any of the three orthogonal planes of the cubic lattice and the charge can circulate round the loop in either direction, making $`6`$ possibilities with the center of the loop at a fixed but randomly selected position. The temperature $`T`$ is then reduced to $`T_1=\alpha T_0`$ with $`\alpha <1`$ whose value is again determined by trial and error and the procedure iterated a large number, $`N`$, of times to reach a lowest temperature $`T_N=\alpha ^NT_00`$. Of course, the system may be trapped in a deep metastable well with barriers too high for the MC passes to overcome so the whole annealing sequence is repeated $`M`$ times from different random initial configurations and the lowest energy out of all the $`NM`$ trials recorded. Again, this does not guarantee that the global minimum energy is found but this method does have a few checks built in. At the crudest level, the best twist condition $`E_{sD}(L)E_0(L)`$ must be obeyed for each sample since $`E_0(L)`$ is, by construction, the lowest energy of the system subject to the BC given by the interactions $`\stackrel{~}{V}=V(\theta _i\theta _jA_{ij})`$ across opposite faces. $`E_{sD}(L)`$ is obtained by $`f_{x1}^0f_{x1}^0+1/2`$ where $`f_{\mu 1}^0`$ is the value of $`f_{\mu 1}`$ which makes the boundary terms of Eqs. (36) and (40) vanish, corresponding to the best twist $`\mathrm{\Delta }_\mu ^0`$. It is clear that $`E_{sD}(L)`$ is the energy of the system containing an extra spin domain wall compared to the system with energy $`E_0(L)`$. As discussed earlier, $`E_0(L)`$ is not necessarily the true GS energy as the system may contain some chiral domain walls. However, by construction, $`E_{sD}(L)`$ is the energy of the system with the same chiral defects and an extra spin defect.
If any sample violates the BT condition, clearly the annealing is not sufficient and one can either increase the number $`MN`$ of annealing attempts or just discard that sample. Increasing $`MN`$ involves a significant increase in CPU time particularly for the larger sizes $`L`$ so the choice depends on the time available. However, even if the BT condition is satisfied for every sample, there is no guarantee that the lowest energies found are true global minima. To improve the chances that true minima are achieved, where possible we did two independent simulations on identical samples with different pseudo random number sequences and, if the same energy minima are found in both simulations, this is defined to be the true minimum. This procedure is very expensive in CPU time and our resources did not permit this last check to be performed for every sample, particularly for the largest systems. This check was done for at least a few randomly selected samples, except for $`L=7`$ in the $`3D`$ gauge glass. Averaging over disorder was performed over as many samples as possible with the aim of making the uncertainty in the mean domain wall energy $`\mathrm{\Delta }E(L)`$ less than $`3\%`$ which requires averaging over at least $`10^3`$ samples. The first set of checks are done on every sample which makes it fairly probable that exact minima $`E_0(L)`$ and $`E_D(L)`$ are obtained. We can assume with some confidence that the uncertainty in $`\mathrm{\Delta }E(L)`$ is purely statistical and $`O(N^{1/2})`$. Averaging over about $`10^3`$ samples leads to an acceptable uncertainty of about $`3\%`$ in $`\mathrm{\Delta }E(L)`$.

### IV.2 $`XY`$ Spin Glass in $`2D`$

This particular system with a Hamiltonian of Eq. (1) in $`2D`$ is not particularly interesting in the sense that there exists no finite temperature transition. It has been known for a long time that $`d_l>2`$ and both spin and chiral stiffness exponents are negative. However, the situation has been controversial due to the possible existence of different stiffness or correlation length exponents for spin and chiral glass order in $`2D`$. Previous estimates of the spin and chiral stiffness exponents are summarized by $`\theta _s2\theta _c0.78`$ based on extensive numerical simulations such as DWRG and finite temperature MC simulations sg-dw-kt ; sg-dw-mg ; sg-dw-k ; sg-by ; sg-rm . However, Ney-Nifle and Hilhorstnh  made a non-rigorous but very plausible conjecture based on analytic considerations that for dimension $`dd_l`$, $`\theta _s=\theta _c0`$. This is supported by exact analytic results on simple models (i) a $`XY`$ spin glass on a ladder lattice nhm  where the common correlation length exponent $`\nu =|\theta |^1=0.5263\mathrm{}`$ and (ii) the $`XY`$ spin glass on a tube lattice tnh  with $`\nu =0.5564\mathrm{}`$. The key disagreement between numerical estimates and analytic theory is in the conjectured equality of $`\theta _s`$ and $`\theta _c`$. If one makes the hypothesis that earlier estimates are in error, then the most likely reason is that the simulations are computing the wrong quantity when estimating one or both of the stiffness exponents. It is extremely unlikely that all the simulations are in error for simple technical reasons as all find the same values of $`\theta _s`$ and $`\theta _c`$, but $`\theta _s\theta _c`$.
To test this hypothesis, we have performed simulations of the $`2DXY`$ spin glass in the CG representation using Eq. (40) to estimate the spin stiffness exponent $`\theta _s`$ using both BT and RT measurements with the results $`\theta _s^{BT}=0.37\pm 0.015`$ and $`\theta _s^{RT}=0.76\pm 0.015`$ xysg-ka . These numbers were obtained from sizes $`L=4,5,6,7,8,10`$ averaging over $`2560`$ samples for $`L8`$ and $`1152`$ samples for $`L=10`$. As expected, the value of $`\theta _s^{RT}`$ agrees with all previous estimates sg-dw-kt ; sg-dw-mg , all of which use the RT measurement in some form. Both the BT and the RT data fit the scaling ansatz of Eq. (4) equally well and some other information is needed to decide which value of $`\theta _s`$, if either, is correct. Both cannot be correct as both are supposed to measure the same quantity. The necessary information is in the chiral stiffness exponent $`\theta _c`$ which we measure by simulating Eq. (44) on a $`2L\times L`$ system. Again, both BT and RT measurements were made with the same range of $`L`$ and the same number of samples as for $`\theta _s`$ with the results $`\theta _c^{BT}=0.37\pm 0.010`$ and $`\theta _c^{RT}=0.37\pm 0.015`$.xysg-ka  At first sight it is surprising that both measurements give the same value for $`\theta _c`$ to within numerical uncertainty while the values of $`\theta _s^{RT}`$ and $`\theta _c^{BT}`$ differ by a factor of $`2`$. Note that the boundary terms $`\sigma _\alpha `$ of Eq. (40) which contain the twist parameter $`f_{\mu 1}`$ vanish in the BT condition and such boundary terms are absent from Eq. (44). Thus, any measurement with a CG representation in which boundary contributions are absent is automatically a BT measurement. Since $`\theta _c=\theta _s^{BT}0.37`$ and $`\theta _c\theta _s^{RT}0.76`$, assuming that the conjecture is correct, we conclude that the BT measurement yields a reasonable estimate of the true $`\theta _s`$ while the commonly used RT measurement yielding $`\theta _s^{RT}\theta _c`$ is not an appropriate method for the small values of $`L`$ accessible at present. Our simulation results are shown in Fig. (4).

We do not understand what, if anything, $`\theta _s^{RT}`$ means despite the apparent excellent fit of $`\mathrm{\Delta }E_s^{RT}(L)`$ to the scaling ansatz as this is the energy difference of the system subject to two random choices of BC and we can see no reason why it should scale as $`L^{\theta _s}`$ over decades of $`L`$.

## V RESULTS

Our investigation of the $`2DXY`$ spin glass establishes that the BT method of measuring domain wall energies can reproduce the one and only available analytic prediction nh  which is relevant for our purposes. This gives some encouragement to venture into areas where no such analytic guide exists. In this Section, we report some new results on the spin stiffness exponent $`\theta _s`$ for the $`XY`$ spin glass in $`3D`$ (V.1), and the gauge glass in both $`2D`$ and $`3D`$ (V.2). We also perform simulations on systems with varying strengths of disorder (V.3) where we study the renormalization group flows for both the coupling constants $`J(L)`$ and disorder strength $`A(L)`$ with increasing length scale $`L`$. All our results are numerical and the exponent $`\theta _s=\theta _s^{BT}`$ is estimated by fitting estimates of $`\mathrm{\Delta }E_s^{BT}(L)`$ to the scaling ansatz of Eq. (4). The system sizes $`L`$ are very small as they are severely constrained by the necessity of controlling the errors in both $`E_0(L)`$ and $`E_D(L)`$, which are estimated by independent simulations. However, the fit to the scaling ansatz is very good despite the very few data points for all cases we have studied and is comparable with the same procedure carried out on the $`2D`$ Ising ferromagnet. In this case, simulations for sizes $`L=2,3,4,5`$ are sufficient to reproduce the exact $`\theta _s=1`$ to high accuracy.

### V.1 $`XY`$ spin glass in $`3D`$

This system has been somewhat controversial for some years and is not yet settled. It has been believed that that $`d_l4`$ for spin glass order, and $`d_l<3`$ for chiral glass order and that, in $`2D`$ and $`3D`$, spin and chiral variables decouple and order separately sg-mc-kt ; sg-dw-kt ; sg-dw-k ; sg-dw-mg ; sg-by ; sg-rm . This allows for the widely accepted scenario that in $`3D`$ spin glass order sets in at $`T_{SG}=0`$ whereas chiral glass order sets in at $`T_{CG}>0`$. This scenario is based on MC simulations at finite $`T`$sg-dw-k ; sg-jy  and on $`T=0`$ defect energy scaling sg-dw-kt ; sg-dw-k  using the RT method. Attempts to show rigorously that $`T_{SG}=0`$ in $`3D`$ana-no  fail if reflection symmetry is brokensg-dw-k ; ana-sy ; ana-on  at finite $`T`$. The first cracks in this widely accepted scenario appeared recently when Maucourt and Grempelsg-dw-mg  published the results of a large scale defect energy scaling study of the $`3D`$ $`XY`$ spin glass model of Eq. (1). They used the RT method, as all previous defect energy scaling studies have done, with sizes $`L12`$ in $`2D`$ and $`L8`$ in $`3D`$. Although the fit to the scaling ansatz is not good due to strong crossover effects and large uncertainties in the large $`L`$ data which was used to estimate $`\theta _{s,c}^{RT}`$,it is clear that both $`\mathrm{\Delta }E_s^{RT}(L)`$and $`\mathrm{\Delta }E_c^{RT}(L)`$ are increasing with $`L`$ which implies that there is both spin and chiral glass order at sufficiently small $`T>0`$. However, as argued above a valid numerical method must yield $`\theta _s=\theta _c`$ in $`2D`$ while they obtain $`\theta _s2\theta _c0.78`$sg-dw-mg  in agreement with other estimates sg-dw-kt ; sg-dw-k . Finite $`T`$ simulations seem to suffer from severe equilibration difficulties sg-dw-k ; sg-jy  which makes any conclusions from them also suspect.
In view of the lack of reliable results for the $`3DXY`$ spin glass, we have done some preliminary simulations on very small systems with $`L=2,3,4,5,6`$ in the CG representation with Eq. (36) where $`f_{\mu 1}=0,1/2`$ in Eq. (37). Following the method outlined in Section IV.2 we estimate $`\theta _s^{BT}=+0.10\pm 0.03`$ while the data for $`\mathrm{\Delta }E_s^{RT}(L)`$ is clearly decreasing with $`L`$ for $`L=2,3,4`$, with roughly the same slope as found by Kawamura sg-dw-k  for the same range of $`L`$. Our value $`\theta _s=0.10\pm 0.03`$ xysg-ka  is to be compared with the recent estimate $`\theta _s^{RT}=+0.052\pm 0.03`$ for $`L=5,6,7`$ sg-dw-mg . Although these values are close numerically and are equal to within error bars, we are more inclined to believe in the former number as we consider $`\theta _s^{BT}`$ to give the true spin stiffness exponent and we speculate that further points will also lie on the scaling ansatz and will significantly reduce the $`30\%`$ uncertainty. The results are shown in Fig. (5).

Unfortunately, we have no estimate for the chiral stiffness exponent $`\theta _c^{BT}`$ for the $`3DXY`$ spin glass as we have been unable to derive the $`3D`$ analogue of Eq. (44) to estimate $`\mathrm{\Delta }E_c^{BT}(L)`$ and $`\mathrm{\Delta }E_c^{RT}(L)`$. However, because the numerical values of $`\theta _c^{BT}`$ and $`\theta _c^{RT}`$ are the same in $`2D`$ and equal to $`\theta _s`$ as required, we are reasonably safe in assuming that existing estimatessg-dw-k ; sg-dw-mg  of $`\theta _c`$ in $`3D`$ are fairly accurate. Even though our spin stiffness exponent $`\theta _s^{BT}+0.1`$ has a rather large uncertainty, it suggests that the lower critical dimension for spin glass order is $`d_l<3`$, as for chiral order. This is to be expected from analytic arguments nh .

### V.2 Gauge Glass in $`2D`$ and $`3D`$

We have also performed simulations on the gauge glass in the CG representation using Eq. (40) in $`2D`$ and Eq. (36) in $`3D`$. The only differences to the spin glass are in the values of $`f_{\mu 1}`$ and $`f_𝐫`$ which can have any value in the interval $`[1/2,1/2)`$. The frustrations $`f_𝐫=_{P(𝐫)}A_{ij}/2\pi `$ are correlated random variables as the $`A_{ij}`$ are the independent random variables. Similarly, in $`3D`$, the frustrations $`𝐟_𝐫`$ of Eq. (26) are correlated random variables with each component in the interval $`[1/2,1/2)`$. The three global frustrations $`f_{\mu 1}`$ can take any value in the same interval. In both spin and gauge glasses, the vorticities $`q_𝐫^\mu `$ and $`q_{\mu 1}`$ are integers. We have not computed $`\theta _c`$ in either $`2D`$ or $`3D`$ for the gauge glass, mainly because we have been unable to obtain the appropriate expression for the Hamiltonian in the CG representation with reflective BC in $`3D`$ and we are unable to understand what information such a simulation would yield. The major reason for doing such a simulation in the $`2D`$ spin glass case was to confirm that our procedure is a useful numerical method to estimate the actual spin stiffness exponent $`\theta _s`$ for both spin and gauge glasses in $`2D`$ and $`3D`$. To compare with earlier work gg-rtyf ; gg-g ; gg-ks ; gg-mg  we have also estimated $`\theta _s^{RT}`$ by keeping $`f_{\mu 1}`$ fixed at some random value during the minimization, exactly as for the spin glass. The results of the simulations gg-ka  in $`2D`$ for system sizes $`L=2,3,4,5,6,8,10`$ are shown in Fig. (6) for $`\mathrm{\Delta }E_s^{BT}(L)`$ and $`\mathrm{\Delta }E_s^{RT}(L)`$.

Averages were performed over about $`10^3`$ samples for each size $`L`$. In this case, all the checks discussed in section IV.1 were performed so we assume the errors are purely statistical. Both fit well to the scaling ansatz of Eq. (4) with very similar errors with the values $`\theta _s^{BT}=0.36\pm 0.010`$ and $`\theta _s^{RT}=0.45\pm 0.015`$. The latter value is consistent with all earlier estimates of $`\theta _s`$ gg-rtyf ; gg-g ; gg-mg  which is not a surprise as all these were done using the RT method in some form. However, as argued in section IV.2, this is not an accurate estimate of $`\theta _s`$ so this number is not the $`T=0`$ spin stiffness exponent. A more accurate estimate of this is $`\theta _s^{BT}=0.36\pm 0.010`$ which is significantly larger than $`\theta _s^{RT}`$, so that the $`2D`$ gauge glass has a longer correlation length $`\xi (T)T^{1/|\theta _s|}`$ than previously thought.
We have also obtained some estimates of $`\theta _s`$ in $`3D`$ by performing simulationsgg-ka  of the gauge glass in the CG representation using Eq. (36) with the distribution of frustrations $`𝐟_𝐫`$ appropriate for this case as determined by taking the $`A_{ij}`$ uniformly distributed in $`(\pi ,+\pi ]`$. The system sizes are $`2L7`$ with disorder averaging over $`10^3`$ samples for $`L5`$, $`300`$ for $`L=6`$ and $`60`$ for the largest size $`L=7`$. The uncertainty in $`\mathrm{\Delta }E_s^{BT}(L)`$ for $`L=7`$ is very large, but this data point is included to check that it is consistent with the behavior deduced from the more reliable data of the smaller sizes. The results are shown in Fig. (7) for $`\mathrm{\Delta }E_s(L)`$ for the unscreened gauge glass using both BT and RT measurements and for a gauge glass in $`3D`$ with screened interactions, $`\lambda <\mathrm{}`$. The BT data fit the scaling form of Eq. (4) very well for sizes $`L6`$ with exponent $`\theta _s^{BT}=+0.31\pm 0.010`$. If the $`L=7`$ point is also included in the fit, we obtain $`\theta _s^{BT}=+0.30\pm 0.015`$. These errors in $`\theta _s`$ come from a naive least squares fit of the data to a straight line and should not be taken too seriously. The $`L=7`$ data is suspect because $`4`$ samples violated the BT condition $`\mathrm{\Delta }E_s^{BT}(L)0`$ out of a total of only $`64`$, despite running highly vectorised code for a few thousand CPU hours on a Cray J90. Time did not permit any further checks for attaining the global energy minimum for $`L=7`$. We cannot be sure that the remaining $`60`$ samples which did not violate the BT condition reached their global energy minima nor that the energies $`E_{sD}`$ are determined sufficiently accurately. The error bar on the $`L=7`$ point in Fig. (7) assumes that the uncertainty in $`\mathrm{\Delta }E_s^{BT}(L=7)`$ is purely statistical and the true uncertainty is probably much larger. At least an order of magnitude more CPU time is needed for sufficient annealing to reach the true minima and to perform the additional simulations with different random number sequences to check that the minimization algorithm is successful. This is just for a single batch of $`64`$ samples and to reduce the uncertainty to $`3\%`$, yet another order of magnitude of CPU time would be needed to average over $`10^3`$ samples. This is totally out of reach with the computing resources available to us. What data we have is entirely consistent with the scaling form of Eq. (4) with $`\theta _s+0.30`$ with no sign of any deviation from this form.
The behavior of $`\mathrm{\Delta }E_s^{RT}(L)`$ is also shown in Fig. (7) for sizes $`L6`$ which is very much like the data obtained by earlier simulations. This clearly does not fit the scaling form of Eq. (4) for these small values of $`L`$, but if one insists on extracting a value of $`\theta _s^{RT}`$ from the data, one obtains consistency with previous estimates gg-g  for the spin stiffness exponent $`\theta _s^{RT}+0.05\pm 0.05`$. As can be seen from Fig. (7), this estimate has no meaning as the data clearly does not obey the scaling ansatz. In fact, as noticed by Maucourt and Grempelgg-mg , $`\mathrm{\Delta }E_s^{RT}(L)`$ seems to start increasing with $`L`$ for $`L>5`$ but, as we argue earlier, this may, or may not, eventually scale as $`\mathrm{\Delta }E_s^{BT}(L)`$ for sufficiently large $`L`$. Speculation along these lines is fruitless until computers which are many orders of magnitude faster become available or until an analytic solution is found. We conclude that, for the unscreened gauge glass in $`3D`$, the spin stiffness exponent $`\theta _s=+0.31\pm 0.010`$, which is considerably larger than earlier estimates and indicates that the lower critical dimension $`d_l<3`$. This is consistent with finite $`T`$ MC results gg-hs ; gg-rtyf ; gg-wy  for the gauge glass in $`3D`$ which indicate that $`T_c=𝒪(J)`$. This value of $`T_c`$ is very difficult to reconcile with the suggestion that $`d_l3`$ from previous DWRG studies gg-rtyf ; gg-g ; gg-ks ; gg-mg  as this implies a very small value of $`T_c/J`$.

We have also studied the effects of screened vortex-vortex interactions on the spin domain wall energy using the BT measurement in $`3D`$. Screening of the Coulomb interaction of vortices is implemented by adding a mass term $`\lambda ^2`$ to the denominator of the Green’s function gg-by

$$G(𝐫)=\frac{1}{L^3}\underset{𝐤0}{}\frac{e^{i𝐤𝐫}1}{62\mathrm{c}\mathrm{o}\mathrm{s}k_x2\mathrm{c}\mathrm{o}\mathrm{s}k_y2\mathrm{c}\mathrm{o}\mathrm{s}k_z+\lambda ^2}$$
(45)
The results are also shown in Fig. (7). We average over $`10^3`$ samples for $`L=2,3,4`$ and $`250`$ for $`L=5`$ for several values of the screening length $`\lambda `$. Screening is clearly a relevant perturbation when $`\lambda `$ is finite and $`\theta _s^{BT}<0`$ but our small sizes do not permit an estimate of the value of $`\theta _s^{BT}`$. For large screening lengths, $`\mathrm{\Delta }E_s^{BT}(L)`$ seems to scale the same as for the unscreened case but we expect that $`\mathrm{\Delta }E_s^{BT}(L)`$ will decrease as $`L^{\theta _s}`$ with $`\theta _s<0`$ at length scales which are beyond our computing power for any $`\lambda <\mathrm{}`$. These results are consistent with those of Bokil and Younggg-by  who studied the question of screening using a RT measurement and with Kisker and Rieger kisker  for very strong screening.

### V.3 Varying disorder strength in $`2D`$ and $`3D`$

We have also performed simulations with various strengths of disorder in the CG representation using Eq. (40) for the $`2D`$ case and Eq. (36) for $`3D`$ where the random bond variables $`A_{ij}`$ are independently uniformly distributed in the range $`[\alpha \pi ,\alpha \pi )`$ with $`0\alpha 1`$ so that $`A_{ij}=0`$ and $`|A_{ij}|=\alpha \pi /2`$. Physical realizations of this model are, e.g., an $`XY`$ magnet with random Dzyaloshinski-Moriya interactions xy-magnet  and Josephson junction arrays with positional disorder gk1 ; gk2  where both the effective coupling constant $`J(L)`$ and the effective disorder strength $`A(L)`$ at length scale $`L`$ play a role. Studies in $`2D`$ xy-magnet ; gk1 ; gk2  suggest that weak disorder $`(\alpha 0)`$ does not affect the existence of an ordered phase at intermediate temperature but there is a re-entrant transition to a disordered phase at low temperature. However, recent analytic xyrdm-nsk ; xyrdm-kn ; xyrdm-s ; xyrdm-cf  and numerical gg-ks  studies show there is an ordered phase for $`T<T_c(\alpha )`$ with $`T_c(\alpha )>0`$ for $`0\alpha <\alpha _c`$.

To study how the disorder strength behaves as the length scale $`L`$ varies, one must, if possible, identify the scaled disorder strength $`A(L)`$ with some measurable quantity. An identification has been proposed in a recent numerical study gg-ks  where the effective disorder strength is defined as $`2\pi A(L)|𝚫^0(L)|`$ with $`A(1)=\alpha /4`$ so that one can follow the flows of both $`J(L)`$ and $`A(L)`$ with increasing length scale $`L`$. With this definition, $`0A(L)1/4`$. $`𝚫^0(L)`$ is the global phase twist minimizing the energy of a system of size $`L`$ for a particular realization of disorder. For two phases with energy $`E_{12}=V(\theta _1\theta _2A_{12})`$, the minimum is at $`\theta _1\theta _2=A_{12}`$ which is satisfied by applying a “global” phase twist of $`A_{12}`$. Hence, follows the definition of $`A(L)`$ as a measure of disorder at scale $`L`$.
Since the numerical study gg-ks  used the Hamiltonian of Eq. (1) in the phase representation, we re-investigate this model in the CG representation. The simulations were performed for $`L=2,3,4,5,6,8,10`$ in $`2D`$, $`L=2,3,4,5`$ in $`3D`$ and averaged over at least $`10^3`$ samples. The results are shown in Fig. (8) in $`2D`$ and Fig. (9) in $`3D`$. We are interested in the stable fixed point values at $`L\mathrm{}`$ $`J^{}`$ and $`A^{}`$ as these determine the nature of the phases.

In $`2D`$, weak disorder $`(\alpha <\alpha _c0.37)`$ seems to be marginal and the system seems to iterate toward a glass phase with quasi-long range order characterized by $`(J^{},A^{})`$ where the fixed point values $`J^{}`$, $`A^{}`$ are finite and depend on $`\alpha `$. This is consistent with recent analytic studies xyrdm-nsk ; xyrdm-kn ; xyrdm-s ; xyrdm-cf . On the other hand, systems with strong disorder, $`\alpha >\alpha _c`$, seem to flow to a disordered fixed point $`(J^{},A^{})=(0,1/4)`$ which corresponds to a non-superconducting glass.
In $`3D`$, for $`\alpha <\alpha _c0.57`$, the system flows to a strong coupling limit $`J^{}=\mathrm{}`$. The disorder strength $`A(L)`$ appears to flow to a finite fixed point value $`A^{}`$ which depends on $`\alpha `$. However this is not conclusive from our simulations as only very small values $`L=2,3,4,5,6`$ are used. This can be interpreted as the zero field version of a Bragg glass phase. For large disorder, $`\alpha >\alpha _c`$, $`(J(L),A(L))`$ seem to iterate to their maximum values of $`(\mathrm{},1/4)`$ corresponding to the gauge glass fixed point. It seems that, in the absence of screening, $`\lambda \mathrm{}`$, there are two glassy superconducting phases at $`T=0`$ which, in an applied magnetic field, correspond to a Bragg glass giamarchi  for $`\alpha <\alpha _c`$ and to a vortex glass fisher:89  for $`\alpha >\alpha _c`$.

## VI Discussion and Summary

In this paper, we re-investigate the possibility of an ordered phase at small but finite temperature $`T`$ by a numerical domain wall renormalization group method in a disordered $`XY`$ model in $`2D`$ and $`3D`$ described by the Hamiltonian of Eq. (1) in the Coulomb gas representation. For the $`\pm J`$ $`XY`$ spin glass in $`3D`$, our simulations yield the spin glass stiffness exponent$`\theta _s^{BT}+0.10`$ which suggests its lower critical dimension is $`d_l<3`$. This value of $`\theta _s^{BT}`$ is very different from existing estimates of the chiral glass stiffness exponent in $`3D`$ $`\theta _c+0.47`$ sg-dw-k  and $`\theta _c+0.56\pm 0.18`$ sg-dw-mg . The difference between $`\theta _s`$ and $`\theta _c`$ seems to support the decoupling of two degrees of freedom in $`3D`$. For the gauge glass, we estimate the stiffness exponent $`\theta _s^{BT}=0.36\pm 0.01`$ in $`2D`$ and $`\theta _s^{BT}=+0.31\pm 0.01`$ in $`3D`$, which are considerably larger than all earlier estimates. The latter value is consistent with $`T_c/J𝒪(1)`$ from finite temperature MC studies gg-hs ; gg-rtyf ; gg-wy  and also strongly suggests $`d_l<3`$. The results for the $`XY`$ spin glass in $`3D`$ are consistent with spin glass order at $`T>0`$ which is in contradiction with all other studies sg-mc-kt ; sg-dw-kt ; sg-dw-k  except one sg-dw-mg .
We also studied the effects of varying the disorder strength. In $`2D`$, our simulations imply that weak disorder is marginal xyrdm-nsk ; xyrdm-kn ; xyrdm-s ; xyrdm-cf  and a system with strong disorder flows to a disordered fixed point. There is no sign of a re-entrant transition in our simulations. In $`3D`$, weak disorder has little effect and the system flows to an ordered phase which is the zero field analogue of a Bragg glass giamarchi . For strong disorder, the system seems to flow to a gauge glass fixed point. The disagreement between the stiffness exponent $`\theta _s^{BT}`$ and previous estimates is because these measure $`\theta _s^{RT}`$ whose meaning is less clear. The quantity $`\mathrm{\Delta }E^{RT}(L)`$ seems more likely to suffer from large corrections to scaling as seen in Fig. (7), especially for the small system sizes $`L`$ which are possible to simulate at present. However, we conjecture that both measurements would coincide if much larger values of $`L`$ could be reached. Since our simulations are also limited to very small sizes $`L`$, it is not possible to draw any definite conclusions from them and more studies are needed to settle these problems in random systems more satisfactorily.
One interesting conclusion we can reach concerns the Bragg and vortex glass states in disordered superconductors in an applied magnetic field. Recently, one of us giardina  studied the model of Eq. (1) in the strong screening limit $`\lambda 0`$ in $`3D`$ and found that two phases exist at $`T=0`$ in the presence of an applied external field. The low field, small disorder phase has a well ordered vortex line lattice as a ground state with stiffness exponent $`\theta _s=+1.0`$, whereas the high field, large disorder ground state is a disordered entangled vortex configuration with $`\theta _s1.0`$ kisker ; giardina . We identify the low field state as a Bragg glass and the high field state as a disordered entangled vortex liquid. In this limit, the evidence is strongly in favor of a direct, disorder or field driven transition from a superconducting Bragg glass to a normal non-superconducting phase. This scenario seems to be favored by recent experiments.
In the absence of screening of the vortex-vortex interactions, the picture which results from this work is somewhat different, although the studies here are all done in zero applied field. One may argue that increasing the disorder is equivalent to increasing the field at fixed disorder. At low field or low disorder, the ground state is a Bragg glass with $`\theta +1.0`$, exactly as with screening. Without screening, the main difference is that the high field, large disorder, phase is a true vortex glass with stiffness exponent $`\theta +0.30`$, as proposed by Fisher et al. fisher:89 . We tentatively conclude that a true superconducting vortex glass phase does not exist in $`3D`$ except in the absence of screening ($`\lambda \mathrm{}`$). Our understanding of the experimental consequences for real systems with mesoscopic penetration depths $`\lambda 𝒪(10^3)`$Å is lacking and may be of some interest.
Computations were performed at the Theoretical Physics Computing Facility at Brown University. JMK thanks A. Vallat and B. Grossman for many discussions about spin glasses, best twists etc.