# Superconductor-insulator transition of Josephson-junction arrays on a honeycomb lattice in a magnetic field

###### Abstract

We study the superconductor to insulator transition at zero temperature in a Josephson-junction array model on a honeycomb lattice with flux quantum per plaquette. The path integral representation of the model corresponds to a (2+1)-dimensional classical model, which is used to investigate the critical behavior by extensive Monte Carlo simulations on large system sizes. In contrast to the model on a square lattice, the transition is found to be first order for and continuous for but in a different universality class. The correlation-length critical exponent is estimated from finite-size scaling of vortex correlations. The estimated universal conductivity at the transition is approximately four times its value for . The results are compared with experimental observations on ultrathin superconducting films with a triangular lattice of nanoholes in a transverse magnetic field.

###### pacs:

74.81.Fa, 73.43.Nq, 74.40.Kb, 74.25.Uv## I Introduction

Superconductor-insulator transitions in Josephson-junction arrays are particularly interesting as physical realizations of a quantum phase transition fazio (); sondhi (); geerligs (); vdzant (); stroud (); doniach (); gk90 (); schon (); eg04 (); cha1 (). These arrays also provide a simple model for studies of phase coherence in inhomogeneous superconductors. Artificial arrays can be fabricated in different geometries both in and dimensions (2D) vdzant () as a lattice of superconducting grains coupled by the Josephson or proximity effect, with well-controlled parameters. As a theoretical model, a Josephson-junction array is closely related to the Bose-Hubbard model, where Cooper pairs interact on a lattice potential, in the limit of large number of bosons per site sondhi (); cha1 () and it also shows interesting analogies with ultracold atoms on optical lattices nature (); micnas ().

If charging effects due to the small capacitance of the grains are sufficiently large, strong quantum fluctuations of the phase of the superconducting order parameter drive the Josephson-junction array into an insulating phase leading to a superconductor-insulator transition for increasing ratio of charging energy to the Josephson coupling. Under an applied magnetic field, frustration effects lead to phase transitions with distinct universality classes, characterized by different critical exponents and universal conductivities at the transition, depending on the geometry of the array and the value of the frustration parameter . The parameter corresponds the number of flux quantum per plaquette. For a Josephson-junction array model on a square lattice, the superconductor-insulator transition is continuous at low-order rational frustrations, and , and the critical behavior have already been investigated in detail numerically cha1 (); cha2 (); cha3 (); stroud08 (), both in relation to experiments geerligs (); vdzant () and theoretical predictions gk90 (); schon (); fisher (). In particular, for the universal conductivity at the transition gk90 (); cha2 () is approximately two times its value for . For , the transition becomes first order cha1 (). Similar results are not available for Josephson-junction arrays on a honeycomb lattice, despite the general interest on quantum phase transitions. In part, this is due to the lack of experimental data for artificial arrays in this geometry. Recently, however, there has been a growing interest in a related system, in the form of an ultra-thin superconducting film with a triangular lattice of nanoholes valles07 (); valles11 (); valles13 (); kopnov (); valles15 (). A simple model for this system consists of a Josephson-junction array on a honeycomb lattice, with the triangular lattice of nanoholes corresponding to the dual lattice, and it has already been used to investigated the thermal resistive transition in absence of charging effects eg13 (). Nevertheless, quantum phase fluctuations should be taken into account since the system undergoes a superconductor to insulator transition for decreasing film thickness. The transition is very sensitive to low applied magnetic fields and is the analog of the transition in a Josephson-junction array for increasing ratio of the charging energy to the Josephson coupling at different frustration parameters. In the insulating phase near the transition, the magnetoresistance oscillates with the applied magnetic field at low temperatures, displaying minima at integer values of and maxima at , corresponding to quantum phase transitions at different critical points.

