{'fraction_non_alphanumeric': 0.03506523765144455, 'fraction_numerical': 0.014655172413793103, 'mean_word_length': 4.406348406600327, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Phase Transition in 3d Heisenberg Spin Glasses with Strong Random Anisotropies, through a Multi-GPU Parallelization

We characterize the phase diagram of anisotropic Heisenberg spin glasses, finding both the spin and the chiral glass transition. We remark the presence of strong finite-size effects on the chiral sector. We find a unique phase transition for the chiral and spin glass sector, in the Universality class of Ising spin glasses. We focus on keeping finite-size effects under control, and we stress that they are important to understand experiments. Thanks to large GPU clusters we have been able to thermalize cubic lattices with up to 64x64x64 spins, over a vast range of temperatures.


I. INTRODUCTION
Spin Glasses (SG) are disordered magnetic alloys. [1][2][3] Their microscopic modelization includes several interactions, such as the RKKY interaction that is invariant over rotations, [4][5][6] and the Dzyaloshinsky-Moriya (DM) interaction that breaks the rotational symmetry. 7,8 Therefore, in theoretical physics SGs are often studied with simplified models that take in account only a few essential characteristics (in particular, quenched disorder and symmetries). 9 The DM interaction, through a spin-orbit coupling with a third spin, causes the interactions between spins in any SG to have a certain degree of random anisotropy. This implies that real SGs are never fully isotropic (this theoretical limit is named Heisenberg SG). In fact, materials are classified according to the degree of anisotropy in their interactions, 10 which turns out to be relevant in their non-equilibrium magnetic response. 11 On one end of the materials' spectrum we find the extremely anisotropic Fe 0.5 Mn 0.5 TiO 3 , which is maybe the best realization of the ideal limit of an Ising SG (Ising SGs correspond to the idealization of uniaxial spins). On the other end, we have very isotropic alloys such as AgMn or CuMn (whose modelization is notoriously difficult, 12 due to the presence of short range spin-density wave ordering [13][14][15] ).
Despite the variety of interactions, already in the early '90s there was general experimental agreement on that SGs undergo a phase transition at sufficiently low temperature. [16][17][18] On the other hand, theoretical work was less advanced, even though one works with extremely simple models. For the Ising SG there were arguments supporting the existence of a phase transition, 19 that were later confirmed numerically. 20,21 In the Heisenberg case, instead, all the attempts carried out during the '80s and '90s failed in finding a phase transition at a finite temperature T SG > 0. [22][23][24][25] In fact, Matsubara et al. argued in 1991 that once a small anisotropic term is added to the Heisenberg Hamiltonian the phase transition becomes visible. 25 This was in agreement with a later domainwall computation. 26 The accepted picture at the time was that the lower critical dimension (i.e. the spatial dimension below which there is no phase transition) lies somewhere between 3D and 4D. 27 However, the story was slightly more complicated. Villain and coworkers made a provocative suggestion hypothesizing that, although maybe there was no spin glass transition, a different order parameter called chirality (or vorticity) could be critical. 28 Chirality is a scalar observable that describes vorticity and alignment between neighboring spins (see below the precise definition in Sect. II B). This idea was elaborated by Kawamura in his 1992 spin-chirality decoupling scenario: in the ideal case of a purely isotropic system the spin and chiral glass order parameters would be decoupled, but the introduction of any small anisotropy would couple them. 29 Kawamura's scenario was apparently consistent with all the observations until 2003, when Lee and Young employed more efficient simulation algorithms and finite-size scaling techniques to show that the spin glass channel is critical also in the fully isotropic model (i.e. the Heisenberg limit). 30 Both order parameters seemed to become positive at the same temperature. Further simulations confirmed the existence of a spin glass phase transition, although uncertainty remains on whether the transition is unique 31,32 or chiralities order at a slightly higher temperature T CG . 33 A parallel issue is measuring the chiral order parameter in experiments. Kawamura proposed in 2003 that the extraordinary Hall resistivity is a simple function of the linear and non-linear chiral-glass (CG) susceptibilities. 34 Experiments based on this proposal observed the chiral transition and measured, for instance, the critical exponent δ. 35 Interestingly enough, the value of δ turned out to be compatible between spin and chiral glass sector. Nonetheless, it was impossible to identify a Universality class despite the critical exponents of these sys-tems had been extensively measured (at least in the SG sector): 10,16,17 the impression was that they change in a continuous way from the Heisenberg to the Ising limit, 36 as we increased the anisotropy.
However, analogy with ferromagnetic materials suggests a different interpretation. Anisotropy would be a relevant parameter in the sense of the Renormalization Group. 37 There should be a new dominant fixed point, and symmetry considerations lead to think it should belong to the Ising-Edwards-Anderson Universality class. Yet, when we add a relevant parameter to the Hamiltonian, there should be some cross-over effects. In other words, one expects that while the correlation length ξ is small, the critical exponents are closer to the Heisenberg-Edwards-Anderson Universality class, and that only for large enough ξ the Universality class reveals its nature.
Notwithstanding, it is very hard, both numerically and experimentally, to prepare a SG with a large correlation length, since one should wait very long times (it has been argued that the waiting time t w required to reach a certain coherence length is proportional to almost its seventh power, see e.g. Refs. 38 and 39). To our knowledge, for this reason, the largest measured correlation lengths are of the order of only one hundred lattice spacings. 11,39 That is a rather small distance to reveal the true Universality class, so it is plausible that experiments will find critical exponents between the two Universality classes.
To further complicate things, in experiments one has to take in account at least two relevant crossovers. The first is the competition, that we just pointed out, between the isotropic and the anisotropic fixed points. It is the one we treat in this paper. The second crossover, that we will not address, is about short versus long range interactions. In fact, the Hamiltonian we treat is short range, but the DM interaction has been shown to be quasi-long-range, in the sense that the interactions are long range, but only until a cut-off distance of the order of some tens of atomic spacings. 40 Aiming to untangle these questions, one of the authors undertook a numerical study of Heisenberg SGs with very weak random anisotropies, 41 but the scenario remained even more foggy, since it was observed that: • The chiral glass critical temperature T CG was significantly higher than T SG , in disagreement with experiments and expectations.
• Apparently, the chiral susceptibility was not divergent at T CG . This is surprising and, apparently, in contrast with experiments. 35 Technically, this lacking divergence appeared as a very large anomalous dimension η CG ∼ 2. 42 • Introducing very weak anisotropies changed dramatically T SG . For example, the T SG found by comparing systems of size L = 6, 12 was about twice its equivalent on the fully isotropic model. This is surprising, since one expects that the critical temperature would change very little from the isotropic case when D is as small as in Ref. 41.
In this paper we will focus on the uniqueness of the phase transition and on the Universality class, proposing that there is a unique transition, belonging to the Ising-Edwards-Anderson (IEA) Universality class. 9 We will also give an interpretation to the results of Ref. 41, showing that the apparent inconsistencies are due to scaling corrections, that we will try to characterize, since we believe them to be fundamental both in the interpretation of numerical simulations and of experiments.
To do all this, we will study numerically the Heisenberg spin glass model with strong random anisotropies, in order to suppress both Finite-Size Effects, and traces of the cross-over from the isotropic limit.
We simulated on the largest lattices to present (up to L = 64), over a wide temperature range. 43 This has been possible thanks to an intense use of graphic accelerators (GPUs) for the computations. We made use of the Tianhe-1A GPU cluster in Tianjin, China, 44 and of the Minotauro GPU cluster in Barcelona. 45 The rest of the paper is organized as follows. In Section II we give an explicit definition of the model, and we introduce the observables we extracted from simulations and analysis. Section III contains details on how we practically conducted the simulations, although much information is relegated to Appendix A, where we also discuss the use of GPUs for spin glasses. On Section IV we recall some Finite-Size scaling concepts we used in our analysis, to find the critical temperatures and exponents (some technical details are given in Section V). Finally, in Section VI we refer the results obtained in this work, and give our conclusions in Section VII.

A. The Model and its symmetries
We study the model introduced by Matsubara et al., 25 which is particularly convenient because of its simplicity. We consider N = L 3 3-dimensional unitary vectors s x = (s 1 x , s 2 x , s 3 x ) on a cubic lattice of linear size L, with periodic boundary conditions. The Hamiltonian is where < · > means the sum goes only over nearest neighbors, and the indices α, β indicate the component of the spins. J xy is the isotropic coupling between sites x and y. D xy is the anisotropy operator: a 3 × 3 symmetric matrix, where the six matrix elements D αβ xy , α ≥ β, are independent random variables.
There is quenched disorder, this means that the time scales of the couplings {J xy , D xy } are infinitely larger than those of our dynamic variables, so we represent them as constant in time random variables, with J xy = D αβ xy = 0, J 2 xy = 1 and (D αβ xy ) 2 = D 2 The overline · · · denotes the averages over the instances of the disorder, while for thermal averages we will use · · · . Each different realization of the couplings {J xy , D xy } is called sample. Independent systems with the same couplings are replicas of the same sample. We use two replicas per sample. Notice that if all the matrix elements D αβ xy are zero we recover the fully isotropic Heisenberg model, with O(3) symmetry. However, if the D αβ xy are non-vanishing, the only remaining symmetry is time-reversal: s x −→ − s x for all the spins in the lattice. Time reversal is an instance of the Z 2 symmetry. This is the symmetry group of the IEA model. 9 Hence, we expect that the Z 2 symmetry will be spontaneously broken in a unique phase transition belonging to the IEA Universality class (see e.g. Ref. 26). Of course, underlying this expectation is the assumption that the anisotropic coupling is a relevant perturbation in the Renormalization Group sense (as it is the case in ferromagnets 37 ). In fact, the infinite-anisotropy limit can be explicitly worked out for a problem with site anisotropy [rather than link anisotropy as in Eq. (1)]: one finds an IEA-like behavior. 46,47 It is widely accepted that the Universality class does not change with the probability distribution of the couplings. 48 We take advantage of this, and choose a bimodal distribution for J xy and D αβ xy , J xy = ±1 and D αβ xy = ±D. These couplings can be stored in a single bit, which is important because we are using GPUs, special hardware devices where memory read/write should be minimized (Appendix A).
We chose the two different values D = 0.5, 1. We want to compare our results with those in Ref. 41, where simulations were done on samples with weak random anisotropies. In that work the D αβ xy did not follow a bimodal distribution, but were uniformly distributed between −0.05 and 0.05. To make proper comparisons we consider the standard deviation of the distribution. For bimodal distributions it is exactly D, in Ref. 41 it is (D 2 ) 1/2 = 1/ √ 1200 ≃ 0.03.

B. The Observables
To define the SG and CG order parameters we use two replicas. The overlap field is q x = s a x · s b x , where a and b are replica indices. Its Fourier Transform at wave vector k isq SG (k) = x q x e ik·x /N . The chirality represents the oriented volume of the parallelepiped we can construct on 3 consecutive spins: where e µ is the unitary vector in the µ direction. The CG overlap is defined similarly to the SG one, as κ x,µ = ζ a x,µ ζ b x,µ . Again a and b indicate the replica. The Fourier Transform of the CG overlap field isq µ CG (k) = x κ x e ik·x /N . We define the wave-vector dependent susceptibilities on the two overlap fields as and from each of them we can compute the correlation length of the related field 37 ξ = 1 2 sin(k min /2) being k min = (2π/L, 0, 0) or permutations. When computing ξ CG , one can choose µ parallel or orthogonal to the wave vector k min . As it was already observed in Ref. 31, there is no apparent difference between the two options, so we averaged over all the values of µ to enhance our statistics.

III. SIMULATION DETAILS AND EQUILIBRATION
We used Monte Carlo dynamics throughout all the work. Previous experience advises to mix several Monte Carlo dynamics. 32,57,58 In fact, our single Monte Carlo step (MCS) consisted of (in successive order): (i) one full lattice sweep with the heat-bath algorithm, (ii) L lattice sweeps of microcanonical overrelaxation algorithm, 59 and (iii) one single Parallel Tempering sweep. 60,61 The combination of the first two, which update one spin at a time, has been shown to be effective in the case of isotropic SGs 75 and other models with frustration. 62,63 Both heatbath and overrelaxation are directly generalized to the anisotropic case. 64 All the simulations were run on NVIDIA Tesla GPUs. Except L = 64, D = 0.5, where we parallelized 45 GPUs, each sample was simulated on a single GPU. The interested reader can find in appendix A details on how they were performed. Table I depicts the relevant simulation parameters. For given L and D, the simulations were all equally long, except for L = 64, D = 0.5, where we extended the simulation of the samples with the longest relaxation times.
To ensure thermalization we made a logarithmic data binning. Each bin had twice the length of the previous, i.e. it contained two times more Monte Carlo Steps (MCS), and had twice the measures. More explicitly, let us call i f the last bin: i f contains the last half of the Monte Carlo time series, i f − 1 the second quarter, i f − 2 the second octave, and so on. This allowed us to create a sequence of values O n (i) , for every observable O, where n indicates the sample, and i identifies the bin, that has length 2 i MCS. A set of samples was considered thermalized if O n (i) − O n (i f ) converged to zero. This test is stricter than merely requesting the convergence of the sequence of O n (i) , because neighboring blocks are statistically correlated, so the fluctuation of their difference is smaller. 65 Physical results were taken only from the last block. We show the simulation parameters for each anisotropy D, and lattice size L. N samples is the number of simulated samples. NT is the number of temperatures that were used in parallel tempering. The temperatures followed a geometric sequence between Tmin and Tmax, and NT was chosen so that the parallel tempering's acceptance was around 15%. N min MCS is the minimum number of MCS for each simulation. The simulation for L = 64, D = 1 was intended only to locate TCG.
Since the L = 64, D = 0.5 samples were the most GPU-consuming, we were more strict with them. To ensure and monitor thermalization, beyond the previous criteria, we measured the integrated autocorrelation time (mixing time) of the random walk in temperatures of each sample. 31 In a thermalized sample, all the replicas stay a significant amount of time at each temperature. We made sure that all the simulations were longer than 10 times this autocorrelation time. The sample-to-sample fluctuations were not extreme, and the autocorrelation times τ spanned between 10000 MCS to 50000 MCS, depending on the sample. Finally, we decided to take measures only over the last 64000 MCS of each simulation.

IV. FINITE-SIZE SCALING
Our simulations were far from the thermodynamic limit, therefore in our analysis we had to take in account finite-size effects. Finite-Size Scaling (FSS) consists in comparing results at different lattice size to characterize the critical point. Specifically, we shall be employing phenomenological-renormalization, also known as the quotients method. 37,66,67 Since FSS applies irrespectively of the considered order parameter, in the current section we will not distinguish between spin and chiral sector. The generic critical temperature will be called T c .
If an observable O diverges at the critical temperature as O ∝ |T − T c | xO , then its thermal average close to the critical point can be expressed like where f O , g O and h O are analytic scaling functions for observable O, while ν is the thermal critical exponent. The exponent ω > 0 is universal, and it expresses the corrections to scaling. The lower dots stand for sub-leading corrections to scaling. Let us name ξ L (T ) the correlation length in a lattice of finite size L, at temperature T . The case O = ξ L (T )/L is of special interest, since ν is the critical exponent for the correlation length. Then, Eq. (6) becomes in this case, up to the leading-order, Therefore, we can identify T c as the temperature where the curves ξ L (T )/L cross for all L for sufficiently large L.
If we let T L,2L be the temperature where ξ L (T )/L crosses ξ 2L (T )/(2L), this regime is reached once the T L,2L has converged. Yet, if ω is small, our lattice sizes may not be large enough, so we will have to take in account the aforementioned corrections to scaling. Including corrections to the order L −2ω , the approach of the crossing temperature T L,2L to the asymptotic value T c can be written as where A and B are non-universal scaling amplitudes.
To compute the critical exponents ν and η we use the quotients' method, taking the quotient of the same observable between different lattice sizes L and 2L. At the temperature T L,2L we get: Again, A xO is a non-universal amplitude, while the dots stand for subleading corrections to scaling. Therefore, if O is the thermal derivative of ξ, we can compute the ν critical exponent through the relation To calculate η we use the susceptibility, as χ ∝ |T −T c | −γ and 2 − η = γ/ν, hence Note that the value of ξ L /L at the crossing tends as well to a Universal quantity:

V. INTERPOLATIONS, EXTRAPOLATIONS AND ERRORS
We have been able to estimate the critical temperature from the crossing of the curves ξ/L at L and 2L, and the exponents ν and η with the method of the quotients, as described in section IV.
To identify the crossing point between the pairs of curves, we used low-order polynomial fits: for each lattice size, we took the four temperatures in the parallel tempering nearest to the crossing point. We fitted these four data points to a linear or quadratic function of the temperature. The obtained results were compatible within one standard deviation (the values reported in this work come from the linear interpolation). In order to calculate ν we needed the derivative of the correlation length at the crossing point. We extracted it by taking the derivative of the polynomial interpolations.
However, there is a difficulty in the calculation of statistical errors: the fits we had to perform came from strongly correlated data (because of the paralleltempering temperature swap). Therefore, to get a proper estimate of the error, we made Jack-knife blocks, fitted separately each block, and calculated the Jack-knife error. 37 The whole mentioned procedure was fluid while T L,2L SG fell in our simulated temperature span. Yet, since T L,2L SG was fairly lower than T L,2L CG , it occurred in four cases that we did not reach low enough temperatures in our simulations to be able to interpolate the crossing, and we had to recur to extrapolations. This happened with D = 1, T 32,64 SG and T 32, 64 CG , and in the lower anisotropy D = 0.5, with T 16,32 SG and T 32, 64 SG . The case of T 32,64 SG (D = 1) and T 16,32 SG (D = 0.5) was not a great issue, because the crossing point was very near to the lowest simulated temperature, so we treated these crossings just like the others.
In the case of T 32,64 SG (D = 0.5), instead, we had to extrapolate at a long distance (see Fig. 1-top, in the next section). Again, we performed the extrapolation through linear in temperature fits. To make the fit of L = 64 more stable, we took in account a progressive number of points (i.e. we fitted to the n lowest temperatures). We increased the number of temperatures, while the crossing temperature was constant. Note that increasing the number of temperatures in the fit results in a smaller statistical error for the crossing-temperature. However, ξ L (T )/L is not a linear function at high T (see Fig. 1). Therefore a tradeoff is needed because, when too high temperatures were included in the fit, the crossing temperature started to change, and we knew that curvature effects were biasing it. Our final extrapolation was obtained from a fit performed on the 10 lowest-temperature points. Unfortunately, this approach was not feasible for the SG susceptibility due to its strongly non-linear behavior. Hence, in the next section we will not give an estimate for η SG (L = 64).
In the case of T 32,64 SG (D = 1), the simulation was not devised to reach that crossing point, and we did not extrapolate data.   (15) TABLE II. For each anisotropy D, and each pair of lattices (L, 2L), we obtain effective size-dependent estimates for TSG, and the universal quantities νSG, ηSG and ξL(TSG)/L. The thermodynamic limit, indicated with L = ∞, is obtained by means of fits to equations (7), (9), (10) and (11) L → ∞ by assuming the IEA Universality class. We took ω IEA = 1.0(1) from Ref. 49, and fitted to Eqs. (9), (10) and (11). In those fits we took in account both the anticorrelation in the data, 70 and the bias arising from the indetermination of the exponent ω IEA . Notice, from table II, that the dependence on L of the data is so weak, that this bias is practically negligible. This situation is different from the one encountered in Ref. 41, where the anisotropy fields were extremely small (D ≃ 0.03). 71 There, the finite-size effects in the SG sector were huge. Overall, the strong consistency of our extrapolations to large L with the IEA exponents shows a posteriori that our assumption was proper.

B. Chiral Glass Transition
In the CG channel (figures 2 and table III) the interpretation is slightly more controversial, since finite-size effects are heavy. For the smaller lattice sizes, T CG is consistently larger than T SG , and ν CG is incompatible with the IEA limit. On the other side, when L is larger, T CG approaches noticeably its SG counterpart, and so does ν CG . We notice that η CG marks the distinction between these two regimes. In fact, when L is small, it is very close to 2. This means that the divergence of χ CG is extremely slow (χ ∼ L 2−η ), 42 revealing we are still far from the asymptotic limit. When L is larger, η CG is consistently smaller, the divergence of χ CG is less suppressed, and we can assume the asymptotic behavior is starting to show up. Consistently with this observation, the value of ξ CG /L at the crossing temperature becomes sizeable [indeed, the second-moment correlation length (4) is well defined only if η < 2, see e.g. Ref. 37].

C. Uniqueness of the transition
Although the SG and CG transitions do not coincide yet with our values of L and D, the critical temperatures, as well as ν, become more and more similar as the linear size of the system increases. Moreover, the decrease of η CG as a function of L has not yet stabilized, so it is likely that the chiral quantities will keep changing with bigger lattice sizes.
As explained in Sect. II A, we expect that the transition should belong to the IEA Universality class. To confirm this expectation, we make the ansatz of a unique transition, of the IEA Universality class, to seek if the two critical temperatures join for L → ∞. Figure 3 (upper half) shows the difference between the critical temperatures as a function of the natural scale for first order corrections to scaling, L −(ωIEA+1/νIEA) [Eq. (7)]. Again, ω IEA and ν IEA are taken from Ref. 49. Not only Fig. 3 (top) reveals a marked increase of the speed of the convergence for L = 64 (to which corresponds the smallest anomalous exponent η CG ), but also, a linear interpolation to infinite volume, taking that point and the previous, extrapolates T SG = T CG within the error. Fig. 4 shows how the SG and CG critical temperatures approach each other with L. Again, T CG gets closer to T SG , and the speed of the approach increases with the lattice size. The points in the intercept represent extrap-  Table II for the extrapolations of T ∞ SG ). The exponents ωIEA and νIEA are taken from Ref. 49. In the upper plot we represent our data, for D = 0.5, 1. The two transitions get closer when we increase L, and the approach appears faster when the lattice size increases. Notice that a linear interpolation between the two largest lattice sizes intercepts the y axis compatibly with a coupling between the two transitions (i.e. TSG = TCG). On the bottom plot we show data from Ref. 41, where much lower anisotropies were considered. Here the scenario is completely different, since the critical temperatures drift apart for large enough L. The horizontal dashed line corresponds to TCG − TSG = 0.
olations to the thermodynamic limit of the T SG . Since the observations are compatible with the ansatz of a unique phase transition, belonging to the IEA universality class, we used the infinite-size limit of T SG to plot the model's phase diagram (Fig. 4, inset). 72

