Conjugate field and fluctuation-dissipation relation for the
dynamic phase transition in the two-dimensional kinetic Ising
model
D.T. Robb,1, 2, ∗ P.A. Rikvold,1, 3, 4 A. Berger,5 and M.A. Novotny6
1School of Computational Science, Florida State
University, Tallahassee, Florida 32306, USA
2Department of Physics, Clarkson University, Potsdam, New York 13699, USA
3Center for Materials Research and Technology and Department of Physics,
Florida State University, Tallahassee, Florida 32306-4350, USA
4National High Magnetic Field Laboratory, Tallahassee, Florida 32310, USA
5San Jose Research Center, Hitachi Global Storage
Technologies, San Jose, California 95120, USA
6Department of Physics and Astronomy and HPC2 Center for Computational Sciences,
Mississippi State University, Mississippi State, Mississippi 39762, USA
(Dated: October 24, 2018)
http://arxiv.org/abs/0704.1123v2
Abstract
The two-dimensional kinetic Ising model, when exposed to an oscillating applied magnetic field,
has been shown to exhibit a nonequilibrium, second-order dynamic phase transition (DPT), whose
order parameter Q is the period-averaged magnetization. It has been established that this DPT falls
in the same universality class as the equilibrium phase transition in the two-dimensional Ising model
in zero applied field. Here we study for the first time the scaling of the dynamic order parameter
with respect to a nonzero, period-averaged, magnetic ‘bias’ field, Hb, for a DPT produced by a
square-wave applied field. We find evidence that the scaling exponent, δd, of Hb at the critical
period of the DPT is equal to the exponent for the critical isotherm, δe, in the equilibrium Ising
model. This implies that Hb is a significant component of the field conjugate to Q. A finite-
size scaling analysis of the dynamic order parameter above the critical period provides further
support for this result. We also demonstrate numerically that, for a range of periods and values
of Hb in the critical region, a fluctuation-dissipation relation (FDR), with an effective temperature
Teff (T, P,H0) depending on the period, and possibly the temperature and field amplitude, holds
for the variables Q and Hb. This FDR justifies the use of the scaled variance of Q as a proxy for
the nonequilibrium susceptibility, ∂〈Q〉/∂Hb, in the critical region.
PACS numbers: 05.70.Ln, 64.60.Ht, 89.75.Da, 75.70.Cn
∗Corresponding author: drobb@clarkson.edu
mailto:drobb@clarkson.edu
I. INTRODUCTION
The dynamic phase transition (DPT) in a ferromagnetic system below its critical temper-
ature was first observed in numerical solutions of a mean-field model exposed to an oscillating
magnetic field [1, 2]. It was then studied further, both in mean-field models [3, 4, 5, 6] and
in kinetic Monte Carlo (KMC) simulations [5, 6, 7, 8, 9, 10]. A review of this early work can
be found in Ref. [11]. More recently, the study of the DPT has expanded to include varying
(and often more physical) model geometries. These include mean-field studies of domain-
wall motion in an anisotropic XY model in one dimension [12, 13, 14], KMC simulations
of a three-dimensional Ising system [15], and KMC simulations of a uniaxially anisotropic
Heisenberg system in an off-axial field [16], an elliptically polarized applied field [17], and
with the effect of a thin-film surface energy [18, 19, 20]. The phenomenon has also been
observed in simulations of CO oxidation under oscillating CO pressure [21, 22]. Further
simulation studies of the DPT in the two-dimensional kinetic Ising model have appeared
[23, 24, 25, 26, 27], as well as analytical studies of the DPT [28, 29, 30, 31].
Here, we concentrate on the DPT in the two-dimensional kinetic Ising model. It was ob-
served in simulations of this model that there exists a singularity at a critical period of the
applied oscillating field [9, 10], and that the critical exponents β and γ (and, with less accu-
racy, ν) are consistent with the universality class of the equilibrium two-dimensional Ising
transition in zero field [23] (β = 1/8, γ = 7/4, and ν = 1). In those studies, the techniques
of finite-size scaling were extended to the study of the dynamic order parameter (Q, defined
in Sec. II) in the non-equilibrium steady state. This provided evidence for a diverging cor-
relation length at a critical value of the period. In particular, because the field conjugate
to Q and a fluctuation-dissipation relation were not known, a susceptibility could not be
measured directly, and the scaled variance X
L = L
2 (〈Q2〉 − 〈Q〉2), where L is the linear
system size, was used as a proxy. An analytical argument, based on the correspondence
of the two-dimensional kinetic Ising model and the continuous, two-dimensional Ginzburg-
Landau model at the equilibrium critical point, provided an effective Hamiltonian for the
non-equilibrium system and confirmed that the DPT is in the Ising universality class [28].
These findings are consistent with earlier symmetry arguments that any continous phase
transition in a stochastic cellular automaton that preserves the Ising up-down symmetry
should be in the equilibrium Ising universality class [32, 33].
Recently, experiments were performed on a [Co(0.4nm)/Pt(0.7nm)]3 multilayer film with
strong uniaxial anisotropy [34], whose equilibrium behavior is known to be Ising-like [35, 36].
The film was exposed to an oscillating (sawtooth) applied field with varying period, in the
presence of constant ‘bias’ magnetic fields Hb of varying strength and sign. (The bias field
is defined explicitly in Sec. II.) The behaviors of the dynamic order parameter and its
variance, as functions of the applied field period and the bias field, provided strong evidence
for the presence of the DPT in this experimental system. The observed behavior of the order
parameter with respect to the bias field supported previous conjectures that the conjugate
field could include the period-averaged magnetic field as an important component, and
stimulated the numerical investigations in this paper.
This paper is organized as follows. In Sec. II, we describe the two-dimensional kinetic
Ising model and our computational methods. In Sec. III, we verify directly the scaling
of the dynamic order parameter with respect to the period-averaged magnetic field at the
critical period, with scaling exponent δd ≈ δe = 15, in agreement with the equilibrium Ising
transition. In Secs. IV and V, we derive the expected asymptotic scaling functions in a finite-
size scaling analysis of the dynamic phase transition with non-zero period-averaged bias field,
and then compare the expected scaling of the dynamic order parameter to our numerical
results. In Sec. VI, we present numerical data assessing the applicability of a fluctuation-
dissipation relation (FDR) to this far-from-equilibrium system. We then in Sec. VII compare
the expected scaling of the susceptibility (of the dynamic order parameter) to our numerical
data, using the results of Sec. VI to reconcile our findings with previous results on the scaling
of the fluctuations of the dynamic order parameter. Finally, we present a summary of our
results in Sec. VIII.
II. COMPUTATIONAL MODEL
In order to facilitate comparison with previous results, we employ the same model and
computational method as in Ref. [23]. Specifically, we perform kinetic Monte Carlo (KMC)
simulations of a two-dimensional periodic square lattice of Ising spins Si, which can take
only the values Si = ±1. The Hamiltonian of the model is
H = −J
〈i,j〉
SiSj −H(t)
Si, (1)
where J > 0 is the ferromagnetic exchange interaction,
〈i,j〉 runs over all nearest-neighbor
pairs,
i runs over all L
2 lattice sites, and H(t) is an oscillating, spatially uniform applied
magnetic field. The form of H(t) is here taken as a square wave with amplitude H0 = 0.3J
and period P , measured in Monte Carlo steps per spin (MCSS). The square-wave form not
only allows for more efficient KMC simulation, but also reduces the critical period and the
finite-size effects for the DPT [23]. Other symmetric field shapes, such as sinusoidal [9, 10]
and sawtooth [34], yield essentially the same results, but with a larger critical period and
with stronger finite-size effects. The Glauber single-spin-flip MC algorithm with updates at
randomly chosen sites is used, in which each attempted spin flip is accepted with probability
W (Si → −Si) =
1 + exp (∆E/T )
, (2)
where ∆E is the energy change that would result from acceptance of the spin flip, and T
is the absolute temperature in energy units (i.e., with Boltzmann’s constant set to unity).
All simulations were performed at T = 0.8Tc, where Tc = 2.269J is the equilibrium critical
temperature of the square-lattice Ising ferromagnet in zero applied field [37].
The system responds to the oscillating field via the time-dependent magnetization per
site,
m(t) =
Si(t). (3)
The dynamic order parameter is defined as the average of m(t) over a given field cycle i [1] :
(i−1)P
m(t)dt. (4)
We define the bias field, so named because it measures the shift (or ‘bias’) of the periodic
field toward either negative or positive field values, as the period-averaged magnetic field,
H(t)dt. (5)
This definition applies generally to any periodic magnetic field H(t). In this paper, the
applied field consists of a square wave with period P superposed with a constant magnetic
field. Applying (5), since the period-average of the square-wave field is zero, the bias field
Hb in this case is simply equal to the superposed constant magnetic field.
III. SCALING WITH RESPECT TO THE BIAS FIELD
In the two-dimensional equilibrium Ising model, the critical isotherm is given (in the
thermodynamic limit, i.e., as L → ∞) as
m (T = Tc, H → 0) ∝ H
1/δe , (6)
where the critical exponent δe = 15 [38]. For finite systems, this relationship breaks down
when the infinite-system correlation length, ξ∞ (T = Tc, H), which diverges as H → 0,
becomes comparable to the linear system size L. The relationship also naturally breaks
down at larger fields away from the critical region. Therefore, a plot of m vs H for a given
system size L will follow the power law (6) for a range of H near the critical value H = 0,
with this range extending to smaller H as L is increased [39].
We can determine directly whether the non-equilibrium system exhibits a similar rela-
tionship,
〈Q〉 (P = Pc, Hb → 0) ∝ H
b (7)
in an analogous way. In Fig. 1, we plot 〈Q〉 vs Hb at the critical value of the period, P = Pc.
In previous work, the reversal time for the magnetization, following instantaneous reversal
of the uniform magnetic field H at (H = 0.3J, T = 0.8Tc), was found as τ = 74.5977 MCSS
[9, 10]. The critical scaled half-period for the square waveform was determined to be Θc =
Pc/ (2τ) = 0.918± 0.005 [23]. This yields Pc = 136.96± 0.75 MCSS, and in our simulations
and analysis in this paper we use Pc = 136.96 MCSS.
A power-law dependence is indeed seen to hold in Fig. 1, within a range which extends
to lower values of Hb as L is increased. We fit the L = 256 data between the points
labeled A and B in Fig. 1, finding a statistically significant fit with power-law exponent
δd = 14.85 ± 0.18. As including points with Hb ≥ 0.01J was found to greatly reduce the
statistical significance of the fit, the value Hb = 0.01J serves as a boundary of the scaling
region at P = Pc. This result is consistent with an exponent δd = δe = 15, suggesting that
the bias field Hb, for these parameters and the square waveform, is the dominant component
of a conjugate field which exhibits the same scaling exponent in the DPT as does the applied
magnetic field in the equilibrium Ising transition.
IV. FINITE-SIZE SCALING ANALYSIS WITH BIAS FIELD
To provide more complete evidence that Hb is the dominant component of the field
conjugate to 〈Q〉, in the next several sections we will demonstrate data collapse onto a two-
parameter finite-size scaling function for the system-size dependent quantity 〈Q〉L at points
(P ≥ Pc, Hb > 0), for lattice sizes L = 90, 128, 180, and (in several cases) L = 256, using
the critical exponents for the equilibrium Ising system. In this section, we briefly review the
theory of finite-size scaling as it applies to this system. We then determine the expected
asymptotic forms of the scaling functions, which are compared in later sections of the paper
to our computational data.
The theory of finite-size scaling [40, 41] states that near a continuous phase transition,
the singular part of the free-energy density for a d-dimensional system of linear size L can
be written as
fL ≈ L
|ǫ|L1/ν , HLβδ/ν
, (8)
where ǫ = (T − Tc)/Tc, H (in units of kBT ) is the field conjugate to the order parameter, ν
is the critical exponent for the correlation length, β is the exponent for the order parameter,
δ is the exponent for the critical isotherm, and Y± are scaling functions above (+) and below
(−) the critical point. This yields for the order parameter at finite L
= L−β/νF0±
|ǫ|L1/ν , HLβδ/ν
, (9)
where the exponent for L in the prefactor is obtained by using the hyperscaling relation
dν = 2 − α and the exponent equality α = 2 − β(δ + 1). Further differentiation yields the
susceptibility,
= Lγ/νG0±
|ǫ|L1/ν , HLβδ/ν
, (10)
where the exponent for L in the prefactor is obtained by using the exponent equality γ =
β(δ − 1).
It has previously been shown analytically that the DPT for a sinusoidal applied field,
which is symmetric under H(t) → −H(t+P/2) and so which can safely be assumed to have
Hc = 0, has an effective Ginzburg-Landau free-energy density in the same universality class
as the equilibrium Ising model [28]. It therefore appears reasonable to write corresponding
scaling functions for the dynamic order parameter 〈Q〉 and its associated susceptibility χ̂,
〈Q〉L = L
−β/νF±
|θ|L1/ν , (Hc/J)L
, (11)
χ̂L = L
γ/νG±
|θ|L1/ν , (Hc/J)L
, (12)
where θ = (P − Pc)/Pc, and Hc is the (as yet unknown) field conjugate to 〈Q〉. In this
paper we express Hc (and Hb) in units of the exchange constant, J , so that the second
scaling parameter in Eqs. 11 and 12 is dimensionless. The specific form, Hc/J , with which
Hc is assumed to enter the second scaling parameter needs more theoretical investigation,
and could conceivably change as the theory of the DPT is further developed. However, this
should not affect our conclusions [42]. Computational results for sinusoidal and square-wave
fields, which both are symmetric underH(t) → −H(t+P/2) and so presumably haveHc = 0,
have previously confirmed the scaling behavior with respect to θ. The exponent values were
determined as γ/ν = 1.74± 0.05, β/ν = 0.126± 0.005, and ν = 0.95± 0.15 [23], consistent
with the exact values for the two-dimensional equilibrium Ising model, γ = 7/4 = 1.75,
β = 1/8 = 0.125, and ν = 1.
We now determine the expected asymptotic forms of the scaling functions F+(y1, y2) and
G+(y1, y2), where we emphasize that the + subscript indicates that the scaling functions refer
to the range P ≥ Pc, and where the scaling parameters are y1 ≡ θL
1/ν and y2 ≡ (Hc/J)L
βδ/ν .
y1 ≫ y2. We expect χ̂L ∼ θ
−γ = Lγ/νy
1 (independent of y2) and 〈Q〉L = χ̂LHc ∼ θ
−γHc ∼
L−β/νy
1 y2, where γ = β(δ − 1) was used to obtain the exponent for L in 〈Q〉L.
y1 ≪ y2. We expect that 〈Q〉L ∼ H
c ∼ L
−β/νy
2 and χ̂L = ∂〈Q〉L/∂Hc ∼ L
(1−δ)/δ
(both independent of y1), where γ = β(δ − 1) was used to obtain the exponent for L in χ̂L.
Thus, the asymptotic forms of the scaling functions are expected to be
F+(y1, y2) ≡ L
β/ν〈Q〉L ∼
1 y2 for y1 ≫ y2
2 for y1 ≪ y2
G+(y1, y2) ≡ L
−γ/ν χ̂L ∼
1 for y1 ≫ y2
(1−δ)/δ
2 for y1 ≪ y2
. (14)
V. COMPARISON OF FIRST SCALING FUNCTION TO COMPUTATIONAL
RESULTS
In Fig. 2(a), using the equilibrium values βe and νe in calculating F+(y1, y2) ≡ L
β/ν〈Q〉L,
we present a plot of the scaling function F+ vs y1 for different values of y2, for lattice
sizes L = 90, 128, and 180. Here and for the remainder of the paper, exponents with the
subscripts ‘d’ and ‘e’ refer to the behavior of the nonequilibrium system (with a dynamic
phase transition) and the equilibrium system, respectively. The scaling function exhibits
a power-law dependence in the regime y1 ≫ y2, which is consistent with Eq. (13). At
progressively larger values of the constant y2, the power-law scaling can be seen to begin at
increasing values of y1, as would be expected. A best-fit line to the final five points of the
L = 180 data at y2 = 3.39 yields an estimate of the scaling exponent −γd = −1.76 ± 0.07
in Eq. (13). This is consistent with the previous results for Hc = 0 cited above [23], and
it supports the hypothesis that γd = γe = 7/4 = 1.75. In Fig. 2(b), we present just the
data for y2 = 3.39, including additional data points at y1 = 280 and 477. The data deviate
from the power-law behavior for L = 90 at y1 > 149, for L = 128 at y1 > 280, and for
L = 180 at y1 > 477. This locates the boundary of the scaling regime (for y2 = 3.39) at
θ = y1/L
1/ν ≈ 2.65 .
In Fig. 3, again using βe and νe in calculating F+, we plot the scaling function F+ vs
y2 at different values of y1, in order to examine the scaling behavior for y1 ≫ y2. For the
constant values y1 = 43.4, 69.7, and 149, power-law scaling can be observed in the regime
y1 ≫ y2. A best-fit line to the y1 = 149 data for the five points from y2 = 3.39 to 84.6 yields
a scaling exponent of 1.01± 0.01, which is consistent with the value of 1 expected from Eq.
(13).
In order to investigate the scaling of F+ in the asymptotic limit y1 ≪ y2, we plot in Fig. 4
the scaling function F+(y1, y2) vs y2 at the critical period P = Pc (i.e., y1 = 0), at lattice
sizes L = 128, 180, and 256. In the range 20 < y2 < 50, power-law scaling is observed for all
three lattice sizes. For y2 > 50, the data deviate from power-law scaling, with the smallest
lattice size deviating first, as expected in a finite-size scaling plot. A fit of the L = 256 data
in the range from y2 = 8.46 to 84.6 produces a scaling exponent 0.0673± 0.0008. Since the
constant factors Lβδ/ν and Lβ/ν do not affect the fit of the scaling exponent, this is the same
exponent found in the fit of 〈Q〉 vs Hb in Fig. 1. The reciprocal of this scaling exponent
is thus δd = 14.85 ± 0.18, which is consistent with the exponent of the critical isotherm,
δe = 15, in the equilibrium Ising model.
The comparison of the second scaling function, G+(y1, y2), to numerical data is more
clearly presented after the relationship of the susceptibility χ̂L and the scaled variance X
has been examined. Therefore, we present in the next section numerical results on the extent
of applicability of an FDR between χ̂L and X
L , before turning in Sec. VII to the second
scaling function.
VI. APPLICABILITY OF A FLUCTUATION-DISSIPATION RELATION
FDRs, such as the Einstein relation, Green-Kubo relations, etc., hold a central place in
equilibrium statistical mechanics. This is essentially a consequence of detailed balance and
the role of the partition function as a moment-generating function, and thus such relations
cannot be readily extended to nonequilibrium steady states. However, it has recently been
shown that certain FDRs can be extended to far-from equilibrium steady states by use of an
effective temperature [43, 44]. Here we will therefore consider whether the nonequilibrium
susceptibility and the scaled variance of the dynamic order parameter can be related as
χ̂L ≡
∂〈Q〉L
L2 (〈Q2〉L − 〈Q〉
, (15)
with an effective temperature Teff , in a way analogous to the equilibrium FDR,
∂〈m〉L
L2 (〈m2〉L − 〈m〉
, (16)
in which T is the temperature. As mentioned in Secs. I and IV, this conjecture motivated
the use in previous work of the scaled variance X
L as a proxy for χ̂L in investigating the
scaling behavior of the nonequilibrium system near its critical period.
To test the extent to which Eq. (15) holds, we computed values of χ̂L and X
L for a range
of periods from P = 140 to 250 MCSS and a range of bias fields from Hb = 0 to an upper
limit between 0.005J and 0.2J . (The bias field necessary to ‘saturate’ the nonequilibrium
system, i.e., to produce values of χ̂L and X
L near zero, increases as the period is increased.)
The computations were perfomed at L = 180. The quantity χ̂L was computed directly as a
numerical derivative:
χ̂L(P,Hb) ≈ (〈Q〉(P,Hb +∆Hb)− 〈Q〉(P,Hb −∆Hb)) /2∆Hb. (17)
The choice of ∆Hb = 0.1Hb was found to produce sufficiently accurate values of the numerical
derivative across the range of bias fields studied. The results for periods P = 140 through
190 MCSS are shown in Fig. 5. A linear relationship is seen to exist between χ̂L and X
L , for
each value of P , over a wide range of χ̂L values. At each period, the dependence becomes
nonlinear below a certain value of χ̂L, as illustrated for periods P = 150, 170, and 190 MCSS
in Fig. 6. Since low values of the susceptibility χ̂L correspond to large values of the bias
field Hb, we interpret this breakdown of linearity as an indication that the FDR in Eq. (15)
holds only in a scaling regime around the critical point, i.e., for a limited range of Hb around
Hb = 0.
The relationship between X
L and χ̂L at the higher periods, P = 220 and P = 250 MCSS,
is more complicated, as shown in Fig. 7. At P = 220 MCSS, following the nonlinear regime
at low χ̂L, there is a linear relationship with slope Teff = (6.27± 0.11)J up to χ̂L ≈ 13 J
followed by a second distinct linear dependence with slope Teff ≈ (4.16 ± 0.29)J above
χ̂L = 13 J
−1. At P = 250 MCSS, the initial nonlinear dependence is again present. Then
the first linear regime has Teff = (6.49 ± 0.07)J up to χ̂L ≈ 13 J
−1, and is followed by a
regime which can be characterized as either linear with very gentle slope (0.27± 0.22)J , or
as an effective ‘saturation’ of X
L past χ̂L = 13 J
In Fig. 8 we plot the best-fit slopes from Figs. 5 and 7, which according to Eq. (15)
represent estimates of Teff , vs the scaling parameter θ = (P − Pc) /Pc. We have included
in the plot the slopes of both linear regimes for the values θ = 0.606 and 0.825 (P = 220
and 250 MCSS). For θ below 0.4 (P ≈ 190 MCSS), Teff increases with θ in a way not
inconsistent with a linear relationship (with slope 2.97J). It may be interesting to note that
an extrapolation of the linear relationship to θ = 0 (P = Pc) yields the value Teff = 3.39J ,
which is significantly higher than the critical temperature, Tc = 2.2619J , of the equilibrium
Ising system. However, one should not put too much emphasis on the numerical values of
Teff [45], as they could easily be changed. For instance, if Hb is only proportional to the full
conjugate field Hc, with a proportionality constant different from unity, this would trivially
change Teff in Eq. (15). The important result, which we have demonstrated to hold in the
critical region, is the linear relationship between X
L and χ̂L.
We can thus characterize the extent of applicability of an FDR to the DPT above the
critical period as follows. For θ < 0.4, an FDR holds outside of a small nonlinear regime at
low χ̂L (high Hb), with an effective temperature Teff which increases approximately linearly
with θ. For θ above 0.4, two linear relationships appear to exist between X
L and χ̂L in
separate regimes, making it impossible to define a unique Teff at a given value of θ. An
understanding of the nonlinear regime, which is present at low χ̂L for all periods examined,
as well as of the complication of the FDR above θ = 0.4, would be highly desirable. We hope
that these numerical results can stimulate the development of, as well as test the accuracy
of, a theoretical description of the non-equilibrium steady states produced in the presence
of non-zero Hb for this DPT.
VII. COMPARISON OF SECOND SCALING FUNCTION TO COMPUTA-
TIONAL RESULTS
We now test the asymptotic scaling forms for G+ in Eq. (14). First, we note that in
performing least-squares fits, one normally requires the goodness-of-fit parameter q, i.e., the
probability that (assuming the fit relationship were true) random error alone could produce
the observed data, to be greater than 10−3 to consider the fit reasonable. Within this section,
however, and in the captions to Figs. 9 through 11, it will be useful for descriptive purposes
to refer to scaling exponents resulting from attempts at least-squares fits with q values below
this acceptable range. We will refer to the results of such unsuccessful fitting attempts as
‘nominal’ scaling exponents, and for clarity will report the value of the parameter q for each
scaling exponent presented in this section.
In Fig. 9, we plot the scaling function G+ vs y1 for y2 = 8.46 at L = 180, using γe and νe
to calculate G+, and evaluating χ̂L numerically according to Eq. (17). In addition, we plot
in the same figure the scaling function GX+ (y1, y2) ≡ X
−γ/ν vs y1, for the same values of
y2 and L. A fit to all four G+ data points yields a scaling exponent −1.60± 0.03 (q = 0.02),
while a fit to the last three G+ data points yields a scaling exponent −1.71±0.05 (q = 0.25).
We will provide evidence in the next paragraph that only the last three G+ data points,
and not the first, satisfy the asymptotic condition y1 ≫ y2. Thus, these data are consistent
with power-law scaling of χ̂L with exponent −γd = −γe = −7/4 = −1.75. Attempts to
fit the GX+ data to all four and the last three data points yield nominal scaling exponents
−1.73 ± 0.01 (q < 10−15) and −1.81 ± 0.02 (q = 2.3 × 10−14), respectively. Thus, while in
Fig. 9 it appears that the power-law relationships with these nominal scaling exponents give
respectable visual fits to the GX+ data, there are variations in the data which, while small,
are larger than the statistical error bars, and which prevent a statistically significant fit. We
will describe the causes of these variations later in this section.
We present in Fig. 10 a plot of G+ and G
+ vs y1 at y2 = 0, again for L = 180, for a
larger range from y1 = 30.3 to 477. With y2 = 0, we expect that y1 = 30.3 (and indeed,
essentially any nonzero value of y1) should satisfy the asymptotic scaling condition y1 ≫ y2.
An attempt to fit to all six G+ data points in Fig. 10 gives a nominal scaling exponent
−1.65 ± 0.03 (q = 3.8 × 10−7), while a fit to the first five points (y1 = 30.3 through 280)
gives a scaling exponent −1.74±0.03 (q = 0.07). Excluding the first data point at y1 = 30.3
has little effect on either fit. This supports the assumption that with y2 = 0, the asymptotic
scaling condition y1 ≫ y2 holds for y1 = 30.3, while for y2 = 8.46, as used in Fig. 9, the
asymptotic scaling condition does not hold for y1 = 30.3. Attempts to fit all of, and the first
five of, the GX+ data points to power-law scaling again yield only nominal scaling exponents
−2.01± 0.01 (q < 10−15) and −2.10± 0.01 (q < 10−15), respectively.
We considered that the low statistical significance of the fits to the GX+ data could be
caused by underestimation of the error bars on X
L . These error bars were calculated by
(i) finding the correlation time in the numerical data series Qi from the simulation, and
sampling data at intervals of twice the correlation time; (ii) dividing this sampled data into
k > 16 groups and calculating the value of X
L within each group; (iii) finding the mean and
standard error of this collection of X
L values. As a check on self-consistency, we performed
several independent calculations of X
L by this method, and found that the standard error
of these values (corrected for small sample size) was comparable to the standard error found
within each calculation. Thus, we have strong evidence that the error bars for X
L (and G
are accurate.
These scaling results can be understood in light of the observations in Sec. VI on the
relationship between X
L and χ̂L. We can reasonably assume that the breakdown in scaling
of G+ past y1 = 280 (P = 350 MCSS, θ = 1.56) in Fig. 10 occurs because this is the boundary
of the critical region. The small variations of the GX+ data for P < 350 MCSS in Fig. 10 from
a scaling relationship with exponent −δd = −1.75 then have three main causes. The first
cause is the multiplication of the accurately scaling function G+ by the θ- and y1-dependent
value Teff , according to Eq. (15). However, such a variation would also occur in an analogous
plot for the scaling of XML ≡ L
2 (〈m2〉 − 〈m〉2) vs y1,e = ǫL
1/ν = ((T − Tc) /Tc)L
1/ν in the
equilibrium Ising model, since the susceptibility χML ≡ ∂〈m〉/∂H scales with exponent −γe,
and XML is related to χ
L by the ǫ-dependent temperature T according to Eq. (16). This
effect is small enough to be neglected in equilibrium critical scaling, and, since the change in
Teff from θ = 0.02 to θ = 0.4 is comparable to the change in T from ǫ = 0.02 to ǫ = 0.4 in the
equilibrium transition, it can also be neglected here. The second cause is the presence of the
nonlinear regimes in the plots of X
L vs χ̂L at low χ̂L, resulting in non-zero X
L -intercepts
in the application of Eq. (15) to Figs. 5 and 7. Because of this, division of the X
L data
by the appropriate Teff values (given in the caption of Fig. 5) does not quite reproduce the
corresponding χ̂L data, and (even below θ = 0.4) does not quite result in scaling consistent
with δd = 1.75 with statistical signficance. The third and most signficant cause of the
variations of the GX+ data is the ‘doubly linear’ behavior observed in Fig. 7 for θ > 0.4,
which prevents identification of a unique Teff in this range.
The assumption that X
L can be used as a proxy for χ̂L is thus fairly well justified
close to the critical period, where Teff varies over a limited range and the more complicated
effects observed at θ > 0.4 are not relevant. This is supported by Fig. 9, where the data
points cluster closely around the line corresponding to power-law scaling with exponent
−1.73 ≈ −γe. However, because of the first two causes just described, there are small
systematic variations in the GX+ data which prevent a statistically significant fit to a pure
power-law relationship as a function of y1.
In order to clarify the relationship of these scaling results to those in previous work, we
also plot in Fig. 11 data of G
+ (y1, y2) ≡ X
−γ/ν ≡ (〈Q2〉 − 〈|Q|〉2)L−γ/ν vs y1 at y2 = 0,
again using the equilibrium values γe and νe to calculate G
+ . This can be directly compared
to Fig. 11(d) in Ref. [23], in which the quantity we call X
L was called X
L . Attempted fits
to all five data points and to the last four data points of G
+ in Fig. 10 produce nominal
scaling exponents −1.60 ± 0.02 and −1.69± 0.02 (both with q < 10−15). The agreement in
Fig. 11(d) of Ref. [23] of the line with slope −7/4 with the data for θ > θc must therefore
be viewed as qualitative. The method used in Ref. [23] to numerically estimate γd, however,
which involves finite-size scaling at the critical period, is fully consistent with the results of
this paper, since at each period with θ < 0.4 we have found that the FDR in Eq. (15) holds
to a very good approximation.
Finally, in Fig. 12, we plot GX+ vs y2 at P = Pc, to study its scaling in the regime y1 ≪ y2.
As just noted, the use of X
L as a proxy for χ̂L is well justified at P = Pc by our results.
Power-law scaling is perhaps suggested in the range 20 < y2 < 50 for L = 180, and it is
clearly obeyed from y2 = 8.42 to 84.2 for L = 256. The scaling exponent was determined as
(1− δd)/δd = −0.914± 0.030, which is consistent with the corresponding equilibrium value
(1− δe)/δe = −14/15 ≈ −0.933.
VIII. CONCLUSIONS AND OUTLOOK
In this article, we have continued the computational study of the dynamic phase transition
(DPT) in the two-dimensional kinetic Ising model exposed to a periodically oscillating field,
which was begun in Refs. [9, 10, 23]. We have established two distinct but related results
about the field conjugate to the dynamic order parameter. First, we have identified the
period-averaged magnetic field, or ‘bias field’, Hb as an important component of the full
conjugate field. This claim is supported by numerical evidence that the dynamic order
parameter and its susceptibility follow critical scaling with respect to Hb. In particular, the
scaling exponent δd of the conjugate field was determined for the first time, and found by
finite-size scaling analysis of large-scale kinetic Monte Carlo simulations to be equal to the
critical-isotherm exponent for the equilibrium Ising transition, δe = 15. Furthermore, in
agreement with previous results [23], the dynamic scaling exponents γd, βd, and νd were also
found to equal their equilibrium Ising counterparts, γe = 7/8, βe = 1/8, and νd = 1.
These results further strengthen previous numerical [9, 10, 23] and analytical [28, 32, 33]
claims that the DPT in a periodically driven two-dimensional kinetic Ising model belongs to
the universality class of the equilibrium two-dimensional Ising model. However, with respect
to the direct applicability of the symmetry arguments of Refs. [32, 33], we caution the reader
that what is claimed in the present paper (as well as in Ref. [28]) is only equivalence of the
phase transitions in the driven kinetic Ising model and the equilibrium Ising model. Outside
the critical region, it is neither clear how closely P −Pc and Hb play the roles of T −Tc and
the ordinary magnetic field, respectively, nor how closely the dynamic order parameter, Q,
corresponds to the average equilibrium magnetization. From our discussion of the FDR in
Sec. VI, it appears likely that one or more of these relations break down outside the critical
region. Much theoretical work remains to be done in this area.
The second main result of this article is that a fluctuation-dissipation relation (FDR),
that is, a proportionality relation between the scaled fluctuations X
L ≡ L
2 (〈Q2〉 − 〈Q〉2)
and the susceptibility χ̂L with a slope we have called Teff , holds for a range of periods above
Pc and for a range of bias fields around Hb = 0. We stress again that we have found the FDR
of Eq. (15) to hold only in the critical region in this nonequilibrium system, in contrast to the
equilibrium FDR of Eq. (16) which follows directly from the partition function, and which
thus holds everywhere. We note that, for the parameters used in our computation at least,
the critical region in which the nonequilibrium FDR holds (P < 190 MCSS) is somewhat
smaller than the critical region in which power-law scaling is obeyed (P < 350 MCSS).
In previous work, when the conjugate field had not been identified, the scaled fluctuations
L were used as a proxy for the (then unknown) quantity χ̂L. The evidence for the FDR
presented here shows this assumption to be fully justified at the critical period (see Fig. 12),
and to be a very good approximation – nearly as good as the use of the scaled fluctuations
as a proxy for the susceptibility in the equilibrium Ising model – in the critical region where
the FDR holds.
There are at least three further computational projects suggested by the progress re-
ported here. The first is to investigate whether the field Hb functions as the conjugate field,
with scaling exponents consistent with the equilibrium Ising transition for periods P < Pc,
below the critical period. In the equilibrium system, the study of critical scaling in nonzero
field for T < Tc is complicated by the long time correlations and strong finite-size effects
which accompany the bimodal distributions of magnetization below Tc. Similar effects would
complicate the investigation of scaling with respect to Hb in the DPT for P < Pc, but the
advanced techniques [46, 47] which make the equilibrium simulations tractable do not ex-
tend obviously to the nonequilibrium case. The second computational project suggested is
to determine the nature of the full conjugate field Hc. The third project would be to study
the FDR at different values of the temperature, T , and the amplitude, H0, of the driving
field. Finally, we remark that it would be very desirable to extend the current understanding
of the theory of nonequilibrium steady states to include the conjugate field Hb, the FDR
found in the critical region, and the scaling of Hb.
Acknowledgments
Research at Florida State and Mississippi State Universities was supported by NSF Grant
No. DMR-0444051, and at Clarkson University by NSF Grant No. DMR-0509104. This
research also used resources of the National Center for Computational Sciences at Oak Ridge
National Laboratory, which is supported by the Office of Science of the U.S. Department of
Energy under Contract No. DE-AC05-00OR22725.
[1] T. Tomé and M. J. de Oliveira, Phys. Rev. A 41, 4251 (1990).
[2] J. F. F. Mendes and E. J. S. Lage, J. Stat. Phys. 64, 653 (1991).
[3] M. F. Zimmer, Phys. Rev. E 47, 3950 (1993).
[4] M. Acharyya and B. K. Chakrabarti, Phys. Rev. B 52, 6550 (1995).
[5] M. Acharyya, Phys. Rev. E 56, 2407 (1997).
[6] M. Acharyya, Phys. Rev. E 58, 179 (1998).
[7] W. S. Lo and R. A. Pelcovits, Phys. Rev. A 42, 7471 (1990).
[8] M. Acharyya, Phys. Rev. E 56, 1234 (1997).
[9] S. W. Sides, P. A. Rikvold, and M. A. Novotny, Phys. Rev. Lett. 81, 834 (1998).
[10] S. W. Sides, P. A. Rikvold, and M. A. Novotny, Phys. Rev. E 59, 2710 (1999).
[11] B. K. Chakrabarti and M. Acharyya, Rev. Mod. Phys. 71, 847 (1999).
[12] T. Yasui, H. Tutu, M. Yamamoto, and H. Fujisaka, Phys. Rev. E 66, 036123 (2002).
[13] N. Fujiwara, H. Tutu, and H. Fujisaka, Phys. Rev. E 70, 066132 (2004).
[14] N. Fujiwara, H. Tutu, and H. Fujisaka, Prog. Theor. Phys. Suppl. 161, 181 (2006).
[15] Y. Z. Shao, W. R. Zhong, and G. M. Lin, Acta Physica Sinica 53, 3165 (2004).
[16] M. Acharyya, Int. J. Mod. Phys. C 14, 49 (2003).
[17] M. Acharyya, Phys. Rev. E 69, 027105 (2004).
[18] H. Jang and M. J. Grimson, Phys. Rev. E 63, 066119 (2003).
[19] H. Jang, M. J. Grimson, and C. K. Hall, Phys. Rev. B 67, 094411 (2003).
[20] H. Jang, M. J. Grimson, and C. K. Hall, Phys. Rev. E 68, 046115 (2003).
[21] E. Machado, G. M. Buend́ıa, P. A. Rikvold, and R. M. Ziff, Phys. Rev. E 71, 016120 (2005).
[22] G. M. Buend́ıa, E. Machado, and P. A. Rikvold, J. Mol. Struct.: THEOCHEM 769, 189
(2006).
[23] G. Korniss, C. J. White, P. A. Rikvold, and M. A. Novotny, Phys. Rev. E 63, 016120 (2001).
[24] G. Korniss, P. A. Rikvold, and M. A. Novotny, Phys. Rev. E 66, 056127 (2002).
[25] A. Chatterjee and B. K. Chakrabarti, Phys. Rev. E 67, 046113 (2003).
[26] A. Chatterjee and B. K. Chakrabarti, Phase Transitions 77, 581 (2004).
[27] M. Acharyya, Phys. Rev. E 59, 218 (1999).
[28] H. Fujisaka, H. Tutu, and P. A. Rikvold, Phys. Rev. E 63, 036109 (2001); 63, 059903(E)
(2001).
[29] H. Tutu and N. Fujiwara, J. Phys. Soc. Jpn. 73, 2680 (2004).
[30] E. Z. Meilikhov, JETP Letters 79, 620 (2004).
[31] S. B. Dutta, Phys. Rev. E 69, 066115 (2004).
[32] G. Grinstein, C. Jayaprakash, and Y. He, Phys. Rev. Lett. 55, 2527 (1985).
[33] K. E. Bassler and B. Schmittmann, Phys. Rev. Lett. 73, 3343 (1994).
[34] D. T. Robb, Y. H. Xu, O. Hellwig, J. McCord, A. Berger, M. A. Novotny, and P. A. Rikvold,
Evidence for a dynamic phase transition in [Co/Pt]3 magnetic multilayers (2007), to be sub-
mitted to Phys. Rev. B.
[35] M. Bander and D. L. Mills, Phys. Rev. B 38, R12015 (1988).
[36] C. Back, C. Würsch, A. Vaterlaus, U. Ramsperger, U. Maier, and D. Pescia, Nature 378, 597
(1995).
[37] L. Onsager, Phys. Rev. 65, 117 (1944).
[38] D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics
(Cambridge University Press, Cambridge, UK, 2000).
[39] K. Binder and D. P. Landau, Phys. Rev. B 30, 1477 (1984).
[40] V. Privman and M. E. Fisher, Phys. Rev. B 30, 322 (1984).
[41] V. Privman, in V. Privman, ed., Finite Size Scaling and Numerical Simulation of Statistical
Systems (World Scientific, Singapore, 1990), pp. 4–98.
[42] Until the theory of the DPT, e.g. in Ref. [28], is extended to include Hc, it is unclear in
what form it should enter the dimensionless scaling parameter, i.e. as Hc/J , Hc/kBT , or
otherwise. Here we have chosen the form Hc/J . Since our numerical data is all taken at a
single temperature, T = 0.8Tc = 1.8152J , a change in this form would at most multiply
values of y2 by a constant factor, and would not affect the scaling relations reported in this
paper.
[43] K. Hayashi and S.-I. Sasa, Phys. Rev. E 69, 066119 (2004).
[44] K. Hayashi and S.-I. Sasa, Phys. Rev. E 71, 046143 (2005).
[45] To avoid misunderstandings, note that it is certainly possible to define other sensible effective
temperatures for an Ising-like model, so that one should be clear about the definition being
used in a given context. For example, the nonequilibrium FDR proposed in Refs. [48, 49]
connects two-time response and correlation functions, and applies to critical ageing in a two-
parameter family of Ising-like models. This family includes many models without detailed
balance (e.g., the noisy voter and majority models), and one model – the equilibrium Ising
model with Glauber dynamics – with detailed balance. In contrast, our Eq. (15) concerns
single-time response and correlations, and the kinetic Ising model driven by an oscillating
field is not included in the two-parameter family of models considered in Refs. [48, 49].
[46] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992).
[47] B. A. Berg and T. Neuhaus, Phys. Lett. B 267, 249 (1991).
[48] C. Godrèche and J. M. Luck, J. Phys. A: Math. Gen. 33, 9141 (2000).
[49] F. Sastre, I. Dornic, and H. Chaté, Phys. Rev. Lett. 91, 267205 (2003).
0.0001 0.001 0.01 0.1
H  (units of   )
L = 90
L = 128
L = 180
L = 256
1/     = 0.0673 (fit)
1/     = 0.0666 (ref.)
FIG. 1: (Color online.) Log-log plot of the dynamic order parameter, 〈Q〉, vs bias field, Hb, at
P = Pc for L = 90, 128, 180, and 256. A least-squares fit to power-law scaling of the L = 256 data,
in the range between the labels A and B above, produced a statistically significant fit with scaling
exponent 1/δd = 0.0673±0.0008 (corresponding to δd = 14.85±0.18). The dotted line corresponds
to the scaling exponent 1/δd = 0.0673. A reference line representing scaling with the equilibrium
Ising exponent, δe = 15 (1/δe = 0.0666), is shown as the dashed line.
1 10 100 1000
y  =   L 
0.001
L = 90
L = 128
L = 180
Slope = -1.76 (fit)
Slope = -1.75 (ref.)
y   = 3.39
y   = 8.46
y   = 16.9
y   = 33.9
y   = 84.6
y   = 169
1 10 100 1000
y  =    L
0.001
L = 90
L = 128
L = 180
Slope = -1.76 (fit)
Slope = -1.75 (ref.)
FIG. 2: (Color online.) Log-log plots of the scaling function F+(y1, y2) vs y1 for lattice sizes
L = 90, 128, and 180. The values of y1 plotted are 4.00, 17.2, 30.3, 43.4, 69.7, 147, and (in (b))
280 and 477. (a) The data are shown for the values of y2 listed on the plot. The best-fit line for
the last five points of the L = 180 data at y2 = 3.39, with slope −1.76 ± 0.07, is included along
with a reference line with the slope −γe = −1.75. (b) The data for y2 = 3.39 with two additional
y1 values illustrates the boundary of the regime of power-law scaling.
1 10 100
0.001
L = 90
L= 128
L = 180
Slope = 1.01 (fit)
Slope = 1.00 (ref.)
ν y  = 0.00
y  = 4.00
y  = 17.1
y  = 30.3
y  = 43.4
y  = 69.7
y  = 149
FIG. 3: (Color online.) Log-log plot of the scaling function F+(y1, y2) vs y2 = (Hb/J)L
βδ/ν for
lattice sizes L = 90, 128, and 180, for the constant values of y1 labeled in the plot. The values of
y2 used are y2 = 3.39, 8.46, 16.9, 33.9, 84.6, and 169. The dotted line represents the best fit to the
first five points of the L = 180 data at y1 = 149, and has a slope of 1.01 ± 0.01. The dashed line
shows the slope value of 1, expected from Eq. (13).
10 100
L = 128
L = 180
L = 256
Slope = 0.0673 (fit)
Slope = 0.0666 (ref.)
FIG. 4: (Color online.) Log-log plot of the scaling function F+(y1, y2) vs y2 = (Hb/J)L
βδ/ν for
lattice sizes L = 128, 180, and 256, at the critical period Pc, where y1 = 0. In the L = 256 data,
near the values y1 = 16, 32, 165, and 335, two closely spaced data points are actually plotted. The
best-fit line to the L = 256 data in the range 8.46 < y2 < 84.6, shown as a dotted line in the
plot, corresponds to a scaling exponent 1/δd = 0.0673 ± 0.0008. A reference line corresponding to
scaling exponent 1/δe = 1/15 = 0.0666 is also shown. These results are in complete agreement
with those shown in Fig. 1.
0 100 200
P = 140 MCSS
P = 150 MCSS
P = 160 MCSS
P = 170 MCSS
P = 190 MCSS
χ̂ (units of     )-1
FIG. 5: (Color online.) The scaled fluctuations X
L of the dynamic order paramater plotted vs its
susceptibility χ̂L to the bias field Hb, calculated at L = 180, for periods P = 140, 150, 160, 170 and
190 MCSS. The quantity χ̂L was calculated using the numerical derivative in Eq. (17). The best-fit
lines shown, whose slopes increase monotically with the period P of the data to which they were fit,
were calculated as (3.239J)χ̂L+10.42, (3.557J)χ̂L+8.735, (3.980J)χ̂L+3.889, (4.232J)χ̂L+5.868
and (4.497J)χ̂L + 5.371, respectively.
0 10 20 30
P = 140 MCSS
P = 150 MCSS
P = 160 MCSS
P = 170 MCSS
P = 190 MCSS
χ̂ (units of     )-1
FIG. 6: (Color online.) Closeup of Fig. 5, showing the relationship of X
L and χ̂L at low values
of χ̂L, which correspond to large values of the bias field Hb. For P = 150, 170 and 190 MCSS,
data have been taken (and are shown) down to very low values of χ̂L, where the breakdown of the
linear relationship between X
L and χ̂L can be clearly seen. The dashed lines are the same best-fit
lines shown in Fig. 5.
0 10 20 30 40
P = 220 MCSS
P = 250 MCSS
χ̂ (units of     )-1
FIG. 7: (Color online.) The scaled fluctuations X
L of the dynamic order paramater plotted vs its
susceptibility χ̂L to the bias fieldHb, calculated for lattice size L = 180, at periods P = 220 and 250
MCSS. At each period, the data were fit (purely phenomenologically) to two linear relationships.
For P = 220 MCSS, the fits were calculated as (6.265J)χ̂L−7.497 at low χ̂L, and (4.161J)χ̂L+19.94
at high χ̂L. For P = 250 MCSS, the fits were (6.485J)χ̂L−9.240 at low χ̂L, and (0.2726J)χ̂L+67.92
at high χ̂L.
0 0.2 0.4 0.6 0.8 1
FIG. 8: (Color online.) The effective temperature Teff , obtained as the slopes of the linear fits to
the data in Figs. 5 and 7, plotted vs θ = (P − Pc) /Pc. For the values θ = 0.606 and 0.825 (P = 220
and 250 MCSS), the slopes of both linear regimes fit in Fig. 7 are plotted as Teff values, using filled
squares and diamonds rather than filled circles. The straight line is a weighted least-squares fit to
the data below θ ≈ 0.4 (P < 190 MCSS), and has a slope of 2.97J .
20 50 100 200
y  =    L
0.001
Slope = -1.73 (nom.)
Slope = -1.81 (nom.)
Slope = -1.75 (ref.)
Slope = -1.60 (fit)
Slope = -1.71 (fit)
FIG. 9: (Color online.) Log-log plots of G+(y1, y2) and G
+ (y1, y2) vs y1, over the range y1 = 30.3
through 149, for y2 = 8.46 at L = 180. The relatively small error bars on each data point can
be seen inside the larger symbols. The solid and dash-dash-dotted lines are fits to all four and
the last three G+ data points, respectively, and correspond to scaling exponents −1.60 ± 0.03 and
−1.71±0.05. The dotted and dash-dotted lines are the result of attempts to fit all four and the last
three GX+ data points, respectively. They correspond to nominal scaling exponents −1.73 ± 0.01
and −1.81± 0.02. The dashed line is a reference line corresponding to scaling exponent −1.75.
20 50 100 200 500
y  =    L
0.0001
0.001
Slope = -2.01 (nom.)
Slope = -2.10 (nom.)
Slope = -1.75 (ref.)
Slope = -1.65 (nom.)
Slope = -1.74 (fit)
FIG. 10: (Color online.) Log-log plots of G+(y1, y2) and G
+ (y1, y2) vs y1, over the range y1 = 30.3
through 477, for y2 = 0 at lattice size L = 180. The relatively small error bars on each data
point can be seen inside the larger symbols. The solid and dash-dash-dotted lines show the result
of attempts to fit all six and the first five G+ data points, and correspond to a nominal scaling
exponent −1.65 ± 0.03, and a statistically significant scaling exponent −1.74 ± 0.03, respectively.
The dotted and dash-dotted lines are the result of attempts to fit all six and the first five GX+ data
points, and correspond to nominal scaling exponents −2.01 ± 0.01 and −2.10± 0.01, respectively.
The dashed line is a reference line corresponding to scaling exponent −1.75.
20 50 100 200
y  =    L
0.001
Slope = -1.59 (nom.)
Slope = -1.70 (nom.)
Slope = -1.75 (ref.)
FIG. 11: (Color online.) Log-log plot of G
+ (y1, y2) vs y1, over the range y1 = 30.3 through
149, for y2 = 0 at lattice size L = 180. The relatively small error bars on each data point can
be seen inside the larger symbols. The solid and dotted lines are the results of attempts to fit
all five and the last four data points with a power-law relationship, and correspond to nominal
scaling exponents of −1.59 ± 0.02 and −1.70 ± 0.03, respectively. The dashed line is a reference
line corresponding to a scaling exponent of −1.75.
10 100
0.001
L = 128
L = 180
L = 256
Slope = -0.914 (fit)
Slope = -0.933 (ref.)
FIG. 12: (Color online.) Log-log plot of GX+ (y1, y2) vs y2 = (Hb/J)L
βδ/ν for lattice sizes L = 128,
180, and 256, at the critical period Pc, where y1 = 0. The best-fit line to the L = 256 data in
the range 8.46 < y2 < 84.6, shown as a dotted line in the plot, corresponds to a scaling exponent
(1 − δd)/δd = −0.914 ± 0.029. A reference line (dashed), corresponding to a scaling exponent
(1− δe)/δe = −14/15 ≈ −0.933, is also shown.
	Introduction
	Computational model 
	Scaling with respect to the bias field 
	Finite-size scaling analysis with bias field
	Comparison of first scaling function to computational results
	Applicability of a fluctuation-dissipation relation
	Comparison of second scaling function to computational results
	Conclusions and outlook
	Acknowledgments
	References

  The two-dimensional kinetic Ising model, when exposed to an oscillating
applied magnetic field, has been shown to exhibit a nonequilibrium,
second-order dynamic phase transition (DPT), whose order parameter Q is the
period-averaged magnetization. It has been established that this DPT falls in
the same universality class as the equilibrium phase transition in the
two-dimensional Ising model in zero applied field. Here we study for the first
time the scaling of the dynamic order parameter with respect to a nonzero,
period-averaged, magnetic `bias' field, H_b, for a DPT produced by a
square-wave applied field. We find evidence that the scaling exponent,
\delta_d, of H_b at the critical period of the DPT is equal to the exponent for
the critical isotherm, \delta_e, in the equilibrium Ising model. This implies
that H_b is a significant component of the field conjugate to Q. A finite-size
scaling analysis of the dynamic order parameter above the critical period
provides further support for this result. We also demonstrate numerically that,
for a range of periods and values of H_b in the critical region, a
fluctuation-dissipation relation (FDR), with an effective temperature
T_{eff}(T, P, H_0) depending on the period, and possibly the temperature and
field amplitude, holds for the variables Q and H_b. This FDR justifies the use
of the scaled variance of Q as a proxy for the nonequilibrium susceptibility,
\partial<Q> / \partial H_b, in the critical region.

Introduction
	Computational model 
	Scaling with respect to the bias field 
	Finite-size scaling analysis with bias field
	Comparison of first scaling function to computational results
	Applicability of a fluctuation-dissipation relation
	Comparison of second scaling function to computational results
	Conclusions and outlook
	Acknowledgments
	References
