Conjugate field and fluctuation-dissipation relation for the dynamic phase transition in the two-dimensional kinetic Ising model

TL;DR

This study investigates the dynamic phase transition in the 2D kinetic Ising model under oscillating magnetic fields, revealing conjugate field relations and a fluctuation-dissipation relation that extends understanding of nonequilibrium critical phenomena.

Contribution

It demonstrates that the bias magnetic field acts as a conjugate field to the dynamic order parameter and establishes a fluctuation-dissipation relation in the nonequilibrium setting.

Findings

The scaling exponent of the bias field matches the equilibrium critical isotherm exponent.

A fluctuation-dissipation relation holds for the dynamic order parameter and bias field.

Finite-size scaling supports the conjugate field relationship.

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, 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.

Equations

$$\mathcal{H}=-J\sum _{⟨i,j⟩}{S}_i{S}_j-H(t)\sum _i{S}_i,$$

$$W\left(S_i\to -S_i\right)=\frac{1}{1+\exp \left(\Delta E/T\right)},$$

$$m(t)=\frac{1}{L^2}\sum _{i=1}^{L^2}{S}_i(t).$$

$$Q_i=\frac{1}{P}{\int }_{(i-1)P}^{iP}m(t)dt.$$

$$H_b=\frac{1}{P}{\int }_0^{P}H(t)dt.$$

$$m\left(T=T_c,H\to 0\right)\propto H^{1/{\delta }_{\mathrm{e}}},$$

$$⟨Q⟩\left(P=P_c,H_b\to 0\right)\propto H_b^{1/{\delta }_{\mathrm{d}}}$$

$$f_L\approx L^{-d}{Y}_{\pm }\left(|ϵ|L^{1/\nu },H{L}^{\beta \delta /\nu }\right),$$

$$m_L=\frac{\partial f}{\partial H}=L^{-\beta /\nu }{\mathcal{F}}_{0\pm }\left(|ϵ|L^{1/\nu },H{L}^{\beta \delta /\nu }\right),$$

$${\chi }_L=\frac{\partial m_L}{\partial H}=L^{\gamma /\nu }{\mathcal{G}}_{0\pm }\left(|ϵ|L^{1/\nu },H{L}^{\beta \delta /\nu }\right),$$

$$⟨Q{⟩}_L=L^{-\beta /\nu }{\mathcal{F}}_{\pm }\left(|\theta |L^{1/\nu },(H_c/J)L^{\beta \delta /\nu }\right),$$

$${\hat{\chi }}_L=L^{\gamma /\nu }{\mathcal{G}}_{\pm }\left(|\theta |L^{1/\nu },(H_c/J)L^{\beta \delta /\nu }\right),$$