D. Comparing with weak anisotropies
Both plots of Fig. 3 show the same observable, for different anisotropies. The top plot depicts our data, in the case of strong anisotropies D = 0.5, 1. The bottom one represents the case of weak anisotropies (D ≃ 0.03), 71 coming from Ref. 41. The behavior is very different between the two cases. For strong anisotropies, the critical temperatures tend to meet as we increase L. That is qualitatively very different from the weak anisotropy case, where their distance increases. We can ask ourselves where this qualitative difference of behavior comes from.
If we compare same system sizes and different D in table III, we notice that finite-size effects are larger (and η closer to two) the smaller the anisotropy. These differences in the finite-size effects are appreciable with a factor 2 change in the anisotropy (from D = 1 to D = 0.5), so it is reasonable that suppressing the anisotropy by a factor 17 or 35 will increase drastically the finite-size effects.
The most economic explanation is then that there is a non-asymptotic effect that disappears with much larger systems or, as we have seen, with larger anisotropies. In other words there is a L * (D) after which T SG and T CG start joining. For D ≃ 0.03, L * is so large that we observe a growing T CG − T SG , while for D ≥ 0.5 we find L * < 8.
Another from the SG transition alone. It had been observed that a very weak perturbation on the symmetry of the isotropic system implied huge changes in the critical temperature, while one would expect that the transition line is smooth.
To solve this dilemma, we take advantage of having strong evidence for the Universality class of the transition. So, we take the data from Ref. 41, and use once again the exponents ν IEA and ω IEA in Ref. 49 to extrapolate the infinite volume limit with second order corrections to scaling (Eq. 7). The fit is good (χ 2 /d.o.f. = 0.70/1), and, as we show in Fig. 5, its L → ∞ extrapolation for the critical temperature is compatible with T SG (D = 0) within one standard deviation. Thus, taming the finite-size effects was enough to make the scenario consistent, and the issue reduces to the fact that finite-size effects are extremely strong when the anisotropy is smaller.