In the absence of magnetic field, , or integer values of , the difference between periodic square and honeycomb lattices is irrelevant to the critical behavior. However, for non integer values of , commensurability effects of the vortex lattice on the pinning centers may lead to transitions which are in different universality classes. The case is of particular interest since the topology of the honeycomb lattice leads to geometrical frustration, which can destroy phase coherence even in the absence of charging effects. For the array on a square or triangular lattice, the dual lattice formed by the plaquette centers is bipartite and a commensurate vortex lattice with density and double degeneracy is possible teitel (); shih85 (). For the honeycomb array, however, the dual lattice is triangular and the resulting vortex lattice is incommensurate, with a macroscopically degenerate ground state shih85 (). Thermal fluctuations leads to unusual phase-coherence and vortex-order transitions as a function of temperature, which has been investigated both analytically korsh04 (); korsh12 () and numerically shih85 (); leeteitel (); reid (); eg12 (); eg13 () but are still not well understood. However, the phase-coherence and vortex-order transitions at zero temperature for increasing quantum phase fluctuations have not been investigated in detail.

In this work, we study the superconductor-insulator transition at zero temperature using a self-charging model of Josephson-junction arrays on a honeycomb lattice, with flux quantum per plaquette. The path integral representation of the model corresponds to a (2+1)-dimensional classical model, which is used to investigate the critical behavior by extensive Monte Carlo simulations on large system sizes. A short account of some preliminary calculations for smaller system sizes and for the phase correlations has been reported in a conference proceedings conf (). In contrast to the results for the square lattice gk90 (); cha1 (); cha2 (); cha3 (), we find that for the honeycomb lattice the transition is first order for and continuous for but in a different universality class. Here, the first order transition is determined from the finite-size behavior of the free-energy barrier between coexistent states and the critical behavior for is determined from the finite-size scaling of vortex correlations. We also obtain numerically the universal conductivity at the transition from a scaling analysis of the phase stiffness cha2 (); cha3 (). The results indicate that the phase-coherence and vortex-order transitions occur at the same critical point with a single divergent correlation length. In particular, the universal conductivity at the transition is found to be about four times its value for . The results are compared with experimental measurements on ultrathin superconducting films with a periodic triangular pattern of nanoholesvalles07 (). The magnetoresistance oscillations and activation energy behavior observed in this system in the insulating phase at low temperatures and low magnetic fields, are consistent with the Josephson-junction array model. We argue that the absence of secondary minima at predicted by the model and the activation-energy behavior with film thickness are due to effects of Josephson-coupling disorder.

## Ii Model and Monte Carlo simulation

We consider a Josephson-junction array model on a honeycomb lattice, as illustrated in Fig. 1, described by the Hamiltonian fazio (); doniach ()

(1) |

The first term in Eq. (1) describes quantum fluctuations induced by the charging energy of a non-neutral superconducting grain located at site , where is the electronic charge. The effective capacitance to the ground of each grain is assumed to be spatially uniform. The second term in (1) is the Josephson-junction coupling between nearest-neighbor grains described by phase variables . In the present calculations we assume a spatially uniform Josephson coupling, . The effect of the magnetic field applied in the perpendicular (-direction) appears through the link variables , where is flux quantum and is the vector potential. Since , the gauge-invariant sum around each elementary hexagonal plaquette of the array is constrained by with . The frustration parameter corresponds to the number of flux quantum per hexagonal plaquette. The properties of the model are periodic in with period and have a reflection symmetry about .

To study the quantum phase transition at zero temperature, it is useful to employ the imaginary-time path-integral formulation of the model sondhi (). In this representation, the 2D quantum model of Eq. (1) maps into a (2+1)D classical statistical mechanics problem. The extra dimension corresponds to the imaginary-time direction. Dividing the time axis into slices , the ground state energy corresponds to the reduced free energy of the classical model per time slice. The classical reduced Hamiltonian can be written as

(3) | |||||

where a re-scaling of the time slices has been performed in order to obtain space-time isotropic couplings. In this equation, label the sites in the imaginary-time direction. The ratio , which drives the superconductor to insulator transition for the model of Eq. (1), corresponds to an effective ”temperature” in the 3D classical model of Eq. (3). The energy gap of the insulating phase is related to the phase correlation length in the time direction by . The classical Hamiltonian of Eq. (3) can be viewed as an XY model on a layered honeycomb lattice, where frustration effects exist only in the honeycomb layers.