$$\mathcal{F}_{+}(y_{1},y_{2})\equiv L^{\beta/\nu}\langle Q\rangle_{L}\sim\left\{\begin{array}[]{lll}y_{1}^{-\gamma}y_{2}&\mbox{for}&y_{1}\gg y_{2}\\ y_{2}^{1/\delta}&\mbox{for}&y_{1}\ll y_{2}\end{array}\right.$$

$$\mathcal{G}_{+}(y_{1},y_{2})\equiv L^{-\gamma/\nu}\hat{\chi}_{L}\sim\left\{\begin{array}[]{lll}y_{1}^{-\gamma}&\mbox{for}&y_{1}\gg y_{2}\\ y_{2}^{(1-\delta)/\delta}&\mbox{for}&y_{1}\ll y_{2}\end{array}\right.\;.$$

$${\hat{\chi }}_L\equiv \frac{\partial ⟨Q{⟩}_L}{\partial H_b}=\frac{L^2\left(⟨Q^2{⟩}_L-⟨Q{⟩}_L^2\right)}{T_{\mathrm{eff}}}\equiv \frac{X_L^Q}{T_{\mathrm{eff}}},$$

$${\chi }_L\equiv \frac{\partial ⟨m{⟩}_L}{\partial H}=\frac{L^2\left(⟨m^2{⟩}_L-⟨m{⟩}_L^2\right)}{T},$$

$${\hat{\chi }}_L(P,H_b)\approx \left(⟨Q⟩(P,H_b+\Delta H_b)-⟨Q⟩(P,H_b-\Delta H_b)\right)/2\Delta H_b.$$

Full text

I Introduction

The dynamic phase transition (DPT) in a ferromagnetic system below its critical temperature was first observed in numerical solutions of a mean-field model exposed to an oscillating magnetic field Tomé and de Oliveira (1990); Mendes and Lage (1991). It was then studied further, both in mean-field models Zimmer (1993); Acharyya and Chakrabarti (1995); Acharyya (1997a, 1998) and in kinetic Monte Carlo (KMC) simulations Lo and Pelcovits (1990); Acharyya (1997b, a, 1998); Sides et al. (1998, 1999). A review of this early work can be found in Ref. Chakrabarti and Acharyya (1999). 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 Yasui et al. (2002); Fujiwara et al. (2004, 2006), KMC simulations of a three-dimensional Ising system Shao et al. (2004), and KMC simulations of a uniaxially anisotropic Heisenberg system in an off-axial field Acharyya (2003), an elliptically polarized applied field Acharyya (2004), and with the effect of a thin-film surface energy Jang and Grimson (2003); Jang et al. (2003a, b). The phenomenon has also been observed in simulations of CO oxidation under oscillating CO pressure Machado et al. (2005); Buendía et al. (2006). Further simulation studies of the DPT in the two-dimensional kinetic Ising model have appeared Korniss et al. (2001, 2002); Chatterjee and Chakrabarti (2003, 2004); Acharyya (1999), as well as analytical studies of the DPT Fujisaka et al. (2001); Tutu and Fujiwara (2004); Meilikhov (2004); Dutta (2004).

Here, we concentrate on the DPT in the two-dimensional kinetic Ising model. It was observed in simulations of this model that there exists a singularity at a critical period of the applied oscillating field Sides et al. (1998, 1999), and that the critical exponents $\beta$ and $\gamma$ (and, with less accuracy, $\nu$) are consistent with the universality class of the equilibrium two-dimensional Ising transition in zero field Korniss et al. (2001) $\left(\beta=1/8,\gamma=7/4,\mathrm{and}~{}\nu=1\right)$. 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 correlation 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}^{Q}=L^{2}\left(\langle Q^{2}\rangle-\langle Q\rangle^{2}\right)$, 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 Fujisaka et al. (2001). 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 Grinstein et al. (1985); Bassler and Schmittmann (1994).

Recently, experiments were performed on a $\left[\mathrm{Co(0.4nm)}/\mathrm{Pt(0.7nm)}\right]_{3}$ multilayer film with strong uniaxial anisotropy Robb et al. (2006), whose equilibrium behavior is known to be Ising-like Bander and Mills (1988); Back et al. (1995). The film was exposed to an oscillating (sawtooth) applied field with varying period, in the presence of constant ‘bias’ magnetic fields $H_{b}$ 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 $\delta_{\mathrm{d}}\approx\delta_{\mathrm{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. Korniss et al. (2001). Specifically, we perform kinetic Monte Carlo (KMC) simulations of a two-dimensional periodic square lattice of Ising spins $S_{i}$, which can take only the values $S_{i}=\pm 1.$ The Hamiltonian of the model is

[TABLE]

where $J>0$ is the ferromagnetic exchange interaction, $\sum_{\langle i,j\rangle}$ runs over all nearest-neighbor pairs, $\sum_{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 $H_{0}=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 Korniss et al. (2001). Other symmetric field shapes, such as sinusoidal Sides et al. (1998, 1999) and sawtooth Robb et al. (2006), 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

[TABLE]

where $\Delta 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.8T_{c}$, where $T_{c}=2.269J$ is the equilibrium critical temperature of the square-lattice Ising ferromagnet in zero applied field Onsager (1944).

The system responds to the oscillating field via the time-dependent magnetization per site,

[TABLE]

The dynamic order parameter is defined as the average of $m(t)$ over a given field cycle $i$ Tomé and de Oliveira (1990) :

[TABLE]

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,

[TABLE]

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 $H_{b}$ 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\rightarrow\infty)$ as

[TABLE]

where the critical exponent $\delta_{\mathrm{e}}=15$ Landau and Binder (2000). For finite systems, this relationship breaks down when the infinite-system correlation length, $\xi_{\infty}\left(T=T_{c},H\right)$, which diverges as $H\rightarrow 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 Binder and Landau (1984).

We can determine directly whether the non-equilibrium system exhibits a similar relationship,

[TABLE]

in an analogous way. In Fig. 1, we plot $\langle Q\rangle$ vs $H_{b}$ at the critical value of the period, $P=P_{c}$. In previous work, the reversal time for the magnetization, following instantaneous reversal of the uniform magnetic field $H$ at $\left(H=0.3J,T=0.8T_{c}\right)$, was found as $\tau=74.5977$ MCSS Sides et al. (1998, 1999). The critical scaled half-period for the square waveform was determined to be $\Theta_{c}=P_{c}/\left(2\tau\right)=0.918\pm 0.005$ Korniss et al. (2001). This yields $P_{c}=136.96\pm 0.75$ MCSS, and in our simulations and analysis in this paper we use $P_{c}=136.96$ MCSS.

A power-law dependence is indeed seen to hold in Fig. 1, within a range which extends to lower values of $H_{b}$ 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 $\delta_{\mathrm{d}}=14.85\pm 0.18$. As including points with $H_{b}\geq 0.01J$ was found to greatly reduce the statistical significance of the fit, the value $H_{b}=0.01J$ serves as a boundary of the scaling region at $P=P_{c}$. This result is consistent with an exponent $\delta_{\mathrm{d}}=\delta_{\mathrm{e}}=15$, suggesting that the bias field $H_{b}$, 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 $H_{b}$ is the dominant component of the field conjugate to $\langle Q\rangle$, in the next several sections we will demonstrate data collapse onto a two-parameter finite-size scaling function for the system-size dependent quantity $\langle Q\rangle_{L}$ at points $\left(P\geq P_{c},H_{b}>0\right)$, 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 Privman and Fisher (1984); Privman (1990) 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

[TABLE]

where $\epsilon=(T-T_{c})/T_{c}$, $H$ (in units of $k_{B}T$) is the field conjugate to the order parameter, $\nu$ is the critical exponent for the correlation length, $\beta$ is the exponent for the order parameter, $\delta$ is the exponent for the critical isotherm, and $Y_{\pm}$ are scaling functions above (+) and below ($-$) the critical point. This yields for the order parameter at finite $L$

[TABLE]

where the exponent for $L$ in the prefactor is obtained by using the hyperscaling relation $d\nu=2-\alpha$ and the exponent equality $\alpha=2-\beta(\delta+1)$. Further differentiation yields the susceptibility,

[TABLE]

where the exponent for $L$ in the prefactor is obtained by using the exponent equality $\gamma=\beta(\delta-1)$.

It has previously been shown analytically that the DPT for a sinusoidal applied field, which is symmetric under $H(t)\rightarrow-H(t+P/2)$ and so which can safely be assumed to have $H_{c}=0$, has an effective Ginzburg-Landau free-energy density in the same universality class as the equilibrium Ising model Fujisaka et al. (2001). It therefore appears reasonable to write corresponding scaling functions for the dynamic order parameter $\langle Q\rangle$ and its associated susceptibility $\hat{\chi}$,

[TABLE]

and

[TABLE]

where $\theta=(P-P_{c})/P_{c}$, and $H_{c}$ is the (as yet unknown) field conjugate to $\langle Q\rangle$. In this paper we express $H_{c}$ (and $H_{b}$) in units of the exchange constant, $J$, so that the second scaling parameter in Eqs. 11 and 12 is dimensionless. The specific form, $H_{c}/J$, with which $H_{c}$ 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 scaling.param.endnote . Computational results for sinusoidal and square-wave fields, which both are symmetric under $H(t)\rightarrow-H(t+P/2)$ and so presumably have $H_{c}=0$, have previously confirmed the scaling behavior with respect to $\theta$. The exponent values were determined as $\gamma/\nu=1.74\pm 0.05$, $\beta/\nu=0.126\pm 0.005$, and $\nu=0.95\pm 0.15$ Korniss et al. (2001), consistent with the exact values for the two-dimensional equilibrium Ising model, $\gamma=7/4=1.75$, $\beta=1/8=0.125$, and $\nu=1$.

We now determine the expected asymptotic forms of the scaling functions $\mathcal{F}_{+}(y_{1},y_{2})$ and $\mathcal{G}_{+}(y_{1},y_{2})$, where we emphasize that the $+$ subscript indicates that the scaling functions refer to the range $P\geq P_{c}$, and where the scaling parameters are $y_{1}\equiv\theta L^{1/\nu}$ and $y_{2}\equiv(H_{c}/J)L^{\beta\delta/\nu}$.

$y_{1}\gg y_{2}$. We expect $\hat{\chi}_{L}\sim\theta^{-\gamma}=L^{\gamma/\nu}y_{1}^{-\gamma}$ (independent of $y_{2}$) and $\langle Q\rangle_{L}=\hat{\chi}_{L}H_{c}\sim\theta^{-\gamma}H_{c}\sim L^{-\beta/\nu}y_{1}^{-\gamma}y_{2}$, where $\gamma=\beta(\delta-1)$ was used to obtain the exponent for $L$ in $\langle Q\rangle_{L}$.

$y_{1}\ll y_{2}$. We expect that $\langle Q\rangle_{L}\sim H_{c}^{1/\delta}\sim L^{-\beta/\nu}y_{2}^{1/\delta}$ and $\hat{\chi}_{L}=\partial\langle Q\rangle_{L}/\partial H_{c}\sim L^{\gamma/\nu}y_{2}^{(1-\delta)/\delta}$ (both independent of $y_{1}$), where $\gamma=\beta(\delta-1)$ was used to obtain the exponent for $L$ in $\hat{\chi}_{L}$.

Thus, the asymptotic forms of the scaling functions are expected to be

[TABLE]

and

[TABLE]

V Comparison of first scaling function to computational results

In Fig. 2(a), using the equilibrium values $\beta_{\mathrm{e}}$ and $\nu_{\mathrm{e}}$ in calculating $\mathcal{F}_{+}(y_{1},y_{2})\equiv L^{\beta/\nu}\langle Q\rangle_{L}$, we present a plot of the scaling function $\mathcal{F}_{+}$ vs $y_{1}$ for different values of $y_{2}$, 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 $y_{1}\gg y_{2}$, which is consistent with Eq. (13). At progressively larger values of the constant $y_{2}$, the power-law scaling can be seen to begin at increasing values of $y_{1}$, as would be expected. A best-fit line to the final five points of the $L=180$ data at $y_{2}=3.39$ yields an estimate of the scaling exponent $-\gamma_{\mathrm{d}}=-1.76\pm 0.07$ in Eq. (13). This is consistent with the previous results for $H_{c}=0$ cited above Korniss et al. (2001), and it supports the hypothesis that $\gamma_{\mathrm{d}}=\gamma_{\mathrm{e}}=7/4=1.75$. In Fig. 2(b), we present just the data for $y_{2}=3.39$, including additional data points at $y_{1}=280$ and $477$. The data deviate from the power-law behavior for $L=90$ at $y_{1}>149$, for $L=128$ at $y_{1}>280$, and for $L=180$ at $y_{1}>477$. This locates the boundary of the scaling regime (for $y_{2}=3.39$) at $\theta=y_{1}/L^{1/\nu}\approx 2.65$ .

In Fig. 3, again using $\beta_{\mathrm{e}}$ and $\nu_{\mathrm{e}}$ in calculating $\mathcal{F}_{+}$, we plot the scaling function $\mathcal{F}_{+}$ vs $y_{2}$ at different values of $y_{1}$, in order to examine the scaling behavior for $y_{1}\gg y_{2}$. For the constant values $y_{1}=43.4,69.7,$ and 149, power-law scaling can be observed in the regime $y_{1}\gg y_{2}$. A best-fit line to the $y_{1}=149$ data for the five points from $y_{2}=3.39$ to $84.6$ yields a scaling exponent of $1.01\pm 0.01$, which is consistent with the value of $1$ expected from Eq. (13).

In order to investigate the scaling of $\mathcal{F}_{+}$ in the asymptotic limit $y_{1}\ll y_{2}$, we plot in Fig. 4 the scaling function $\mathcal{F}_{+}(y_{1},y_{2})$ vs $y_{2}$ at the critical period $P=P_{c}$ (i.e., $y_{1}=0$), at lattice sizes $L=128,180,$ and 256. In the range $20<y_{2}<50$, power-law scaling is observed for all three lattice sizes. For $y_{2}>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 $y_{2}=8.46$ to $84.6$ produces a scaling exponent $0.0673\pm 0.0008$. Since the constant factors $L^{\beta\delta/\nu}$ and $L^{\beta/\nu}$ do not affect the fit of the scaling exponent, this is the same exponent found in the fit of $\langle Q\rangle$ vs $H_{b}$ in Fig. 1. The reciprocal of this scaling exponent is thus $\delta_{\mathrm{d}}=14.85\pm 0.18$, which is consistent with the exponent of the critical isotherm, $\delta_{\mathrm{e}}=15$, in the equilibrium Ising model.

The comparison of the second scaling function, $\mathcal{G}_{+}(y_{1},y_{2}$), to numerical data is more clearly presented after the relationship of the susceptibility $\hat{\chi}_{L}$ and the scaled variance $X_{L}^{Q}$ has been examined. Therefore, we present in the next section numerical results on the extent of applicability of an FDR between $\hat{\chi}_{L}$ and $X_{L}^{Q}$, 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 Hayashi and Sasa (2004, 2005). Here we will therefore consider whether the nonequilibrium susceptibility and the scaled variance of the dynamic order parameter can be related as

[TABLE]

with an effective temperature $T_{\mathrm{eff}}$, in a way analogous to the equilibrium FDR,

[TABLE]

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}^{Q}$ as a proxy for $\hat{\chi}_{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 $\hat{\chi}_{L}$ and $X_{L}^{Q}$ for a range of periods from $P=140$ to 250 MCSS and a range of bias fields from $H_{b}=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 $\hat{\chi}_{L}$ and $X_{L}^{Q}$ near zero, increases as the period is increased.) The computations were perfomed at $L=180$. The quantity $\hat{\chi}_{L}$ was computed directly as a numerical derivative:

[TABLE]

The choice of $\Delta H_{b}=0.1H_{b}$ 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 $\hat{\chi}_{L}$ and $X_{L}^{Q}$, for each value of $P$, over a wide range of $\hat{\chi}_{L}$ values. At each period, the dependence becomes nonlinear below a certain value of $\hat{\chi}_{L}$, as illustrated for periods $P=150$, 170, and 190 MCSS in Fig. 6. Since low values of the susceptibility $\hat{\chi}_{L}$ correspond to large values of the bias field $H_{b}$, 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 $H_{b}$ around $H_{b}=0$.

The relationship between $X_{L}^{Q}$ and $\hat{\chi}_{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 $\hat{\chi}_{L}$, there is a linear relationship with slope $T_{\mathrm{eff}}=(6.27\pm 0.11)J$ up to $\hat{\chi}_{L}\approx 13~{}J^{-1}$, followed by a second distinct linear dependence with slope $T_{\mathrm{eff}}\approx(4.16\pm 0.29)J$ above $\hat{\chi}_{L}=13~{}J^{-1}$. At $P=250$ MCSS, the initial nonlinear dependence is again present. Then the first linear regime has $T_{\mathrm{eff}}=(6.49\pm 0.07)J$ up to $\hat{\chi}_{L}\approx 13~{}J^{-1}$, and is followed by a regime which can be characterized as either linear with very gentle slope $(0.27\pm 0.22)J$, or as an effective ‘saturation’ of $X_{L}^{Q}$ past $\hat{\chi}_{L}=13~{}J^{-1}$.

In Fig. 8 we plot the best-fit slopes from Figs. 5 and 7, which according to Eq. (15) represent estimates of $T_{\mathrm{eff}}$, vs the scaling parameter $\theta=\left(P-P_{c}\right)/P_{c}$. We have included in the plot the slopes of both linear regimes for the values $\theta=0.606$ and 0.825 ($P=220$ and 250 MCSS). For $\theta$ below 0.4 ($P\approx 190$ MCSS), $T_{\mathrm{eff}}$ increases with $\theta$ 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 $\theta=0$ ($P=P_{c}$) yields the value $T_{\mathrm{eff}}=3.39J$, which is significantly higher than the critical temperature, $T_{c}=2.2619J$, of the equilibrium Ising system. However, one should not put too much emphasis on the numerical values of $T_{\mathrm{eff}}$ eff.temp.endnote , as they could easily be changed. For instance, if $H_{b}$ is only proportional to the full conjugate field $H_{c}$, with a proportionality constant different from unity, this would trivially change $T_{\mathrm{eff}}$ in Eq. (15). The important result, which we have demonstrated to hold in the critical region, is the linear relationship between $X_{L}^{Q}$ and $\hat{\chi}_{L}$.

We can thus characterize the extent of applicability of an FDR to the DPT above the critical period as follows. For $\theta<0.4$, an FDR holds outside of a small nonlinear regime at low $\hat{\chi}_{L}$ (high $H_{b}$), with an effective temperature $T_{\mathrm{eff}}$ which increases approximately linearly with $\theta$. For $\theta$ above 0.4, two linear relationships appear to exist between $X_{L}^{Q}$ and $\hat{\chi}_{L}$ in separate regimes, making it impossible to define a unique $T_{\mathrm{eff}}$ at a given value of $\theta$. An understanding of the nonlinear regime, which is present at low $\hat{\chi}_{L}$ for all periods examined, as well as of the complication of the FDR above $\theta=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 $H_{b}$ for this DPT.

VII Comparison of second scaling function to computational results

We now test the asymptotic scaling forms for $\mathcal{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 $\mathcal{G}_{+}$ vs $y_{1}$ for $y_{2}=8.46$ at $L=180$, using $\gamma_{\mathrm{e}}$ and $\nu_{\mathrm{e}}$ to calculate $\mathcal{G}_{+}$, and evaluating $\hat{\chi}_{L}$ numerically according to Eq. (17). In addition, we plot in the same figure the scaling function $\mathcal{G}^{X}_{+}(y_{1},y_{2})\equiv X_{L}^{Q}L^{-\gamma/\nu}$ vs $y_{1}$, for the same values of $y_{2}$ and $L$. A fit to all four $\mathcal{G}_{+}$ data points yields a scaling exponent $-1.60\pm 0.03$ ($q=0.02$), while a fit to the last three $\mathcal{G}_{+}$ data points yields a scaling exponent $-1.71\pm 0.05$ ($q=0.25$). We will provide evidence in the next paragraph that only the last three $\mathcal{G}_{+}$ data points, and not the first, satisfy the asymptotic condition $y_{1}\gg y_{2}$. Thus, these data are consistent with power-law scaling of $\hat{\chi}_{L}$ with exponent $-\gamma_{\mathrm{d}}=-\gamma_{\mathrm{e}}=-7/4=-1.75$. Attempts to fit the $\mathcal{G}^{X}_{+}$ data to all four and the last three data points yield nominal scaling exponents $-1.73\pm 0.01$ ($q<10^{-15}$) and $-1.81\pm 0.02$ ($q=2.3\times 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 $\mathcal{G}^{X}_{+}$ 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 $\mathcal{G}_{+}$ and $\mathcal{G}^{X}_{+}$ vs $y_{1}$ at $y_{2}=0$, again for $L=180$, for a larger range from $y_{1}=30.3$ to 477. With $y_{2}=0$, we expect that $y_{1}=30.3$ (and indeed, essentially any nonzero value of $y_{1}$) should satisfy the asymptotic scaling condition $y_{1}\gg y_{2}$. An attempt to fit to all six $\mathcal{G}_{+}$ data points in Fig. 10 gives a nominal scaling exponent $-1.65\pm 0.03$ ($q=3.8\times 10^{-7}$), while a fit to the first five points ($y_{1}=30.3$ through 280) gives a scaling exponent $-1.74\pm 0.03$ ($q=0.07$). Excluding the first data point at $y_{1}=30.3$ has little effect on either fit. This supports the assumption that with $y_{2}=0$, the asymptotic scaling condition $y_{1}\gg y_{2}$ holds for $y_{1}=30.3$, while for $y_{2}=8.46$, as used in Fig. 9, the asymptotic scaling condition does not hold for $y_{1}=30.3$. Attempts to fit all of, and the first five of, the $\mathcal{G}^{X}_{+}$ data points to power-law scaling again yield only nominal scaling exponents $-2.01\pm 0.01$ ($q<10^{-15}$) and $-2.10\pm 0.01$ ($q<10^{-15}$), respectively.

We considered that the low statistical significance of the fits to the $\mathcal{G}^{X}_{+}$ data could be caused by underestimation of the error bars on $X_{L}^{Q}$. These error bars were calculated by (i) finding the correlation time in the numerical data series $Q_{i}$ 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}^{Q}$ within each group; (iii) finding the mean and standard error of this collection of $X_{L}^{Q}$ values. As a check on self-consistency, we performed several independent calculations of $X_{L}^{Q}$ 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}^{Q}$ (and $\mathcal{G}^{X}_{+}$) are accurate.

These scaling results can be understood in light of the observations in Sec. VI on the relationship between $X_{L}^{Q}$ and $\hat{\chi}_{L}$. We can reasonably assume that the breakdown in scaling of $\mathcal{G}_{+}$ past $y_{1}=280$ ($P=350$ MCSS, $\theta=1.56$) in Fig. 10 occurs because this is the boundary of the critical region. The small variations of the $\mathcal{G}^{X}_{+}$ data for $P<350$ MCSS in Fig. 10 from a scaling relationship with exponent $-\delta_{\mathrm{d}}=-1.75$ then have three main causes. The first cause is the multiplication of the accurately scaling function $\mathcal{G}_{+}$ by the $\theta$- and $y_{1}$-dependent value $T_{\mathrm{eff}}$, according to Eq. (15). However, such a variation would also occur in an analogous plot for the scaling of $X_{L}^{M}\equiv L^{2}\left(\langle m^{2}\rangle-\langle m\rangle^{2}\right)$ vs $y_{1,\mathrm{e}}=\epsilon L^{1/\nu}=\left(\left(T-T_{c}\right)/T_{c}\right)L^{1/\nu}$ in the equilibrium Ising model, since the susceptibility $\chi_{L}^{M}\equiv\partial\langle m\rangle/\partial H$ scales with exponent $-\gamma_{\mathrm{e}}$, and $X_{L}^{M}$ is related to $\chi_{L}^{M}$ by the $\epsilon$-dependent temperature $T$ according to Eq. (16). This effect is small enough to be neglected in equilibrium critical scaling, and, since the change in $T_{\mathrm{eff}}$ from $\theta=0.02$ to $\theta=0.4$ is comparable to the change in $T$ from $\epsilon=0.02$ to $\epsilon=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}^{Q}$ vs $\hat{\chi}_{L}$ at low $\hat{\chi}_{L}$, resulting in non-zero $X_{L}^{Q}$-intercepts in the application of Eq. (15) to Figs. 5 and 7. Because of this, division of the $X_{L}^{Q}$ data by the appropriate $T_{\mathrm{eff}}$ values (given in the caption of Fig. 5) does not quite reproduce the corresponding $\hat{\chi}_{L}$ data, and (even below $\theta=0.4$) does not quite result in scaling consistent with $\delta_{\mathrm{d}}=1.75$ with statistical signficance. The third and most signficant cause of the variations of the $\mathcal{G}^{X}_{+}$ data is the ‘doubly linear’ behavior observed in Fig. 7 for $\theta>0.4$, which prevents identification of a unique $T_{\mathrm{eff}}$ in this range.

The assumption that $X_{L}^{Q}$ can be used as a proxy for $\hat{\chi}_{L}$ is thus fairly well justified close to the critical period, where $T_{\mathrm{eff}}$ varies over a limited range and the more complicated effects observed at $\theta>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\approx-\gamma_{\mathrm{e}}$. However, because of the first two causes just described, there are small systematic variations in the $\mathcal{G}^{X}_{+}$ data which prevent a statistically significant fit to a pure power-law relationship as a function of $y_{1}$.

In order to clarify the relationship of these scaling results to those in previous work, we also plot in Fig. 11 data of $\mathcal{G}^{|X|}_{+}(y_{1},y_{2})\equiv X_{L}^{|Q|}L^{-\gamma/\nu}\equiv\left(\langle Q^{2}\rangle-\langle\left|Q\right|\rangle^{2}\right)L^{-\gamma/\nu}$ vs $y_{1}$ at $y_{2}=0$, again using the equilibrium values $\gamma_{\mathrm{e}}$ and $\nu_{\mathrm{e}}$ to calculate $\mathcal{G}^{|X|}_{+}$. This can be directly compared to Fig. 11(d) in Ref. Korniss et al. (2001), in which the quantity we call $X_{L}^{|Q|}$ was called $X_{L}^{Q}$. Attempted fits to all five data points and to the last four data points of $\mathcal{G}^{|X|}_{+}$ in Fig. 10 produce nominal scaling exponents $-1.60\pm 0.02$ and $-1.69\pm 0.02$ (both with $q<10^{-15}$). The agreement in Fig. 11(d) of Ref. Korniss et al. (2001) of the line with slope $-7/4$ with the data for $\theta>\theta_{c}$ must therefore be viewed as qualitative. The method used in Ref. Korniss et al. (2001) to numerically estimate $\gamma_{\mathrm{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 $\theta<0.4$ we have found that the FDR in Eq. (15) holds to a very good approximation.

Finally, in Fig. 12, we plot $\mathcal{G}^{X}_{+}$ vs $y_{2}$ at $P=P_{c}$, to study its scaling in the regime $y_{1}\ll y_{2}$. As just noted, the use of $X_{L}^{Q}$ as a proxy for $\hat{\chi}_{L}$ is well justified at $P=P_{c}$ by our results. Power-law scaling is perhaps suggested in the range $20<y_{2}<50$ for $L=180$, and it is clearly obeyed from $y_{2}=8.42$ to $84.2$ for $L=256$. The scaling exponent was determined as $(1-\delta_{\mathrm{d}})/\delta_{\mathrm{d}}=-0.914\pm 0.030$, which is consistent with the corresponding equilibrium value $(1-\delta_{\mathrm{e}})/\delta_{\mathrm{e}}=-14/15\approx-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. Sides et al. (1998, 1999); Korniss et al. (2001). 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’, $H_{b}$ 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 $H_{b}$. In particular, the scaling exponent $\delta_{\rm 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, $\delta_{\rm e}=15$. Furthermore, in agreement with previous results Korniss et al. (2001), the dynamic scaling exponents $\gamma_{\rm d}$, $\beta_{\rm d}$, and $\nu_{\rm d}$ were also found to equal their equilibrium Ising counterparts, $\gamma_{\rm e}=7/8$, $\beta_{\rm e}=1/8$, and $\nu_{\rm d}=1$.

These results further strengthen previous numerical Sides et al. (1998, 1999); Korniss et al. (2001) and analytical Fujisaka et al. (2001); Grinstein et al. (1985); Bassler and Schmittmann (1994) 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. Grinstein et al. (1985); Bassler and Schmittmann (1994), we caution the reader that what is claimed in the present paper (as well as in Ref. Fujisaka et al. (2001)) 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-P_{c}$ and $H_{b}$ play the roles of $T-T_{c}$ 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}^{Q}\equiv L^{2}\left(\langle Q^{2}\rangle-\langle Q\rangle^{2}\right)$ and the susceptibility $\hat{\chi}_{L}$ with a slope we have called $T_{\mathrm{eff}}$, holds for a range of periods above $P_{c}$ and for a range of bias fields around $H_{b}=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 $X_{L}^{Q}$ were used as a proxy for the (then unknown) quantity $\hat{\chi}_{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 reported here. The first is to investigate whether the field $H_{b}$ functions as the conjugate field, with scaling exponents consistent with the equilibrium Ising transition for periods $P<P_{c}$, below the critical period. In the equilibrium system, the study of critical scaling in nonzero field for $T<T_{c}$ is complicated by the long time correlations and strong finite-size effects which accompany the bimodal distributions of magnetization below $T_{c}$. Similar effects would complicate the investigation of scaling with respect to $H_{b}$ in the DPT for $P<P_{c}$, but the advanced techniques Berg and Neuhaus (1992, 1991) which make the equilibrium simulations tractable do not extend obviously to the nonequilibrium case. The second computational project suggested is to determine the nature of the full conjugate field $H_{c}$. The third project would be to study the FDR at different values of the temperature, $T$, and the amplitude, $H_{0}$, 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 $H_{b}$, the FDR found in the critical region, and the scaling of $H_{b}$.

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.