VII. CONCLUSIONS
We performed a numerical study of the critical behavior of Heisenberg Spin Glasses with strong bimodal anisotropies. Our aim was to clarify the role of scalingcorrections, as well as the crossover effects between the Heisenberg and Ising Universality classes, to be expected when the anisotropic interactions are present. In fact, we show that anisotropic interactions are a relevant perturbation in the Renormalization Group sense: no matter how small the anisotropy, the asymptotic critical exponents are those of the Ising-Edwards-Anderson model. However, a fairly large correlation length maybe needed to reach the asymptotic regime. This observation is relevant for the interpretation of both numerical simulations, 41 and experiments. 10 It is then clear that large system sizes are needed to make progress, something that calls for extraordinary simulation methods. Therefore, we performed single-GPU and multi-GPU simulations to thermalize lattices up to L = 64 at low temperatures. As side benefit, our work provides a proof-of-concept for GPU and multi-GPU massive simulation of spin-glasses with continuous degrees of freedom. This topic is elaborated further in Appendix A.
We performed a finite-size scaling analysis based on phenomenological renormalization. 66,67 We imposed scale-invariance on the second-moment correlation length in units of the system size, ξ L /L. We followed this approach for both the chiral and spin glass order parameters.
Our results for the spin-glass sector were crystal clear: all the indicators of the Universality class were compatible with their counterparts in the Ising-Edwards-Anderson model. On the other hand, in the chiral sector scaling-corrections were annoyingly large, despite they decrease upon increasing the magnitude of the anisotropic interactions.
Regarding the coupling of chiral and spin glass transition, our numerical results seem to indicate that the two phase-transitions take place at the same temperature (i.e. T CG = T SG ). However, it is important to stress that we need our very largest lattices to observe this trend. Nevertheless, what we see is in agreement with both Kawamura's prediction and experiments, where the phase transitions are apparently coupled, and the chiral glass susceptibility is divergent. 35 Moreover, we were able to rationalize the numerical results in Ref. 41 with corrections to scaling, by assuming the Ising-Edwards-Anderson Universality class.
We remark that there are strong analogies between the interpretation of numerical and experimental data. In both cases, there is a relevant length-scale (the correlation length for experiments, the system size for simulations). If that length is large enough, the asymptotic Ising-Edwards-Anderson Universality class should be observed. Otherwise, intermediate results between Heisenberg and Ising are to be expected, and indeed appear. 10 The difficulty in reaching the asymptotic regime lies on time: the time growth of the correlation length is remarkably slow (ξ(t w ) ∼ t 1/z w with z ≈ 7, 38,39 where t w is the waiting time). Indeed, the current experimental record is around ξ ∼ 100 lattice spacings, 11,39 pretty far from the thermodynamic limit. 73 Hence attention should shift to the study of the intermediate cross-over regime.
An intriguing possibility appears: one could envisage an experimental study of the crossover effects as a function of the waiting time. In fact, t w varies some four orders of magnitude in current experiments, 74 which should result in a factor 4 variation of ξ(t w ). The appendix is structured as follows. The specific algorithms that we have used are explained in Sect. A 1 with no reference to their implementation. However, implementation is crucial: our simulations are so demanding that we have used special hardware described in Sect. A 2. This special hardware speeds up the simulations thanks to parallelization, so in Sect. A 3 we give some brief details about it. Finally, we address in Sect. A 4 some issues regarding the generation of pseudorandom numbers.

