# Simulation of Magnetization Switching in Biaxial Single-Domain Ferromagnetic Particles

## I Introduction

Biaxial media have recently attracted much research interest as potential materials for applications in high-density magnetic recording. This is a consequence of their good thermal stability, high coercivity and coercivity squareness, and low medium noise -. With the current rapid increase in recording density, the size of a magnetic particle used to store a single bit may soon become so small that it can only contain one single magnetic domain. Understanding the dynamics of magnetization switching in individual single-domain particles is therefore essential. Biaxial media may be obtained, for example, by sputtering a Co-based alloy onto a Cr substrate , by fabricating perovskite superlattices consisting of different ferromagnetic metallic oxides , or by epitaxially growing Fe on GaAs(001) . By applying lithographic patterning techniques to precisely control the particle shape, a uniaxial medium may also exhibit biaxial properties due to competition between the magneto-crystalline and shape anisotropies . The two easy axes in biaxial media are often perpendicular to each other, and they may have the same or different magnetic anisotropies. In thin films, these two axes can both be in-plane, or one of them can be out-of-plane, often perpendicular to the plane of the film. New techniques, such as nanolithography and ultra-high resolution scanning microscopies, have recently made the synthesis and experimental observation of isolated single-domain particles possible. For instance, the magnetization switching of individual barium ferrite nanoscale particles has been observed using Magnetic Force Microscopy (MFM), and a maximum switching field was observed for particles of diameter near 55 nm . For other particles, the angular dependence of the switching (or coercive) field has been measured , the switching dynamics investigated , and the switching statistics measured . In Refs.  and , arrays of single-domain Co dots were fabricated, and the magnetization switching was observed by MFM. Depending on the shape of the Co dot, both in-plane and out-of-plane magnetic moments were observed.

A powerful non-perturbative theoretical method to study models of magnetization switching dynamics is Monte Carlo simulation. One standard model used to study uniaxial magnets is the kinetic Ising model . However, the Ising model can only be used as a reasonably realistic representation of extremely anisotropic uniaxial magnets, and a different model is required to capture the unique properties of biaxial materials. In this work we study a 4-state clock model with two perpendicular axes. In the most symmetric case there is no applied field and four states are equivalent. Suzuki has found a mapping between this symmetric 4-state clock model and two superimposed Ising models . When a field is applied, this symmetry is lost and o the results of extensive simulations of the switching dynamics. In order to simplify the analysis of the simulation results, we extend Suzuki’s mapping to non-zero applied field and to the dynamics.

The rest of this paper is organized as follows. In Sec. II we describe the model used in the simulations. In Sec. III we briefly review the aspects of droplet theory used to analyze the numerical results. In Sec. IV and Sec. V we present and analyze our numerical results for lifetime and switching field, respectively. Finally, conclusions are drawn in Sec. VI.

## II Model Description

The 4-state clock model is defined by the microscopic Hamiltonian

$$=J\underset{i,j}{}\stackrel{}{\mu }_i\stackrel{}{\mu }_j\stackrel{}{H}\underset{i}{}\stackrel{}{\mu }_iA\underset{i}{}(\stackrel{}{\mu }_i\stackrel{}{e}_0)^2.$$
(1)   
The unit vector $`\stackrel{}{\mu }_i`$ is the orientation of the spin at site $`i`$, which can be in one of 4 directions: $`\stackrel{}{e}_0`$, $`\stackrel{}{e}_1`$, $`\stackrel{}{e}_2`$, and $`\stackrel{}{e}_3`$, as shown in Fig. 1. The angle between the applied magnetic field $`\stackrel{}{H}`$ and $`\stackrel{}{e}_0`$ is $`\theta `$. The dimensionless system magnetization is given by