MC simulations are carried out using the 3D classical Hamiltonian of Eq. (3). The parallel tempering method nemoto () is used in the simulations with periodic boundary conditions. The method is implemented using a number of replicas of the system, each one with a different coupling within a range around the critical point. These replicas are simulated in parallel and the corresponding configurations are exchanged with a probability distribution satisfying detailed balance. For convenience, the honeycomb lattice is defined on a rectangular geometry (Fig. 1) with linear size given by a dimensionless length . In terms of , the linear size in the and directions correspond to and , respectively. We choose a gauge where , on the (tilted) bonds along the rows in the direction numbered by the integer and on the bonds along the direction. For the finite-size scaling analysis, calculations are done for different values of with the linear size in the time direction set equal to linear size in the spatial direction . This choice assumes implicitly that the dynamic critical exponent . More generally, a quantum phase transition shows intrinsic anisotropic scaling, with diverging correlation lengths and in the spatial and time directions sondhi (), respectively. They are related by the dynamic critical exponent as . The scaling behavior discussed in the next Section is indeed consistent with . In particular, for the unfrustrated cased, , the transition is in the universality class of the 3D classical XY model, for which it is known that . For a Josephson-junction array on square lattice cha1 (); cha2 (), it has already been verified that is consistent with the scaling behavior in absence of disorder.

## Iii Numerical results and discussion

We first consider the histogram of the probability distribution for the energy density near the transition. The structure and finite-size dependence of provide information on the nature of the transition, if it is first order or continuous LeeK (). For (Inset of Fig. 2), the restricted free energy defined as shows a double minima structure. This suggests that there are two coexisting phases, separated by a free-energy barrier , where corresponds to one of the minima of and corresponds to the maxima between them. Fig. 2 shows that increases with system sizes indicating that the transition is first order. On the other hand, the double minima feature is absent for even for larger system sizes, which is consistent with a continuous phase transition. This is in sharp contrast with the results for the same model on a square lattice cha1 (); cha2 (), where the transition becomes first order at much smaller frustration . For , the transition is also continuous and in the universality class of the classical XY model in three dimensions, where the dynamic exponent is and the correlation length exponent justin () is . However, for the critical properties are not known in detail. Some results were obtained previously for the phase-coherence transition conf (). Here,, we consider the scaling behavior of the vortex transition and determine the corresponding critical exponents. This allows us to consider the interesting question of the sequence of these transitions, if they occur separately or at the same critical point as a function of the coupling .

To investigate the vortex-order transition, we study the scaling behavior of the finite-size correlation length given by cooper ()

(4) |

Here is the Fourier transform of the vortex correlation function and is the smallest nonzero wave vector. This definition of finite-size correlation length corresponds to a finite-difference approximation to the infinite system correlation length , taking into account the lattice periodicity. If the phase transition is continuous, should satisfy the scaling form cooper ()

(5) |

where is a scaling function. According to this scaling form, curves of as a function of , should cross at the same critical coupling , for different system sizes . Moreover, a scaling plot of sufficiently close to for different should collapse on to the same curve.

For , each honeycomb layer of the 3D classical XY model in Eq. (3) is fully frustrated with a highly degenerate ground state shih85 (). This leads to low-energy states for the layered model, which are macroscopically degenerate. In this case, the choice of an order parameter to describe the correlation function is not obvious since the ordered phase may correspond to an aperiodic or glassy-like pattern. Here we find convenient in the numerical simulations to define the correlation function in terms of an overlap order parameter bhat (). The overlap order parameter of vortex variables is given by , where is the net vorticity at site of the dual lattice (Fig. 1). The vorticity is defined as , where the summation is taken around the corresponding elementary hexagonal plaquette and the gauge-invariant phase difference is restricted to the interval . For the fully frustrated case, , the vortex variables. or chirality variables, are Ising-like variables measuring the direction of the circulating current in an elementary plaquette. The vortex correlation function in the spatial direction is obtained as