Simulation algorithms
As explained in Sect. III, we used a blend of several Monte Carlo dynamics. Specifically, our single Monte Carlo step (MCS) consisted of (in successive order): • 1 full lattice sweep with the Heat-Bath algorithm, • L lattice sweeps of microcanonical overrelaxation algorithm, • 1 Parallel Tempering sweep. 60,61 Heat-bath by itself would provide correct (but inefficient) dynamics. It actually mimics the natural evolution followed by real spin glasses (that never reach equilibrium near or below the critical temperature). For this reason we enhance it with two more algorithms. However, heat-bath does play a crucial role, since it is irreducible (i.e. the full configuration space is reachable, at least in principle), at variance with overrelaxation, which keeps the total energy constant, and parallel-tempering, which changes the temperature but not the spin configuration. Crucial to perform the heat-bath and overrelaxation dynamics is a factorization property of the Boltzmann weight for the Hamiltonian (1). The conditional probability-density for spin s x , given the rest of the spins of the lattice is where h x is the local field produced by the lattice nearestneighbors of spin s x (its precise definition is given in footnote 64).
In the heat-bath update, a new orientation for spin s x is drawn from the conditional probability (A1), see Ref. 37 for instance.
The overrelaxation update is deterministic. Given a spin s x and its local field, we change the spin as much as possible while keeping the energy constant: Contrarily to heat-bath, the order in which the spins are updated is important in overrelaxation. Accessing the lattice randomly increases the autocorrelation time in a substantial way. On the other hand, a sequential update generates a microcanonic wave that sweeps the lattice. The resulting change in the configuration space is significantly larger. A similar microcanonic wave is generated with other types of deterministic lattice sweeps. For instance, one could partition the lattice in a checker-board way and first update all spins in the black sublattice, updating the white spins only afterwards. The combination of heat-bath and overrelaxation has been shown to be effective in the case of isotropic spin glasses 75 and other models with frustration. 62,63 However, if one is interested on very low temperatures or large systems, parallel tempering is often useful. For each sample we simulate N T different copies of the system, each of them at one of the temperatures T 1 < T 2 < . . . < T NT . A parallel tempering update consists in proposing, as configuration change, a swap between configurations at neighboring temperatures. The exchange has the Metropolis acceptance. Evidently, the acceptance is higher if the temperatures T i are closer to each other, since the energy of the configurations will be similar. Notice that exchanging configurations is equivalent to exchange temperatures, so the data transfer is reduced to a single number.