$$\stackrel{}{m}=L^2\underset{i}{}\stackrel{}{\mu }_i.$$
(2)   
We consider a two-dimensional square lattice with $`L`$ spins on each side, and periodic boundary conditions are applied in both directions to suppress boundary effects. (Boundary effects for the Ising model are discussed in Ref. .) The lattice constant is set to unity. In Eq. (1), $`_i`$ runs over all $`L^2`$ sites, and $`_{i,j}`$ runs over all nearest-neighbor pairs on the lattice. The first term in Eq. (1) is the energy due to exchange interactions between nearest-neighbor spins, and $`J>0`$ is the ferromagnetic exchange interaction. The second term is the energy, $`L^2\stackrel{}{H}\stackrel{}{m}`$ due to the applied field $`\stackrel{}{H}`$. The third term represents the energy due to a different preference for one of the two axes, $`\stackrel{}{e}_0`$$`\stackrel{}{e}_2`$ and $`\stackrel{}{e}_1`$$`\stackrel{}{e}_3`$. If $`A=0`$, the model has 4-fold symmetry for $`\stackrel{}{H}=0`$, and the statics of this model is equivalent to two decoupled Ising models . In Appendix A, this mapping is extended to non-zero fields and to the dynamics. A nonzero value of $`A`$ destroys this 4-fold symmetry. Therefore, throughout this paper we refer to the model with $`A0`$ as asymmetric. (In the statistical mechanics literature, clock models with $`A0`$ are often referred to as “anisotropic.” However, in this paper we use the terms “anisotropy” and “anisotropic” as they are customarily understood in the literature on magnetic materials.) For $`A>0`$ and $`\stackrel{}{H}=0`$, the $`\stackrel{}{e}_0`$$`\stackrel{}{e}_2`$ axis is favorable energetically, while the $`\stackrel{}{e}_1`$$`\stackrel{}{e}_3`$ axis is preferred for $`A<0`$.

The critical temperature, $`T_\mathrm{c}`$, of the model with $`A=0`$ and $`\stackrel{}{H}=0`$ is known exactly. Its value is $`k_\mathrm{B}T_\mathrm{c}/J=1/\mathrm{ln}(1+\sqrt{2})1.135`$ , exactly that of a two-dimensional square-lattice Ising model with interaction constant $`J/2`$. Here $`k_\mathrm{B}`$ is Boltzmann’s constant. The simulation temperature is chosen to be $`T=0.8T_\mathrm{c}`$. This simulation temperature is high enough to get a reasonable Monte Carlo acceptance rate, while it is low enough to avoid complications near the critical point such as critical slowing-down and finite-size effects . Based on studies of the kinetic Ising model , simulation results at other moderately low temperatures should be qualitatively similar to those obtained here.

The initial configuration of the system is ordered, with $`\stackrel{}{\mu }_i=\stackrel{}{e}_0`$ for all spins. This may be achieved by applying a strong field in the $`\stackrel{}{e}_0`$ direction for time $`t<0`$. Then at $`t=0`$ the external magnetic field is changed to $`\stackrel{}{H}`$ (Fig. 1). During the simulation, first one spin is selected at random, then a new candidate state is chosen. At low temperatures, large orientational changes of the spin ($`\stackrel{}{e}_0\stackrel{}{e}_2`$ and $`\stackrel{}{e}_1\stackrel{}{e}_3`$) are very unlikely . Consequently, the candidate state is obtained by rotating the spin by $`90^{}`$ or $`90^{}`$. The a priori probabilities of clockwise and anti-clockwise rotations are each set to $`1/2`$.

The probability that the spin flips to the chosen state is given by Glauber dynamics 

$$W=\frac{\mathrm{exp}(\beta \mathrm{\Delta }E)}{1+\mathrm{exp}(\beta \mathrm{\Delta }E)},$$
(3)   
where $`\mathrm{\Delta }E`$ gives the change in the energy of the system that would result if the spin flip were accepted, and $`\beta =1/k_\mathrm{B}T`$. The Glauber dynamic is chosen because it can be derived from quantum mechanics under certain stringent conditions . Here time is reported in the unit of one Monte Carlo steps per spin (MCSS).

The lifetime is defined as the average time, $`\tau `$, needed to satisfy a specific stopping criterion (discussed below), where $`\tau `$ is the first-passage time to the stopping criterion for an individual Monte Carlo run. The average is taken over at least 1000 independent Monte Carlo runs in order to obtain good statistics. In addition to the mean, we also measure the relative standard deviation

$$r=\sqrt{\tau ^2\tau ^2}/\tau .$$
(4)   
The switching field is measured by observing whether on average the system magnetization has reversed after a fixed waiting time, here chosen as $`t_\mathrm{w}=200`$ MCSS. As Ising model studies have demonstrated , the switching fields for different waiting times are qualitatively similar.

To measure the lifetime of a metastable state, a criterion for escape from the metastable state must be chosen. In studies of the kinetic Ising model, the usual criterion is that an escape is registered and the simulation stopped the first time the magnetization reaches a fixed cutoff value, $`m_{\mathrm{stop}}`$. (In other words: the stopping time is the first-passage time to $`m_{\mathrm{stop}}`$.) Results for different values of $`m_{\mathrm{stop}}`$ are significantly different only in the deterministic regime, and even in this regime the results are qualitatively similar .

