Critical Dynamics of Anisotropic Antiferromagnets in an External Field
CCritical Dynamics of Anisotropic Antiferromagnets in an External Field
Riya Nandi ∗ and Uwe C. T¨auber † Department of Physics and Center for Soft Matter and Biological Physics,Virginia Tech, Blacksburg, VA 24061, USA (Dated: October 26, 2020)We numerically investigate the non-equilibrium critical dynamics in three-dimensional anisotropicantiferromagnets in the presence of an external magnetic field. The phase diagram of this systemexhibits two critical lines that meet at a bicritical point. The non-conserved components of thestaggered magnetization order parameter couple dynamically to the conserved component of themagnetization density along the direction of the external field. Employing a hybrid computationalalgorithm that combines reversible spin precession with relaxational Monte Carlo updates, we studythe aging scaling dynamics for the model C critical line, identifying the critical initial slip, autocor-relation, and aging exponents for both the order parameter and conserved field, thus also verifyingthe dynamic critical exponent. We further probe the model F critical line by investigating thesystem size dependence of the characteristic spin wave frequencies near criticality, and measure thedynamic critical exponents for the order parameter including its aging scaling at the bicritical point.
I. INTRODUCTION
Classical anisotropic antiferromagnets governed byHeisenberg spin exchange in an external magnetic fieldexhibit a rich phase diagram with multiple thermody-namic ground states separated by continuous and dis-continuous transition lines that meet at a multicriti-cal point [1, 2]. This paradigmatic model system de-scribes various magnetic compounds including MnF [3],GdAlO [4], NiCl · O [5], and thus has been ex-tensively investigated theoretically as well as experimen-tally. Early renormalization group [6], Monte Carlo sim-ulation [7], and high-temperature expansion [8] studieshave systematically explored its complex phase diagramand characterized the properties of the different orderedphases and determined the universality classes for its crit-ical transition lines.The anisotropy term along one of the crystal axesbreaks the rotational symmetry of the Heisenberg anti-ferromagnet enforcing an antiparallel spin ordering alongthat axis in the low-temperature, low-field ground state,i.e., the Ising antiferromagnetic phase (AF), c.f. Fig. 1.As the external field strength is tuned up while keepingthe temperature low, the ground state switches, via a dis-continuous (first-order) transition, to a spin-flop phase(SF). At higher values of either temperature or magneticfield, the system becomes paramagnetic (PM). The phasetransitions between the PM and the AF and SF statesare both continuous (second-order). While the associ-ated static critical properties of the system are charac-terized solely by the symmetry and the dimensionality ofthe anisotropic Heisenberg Hamiltonian, the dynamics inthe vicinity of the phase transitions driving it from thedisordered phase to the ordered AF and SF phases, re-spectively, are distinctly different and crucially depend on ∗ [email protected] † [email protected] the microscopic, reversible dynamical couplings betweenthe corresponding order parameters and the conservedmagnetization components.Indeed, the transition from the antiferromagnetic tothe paramagnetic phase is described by the dynamic crit-ical behavior of model C, while the spin-flop to param-agnetic phase transition belongs to the dynamical uni-versality class of model F; here we invoke the classifi-cation introduced by Halperin and Hohenberg in 1977 intheir comprehensive early review of dynamic critical phe-nomena [9]. The presence of a non-ordering field shiftsthe non-universal parameters of the system such as thecritical temperature, but it does not change the natureof the critical point, i.e, the universal scaling exponentscharacterizing the critical power law divergences remainunaltered. Yet an intriguing distinct physical scenario re-sults in the vicinity of the multicritical point , where bothcritical lines as well as the discontinuous phase transi-tion meet. At the multicritical point, all three differentphases compete for the lowest free energy configuration;thus the system’s symmetry at this special point is higherthan in the adjacent parameter space. Hence, the dynam-ical properties in the vicinity of such a point are expectedto be characterized by new exponents which are distinctfrom those of the conjoining critical lines.The nature of the multicritical point for anisotropicantiferromagnets subject to an external magnetic field(oriented along the z direction) has been somewhat con-troversial. This special point in parameter space is char-acterized by two coupled order parameter fields with O ( n ⊥ ) (cid:76) O ( n (cid:107) ) symmetry; it displays long-range orderboth along the magnetic field z axis and in the perpendic-ular xy plane. A series of papers by Folk, Holovatch, andMoser employed the renormalization group to investigatethe static and dynamic critical behavior and the stabil-ity of the associated fixed points of the system [10–13].They predicted the emergence of either bicritical behav-ior associated with a Heisenberg fixed point or tetracrit-ical behavior, in turn associated with a biconical or adecoupled renormalization group fixed point. An analy- a r X i v : . [ c ond - m a t . s t a t - m ec h ] O c t sis to higher order resulted in the tetracritical point witha biconical phase to represent the stable fixed point [14].A recent six-loop perturbative dimensional (cid:15) expansionof the three-dimensional n -vector model has shown thatthe for n = 3, the obtained critical exponents are veryclose to the Heisenberg fixed point values even thoughthat is not the stable fixed point in the asymptotic in-frared limit [15]. Thus, the observed critical behaviorin experiments and simulations may be indistinguishablefrom the isotropic Heisenberg fixed point scaling. Thishas indeed been observed to be true in a series of detailedMonte Carlo simulations analyzing order parameter sus-ceptibilities, the Binder cumulant, and associated prob-ability distributions, which reported that the nature ofthe multicritical point is in fact bicritical with Heisenbergsymmetry [16, 17].While the static critical behavior and stationary crit-ical dynamics of this system has been investigated com-prehensively in the literature, we are not aware of previ-ous computational work addressing the non-equilibriumcritical relaxation of the anisotropic Heisenberg antiferro-magnet in an external magnetic field near either contin-uous phase transition line, nor at the multicritical point.To this end, the system is initially prepared in at a dis-ordered configuration and then quenched precisely to itscritical point such that the dynamics algebraically slowlyevolves towards the asymptotic stationary state. Dur-ing this early non-equilibrium relaxation time window,the system retains the memory of its initial state, andthus manifests broken time translation invariance. Bystudying the ensuing aging scaling behavior of two-timequantities at these early times, one may fully character-ize the dynamics near the distinct critical points and thecorresponding universality classes [18, 19].In this work, we utilize a hybrid simulation methodcombining relaxational Monte Carlo kinetics and re-versible spin precession processes [20] to explore the agingscaling behavior of the AF-to-PM model C critical line,and in the vicinity of the multicritical point. We mea-sure the critical aging scaling, autocorrelation decay, andinitial slip exponents for the model C universality class,and determine the associated dynamic critical exponent.We also investigate the non-equilibrium dynamics of theconserved magnetization component. Furthermore, weutilize Fourier spectral analysis for the spin wave excita-tions in the xy plane in the spin-flop ordered phase toverify the dynamic critical exponent for the SF-to-PMmodel F critical line. Finally, we investigate the criticalorder parameter dynamics at the bicritical point. II. MODEL AND SIMULATION METHOD
The Hamiltonian of a three-dimensional antiferromag-net with anisotropic Heisenberg exchange interactions in
FIG. 1. Schematic H ext - T phase diagram of an anisotropicantiferromagnet in the presence of an external field showingtwo ordered phases separated by a first-order transition line.The antiferromagnetic and spin-flopped ordered phases areseparated from the disordered paramagnetic phase by lines ofcontinuous phase transitions that meet at a bicritical point. an external magnetic field is given by H = J N (cid:88) (cid:104) ij (cid:105) (cid:2) ∆( S xi S xj + S yi S yj ) + S zi S zj (cid:3) − H ext N (cid:88) i =0 S zi , (1)where S ix , S iy , S iz represent the components of the three-dimensional spin vector (cid:126)S i at the i th site of a simplecubic lattice of linear extension L with total spin num-ber N = L . The magnitude of the spins are fixed tounit magnitude, S i x + S i y + S i z = 1. J > z axisbetween nearest-neighbor spin pairs (cid:104) ij (cid:105) ; we set J = +1,i.e., measure temperature in units of J/k B and the exter-nal field in units of J . The uniaxial anisotropy 0 < ∆ < z axis such that the spins would orderanti-parallel along this direction in the absence of an ex-ternal field [16]; we choose ∆ = 0 .
8. The presence of thisanisotropy explicitly breaks the O (3) rotational symme-try of the Hamiltonian in spin space and splits it into twosubspaces of dimensions n (cid:107) = 1 and n ⊥ = 2. Thus itsstatic critical properties are governed by the universalityclass with O (1) (cid:76) O (2) symmetry with associated Fisherexponent η ≈ .
04 [10].Applying an external magnetic field H ext (cid:54) = 0 along the z axis forces the uniaxially aligned spins to flop over intothe xy plane beyond some critical field strength H c ext .Thus, at low temperature T and external field H ext , i.e.in the AF phase, the z component of the staggered mag-netization φ (cid:107) = L (cid:88) i,j,k =0 ( − p S zi,j,k (2)represents the non-conserved order parameter for the sys-tem; here the indices ( i, j, k ) denote the three spatial di-rections, and p = i + j + k ensures that the sum ex-tends over the differences between every alternate spinin the lattice. Upon increasing H ext , the value of φ (cid:107) is diminished, and instead the staggered magnetizationcomponents in the xy plane perpendicular to the appliedfield become appreciable. Thus in the SF phase, an ap-propriate, also non-conserved, order parameter is a two-component vector (cid:126)φ ⊥ = ( φ x , φ y ) with magnitude φ ⊥ = (cid:118)(cid:117)(cid:117)(cid:116)(cid:18) L (cid:88) i,j,k =0 ( − p S xi,j,k (cid:19) + (cid:18) L (cid:88) i,j,k =0 ( − p S yi,j,k (cid:19) . (3)It is important to note that the Hamiltonian (1) con-serves the z component of the total magnetization M z = N (cid:88) i,j,k =0 S zi,j,k (4)under the dynamics, {H , M z } = 0; here the Pois-son bracket constitutes the classical counterpart of thequantum-mechanical commutation relation between thespin angular momentum and the Hamiltonian. The dy-namical mode couplings between conserved magnetiza-tion fluctuations and the order parameter componentsdecisively influence the antiferromagnet’s critical dynam-ics [9, 19, 21, 22]. Indeed, in addition to the irreversible,relaxational terms arising from the static couplings in theHamiltonian, one must account for the reversible kineticscaused by the underlying microscopic dynamics betweenthe order parameter and any conserved modes. At zerotemperature, the microscopic equations of motion obeyedby each spin variable are dS αi ( t ) /dt = {H , S αi ( t ) } , wherethe spin vector components satisfy the standard angularmomentum Poisson brackets { S αi , S βj } = (cid:80) γ (cid:15) αβγ S γi δ ij with the fully antisymmetric unit tensor (cid:15) αβγ = ± n (cid:107) = 1 subspace { M z , φ (cid:107) } = 0; thus there isno reversible coupling term. However, in the n ⊥ = 2subspace, the non-conserved vector order parameter (cid:126)φ ⊥ couples reversibly to the conserved magnetization (4) { M z , φ α } = (cid:15) αβz φ β , (5)where α, β ∈ { x, y } . This non-vanishing mode couplinggives rise to the following deterministic equations of mo-tion of the microscopic spin components at T = 0 [21, 23] d(cid:126)S i ( t ) dt = (cid:126)S i ( t ) × ∂ H ∂ (cid:126)S i ( t ) , (6)which describe precession of the unit spin vector in thelocal effective field.To simulate the dynamics of this system at finitetemperature, one needs to implement relaxation termsas well as the reversible microscopic equations of mo-tion. For convenience, we work with the two angu-lar degrees of freedom ϑ and ϕ that are related tothe unit vector spin components through ( S x , S y , S z ) = (sin ϑ cos ϕ, sin ϑ sin ϕ, cos ϑ ). To respect the underlyingconservation property, the azimuthal angles ϑ are up-dated by means of Kawasaki Monte Carlo kinetics wheretwo randomly picked neighboring spins exchange their ϑ values following the standard Metropolis rules [24]. Incontrast, the polar angles ϕ are evolved using Glauberdynamics where the spin component at the selected lat-tice site is subjected to a finite rotation with againMetropolis updates [25]. These Monte Carlo updatesteps of the spin configurations are alternated with afourth-order Runge-Kutta integration of the equationsof motion (6) [20]. For our simulation, the integrationwas performed in parallel on all spins over discrete timeincrements ∆ t = 0 . /J with each integration step sep-arated by 10 Monte Carlo sweeps over the entire lattice.We determined this combination to be optimal in main-taining the conservation laws within the truncation errorbounds of the numerical integration scheme. III. MODEL C DYNAMICAL SCALING
The dynamical universality class conventionally la-beled as model C describes the pure relaxation dynamicsof a non-conserved n -component critical order parameterfield, coupled to a conserved density [9, 19, 22, 26]. Inthe present study, the low-temperature, low-field groundstate is an Ising antiferromagnet; hence n = 1 with the z -component of the magnetization density constitutingthe conserved scalar field m . A two-loop renormalizationgroup calculation demonstrated that the scalar modelC ( n = 1) is governed by a strong-scaling fixed pointwith both the order parameter relaxation and the con-served density diffusion scaling with the same anoma-lous exponent α/ν [27]. This results in the dynamiccritical exponent z z = z m = 2 + α/ν ≈ . ν ≈ .
72 [28, 29] describes the algebraic divergence ofthe correlation length ξ ∼ | τ | − ν ( τ ∼ T − T c ), and α isthe specific heat critical exponent, C ∼ | τ | − α , which canbe obtained from the hyperscaling relation α = 2 − dν in d dimensions.The non-equilibrium relaxation of model C was in-vestigated by Oerding and Janssen using the dynamicrenormalization group approach [30]. Following a crit-ical quench, the two-time order parameter correlationfunction relating two space-time points at distance r andtimes s < t satisfies the scaling law C ( t, s, r, τ ) = r − ( d − η ) ( t/s ) θ − ˆ C ( r/ξ, t/ξ z ) , (7)where θ is the initial slip exponent representing a new in-dependent universal exponent for purely dissipative sys-tems with non-conserved order parameter [31]. It alsodescribes the power law growth of the order parameterin the early-time universal regime which sets in right af-ter the microscopic time during the non-equilibrium re-laxation process. At the critical temperature T = T c ( τ = 0), the two-time autocorrelation function ( r = 0) . . . .
11 1000 3000 s . C z S M ( t , s ) t / s s = 40 s = 80 s = 120 s = 160 s = 200 C z S M ( t , s ) t − s s = 20012040 (a) . . . . . .
91 0 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 .
07 0 . b ∞ = 0 . ± . b ( L ) / L extrapolated (b) FIG. 2. (a) Aging scaling plots for the two-time spin autocorrelation function C z SM ( t, s ) of the Ising antiferromagnetic orderparameter on a simple cubic lattice of linear system size L = 40 with periodic boundary conditions. The system is quenchedfrom an initially disordered configuration to the critical point at T c = 1 . H c ext = 3 .
0. Double-logarithmic rescaled graphs fordifferent waiting times s collapse with the scaling exponent b = 0 .
6. The inset shows the autocorrelation plots as a functionof t − s (in simulation time steps) for s = 200 , ,
40 (top to bottom), demonstrating broken time translation invariance.Statistical errors are indicated in the graph for the shortest waiting time s = 40. (b) Finite-size extrapolation analysis for theaging exponent b ( L ) plotted vs. 1 /L for six different system sizes L = 16 . . .
40. A linear extrapolation to the infinite systemsize limit L → ∞ yields b ∞ = 0 . ± . assumes the simple-aging scaling form [18], C ( t, s ) ∼ s − b ( t/s ) − λ/z , (8)with the static and dynamic exponents related to the scaling collapse exponent b via b = ( d − η ) /z , (9)and to the autocorrelation exponent λ according to λ = d − η + z (1 − θ ) = z (1 + b − θ ) . (10)For model C with n = 1, a second-order perturbativerenormalization calculation predicts θ ≈ .
27 in threedimensions, if one boldly extrapolates the dimensionalexpansion in (cid:15) = 4 − (cid:15) to (cid:15) = 1 [30].To numerically study the critical relaxation and theaging scaling regime we initialized the system in a disor-dered spin orientation configuration corresponding to avery high temperature, and subsequently performed crit-ical quenches to a point ( T c = 1 . H c ext = 3 .
0) on themodel C critical line. As shown in Fig. 2(a), we obtainan aging scaling window for waiting times s = 120 . . . C z SM ( t, s ) is notsimply a function of the time difference t − s as it wouldbe in the stationary limit, but evolves differently for thedistinct waiting times s in this temporal window. By col-lapsing the data of C z SM ( t, s ) for several waiting times s plotted as a function of the time ratio t/s in accordancewith Eq. (8), one can obtain the collapse exponent b forwhich we find b ≈ . L = 40.The data collapse is noticeably improved for both laterwaiting times s and larger observation times t . This is . . . − − . − . − . − . −
115 20 25 30 35 40 45 h λ/ z z i = . ± . C z S M ( t , s ) t / s L = 20 L = 24 L = 30 L = 36 L = 40 λ / z z L FIG. 3. Double-logarithmic plots of the two-time staggeredmagnetization autocorrelation function vs. t/s for differentlinear system sizes L = 20 , . . . ,
40 (top to bottom) taken atan early waiting time s = 40 exhibit power law decays in thelong-time limit t (cid:29) s . The mean value of the autocorrelationexponent is determined to be (cid:104) λ/z z (cid:105) = 1 . ± . . / z m ∼ . C m ( t ) t L = 16 L = 20 L = 24 L = 30 L = 36 L = 40 (a) − . − . − . − . − . − . − . − .
26 0 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 .
07 0 . / z m | ∞ = . − / z m / L extrapolated (b) FIG. 4. (a) Double-logarithmic plots of the autocorrelation function for the conserved magnetization component M z vs. time t for different linear system sizes L display an intermediate regime governed by algebraic decay. The system was quenched tothe critical point at T c = 1 . H c ext = 3 .
0. (b) Finite-size extrapolation analysis for the decay exponent plotted vs. 1 /L for sixdifferent system sizes L = 16 . . .
40. A linear extrapolation to the infinite system size limit L → ∞ yields 1 /z m = 0 . expected since the simple-aging scaling form (8) is sup-posed to hold only for sufficiently large t (cid:29) s and s .However, in our finite simulation domain, data forlarge times are inevitably hindered by finite-size effects.Thus, in order to better estimate the asymptotic col-lapse exponent, we perform a systematic finite-size ex-trapolation analysis by plotting b vs. 1 /L for systemsizes L = 16 , , , , ,
40, c.f. Fig. 2(b). Linearextrapolation to infinite system size L → ∞ leads to b ∞ = 0 . ± .
02. We then obtain from Eq. (9) thedynamic exponent for the order parameter in model C in d = 3 dimensions, z z = (1 + η ) /b = 2 . ± .
09. Thisresult agrees well within our errors with the theoreticallypredicted value z z ≈ . λ/z z can be extracted from the power law tails of C z SM ( t, s ), apparent in double-logarithmic plots vs. t/s for times t (cid:29) s . Fig. 3 displays the data for five differ-ent system sizes, from which we obtain the mean value (cid:104) λ/z z (cid:105) = 1 . ± .
06. Using Eq. (10), one may obtain theinitial slip exponent θ = 0 . ± .
08. This value shows asimilar trend as the theoretical prediction, but unsurpris-ingly, differs in magnitude by about a factor 2 from thenaive extrapolation of the second-order results of the per-turbative (cid:15) expansion about the mean-field value θ = 0.We further explore the non-equilibrium dynamics ofthe conserved magnetization component M z which is re-versibly coupled to the ordered parameter. This non-critical conserved field undergoes diffusive relaxationwith correlations C m ( (cid:126)q, ω ) ∼ q − ˆ C ( ω/q z ) or equivalently C m ( (cid:126)r, t ) ∼ r − ( d − ˜ C ( r/t z ). The asymptotic long-timescaling form for the temporal magnetization correlationfunction at the critical temperature thus becomes C m ( t ) ∼ t − ( d − /z . (11) In the data depicted in Fig. 4(a), we discern an intermedi-ate power law region in the decay of the spin autocorrela-tion function before the graphs fall off exponentially dueto finite-size effects. As expected, this algebraic regimebecomes more prominent upon increasing the linear sys-tem size L . We again perform a systematic finite-sizeextrapolation to obtain the asymptotic value of the de-cay exponent and find (1 /z m ) ∞ ≈ . z m ≈ . ± .
01. This value however differs from the . .
11 1 10 100 ∼ . ± . .
11 2500 5000 7500 s − . C m ( t , s ) t / s s = 80 s = 120 s = 160 s = 200 C m ( t , s ) t − s s = 20012040 FIG. 5. Aging scaling plots for the two-time autocorrelationfunction C m ( t, s ) of the conserved magnetization for linearsystem size L = 40. Double-logarithmic rescaled plots for dif-ferent waiting times s collapses with exponent b m ≈ − . ≈ . ± .
1. The inset shows the autocor-relation plots as a function of t − s for s = 200 , ,
40 (top tobottom), demonstrating broken time translation invariance. . . − . .
05 2000 6000 10000 L = 30 C xy S M ( ~ q , ω ) ω L = 30 L = 24 L = 20 C xy S M ( ~ r , t ) t (a) τ r e l Lz xy = . ± . (b) FIG. 6. (a) Fourier spectrum of the staggered magnetization correlation function in the xy plane for different linear systemsizes L = 20 , ,
30 for a non-zero wave vector (cid:126)q vs. frequency ω . The inset shows propagating spin waves in the spin-spincorrelation function with spatial separation (cid:126)r as a function of time. (b) Double-logarithmic plot of the relaxation time τ rel vs.linear system sizes L = 16 . . .
32. Fitting the data for large
L >
20 yields the dynamic critical exponent z = 1 . ± . order parameter dynamic exponent z z ≈ . ± .
1, in-dicating that likely the time-scale ratio between the stag-gered magnetization relaxation and magnetization diffu-sion still has not reached its asymptotic fixed point, andthe strong dynamic scaling hypothesis cannot be vali-dated. In this context, we direct the readers towardsprevious work by Koch and Dohm discussing the effectof finite system size on the relaxation and diffusion timescales of model C [32].One can also obtain the aging scaling data from thetwo-time autocorrelation function for the conserved mag-netization, c.f. Fig. 5. We note though that for large wait-ing times we observe another power-law region at latertimes which is distinctly different from the previouslyobtained algebraic decay in the intermediate relaxationregime of the single-time autocorrelation function. It isat these later times that the rescaled plots for differentwaiting times collapse with an exponent b m ∼ − . . ± .
01. Earlier analyses ofconserved spin systems have predicted two regimes withdifferent power laws, with a new length scale governingthe crossover between both algebraic regimes [33, 34]. Ul-timately in the long-time limit, however, the decay of theautocorrelations is determined by only one length scale,independent of the waiting time s . In a similar vein,the non-critical conserved magnetization here displaysthe signature of two distinct scaling regimes. Yet unlikein the conserved spin systems, we observe an early relax-ation regime with a faster power law decay, prominentin the single-time autocorrelation plots in Fig. 4, whichsubsequently crosses over to a slower algebraic decay un-til ultimately finite-size effects dominate. Moreover, thenegative value of the exponent b m suggests the presenceof long-lived metastable states. The precise nature of thecrossover scaling for the conserved magnetization in our system thus remains open for future investigation. IV. MODEL F DYNAMICAL SCALING
The continuous phase transition between the spin flopand the paramagnetic phases is described by the dy-namic universality class model F [9, 22]. Also knownas the “asymmetric planar spin model” [35], this uni-versality class describes the critical dynamics of a two-component vector order parameter coupled reversibly toa conserved scalar density in the presence of an exter-nal Z symmetry-breaking field. The only other knownand prominent physical system described by model Fis the normal- to superfluid transition in He [36]. Inanisotropic antiferromagnets, the non-conserved com-ponents of the planar staggered magnetization φ x and φ y couple reversibly through the non-vanishing Poissonbrackets to the conserved magnetization component M z ,resulting in the precession motion (5) of the spin vectorsaround a local field produced by their exchange interac-tion with their nearest neighbors and the external field.One may view the conserved magnetization compo-nents acting as the infinitesimal rotation generators forthe order parameter components, resulting in propa-gating spin waves in the ordered phase [19, 37]. Thespin wave damping Γ c decreases as the critical temper-ature is approached, with the associated relaxation time τ rel = 1 / Γ c diverging at T c . Hence procuring the agingscaling data by probing the conventional two-time corre-lations turns out not to be a viable approach at the modelF critical line. Moreover, owing to the reversible modecouplings between the conserved magnetization and thenon-conserved order parameter components, the initial . .
11 500 1000 1500 s . C z S M ( t , s ) t / s s = 40 s = 80 s = 120 s = 160 s = 200 C z S M ( t , s ) t s = 2008040 (a) . . . . . . . .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 .
07 0 . b ∞ = 0 . ± . b ( L ) / L extrapolated (b) FIG. 7. (a) Aging scaling plots for two-time spin autocorrelation function C z SM ( t, s ) of the Ising antiferromagnetic orderparameter on a simple cubic lattice of linear system size L = 30 with periodic boundary conditions. The system is quenchedfrom an initially disordered configuration to the bicritical point at T c = 1 . H c ext = 3 . s collapse with the scaling exponent b ≈ .
75. The inset shows the autocorrelation plotsas a function of t − s for s = 200 , ,
40 (top to bottom), demonstrating broken time translation invariance. (b) Finite-sizeextrapolation analysis for the aging exponent b plotted vs. 1 /L for five different linear system sizes L = 16 . . .
36. A linearextrapolation using the data from the four largest systems to the infinite system size limit L → ∞ yields b ∞ = 0 . ± . slip exponent θ and hence the autocorrelation exponent λ are expected to be non-universal in this case; specifically,these exponents should depend on the initial distributionof the magnitudes of the conserved modes [20, 38].However, one can extract the dynamic exponent fromthe temporal evolution of the stationary correlation func-tion in the vicinity of the critical parameters in the or-dered phase. Near the critical temperature, the spin waveoscillations have an exponentially decreasing amplitude ∼ e − Γ c t ∼ e − t/τ rel . In a finite system near T c , the sta-tionary relaxation time diverges with linear system sizewith the dynamic critical exponent z characterizing its critical slowing down : τ rel ∼ L z . We have obtained therelaxation time via measuring the half-peak width Γ c ofthe Fourier transform of the spin-spin correlation func-tion [39], C xy SM ( (cid:126)q, ω ) = (cid:82) C xy SM ( (cid:126)r, t ) e iωt dt , see Fig. 6(a).The asymptotic value of the dynamic exponent for modelF is known exactly from the dynamic renormalizationgroup, z xy = d/ d ≤ L the best fit line gives z xy ≈ .
46, within our error bars in agreement with thetheoretical predition 1 .
5, c.f. Fig. 6(b). As one wouldexpect, with larger system sizes z xy tends towards theasymptotic value. V. BICRITICAL DYNAMICAL SCALING
The two continuous phase transition lines describedby models C and F meet at a bicritical point which is described by a different dynamical universality class.In their field-theoretical analysis, Folk, Holovatch, andMoser found that irrespective of whether the static be-havior of the system is described by the Heisenberg orbiconical renormalization group fixed point, the paral-lel and perpendicular order parameter components scale . . . − . − . − . − . − . − . − . − . −
115 20 25 30 35 40 h λ/ z z i = . ± . C z S M ( t , s ) t / s L = 16 L = 20 L = 24 L = 30 L = 36 λ / z z L FIG. 8. Double-logarithmic plots of the two-time staggeredmagnetization autocorrelation function vs. t/s for differentlinear system sizes L = 16 , . . . ,
36 (top to bottom) taken atan early waiting time s = 40 exhibit power law decays in thelong-time limit t (cid:29) s . The mean value of the autocorrelationexponent is determined to be (cid:104) λ/z z (cid:105) = 1 . ± . similarly in time with dynamic critical exponents z z ∼ z xy ≈ .
003 and z m ≈ .
542 in the asymptotic limit [13].However, strong non-asymptotic effects originating fromthe mode coupling terms in the vicinity of the bicriticalpoint lead to very different crossover dynamical expo-nents which exhibit weak dynamical scaling.Similar to the model C analysis, we obtain a dynamicaging scaling window for waiting times s = 80 . . . L = 30 with aging exponent b ≈ .
75. A subsequent system size extrapolation yieldsthe asymptotic value b ∞ = 0 . ± .
03. Using Eq. (9), onecan then infer the dynamic critical exponent z z = 1 . ± .
15 which is in agreement with the theoretical predictionwithin our error bars. From the mean value of (cid:104) λ/z z (cid:105) =1 . ± .
03 over five different system sizes (c.f. Fig. 8),we also obtain the bicritical initial slip exponent θ =0 . ± .
05 for the staggered magnetization along the z axis. VI. CONCLUSION
We have utilized a hybrid numerical method that in-corporates reversible spin precession dynamics through adeterministic integration scheme with relaxational MonteCarlo kinetics to investigate the both the stationary crit-ical dynamics and the non-equilibrium critical relaxationin three-dimensional anisotropic antiferromagnets in anexternal magnetic field. From the aging scaling data ofthe order parameter spin autocorrelation function at themodel C critical line, we obtained the aging and auto-correlation exponents. A systematic finite-size extrapo-lation analysis allowed for the extraction of the asymp-totic value of the aging collapse exponent b ≈ . z z ≈ . θ ≈ .
14. Additionally, we extract the dynamic exponentfor the conserved magnetization z m ≈ .
62 and observetwo distinct time scales in the decay of its two-time spinautocorrelation function.In the vicinity of the model F critical line, the pres-ence of spin waves hindered the aging scaling analysis.However, from the Fourier transform analysis of the spinwaves we obtained the critical relaxation times that in-crease algebraically with system size, with the dynamiccritical exponent z xy ≈ .
46, which agrees with the the-oretical value z xy = 3 / z z ≈ .
962 is different from the corresponding values atboth the model C and model F critical lines, contrast-ing the nature of the dynamical universality class at thebicritical point.
ACKNOWLEDGMENTS
We would like to thank Michel Pleimling for fruitfuldiscussions and valuable suggestions that aided us in thisproject. Research was sponsored by the U.S. Army Re-search Office and was accomplished under Grant Num-ber W911NF-17-1-0156. The views and conclusions con-tained in this document are those of the authors andshould not be interpreted as representing the official poli-cies, either expressed or implied, of the Army ResearchOffice or the U.S. Government. The U.S. Government isauthorized to reproduce and distribute reprints for Gov-ernment purposes notwithstanding any copyright nota-tion herein. [1] J.M. Kosterlitz, D.R. Nelson, M.E. Fisher, Phys. Rev. B , 412 (1976)[2] K.S. Liu, M.E. Fisher, J. Low Temp. Phys. , 655 (1973)[3] Y. Shapira, S. Foner, Phys. Rev. B , 3083 (1970)[4] H. Rohrer, C. Gerber, Phys. Rev. Lett. , 909 (1977)[5] N.F. Oliveira Jr, A. Paduan Filho, S.R. Salinas, C.C.Becerra, Phys. Rev. B , 6165 (1978)[6] M.E. Fisher, D.R. Nelson, Phys. Rev. Lett. , 1350(1974)[7] D.P. Landau, K. Binder, Phys. Rev. B , 2328 (1978)[8] O.G. Mouritsen, E.K. Hansen, S.J.K. Jensen, Phys. Rev.B , 3256 (1980)[9] P.C. Hohenberg, B.I. Halperin, Rev. Mod. Phys. , 435(1977)[10] R. Folk, Y. Holovatch, G. Moser, Phys. Rev. E ,041124 (2008)[11] R. Folk, Y. Holovatch, G. Moser, Phys. Rev. E ,041125 (2008) [12] R. Folk, Y. Holovatch, G. Moser, Phys. Rev. E ,031109 (2009)[13] R. Folk, Y. Holovatch, G. Moser, Phys. Rev. E ,021143 (2012)[14] P. Calabrese, A. Pelissetto, E. Vicari, Phys. Rev. B ,054505 (2003)[15] L.T. Adzhemyan, E.V. Ivanova, M.V. Kompaniets,A. Kudlis, A.I. Sokolov, Nucl. Phys. B. , 332 (2019)[16] W. Selke, Phys. Rev. E , 042102 (2011)[17] S.H. Tsai, S. Hu, D.P. Landau, J. Phys.: Conf. Ser. ,012005 (2014)[18] M. Henkel, M. Pleimling, Non-equilibrium phase transi-tions , vol. 2 (Springer, 2010)[19] U.C. T¨auber,
Critical dynamics: a field theory ap-proach to equilibrium and non-equilibrium scaling behav-ior (Cambridge University Press, 2014)[20] R. Nandi, U.C. T¨auber, Phys. Rev. B , 064417 (2019) [21] R. Freedman, G.F. Mazenko, Phys. Rev. B , 4967(1976)[22] R. Folk, G. Moser, J. Phys. A: Math. Gen. , R207(2006)[23] S. Ma, G.F. Mazenko, Phys. Rev. Lett. , 1383 (1974)[24] K. Kawasaki, Phys. Rev. , 224 (1966)[25] R.J. Glauber, J. Math. Phys. , 294 (1963)[26] V.K. Akkineni, U.C. T¨auber, Phys. Rev. E , 036113(2004)[27] R. Folk, G. Moser, Phys. Rev. E , 036101 (2004)[28] M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi,E. Vicari, Phys. Rev. B , 144520 (2002)[29] D.P. Landau, Physica A , 41 (1994)[30] K. Oerding, H.K. Janssen, J. Phys. A: Math. Gen. ,3369 (1993)[31] H.K. Janssen, B. Schaub, B. Schmittmann, Z. Phys. B –Cond. Matt. , 539 (1989) [32] W. Koch, V. Dohm, Phys. Rev. E , R1179 (1998)[33] C. Sire, Phys. Rev. Lett. , 130602 (2004)[34] C. Godr`eche, F. Krzaka(cid:32)la, F. Ricci-Tersenghi, J. Stat.Mech. , P04007 (2004)[35] B.I. Halperin, P.C. Hohenberg, E.D. Siggia, Phys. Rev.B , 1299 (1976)[36] B.D. Josephson, Phys. Lett. , 608 (1966)[37] B.I. Halperin, P.C. Hohenberg, Phys. Rev. , 898(1969)[38] K. Oerding, H.K. Janssen, J. Phys. A: Math. Gen. ,5295 (1993)[39] S. Chen, U.C. T¨auber, Phys. Biol. , 025005 (2016)[40] J.D. Gunton, K. Kawasaki, Prog. Theor. Phys.56