Hardware features
The GPUs we used were of the Tesla generation, produced by NVIDIA, with an SIMD architecture (Single Instruction, Multiple Data), 76 optimized for the parallel processing of large amounts of double precision data.
We had access to Tesla M2050 GPUs in the Tianhe-1A supercomputer in Tianjin, 44 China, and Tesla M2090 GPUs on the Minotauro cluster 45 in Barcelona, Spain. Despite the extremely high performances claimed by NVIDIA (e.g. 665 Gflops in double precision in the case of the M2090 GPUs), it is practically impossible to reach that limit, because the major bottleneck does not reside in the computing speed, but in the memory access. Yet GPUs keep being a valid tool to simulate on spin glasses, as they typically allow the same function to be launched concurrently on thousands of threads. This is exactly what we need, since we can update simultaneously different replicas, and also non-neighboring spins within the same replica, because the interactions are only between nearest neighbors (see Sect. A 3).
More details on the specific hardware and codes will be given in Ref. 77.

Parallelization
Our update-schemes support two levels of parallelism. Heat-bath and overrelaxation are parallelized within a single lattice. On the other hand, parallel tempering concerns 2N T independent lattices (two replicas, see Sect. II B, at N T temperatures). Clearly, spins in different lattices can be updated simultaneously (between temperature swaps). For small system sizes, the 2N T lattices can be updated efficiently within a single GPU. Yet, for L = 64 we have found it convenient to speed up by employing N T GPUs, each of them simulating two lattices.