(6) |

where is total number of dual sites. An analogous expression is used for the correlation function in the time direction. Eq. (4) is used to obtain the finite-size correlation lengths in the spatial and time directions from the corresponding correlations functions. Note that the correlation length defined in terms of the overlap order parameter may have a different magnitude from the one defined in terms of a single copy. However, they should display the same critical behavior near the transition op ().

In Figs. 3 and 4 we show the finite-size behavior of the vortex correlation length in the time and spatial direction, and , scaled by the system size as function of . The curves for the largest systems cross at the same point, providing evidence of a continuous phase transition 2d (). In the insets of Figs. 3 and 4, we show a scaling plot of the data which verifies the scaling form of Eq. (5) and gives the estimates and from , and and from . The estimates of and are in reasonable agreement, within the errorbars, with the results obtained from the phase correlation length conf (), and , suggesting that the superconductor-insulator transition is accompanied by a vortex disordering transition and described by a single divergent length scale.

In addition to different critical exponents, the superconductor-insulator transition for on a honeycomb lattice, is characterized by a universal conductivity at the critical point, which is significantly different from the square-lattice case. To determine its value we follow the scaling method described by Cha et al. cha2 (); cha3 (). The conductivity is given by the Kubo formula

(7) |

where is the quantum of conductance and is a frequency dependent phase stiffness evaluated at the finite frequency , with an integer. The phase stiffness in the direction is given by

(10) | |||||

where , is the total number of sites in each layer,

(11) |

and is a unit vector between nearest neighbors sites from to . At the transition, vanishes linearly with frequency and assumes a universal value , which can be extracted from its frequency and finite-size dependence cha2 ()

(12) |

The parameter is determined from the best data collapse of the frequency dependent curves for different systems sizes in a plot of versus . The universal conductivity is obtained from the intercept of these curves with the line . From the scaling behavior of the conductivity for shown in Fig. 5 we obtain . In Fig. 6 we show the corresponding behavior for , which gives . As would be expected, the value of the universal conductivity for on a honeycomb lattice agrees with the known result for the square lattice cha3 (), , since the transitions are in the same universality class. However, the corresponding value for is significantly different. Our estimate of the universal conductivity for on a honeycomb lattice is about four times its value for while for the square lattice it is approximately cha2 () two times. It is a clear evidence that, for , these transitions belong to different universality classes.

The critical couplings obtained numerically for different values of frustration are shown in Fig. 7. For comparison, we also show in the inset the critical coupling obtained for the model defined on a square lattice. The lines connecting the data are just guide to the eyes. varies nonmonotonically with a large maximum at integer values of and a secondary maximum at . This is quite different from the behavior for the square lattice, where a secondary maximum occurs at . The sharp variation of the critical coupling with frustration shows up in the behavior of the energy gap with frustration, near the transition. The phase correlation length is obtained from Eq. (4) using the corresponding correlation function in terms of the phase variables , analogous to Eq. (6). Since , for a system with , the energy gap displays oscillations with deep minima at integer values and maxima at , as shown in Fig. 8. A secondary minima at is also observed. In contrast, as shown in the Inset, for the square lattice a secondary minima is expected at . For a continuous superconductor-insulator transition, the gap of the insulating phase vanishes as fisher () , with the exponent . Given the above numerical estimates for the critical exponents and , the energy gap should vanish near the transition with a power-law exponent at zero magnetic field and integer values of , and with at . This sharp dependence on frustration determines the magnetoresistance behavior of the Josephson-junction array in the insulating phase near the transition. In fact, at sufficient low temperatures, the resistance is thermally activated , and the activation energy corresponds to the gap of the insulating phase. Therefore the magnetoresistance should display oscillations with deep minima at integer values of frustration, , and maxima at . Secondary minima at are also expected. In contrast, for the square lattice secondary minima are expected at .