In some experimental studies of magnetization switching, the sign of the magnetization along a particular direction is much easier to measure than its magnitude . The stopping criterion most directly relevant to these experiments is therefore to use $`m_{\mathrm{stop}}=0`$. For the Ising model, this is obtained by stopping the simulation when half of the spins have entered the stable state (and half have left the metastable state). For other experimental situations, a different criterion might be easier to measure. However the stopping criterion chosen should not have a large effect on the results.

In order to completely map the the 4-state clock model onto two superimposed Ising systems (see Appendix A), the stopping criterion for the 4-state clock model must map onto the stopping criteria for the two individual Ising systems. If the stopping criteria for both Ising systems are given by $`m_{\mathrm{stop}}=0`$, two different stopping criteria for the 4-state clock model can be naturally envisioned. Denoting these by $`S_{}`$ and $`S_{}`$, they are given in terms of the first passage to the stopping parameter $`\mu _{\mathrm{stop}}`$ by

$`S_{}\mathrm{has}`$ $`\mu _{\mathrm{stop}}=m_{I0}m_{I1}`$ (5)   
$`S_{}\mathrm{has}`$ $`\mu _{\mathrm{stop}}=1\eta (m_{I0})\eta (m_{I1}),`$ (6)   
where each $`m_{I\alpha }`$ for $`\alpha =0,1`$ is the time-dependent magnetization of one of the Ising systems and $`\eta (x)`$ is the unit step function. The criterion $`S_{}`$ is to stop when the first Ising system reaches zero magnetization, while for $`S_{}`$ the simulation is continued until the last of the two Ising systems reaches zero magnetization.

The mapping between the clock model and the two Ising models illustrates that in order for a clock spin to escape from the initial state, only one of the two Ising spins on the site needs to flip, while both Ising spins on the same site must flip in order for the clock spin to enter the stable state. Therefore, for the clock model the $`S_{}`$ stopping criterion implies that $`\frac{1}{2}`$ of the clock spins have escaped from the initial state. The $`S_{}`$ stopping criterion means a certain fraction ($`\frac{1}{2}`$) of the clock spins have entered the stable state. For the symmetric clock model, the two Ising models are decoupled, and their switching processes are independent of each other. Therefore, for $`S_{}`$ on average $`1/4`$ of the clock spins will have entered the stable state at the time the simulation is stopped. In this paper, we use only the $`S_{}`$ stopping criterion since it more closely models an experimental stopping criterion of zero magnetization. For the asymmetric model, we also use this stopping criterion in order to observe the effects of $`A`$ and $`\stackrel{}{H}`$ consistently. In this case the average fraction of clock spins in the stable state at the end of the simulation is between $`1/4`$ and $`1/2`$, depending on $`A`$. A more detailed discussion of the quantitative differences between different stopping criteria is given in Ref. .

## III Droplet Theory

During magnetization switching in single-domain particles, small “droplets” of the stable phase are continually created and destroyed by thermal fluctuations . The fate of a droplet is determined by the competition between its total surface free energy, which is proportional to its surface area, and the coupling of its magnetization with the applied field, which is proportional to its volume. If the droplet is small, further growth produces an increase of the free energy due to the surface tension. Therefore, the droplet will be destroyed with high probability. If the droplet is sufficiently large, this free-energy penalty is offset by the benefit obtained from increasing the number of spins parallel to the applied magnetic field. Therefore the droplet is very likely to grow further. For a droplet with a critical radius $`R_c`$, the tendencies towards growth and shrinkage balance, and the droplet will grow or shrink with equal probability. By comparing $`R_c`$ with the system size $`L`$, the lattice spacing $`a`$, and the average distance between supercritical droplets $`R_0`$, the magnetization switching dynamics is divided into four regimes :

* The Coexistence (CE) regime $`[a<L<R_c]`$, in which the switching is governed by subcritical fluctuations on the scale of the system size.
* The Single-droplet (SD) regime $`[a<R_c<L<R_0]`$, characterized by switching via a single critical droplet. Once a supercritical droplet is nucleated, it will most likely grow to occupy the entire system before another droplet is nucleated.
* The Multi-droplet (MD) regime $`[a<R_c<R_0<L]`$, characterized by switching via a finite density of critical droplets.
* The Strong-field (SF) regime $`[aR_cR_0<L]`$, in which the concept of a droplet is no longer applicable.