a. Single-GPU
Our parallelization scheme was not very different from the one described extensively in previous works such as Refs. 78 and 79, so we limit ourselves to remark that we used binary couplings in order to be able to store a full coupling in a single byte. Also, due to the fact that the lattice positions were evaluated with bitwise operations, and to our coalesced memory-reading scheme, 77 our program was mostly efficient when the size of the lattice was a power of two, so we favored simulations on those sizes.

b. Multi-GPU
For L = 64 and D = 0.5 the relaxation times were too long to be able to thermalize on a single GPU. Therefore, we prepared a code that mixed CUDA and MPI, in order to be able to concentrate a major computing capability on a single sample. We took advantage of the two levels of parallelization that our update algorithms allow. We used N GPU = N T = 45 GPUs, each updating only two independent lattices with the same couplings, but not necessarily with the same temperature. At the level of the single GPU, the way we swept the lattice with heat-bath and overrelaxation was similar to the single-GPU version. Yet, we had to arrange it in order to get the same thread occupancy as in the single-GPU version. Our choice has been to divide the lattice in rows of 8 spins along the x axis. Non-neighboring rows were updated at the same time. A side advantage of this scheme was that we could use for it the same type of coalesced memory reading that we developed for the single-GPU lattice sweeping.
This arrangement resulted in an extremely small overhead when passing from the single to the multiple-GPU algorithm. We were also favored by other factors. Parallel tempering only requires the exchange with the master of a double precision number. Also, the long correlation times allow to take measurements with low frequency. As a consequence of all this, we obtained a linear scaling of the computing time with the number of GPUs, N T (Fig.  6).