We now compare the expected behavior of the magnetoresistance with experimental observations on ultrathin superconducting films with a triangular lattice of nanoholesvalles07 () at low applied magnetic fields. In the regime where phase fluctuations of the superconducting order parameter are more important than amplitude fluctuations, cha1 (); giroud (); emery () this system can be described by an array of superconducting ”grains” coupled by Josephson junctions in a suitable geometry. The simplest model consists of a Josephson-junction array on a honeycomb lattice eg13 (), with the triangular lattice of nanoholes corresponding to the dual lattice (Fig. 1), which act as vortex pinning centers. The number of flux quantum per unit cell of the triangular nanohole lattice corresponds to the frustration parameter of the honeycomb array. The superconductor to insulator transition is tuned by decreasing the film thickness, which corresponds to increasing the coupling in the Josesphson-junction array model. In fact, measurements of the resistance in the insulating phase near the transition show magnetoresistance oscillations with the applied magnetic field, displaying minima at integer values of and maxima at , as expected. However, although the resistance is thermally activated, the activation energy increases roughly linearly with the deviation from the critical value of the film thickness both for and , which corresponds to an exponent . The secondary minimum at is not observed. This suggests that effects of quenched disorder introduced by the fabrication process may be significant. Disorder in the weak links between superconducting ”grains” can arise from inhomogeneities in the film thickness induced by the substrate valles11 (). In fact, it has been shown numerically eg13 () that for the model of Eq. (1) in absence of charging effects, , disorder in Josephson-junction couplings smooths out the secondary minima at near the thermally induced resistive transition. This kind of disorder should also affect the critical couplings and and critical exponents and for the superconductor-insulator transition.

It turns out to be unfeasible to determine the effects of disorder with the present Monte Carlo method due to the long computer time required for averaging over many realization of disorder and large system sizes. As an alternative, we have performed additional simulations with a driven MC dynamics method eg04 (). In this method, the layered honeycomb model of Eq. (3) is viewed as 3D superconductor and the corresponding ”current-voltage” scaling near the transition is used to determine the critical coupling and critical exponent wengel (). The numerical results egunp () show that indeed increasing disorder in washes out the secondary minima at and lead to an exponent , consistent with the experimental observation.

Geometrical disorder due to spatial irregularities of the system also leads to randomness in the Josephson-coupling but it has a more significant effect for increasing frustration. In fact, weak positional disorder of the grains or weak disorder in the plaquette areas gk86 (), leads to disorder in the magnetic flux per plaquette which increases with the field. Besides changing the universality class of the superconductor-insulator transition, it limits the number of oscillations in the magnetoresistance. This has already been observed in experiments on artificial Josephson-junction arrays with controlled amount of positional disorderbenz88 (), near the thermal resistive transition and in absence of charging effects. Very recently valles15 (), it has also been demonstrated in superconducting films with a pattern of nanoholes with controlled amount of positional disorder near the superconductor-insulator transition.

## Iv Summary and Conclusions