Since the symmetric 4-state clock model is equivalent to two decoupled Ising models, we can simplify the analysis of the clock model based on extensive studies of kinetic Ising models . Here a brief review of the droplet theory for the Ising model is given. For more details, see the references listed above.

In the SD regime, the lifetime is given by

$$\tau _\mathrm{I}[L^d\mathrm{\Gamma }(T,H_\mathrm{I})]^1,$$
(7)   
where $`H_\mathrm{I}`$ is the field in the Ising Hamiltonian, and the nucleation rate per unit volume $`\mathrm{\Gamma }(T,H_\mathrm{I})`$ is 

$$\mathrm{\Gamma }(T,H_\mathrm{I})A(T)|H_\mathrm{I}|^K\mathrm{exp}\left[\frac{\beta \mathrm{\Xi }(T)}{|H_\mathrm{I}|^{d1}}(1+O(H_\mathrm{I}^2))\right].$$
(8)   
Here $`A(T)`$ is a non-universal prefactor, and $`K`$ is believed to be 3 for two dimensional ($`d=2`$) Ising systems . At $`T=0.8T_c`$, $`\mathrm{\Xi }(T)/J_\mathrm{I}^20.92`$ , where $`J_\mathrm{I}`$ is the exchange interaction in the Ising Hamiltonian.

In the MD regime, the system can be approximately partitioned into $`(L/R_0)^d1`$ cells, each of which decays via an independent Poisson process with rate $`R_0^d\mathrm{\Gamma }`$. Therefore, the volume fraction is self-averaging, and the average lifetime $`\tau _\mathrm{I}`$, defined as the average first-passage time to zero magnetization, is given by the well-known Avrami theory of metastable decay as 

$$\tau _\mathrm{I}\left[|H_\mathrm{I}|^d\mathrm{\Gamma }(T,H_\mathrm{I})\right]^{1/(d+1)}.$$
(9)   
By taking the derivative of the logarithm of $`\tau _\mathrm{I}`$ with respect to $`|J_\mathrm{I}/H_\mathrm{I}|^{d1}`$, one obtains

$$\frac{\mathrm{\Lambda }_{\mathrm{Ieff}}}{J_\mathrm{I}^{d1}}\frac{d\mathrm{ln}\tau _\mathrm{I}}{d(|J_\mathrm{I}/H_\mathrm{I}|^{d1})}=\lambda |H_\mathrm{I}/J_\mathrm{I}|^{d1}+\frac{\mathrm{\Lambda }}{J_\mathrm{I}^{d1}},$$
(10)   
with $`\lambda =K/(d1)`$ and $`\mathrm{\Lambda }=\beta \mathrm{\Xi }(T)`$ in the SD regime, while $`\lambda =(K+d)/(d^21)`$ and $`\mathrm{\Lambda }=\beta \mathrm{\Xi }(T)/(d+1)`$ in the MD regime .

For the symmetric 4-state clock model, if the applied field is along the $`\stackrel{}{e}_2`$ direction, the applied field $`H_\mathrm{I}`$ and the exchange interaction $`J_\mathrm{I}`$ of the Ising model are related to the corresponding parameters for the clock model by $`H_\mathrm{I}=H/2`$ and $`J_\mathrm{I}=J/2`$ (Table I). By replacing the parameters for the Ising model by the corresponding ones for the clock model in the above analysis, we get the dimensionless value $`\mathrm{\Lambda }_{\mathrm{eff}}/J^{d1}`$ for the clock model, which is

$$\frac{\mathrm{\Lambda }_{\mathrm{eff}}}{J^{d1}}=\frac{d\mathrm{ln}\tau }{d(|J/H|^{d1})}=\lambda |H/J|^{d1}+\frac{\mathrm{\Lambda }}{(J/2)^{d1}},$$
(11)   
Since $`J_\mathrm{I}=1/2J`$, comparing Eq. (10) and Eq. (11), we have

$$\mathrm{\Lambda }_{\mathrm{eff}}/J^{d1}=\mathrm{\Lambda }_{\mathrm{Ieff}}/J_\mathrm{I}^{d1}.$$
(12)   
Therefore we can compare the simulation results for the symmetric 4-state clock model directly with the theoretical predictions for the Ising models. For the asymmetric clock model, we calculate the same derivative.

## IV Lifetime Measurements

The statistical properties of the lifetime depend on the switching regime of the system. For an Ising system in the SD regime, the probability density of the lifetime is exponential 