Pseudo-Random Number Generator
Pseudo-random number generators (PRNGs) are a critical issue in the implementation of stochastic algorithms, 80 but even more in cases like ours, where each of the N threads threads had to carry its own PRNG, and we had a large number of them acting in parallel on the same lattice. This became a major problem especially in the simulations with MPI, where a huge number of PRNGs was concentrated on only two lattices. It was crucial to guarantee the statistical independence of the N threads pseudo-random sequences. We consider three different aspects: (a) the PRNG that each thread uses, (b) the initialization of the generators and (c) our tests on the generators.

a. The generator
We resorted to a linear combination of Parisi-Rapuano with congruential generators. 81 With the Parisi-Rapuano sequence, 82 the n th pseudorandom number P n is generated through the following relations: y n = (y n−24 + y n−55 ) mod 2 64 (A3) P n = y n XOR y n−61 , where XOR is the exclusive OR logic operator, and y i are 64-bit unsigned integers. Although some pathologies have been found in the 32-bit Parisi-Rapuano PRNG, 83 it looks like its 64-bit version is solid. 84 On the other side, we used a 64-bit congruential generator, where the n th element of the sequence, C n , was given by: 80,85 C n = (C n−1 ×3202034522624059733+1) mod 2 64 . (A4) Also this generator is not reliable when used alone. 81,86 The final pseudo-random number R n was obtained by summing P n and C n : R n = (P n + C n ) mod 2 64 . (A5)