We have studied the superconductor to insulator transition in a self-charging model of Josephson-junction arrays on a honeycomb lattice with flux quantum per plaquette. From Monte Carlo simulations in the path-integral representation of the model, we found that for the superconductor to insulator transition is first order. For , the transition is continuous and the correlation length exponent and universal conductivity at the transition were estimated from finite-size scaling. This is in contrast to the known results for the square lattice gk90 (); cha1 (); cha2 (); cha3 (), for which the transition is also continuous for . Moreover, the critical behavior for is in a different universality class. In particular, the universal conductivity is found to be about four times its value for . As for the square-lattice case vdzant (), it should be interesting to experimentally test this prediction on artificially fabricated Josephson-junction arrays on a honeycomb lattice with controlled parameters. It is interesting to note that the nature of the transitions for and on a honeycomb lattice is suggested from the ground state properties, in the absence of charging effects . For , a vortex pattern with density commensurate with the triangular lattice (dual lattice), with three-fold degeneracy is the ground state shih85 (). If one neglects the coupling to phase variables, the vortex order should then be described by a three-state Potts model, which has a first-order transition in dimensions janke (). Since the quantum phase transition can be described by a (2+1)-dimensional classical model, one thus should expect a first order transition for . On the other hand, for , similar arguments suggest that the vortex order should be described by an antiferromagnetic Ising model on a triangular lattice eg12 (). This Ising model is geometrically frustrated and has a highly degenerate ground state but shows a continuous phase transition in dimensions, as a layered triangular lattice vojta (). In the present case for , however, the transition is continuous but in a different universality class. The behavior of the magnetoresistance oscillations observed experimentally in superconducting films with a triangular lattice of nanoholes valles07 (); valles11 (); kopnov () are qualitatively consistent with the predictions from the model. We argued that the absence of secondary minima at predicted by the model and the approximately linear behavior of the activation energy with film thickness are due to effects of Josephson-coupling disorder. The same system fabricated to be uniformly thick valles13 () is not described by the present model, which assumes superconducting ”grains” and weak links on a length scale of nanohole size and should belong to a different universality. The case is of particular interest since geometrical frustration combined with thermal fluctuations leads to an unusual phase transition as a function of temperature eg13 (); eg12 (); korsh04 (); korsh12 (). It should be of interest to investigate the effects of geometrical disorder in the quantum transition of this system gk86 (); stroud08 (); valles15 ().

###### Acknowledgements.