$$p_\mathrm{I}(\tau )=\tau _\mathrm{I}^1\mathrm{exp}(\tau /\tau _\mathrm{I}),$$
(13)   
where $`\tau _\mathrm{I}`$ is the mean lifetime. As discussed in Sec. II, for the 4-state clock model, magnetization switching, as defined by the stopping criterion used in this paper, implies that both constituent Ising lattices must have reached zero magnetization. Therefore, the lifetime of the clock model is the longer of the lifetimes of the two Ising models. If the applied field is in the $`\stackrel{}{e}_2`$ direction, the two Ising systems are identical, and the probability density of the lifetime for the clock model can be simply written as

$$p_\mathrm{c}(\tau )=\tau _\mathrm{I}^1\mathrm{exp}(\tau /2\tau _\mathrm{I})[1\mathrm{exp}(\tau /2\tau _\mathrm{I})].$$
(14)   
The mean and relative standard deviation of this distribution are $`\tau =3\tau _\mathrm{I}`$ and $`r=\sqrt{5}/3`$, respectively. The derivation of $`p_\mathrm{c}(\tau )`$ for the general case that $`\stackrel{}{H}`$ is not parallel to $`\stackrel{}{e}_2`$ is detailed in Appendix B with the general result given as Eq. (27).

In the MD regime, since the switching process of the Ising model is deterministic, based on the mapping between the Ising model and the 4-state clock model, the switching dynamics of the 4-state clock model is also deterministic. By “deterministic” we mean that the lifetime distribution is characterized by a very small relative standard deviation, $`r1`$. Figure 2 and Fig. 3 show $`\tau `$ and $`r`$ as functions of the inverse applied field, $`J/|H|`$, for different lattice sizes. In these measurements, the applied field $`\stackrel{}{H}`$ is always parallel to $`\stackrel{}{e_2}`$. In the panels showing $`r`$, the horizontal solid and dashed lines correspond to $`r=1`$ and $`r=\sqrt{5}/3`$, respectively. As predicted, for the isotropic clock model, $`r`$ approaches $`\sqrt{5}/3`$ as $`|H|`$ decreases. Also shown in these figures are results for the asymmetric clock model. When the easy axis is along $`\stackrel{}{e_0}`$$`\stackrel{}{e}_2`$ ($`A>0`$), the lifetime is longer than that of the symmetric model with the same applied field because a stronger field is needed to flip the spins from the metastable state $`\stackrel{}{e}_0`$ to the intermediate states, $`\stackrel{}{e}_1`$ and $`\stackrel{}{e}_3`$. When the easy axis is along $`\stackrel{}{e_1}`$$`\stackrel{}{e}_3`$ ($`A<0`$), however, the lifetime is shorter than that for the symmetric model in the MD regime, and longer in the SD regime. The lifetime distribution is wider in the SD regime, and $`r`$ has a system size-dependent maximum for $`A<0`$. In either case, as $`|H|`$ decreases, the asymmetry dominates, and the two states along the easy axis are strongly favored over the other two states. Therefore for weak $`\stackrel{}{H}`$ we expect the asymmetric clock model to behave essentially like a single Ising model, so that $`r`$ will approach unity.

In order to understand the switching process for the asymmetric model better, we show selected snapshots of the model during the simulations in Figs. 4-6. The darkest regions represent the initial state, the lightest ones represent the stable state, and the two intermediate shades of gray represent the intermediate states. The applied field is chosen so that $`r`$ is near its maximum for $`A<0`$ (Fig.2). For $`A>0`$ (Fig. 4), there is only one channel for the model to switch the magnetization: first a droplet of one intermediate state is nucleated and grows quickly. Then a droplet of the stable state nucleates and grows on this intermediate background. This process is similar to single-droplet decay in the Ising model, except that there is no intermediate phase in the Ising model. For $`A<0`$, however, two channels exist. In the slow channel (Fig. 5), the system first enters a uniform intermediate state, from which it takes a long time to nucleate a stable droplet. In the fast channel (Fig. 6), droplets in different intermediate states are formed simultaneously. As they meet, a stable droplet is quickly nucleated on their common boundary to decrease the free energy. Because of the existence of these two channels, the relative standard deviation of the lifetime becomes very large, exceeding unity, which is the maximum value for a single Ising model.