b. Initializing the generators
We have found that problems arise if special care is not devoted to the initialization of the random numbers. This is particularly important in the case of multiple GPUs where N threads = 32768 threads concurrently update the spins in only two lattices.
We decided to use one seed per node. This seed was used to initialize a 64-bit Congruential PRNG, Eq. (A4). We employed it to initialize the state vector of 24-bit Luescher PRNG. 87 The 24-bits words were obtained from three consecutive congruential calls (we kept the most significant byte from each call). As for the Luescher generator, we employed the full luxury version, which is fireproof but slow. We took the 8 most significant bits from each Luescher call to fill up the state vector of the 64-bit PRNGs in Eq. (A5). We were probably excessively cautious, given the high quality of the full-luxury generator, but initialization takes only a small fraction of the total computing time.

c. Tests
We tested with success our random sequences through the whole battery of tests proposed by Marsaglia in Ref. 88. To be sure the sequences were reliable also with concurrent threads, we also generated N threads sequences and tested them horizontally, i.e. taking first the first number of each sequence, then the second, and so on.
Also, we made simulations with ferromagnetic couplings demanding the energies to be equal, up to the 7 th significant digit, to those obtained with an independent CPU program.
Finally, it has been pointed out that local Schwinger-Dyson relations (see e.g. Ref. 89) can be useful to assess the quality of PRNGs. 83 The relevant identity here is We averaged it over all the sites in the lattice, in order to obtain a more stringent test for the simulations.