The author thanks J. M. Valles Jr. for helpful discussions. This work was supported by São Paulo Research Foundation (FAPESP, Grant # 2014/15372-3 ) and computer facilities from CENAPAD-SP. Author contribution statement The sole author had responsibility for all parts of the manuscript.## References

## References

- (1) R. Fazio and H. van der Zant, Phys. Rep. 355, 235 (2001).
- (2) S. L. Sondhi, S. M. Girvin, J.P. Carini, and D. Sahar, Rev. Mod. Phys. 69, 315 (1997).
- (3) L. J. Geerligs, M. Peters, L. E. M. de Groot, A. Verbruggen, and J. E. Mooij, Phys. Rev. Lett. 63, 326 (1989).
- (4) H. S. J. van der Zant, L. J. Geerligs, and J. E. Mooij, Europhys. Lett. 19, 541 (1992).
- (5) W. A. Al-Saidi and D. Stroud, Physica C 402, 216 (2004).
- (6) R. M. Bradley and S. Doniach, Phys. Rev. B 30, 1138 (1984).
- (7) E. Granato and J. M. Kosterlitz, Phys. Rev. Lett. 65, 1267 (1990).
- (8) R. Fazio and G. Schön, Phys. Rev. B 43, 5307 (1991).
- (9) E. Granato, Phys. Rev. B 72, 104521 (2005).
- (10) H. Lee and M.-C. Cha, Phys. Rev. B 65, 172505 (2002).
- (11) M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Nature Physics 10, 588593 (2014).
- (12) A.S. Sajna, T.P. Plak, and R. Micnas, Phys. Rev. A 89, 023631 (2014).
- (13) M.-C. Cha and S. M. Girvin, Phys. Rev. B 49, 9794 (1994).
- (14) M.-C. Cha, M. P. A. Fisher, S. M. Girvin, M. Wallin, and A. P. Young, Phys. Rev. B 44, 6883 (1991).
- (15) K. Kim and D. Stroud, Phys. Rev. B 78, 174517 (2008).
- (16) M. P. A. Fisher, G. Grinstein, and S. M. Girvin, Phys. Rev. Lett. 64 587 (1990).
- (17) E. Granato, Phys. Rev. B 87, 094517 (2013).
- (18) M. D. Stewart Jr., Aijun Yin, J. M. Xu, and J. M. Valles Jr., Science 318, 1273 (2007) .
- (19) S. M. Hollen, H. Q. Nguyen, E. Rudisaile, M. D. Stewart, Jr., J. Shainline, J. M. Xu, and J. M. Valles, Jr., Phys. Rev. B 84, 064528 (2011).
- (20) S. M. Hollen, G. E. Fernandes, J. M. Xu, and J. M. Valles, Phys. Rev. B 87, 054512 (2013).
- (21) G. Kopnov, O. Cohen, M. Ovadia, K. Hong Lee, C.C. Wong, and D. Shahar, Phys. Rev. Lett. 109, 167002 (2012).
- (22) H. Q. Nguyen, S. M. Hollen, J. M. Valles Jr., J. Shainline anJ.M. Xu, Phys. Rev. B 92, 140501 (2015).
- (23) S. Teitel and C. Jayaprakash, Phys. Rev. Lett. 51, 1999 (1983); Phys. Rev. B 27, 598 (1983).
- (24) W.Y. Shih and D. Stroud, Phys. Rev. B 32, 158 (1985); Phys. Rev. B 30, 6774 (1984).
- (25) S.E. Korshunov and B. Douçot, Phys. Rev. Lett. 93, 097003 (2004).
- (26) S.E. Korshunov, Phys. Rev. B 85, 134526 (2012).
- (27) J.R. Lee and S. Teitel, Phys. Rev. Lett. 66, 2100 (1991).
- (28) R. W. Reid, S.K. Bose and B. Mitrović, J. Phys.: Condens. Matter. 9, 7141 (1997).
- (29) E. Granato, Phys. Rev. B 85, 054508 (2012).
- (30) E. Granato, Proceedings of the 27th International Conference on Low Temperature Physics (LT27) , August 2014, Buenos Aires, Argentina, J. Phys.: Conference Series 568, 022017 (2014).
- (31) J. Lee and J.M. Kosterlitz, Phys. Rev. Lett. 65, 137 (1990).
- (32) K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).
- (33) R. Guida and J. Zinn-Justin, J. Phys. A 31, 8103 (1998).
- (34) H. G. Ballesteros, A. Cruz, L. A. Fenández, V. Martin-Mayor, J. Pech, J. J. Ruiz-Lorenzo, A. Tarancón, P. Téllez, C. L. Ullod, and C. Ungil, Phys. Rev. B 62, 14237 (2000)
- (35) R.N. Bhatt and A.P. Young, Phys. Rev. B 37, 5606 (1988).
- (36) We have performed additional calculations using as an order parameter. The wavevector was determined numerically from the Fourier transform of for low energy states, corresponding to a maximum in . Indeed, the results obtained for the critical coupling and the correlation length exponent agree with those of obtained from the overlap order parameter within the estimated errorbars.
- (37) Similar calculations of the finite-size correlation length for different system sizes as a function of temperature for the model of Eq. (1) without charging effects, , shows no crossing point or merging at finite temperature, suggesting that the critical temperature of the thermal transition vanishes eg12 ().
- (38) M. Giroud, O. Buisson, Y.Y. Wang and B. Pannetier, J. Low Temp. Phys. 87, 683 (1992).
- (39) V.J. Emery and S.A. Kivelson, Nature (London) 374, 434 (1995).
- (40) E. Granato, Phys. Rev. B 69, 144203 (2004).
- (41) C. Wengel and A.P. Young, Phys. Rev. B 56, 5918 (1997).
- (42) E. Granato, unpublished.
- (43) E. Granato and J.M. Kosterlitz, Phys. Rev. B 33, 6533 (1986); Phys. Rev. Lett. 62, 823 (1989).
- (44) M.G. Forrester, Hu Jong Lee, M. Tinkham, and C.J. Lobb, Phys. Rev. B 37, 5966 (1988); S.P. Benz, M.G. Forrester, M. Tinkham, and C.J. Lobb, Phys. Rev. B 38, 2869 (1988).
- (45) S. J. Knak Jensen and O. G. Mouritsen, Phys. Rev. Lett. 43, 1736 (1979).
- (46) S. V. Isakov and R. Moessner, Phys. Rev. B 68, 104409 (2003).