Figure 7 shows the derivatives of the logarithm of the lifetime with respect to $`J/|H|`$, $`\mathrm{\Lambda }_{\mathrm{eff}}`$ of Eq. (10), for different system sizes. In this figure, the straight lines are the theoretical results obtained from Eq. (10). Indeed, for the symmetric clock model, our numerical results show good agreement with the theory. Note that there are no adjustable parameters used in obtaining the theoretical results. In the asymmetric cases, when the applied field is strong the effect of the asymmetry is small. For weak fields, the asymmetry dominates the switching process, and the results for the clock model cannot be mapped onto those for the Ising model. Unfortunately, we could not extend the simulations deep into this regime, due to the prohibitive computer resources required: it takes at least two months to simulate 1000 individual realizations for $`L=100`$ and $`J/|H|>12.5`$ (Fig. 2) on a single node of an IBM sp2 computer.

## V Switching-Field Measurements

We define the switching field as follows. We apply a field $`\stackrel{}{H}`$ to the system in the initial state for a chosen waiting time $`t_\mathrm{w}=200`$ MCSS. If the magnetization has switched $`N/2`$ times in $`N`$ simulations, then the magnitude of the applied field is called the switching field, $`H_{\mathrm{sw}}(t_\mathrm{w})`$. We choose this definition in order to simulate qualitatively the experiment discussed in Ref. . The relationship between the switching field and the coercivity is that in coercivity measurements, the applied field is changed sinusoidally or linearly , while in our simulation the applied field is reversed abruptly.

First we compare the switching field for the symmetric clock model and the Ising model. The switching field as a function of the system size $`L`$ is illustrated in Fig. 8. Since the time scale of the clock model is twice that for the Ising model, the switching field for the Ising model is measured with waiting time $`t_\mathrm{w}=100`$ MCSS. In the MD regime, the lifetime for the symmetric clock model is twice that for the Ising model, thus the switching fields for the clock model and the Ising model are the same. In the stochastic regime, from Eq. (13) and Eq. (14), the lifetime for the clock model is three times that of the Ising model. Therefore, the switching field for the clock model is larger than that for the Ising model. Here we emphasize that the maximum shown in Fig. 8 is the result of different switching mechanisms dominating the decay in the CE, SD and MD regimes. It is not the effect of demagnetization since there are no dipole-dipole interactions in our Hamiltonian. For more detailed discussions of the dependence of $`H_{\mathrm{sw}}`$ on $`L`$, $`t_\mathrm{w}`$ and $`T`$ in the Ising model, see Ref. .

Next we consider the effects of the asymmetry $`A`$ on the switching field. The results are shown in Fig. 9. Again, $`\stackrel{}{H}`$ is in the $`\stackrel{}{e_2}`$ direction. When the easy axis is along $`\stackrel{}{e}_0`$$`\stackrel{}{e}_2`$ ($`A>0`$), increasing $`A`$ leads to a monotonic increase of the switching field. When $`A<0`$, on the other hand, making $`A`$ more negative decreases the switching field in the MD regime. However, in the SD regime different $`A`$ shows little effect on the switching field.

The effect of the direction of the applied field on the switching field is shown in Fig. 10. Here $`H_{\mathrm{sw}}`$ is shown as a function of $`L`$ for different $`\theta `$ at $`A=0`$. The symmetric clock model is equivalent to two decoupled Ising models, and the field components acting on these two Ising systems change with $`\theta `$. Thus the switching behaviors of the two Ising systems are different, and the switching of the clock model is determined by the slower one. Therefore $`H_{\mathrm{sw}}`$ increases as $`\stackrel{}{H}`$ deviates further away from the $`\stackrel{}{e_2}`$ direction.

Finally we have measured the $`\theta `$-dependence of $`H_{\mathrm{sw}}`$ for two system sizes ($`L=8`$ and $`L=40`$) with different values of $`A`$. The results are shown in Fig. 11. When the direction of $`\stackrel{}{H}`$ is close to $`\stackrel{}{e}_2`$, the energetic preference for the two intermediate states is similar and mainly determined by the asymmetry. Therefore $`H_{\mathrm{sw}}`$ is governed by $`A`$. As $`\theta `$ approaches $`135^{}`$, $`\stackrel{}{H}`$ increases the stability of one intermediate state and decreases that of the stable state. When $`A<0`$, this tendency is further amplified, therefore the switching field increases. For $`A>0`$, $`\stackrel{}{H}`$ makes the flip from $`\stackrel{}{e}_0`$ to the favored intermediate state easier, whereafter $`A`$ acts like an applied field in flipping the spin to $`\stackrel{}{e}_2`$. Thus $`H_{\mathrm{sw}}`$ decreases compared with the isotropic model.

## VI Conclusion

In this work we have studied magnetization switching in biaxial magnetic media by large-scale Monte Carlo simulation and theoretical analysis of a kinetic 4-state clock model. The model reduces the infinite anisotropy in the kinetic Ising model by introducing two intermediate magnetization directions between the magnetization states parallel and antiparallel to the applied field. If the two perpendicular axes are equivalent, the switching process is equivalent to that of two decoupled Ising models. With asymmetry, the clock model shows a more complicated switching behavior. When the axis along the initial state is favored ($`A>0`$), the model becomes Ising-like. On the other hand, if the other easy axis is favored, the system can switch its magnetization via two different channels, which can lead to a wider distribution of the lifetimes in the SD regime. Another factor affecting the switching process is the direction of the applied field. We studied the direction dependence of the switching field in the isotropic model. With the applied field deviating further from the stable state, a stronger switching field is needed. The different switching mechanisms in the CE, SD, and MD regimes lead to a maximum switching field for a specific system size, in agreement with experimental results. We conclude that the 4-state clock model is a simple but informative platform to study the magnetization switching in biaxial magnetic particles.

## Acknowledgment

The authors wish to thank Dr. Howard L. Richards for proposing this study. We also thank Dr. S. M. Durbin, Dr. S. W. Sides, Dr. G. Korniss, Dr. H. L. Richards, Dr. M. Kolesik and Dr. G. Brown for useful comments. This work is supported in part by the US National Science Foundation Grant No. DMR-9520325, and by Florida State University through the Supercomputer Computations Research Institute (US Department of Energy Contract No. DE-FC05-85ER25000) and the Center for Materials Research and Technology.

## Appendix A

In this appendix, the dynamical mapping between the Ising model and the symmetric 4-state clock model is derived.

First we define two unit vectors, $`\stackrel{}{s}_0`$ and $`\stackrel{}{s}_1`$, such that the angle between $`\stackrel{}{s}_0`$ ($`\stackrel{}{s}_1`$) and $`\stackrel{}{e}_0`$ is $`45^{}`$ ($`+45^{}`$) (Fig. 12). Let $`\stackrel{}{\mu }_i`$ denote the direction of the clock spin at site $`i`$, which can be $`\stackrel{}{e}_0`$, $`\stackrel{}{e}_1`$, $`\stackrel{}{e}_2`$, or $`\stackrel{}{e}_3`$. The unit vector $`\stackrel{}{\mu }_i`$ can be expressed in terms of $`\stackrel{}{s}_0`$ and $`\stackrel{}{s}_1`$ as

$$\stackrel{}{\mu }_i=\frac{1}{\sqrt{2}}(\omega _{0i}\stackrel{}{s}_0+\omega _{1i}\stackrel{}{s}_1),$$
(15)   
where $`\omega _{0i}`$ and $`\omega _{1i}`$ can take the values $`\pm 1`$. The applied field can also be decomposed as

$$\stackrel{}{H}=H_0\stackrel{}{s}_0+H_1\stackrel{}{s}_1,$$
(16)   
where

$`H_0`$ $`=`$ $`|\stackrel{}{H}|\mathrm{cos}(45^{}+\theta ),`$ (17)   
$`H_1`$ $`=`$ $`|\stackrel{}{H}|\mathrm{sin}(45^{}+\theta ).`$ (18)   
Substituting Eqs. (15)-(17) into Eq. (1), we get

$``$ $`=`$ $`_0+_1,`$ (19)   
$`_0`$ $`=`$ $`{\displaystyle \frac{J}{2}}{\displaystyle \underset{i,j}{}}\omega _{0i}\omega _{0j}{\displaystyle \frac{H_0}{\sqrt{2}}}{\displaystyle \underset{i}{}}\omega _{0i},`$ (20)   
$`_1`$ $`=`$ $`{\displaystyle \frac{J}{2}}{\displaystyle \underset{i,j}{}}\omega _{1i}\omega _{1j}{\displaystyle \frac{H_1}{\sqrt{2}}}{\displaystyle \underset{i}{}}\omega _{1i}.`$ (21)   
From these equations, it is clear that the Hamiltonian of the symmetric 4-state clock model is equivalent to the sum of two decoupled Ising Hamiltonians with exchange interactions $`J/2`$. The applied field in each Ising Hamiltonian equals the field component times $`1/\sqrt{2}`$. Therefore, the 4-state clock model in equilibrium can be treated as two superimposed Ising systems, with two Ising spins at each site .

The next step is to construct a mapping of the spin-flip dynamics between the 4-state clock model and the Ising models. Since each spin can only flip to its neighboring state, this can be done by choosing one of the two Ising systems randomly and flipping the corresponding Ising spin. Because each time only one Ising lattice is selected, the change of the energy due to the spin flip is obtained by updating either $`_0`$ or $`_1`$ in Eq. (19).

With the above analysis, it follows that the dynamics of the symmetric 4-state clock model is equivalent to that of two decoupled Ising models. Not only are the spin-flip mechanisms interchangeable, but also only one Ising Hamiltonian in Eq. (19) is affected in each flip because $`_0`$ and $`_1`$ are decoupled. Since each clock spin can be thought to be composed of two Ising spins, during the simulation the time scale for the 4-state clock model is twice that for the Ising model. The mapping of the parameters between the symmetric 4-state clock model and Ising models is summarized in Table I. In this table, $`t`$ is the time measured in Monte Carlo steps per spin (MCSS).

## Appendix B

In this appendix, the lifetime distribution for the isotropic 4-state clock model is derived. We denote the lifetimes of the two Ising models by $`\tau _0`$ and $`\tau _1`$, and the lifetime of the 4-state clock model by $`\tau _\mathrm{c}`$. With the stopping criterion used in this work, we have

$$\tau _\mathrm{c}=2\mathrm{max}(\tau _0,\tau _1).$$
(22)   
Thus, for a given value $`\tau `$ of the lifetime, the probability that $`\tau _\mathrm{c}<\tau `$ can be written as

$$P(\tau _\mathrm{c}<\tau )=P(\mathrm{max}(\tau _0,\tau _1)<\tau /2).$$
(23)   
Here the factor 2 comes from the mapping of the time scale between the clock and the Ising models given in Table I. Since the two Ising models are decoupled, the lifetime distributions are independent. Therefore

$$P(\tau _\mathrm{c}<\tau )=P(\tau _0<\tau /2)P(\tau _1<\tau /2).$$
(24)   
In the SD regime, the lifetime of each Ising model is exponentially distributed,

$$P(\tau _k<\tau )=1\mathrm{exp}(\tau /\tau _k),$$
(25)   
where $`k=0,1`$. Substituting Eq. (25) into Eq. (24) and differentiating with respect to $`\tau `$, we obtain the probability density of the lifetime of the clock model as

$`p_\mathrm{c}(\tau )`$ $`=`$ $`{\displaystyle \frac{\mathrm{exp}(\tau /2\tau _0)\left[1\mathrm{exp}(\tau /2\tau _1)\right]}{2\tau _0}}`$ (27)   
$`+{\displaystyle \frac{\mathrm{exp}(\tau /2\tau _1)\left[1\mathrm{exp}(\tau /2\tau _0)\right]}{2\tau _1}}.`$   
If $`\tau _0=\tau _1=\tau _\mathrm{I}`$, we get Eq. (14).

In the MD regime, since the magnetic switching is deterministic, the lifetime probability density of the Ising model is well approximated by a $`\delta `$-function. Thus for the Ising model,

$$P(\tau _\mathrm{I}<\tau )\eta (\tau \tau _\mathrm{I}),$$
(28)   
where $`\eta `$ is the unit step function. Following the same procedure as above we get the lifetime distribution for the clock model,

$$p_\mathrm{c}(\tau )\delta (\tau 2\mathrm{max}[\tau _0,\tau _1]).$$
(29)   
Finally we consider the case that the field is applied in such a direction that one Ising system is in the MD regime with average lifetime $`\tau _{\mathrm{MD}}`$, and the other in the SD region with the average lifetime $`\tau _{\mathrm{SD}}`$. Starting from Eq. (24), we then have

$$P_\mathrm{c}(\tau _\mathrm{c}<\tau )=\eta (\tau 2\tau _{\mathrm{MD}})[1\mathrm{exp}(\tau /2\tau _{\mathrm{SD}})],$$
(30)   
which yields

$`p_\mathrm{c}(\tau )`$ $`=`$ $`[1\mathrm{exp}(\tau _{\mathrm{MD}}/2\tau _{\mathrm{SD}})]\delta (\tau 2\tau _{\mathrm{MD}})`$ (32)   
$`+{\displaystyle \frac{\mathrm{exp}(\tau /2\tau _{\mathrm{SD}})}{2\tau _{\mathrm{SD}}}}\eta (\tau 2\tau _{\mathrm{MD}}).`$   
A more detailed calculation for the MD cases would replace the $`\delta `$-functions by Gaussians and the step functions by the corresponding error functions , with no essential change from the approximate results given here.