{'fraction_non_alphanumeric': 0.02972764492968611, 'fraction_numerical': 0.009968551593188156, 'mean_word_length': 4.24155090192826, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hysteresis in Random Field XY and Heisenberg Models: Mean Field Theory and Simulations at Zero Temperature

We examine zero temperature hysteresis in random field XY and Heisenberg models in the zero frequency limit of a cyclic driving field. Exact expressions for hysteresis loops are obtained in the mean field approximation. These show rather unusual features. We also perform simulations of the two models on a simple cubic lattice and compare them with the predictions of the mean field theory.


I. INTRODUCTION
Random field XY and Heisenberg models provide a simple framework for exploring the effects of quenched disorder in classical systems of continuous symmetry [1,2]. These models and their variants have helped in understanding a wide range of phenomena including random pinning of spin and charge density waves in metals [3][4][5], vortex lattices in disordered type-II superconductors [6,7], liquid crystals in porous media [8][9][10][11], and disordered ferromagnets [12][13][14][15][16][17]. Initially the models were used to understand the effect of disorder on equilibrium properties of materials. However, with shifting trends in statistical mechanics towards non-equilibrium phenomena, the same models have been supplemented with a simple relaxation dynamics and adapted to study non-equilibrium behavior of systems including their response to a driving field. In the present paper we examine zero temperature hysteresis in these models when the frequency of the cyclic driving field goes to zero, and provide an exact solution of the hysteresis loop in the mean field limit. The zero temperature dynamics is deterministic and therefore simpler to analyze theoretically. But this is not the only reason for using it. It is also meaningful for describing disorder-driven hysteresis in real materials at a finite temperature. Materials with quenched disorder are characterized by a large number of metastable states separated from each other by energy barriers that are much larger than the thermal energy of the system. These findings are based on extensive studies of spinglasses and other random field systems [2] but are intuitive as well. Intuition tells us that if the disorder remains frozen over experimental time scales, thermal energy must be smaller than the barriers due to disorder. Similarly metastable states must be aplenty because most of these have an apparently random configuration. A random configuration does not bring to mind any specific configuration but rather a large number of possible configurations.
Zero temperature hysteresis in XY and Heisenberg models in the zero frequency limit of driving field has been studied by Rava da Silveira and Mehran Kardar [13] as well. They use a slightly different variant of the model than the one studied here. The random field in their model has a Gaussian distribution centered at zero. In our model, the random fields are in the form of randomly oriented unit vectors. We determine the hysteretic response of the system to a changing field by solving the equations of motion directly for a given initial condition. Silveira and Kardar take an indirect approach. They recast the equations of motion into a path integral. The path integral is a sum over all paths of the exponential of an action. It includes paths corresponding to different initial conditions. Silveira and Kardar employ a method for extracting the physically relevant hysteretic path from amongst multiple solutions. We refer the reader to reference [13] for details. The main object of their study is to examine critical points in the hysteretic response of a system. They focus on a point on the hysteresis loop where the susceptibility of the system diverges. If there is such a point on one half of the hysteresis loop, say in increasing applied field, there is also a symmetrically placed point on the other half of the loop corresponding to decreasing field. These points are called non-equilibrium critical points because they are characterized by a diverging correlation length, and show scaling of various quantities and universality of critical exponents that is reminiscent of equilibrium critical point phenomena. Sethna et al [18][19][20] studied the non-equilibrium critical points on the hysteresis loop in the random field Ising model with a Gaussian distribution of random fields. Silveira and Kardar [13] ask the question if the universality class of critical hysteresis studied by Sethna et al would change if we go from Ising spins to vector spins. They find no change in the case when the critical point occurs at a nonzero value of either the applied field or the magnetization. However, if the critical point were to occur when both the applied field and the magnetization vanish, all components of the order parameter may become critical simultaneously. In this case the critical point would have full rotational symmetry of vector spins with a new set of critical exponents.
We focus on the shape of hysteresis loop rather than the critical points on it. The shape of hysteresis is not a universal object like a set of critical exponents, but nonetheless it is of practical interest. An exact calculation of hysteresis loop also determines if there are first order jumps or critical points on the loop. The calculations presented here bring out two rather unexpected but interesting results. In the random field XY model, there is a window in the value of the ferromagnetic coupling parameter where the hysteresis loop splits into two small loops at large values of the cyclic field but there is no hysteresis at small values of the field. This prediction of the mean field theory is also seen qualitatively in our simulations of the model on simple cubic lattices with nearest neighbor interactions. Similar shapes have been observed earlier in the random field Blume-Emery-Griffiths model for martensitic transitions [21], and other theoretical models and experiments [12,[15][16][17] but they do not appear to be known very widely. The other point is that our mean field theory predicts a different kind of phase transition in random field XY model than in random field Heisenberg model. This is somewhat surprising at first sight because the two models have the same critical behavior in the mean field limit of the renormalization group theory. However, it is understandable if we keep in mind that our model has a different distribution of the random field than the one used in reference [13]. We shall return to this issue after presenting our results.
The present paper is organized as follows. In section II, we explain the model based on n-component unit vector spins, and zero temperature dynamics. A discrete-time equation of motion for the magnetization along the applied field is also set up in this section. Section III is devoted to the analysis of the shapes of hysteresis curves. It has four sub-sections. Sub-sections A and B are devoted to the determination of hysteresis loops in the mean field theory for XY and Heisenberg models respectively. Sub-section C takes a closer look at the nature of criticality in the mean field theory of hysteresis. Sub-section D presents the results of numerical simulations of the model on simple cubic lattices and comparison of these results with the predictions of the mean field theory. Section IV contains some concluding remarks.

We consider the Hamiltonian
Here S i and h i are n-component unit vectors located at site-i (i = 1, 2, . . . N ) of a d-dimensional lattice. In the context of magnetic systems, { S i } are classical spins, { h i } a set of on-site random fields, and h is a uniform applied field of magnitude |h|. We focus on n = 2 (XY spins), and n = 3 (Heisenberg spins) since the case n = 1 (random field Ising model) has been studied rather thoroughly [19] albeit for a Gaussian distribution of the random field. The summation over j on the right-hand-side (rhs) is restricted over the nearest neighbors of site-i. The first and the third terms on the rhs promote uniform order in the system: J (J > 0) is ferromagnetic exchange interaction that aligns nearest neighbors parallel to each other; the external field aligns each spin S i along h. The second term on the rhs disorders the system by attempting to align each spin S i in a random direction h i . The random fields { h i ; | h i | = 1} are quenched, i.e. they do not evolve in time. The spins { S i (t)} are the dynamical degrees of freedom.
At zero temperature, an initial configuration { S i (0)} evolves in time so as to lower the energy of the system. The evolution ends when each S i (t) is aligned along the local field f i (t) at that site. Let { S * i } denote a configuration at the termination of the zero-temperature single-spin-flip dynamics. We call it a fixed point configuration because it remains unchanged under the dynamics. In the absence of disorder, the fixed point has all spins parallel to each other irrespective of the starting point { S i (0)}. This corresponds to the lowest energy of the system. In the presence of random fields { h i }, the fixed point becomes rather non trivial on two accounts. First, it may and generically does lose its translational symmetry. Second, it is no longer independent of the starting point. There is now a large set of fixed points each with its domain of attraction from where it can be reached. Each of these fixed points is a local minimum of energy. A local minimum is a stable state at zero temperature because there is no mechanism of escape from it unless the applied field is jacked up sufficiently. It would correspond to a metastable state under finite temperature dynamics if the thermal energy is smaller than the barriers of disorder, but we consider zero temperature dynamics only. In equilibrium problems with quenched disorder, one needs to know the lowest of the local minima. This is a difficult task analytically or computationally. Fortunately, the problem of hysteresis does not require the knowledge of the global minimum. Hysteresis is determined by the sequence of local minima visited by the system as it tries to follow a changing field. Our object is to determine this sequence as the applied field is cycled from −∞ to +∞ and back to −∞ in small steps. At each step, we allow the zero temperature dynamics as much time as it requires to come to a fixed point.
We obtain the local minima by using a discrete time dynamics that progressively lowers the energy of the system. The dynamics transforms a spin configuration { S i (t)} at time t into a lower energy configuration { S i (t + 1)} at time t + 1. The fixed point of this iterative procedure corresponds to a local minimum of the energy of the system. Our dynamics can be stated in a simple form if we rewrite equation (1) in terms of a local effective field f i (t) at site-i, The dynamics is given by the equation: At each site, a new spin S i (t + 1) is obtained that points in the direction of the local field f i (t) at that site. The denominator in equation (3) ensures that the new spin S i (t + 1) has unit length; S i (t + 1) is therefore a rotated form of S i (t). The rotation lowers the energy of each spin, and therefore that of the entire system. However, after the spins are rotated the local field changes as well. Thus the rotated spin S i (t + 1) is generally not aligned along the new local field f i (t + 1) at site i. We can reduce the energy of the system further by repeating the dynamics. Indeed, we start with a random initial configuration { S i (0)} and subject it to repeated applications of equation (3) until an attractor of the dynamics is reached. In discrete time dynamics, an attractor could be in principle a fixed point or a limit cycle. However in our analysis as well as simulations, we find the dynamics always reaches a fixed point configuration { S * i }. The fixed point configuration corresponds to a local minimum of energy. The initial configuration { S i (0)} and the configurations along the path to the fixed point lie in the domain of attraction of the fixed point.
For simplicity, we characterize each configuration of spins by a single parameter that measures the magnetization of the system along the applied field h. We assume that the applied field h is along the x-axis. The equations of motion for the magnetization of Heisenberg and XY spins are quite similar. We first consider the case of XY spins. In this case, S i (t) and h i can be completely specified by the angles θ i (t) and α i (t) that they make with the x-axis. The x-component of equation (3) gives, The above equation is rather difficult to solve analytically except in the mean field limit when a site-i interacts with every other site-j of the system (j = i) with strength J = J 0 /N . Let S x i and S y i be the components of XY spin S i along the x and y axes respectively. We look for a solution of equation (4) in the case when the spins may be ordered along the x-axis, but there is no global ordering in the system in the y direction. We write, The above equation defines a time dependent order parameter cos θ(t), or equivalently a magnetization m(t) = cos θ(t) as the average value of the component of S i (t) along the applied field h. We shall mostly use the notation m(t), but keep cos θ(t) for occasional use when convenient to do so.
Substituting from equation (5) into equation (4) we get, Equation (6) has a nice geometrical interpretation suggested by Mirollo and Strogatz [22] who analyzed the fixed point equations for the XY model for h = 0 rather than the time dependent equation for m(t). Note that the quantity J 0 m(t) + h is the mean field trying to align S i (t) along the x-axis. The mean field has the same value at each site. The random field at each site has a component equal to cos α i that (depending upon the sign of cos α i ) supports or opposes the alignment of S i (t) along the x-axis. A geometrical relationship between the angles θ i (t), α i , and the mean field at time t is illustrated by figure (1) which shows two unit vectors separated from each other by a distance J 0 m(t) + h along the x-axis, and making angles θ i (t + 1) and α i respectively with the x-axis. From the geometry of figure (1), we may write Also, a well known identity relating the sines of the angles of a triangle to the its sides gives, Equation (6) is the most convenient form for studying the evolution of the order parameter m(t) but equations (7) and (8) are useful to get a geometrical picture of the spin configuration of the system. For example, in the limit J 0 m(t) + h → 0, equation (8) gives θ i = α i as may be expected. If J 0 + m(t) = 1, equation (8) gives θ i (t + 1) = α i /2. This is expected as well. In this case the mean field as well as the random field have unit magnitude. One acts along the x-axis and the other makes at an angle α i with the x-axis. Therefore the resultant field aligns the spin at an angle α i /2 with the x-axis.
We obtain a recursion relation for m(t + 1) = cos θ(t + 1) by averaging equation (6) over all sites, A similar mean field equation is obtained for the Heisenberg model. A Heisenberg spin S i (t) may be specified by an azimuthal angle φ i (t) that the spin makes from a fixed axis (say the y-axis) in the yz plane and the polar angle θ i (t) that it makes with the x-axis. The random field h i is also to be specified by a polar angle α i (t), and an azimuthal angle ψ i (t). As in the case of the XY model, we assume that the field h is applied in the x-direction, and any global order in the system lies along the x-direction only.
This gives us the following mean field equation for the Heisenberg model analogous to equation (9) for the XY model:

III. HYSTERESIS
We use the dynamics described above to obtain magnetization curves in a slowly varying cyclic field. The field is increased from h = −∞ to h = ∞ and then decreased to h = −∞ so very slowly that the system has sufficient time to settle into a local minimum of energy at each point. In practice we start with a large negative field when the stable configuration of the system has all spins aligned along the negative x-axis, and then increase the field in small steps till all spins point along the positive x-axis. At each step, the field is held fixed while the system relaxes to a fixed point configuration under the dynamics considered above. This yields a line of fixed points. The graph of magnetization of fixed point configurations versus the applied field gives the magnetization curve in increasing field. Magnetization in decreasing field is obtained similarly. If the magnetization in decreasing field follows a different path than the one in increasing field, the system is said to show hysteresis i.e. history-dependent effects.
We wish to know if the system characterized by Hamiltonian (1) exhibits hysteresis, and if so what is the shape of the hysteresis loop. Another question of interest is whether there is a critical value of disorder that qualitatively separates the hysteretic response of weakly disordered systems from that of strongly disordered systems. The meaning of critical disorder in this context is best explained by a reference to earlier work of Sethna et al [18] on disorder-driven hysteresis in the random field Ising model. They consider a Hamiltonian similar to (1) but in their case the spins and the fields h and h i are scalar quantities; spins take the values ±1, and h i is a random variable chosen from a Gaussian distribution centered at zero and having variance equal to σ 2 . Their results are based on a combination of numerical simulations and analysis, but are quite intuitive as well. These may be summarized as follows. In the limit σ → 0, as the applied field is increased from h = −∞ to h = ∞, each spin and therefore the magnetization per site flips up from −1 to +1 at h = zJ where z is the number of nearest neighbors on the lattice. As σ is increased, the size of the jump in the magnetization decreases and eventually vanishes at h = h c if σ = σ c . For σ > σ c , the magnetization becomes a smooth function of the applied field. The point {h = h c , σ = σ c } is a non-equilibrium critical point characterized by diverging correlation length and scaling laws reminiscent of equilibrium critical phenomena. The parameter σ measures the width of the random field distribution and therefore the amount of disorder in the system. The disorder is said to be critical if σ = σ c . The non-equilibrium critical point may also be studied by fixing the disorder in the system, say by setting σ = 1 and tuning the exchange interaction J and the applied field h to the critical point {h = h c , J = J c }. Now the magnetization curves in increasing and decreasing fields would be smooth for J < J c , but discontinuous for J > J c . The size of the discontinuity would go to zero as J approaches J c from above.
The question is if there is a critical value J c as we go from scalar to vector spins? In the random field Ising model the spins have the value +1 or −1. Therefore the boundaries between domains of positive and negative magnetization are sharp. The width of the domain wall is equal to the distance between nearest neighbors on the lattice. Vector spins can continuously change their orientation from one domain to another over arbitrarily thick domain walls. Phase transition in a system depends on the balance between energy gained by forming a large domain, and energy lost in having to protect it by a domain wall. The energetics of this competition in continuous spins is very different from that in Ising spins [1]. It shows that continuous spins in random fields cannot acquire a spontaneous long range order below four dimensions, while the lower critical dimension for Ising spins is two. This means that the critical hysteresis observed in the random field Ising model in three dimensions may disappear when we go over to vector spins. Although the focus of our work is on the shapes of hysteresis rather than criticality, we shall return to this point after presenting our results.
It is useful to have a brief preview of our results before getting into details. It also gives us an opportunity to mention some unusual aspects of hysteresis in continuous spin systems. In our model, the disorder has a fixed magnitude and sets the energy scale of the system. The behavior of the model is therefore determined by the parameter J. If J = 0, the spins decouple and there can be no hysteresis in the zero frequency limit of the driving field. We find that the behavior of the model for small values of J is qualitatively similar to the behavior for J = 0. This is true in the mean field analysis as well as numerical simulations of the model on a lattice with nearest neighbor interactions. For large values of J we may expect hysteresis as well as jumps in the magnetization. The basis for this expectation is the following. Large J means relatively weak disorder. Thus the spins are mostly aligned parallel to each other. As the applied field is swept from h = −∞ to h = ∞, we expect the majority of spins to reverse their direction at a critical field h = h c . The field h c is determined by the energy required to flip the least stable spin in the system that triggers a large avalanche of flipped spins. For discrete Ising spins with z nearest neighbors, h c is of the order of zJ in the limit of weak disorder. However, in the case of continuous spins, the least stable spin can reverse itself by rotating smoothly along with its neighbors. In other words, the energy barrier for magnetization reversal may be zero for continuous spins in the strong coupling limit just as it is in the weak coupling limit. We find that the mean field theory predicts a non-zero value for h c but simulations based on short range interactions on a lattice indicate h c → 0 in the limit J → ∞.
Hysteresis in continuous spin systems at intermediate values of J where order and disorder compete with each other has several unusual features. Normally if a system shows hysteresis, the magnetization curves in increasing and decreasing fields are separated by the widest margin in the middle at h = 0. We find that there is a range of J values where the magnetization curves for the XY model in the mean field approximation overlap each other in the middle but split from each other as we go away from h = 0 in either direction. Numerical simulations of the XY model show a qualitatively similar behavior although there are significant differences between simulations and the predictions of the mean field theory. Broadly speaking, discontinuities in the magnetization curves predicted by the mean field theory appear to be absent in simulations. The mean field theory of hysteresis in the Heisenberg model has an unusual feature as well. Usually the mean field solution is determined by the intersection of a straight line with an S-shaped curve. In this case the mid portion of the S-shaped curve is a straight line itself. This gives rise to some interesting effects that are seen in corresponding simulations as well. In the following, we examine these issues in detail.

A. XY model
It is instructive to look at the mean field dynamics of the XY model numerically before presenting the analytic solution. Let us set the applied field equal to zero (h = 0), start with an arbitrary initial state characterized by magnetization m 0 , and iterate equation (9) until a fixed point is reached. The results are shown in figure (2). We find two critical values of J: J c1 ≈ 1.489, and J c2 = 2. These values characterize discontinuities in the fixed point behavior in increasing and decreasing J respectively as described below.
The blue curve in figure (2) shows magnetization of fixed points of equation (9) for increasing J. We start with J = 0, and increase J in small steps of ∆J. At each value of J, the magnetization m(J − ∆J) of the previous fixed point is used as an starting point for the iteration of equations. The precise value of ∆J is unimportant. We have chosen a value of ∆J that is small enough so that the line of fixed points appears as a continuous curve on the scale of figure (2). For increasing J, the fixed point magnetization m is zero in the range 0 ≤ J < J c2 . At J = J c2 , it jumps to m ≈ .92, and follows the blue curve as J is increased further. The return path in decreasing J is identical with the blue curve up to J ≤ J c2 , but there is no discontinuity in the return path at J = J c2 . It continues smoothly along the green curve up to J = J c1 at which point it jumps down to zero and remains zero for 0 ≤ J < J c1 . The red curve shows a set of unstable fixed points in the range J c1 ≤ J ≤ J c2 . An unstable fixed point is not realized under iterations of equation (9) because its domain of attraction is zero. However, for a fixed J, the magnetization m u of the unstable fixed point separates the domains of attraction of the two stable fixed points at the same value of J. If the magnetization of the starting state m 0 is less than m u , the equations iterate to the fixed point associated with increasing J. If m 0 > m u the equations iterate to the corresponding fixed point for decreasing J. The reason for the existence of two stable and one unstable fixed point in the range J c1 ≤ J ≤ J c2 may be understood analytically as follows. Let us define, The quantity f (u t ) can be written in terms of complete elliptic integrals of the first and second kinds [22]: The red curve in figure (3) shows a graph of f (u) vs u. It is a continuous, increasing function of u in the range 0 ≤ u ≤ ∞. Some special values are: f (0) = 0; f (1) = 2 π ; f (∞) = 1; f ′ (0) = 1 2 ; lim u→1 f ′ (u) = ∞. Using equation (12), equation (9) may be rewritten as, Fixed points of equation (14) are the roots of the equation, The roots of equation (15) are determined by the intersection of the curve f (u) with the straight line (u−h)/J 0 . For h = 0, u = 0 is always a root because f (0) = 0, and the straight line u/J 0 passes through origin. However, equation (15) may have up to three more roots because f (u) is an S-shaped curve, and a straight line with appropriate slope may cut it at three points. Consider straight lines passing through the origin and having decreasing slopes i.e. lines u/J 0 with increasing J 0 . Let J c1 and J c2 be the smallest and the largest values of J 0 respectively at which the line u/J 0 meets the S-shaped curve f (u) tangentially as shown in figure (3). As mentioned at the beginning of this section, J c1 ≈ 1.489 and J c2 = 2. There are three non-zero roots of equation (15) in the range J c1 ≤ J 0 ≤ J c2 , and only one non-zero root for J 0 > J c2 . Which of these roots is actually realized by the dynamics is determined by the starting point u 0 used in iterating equation (14). The stability of a root can be checked by analyzing equation (14) in the neighborhood of its fixed point [22]. However, the full equation is necessary to determine the domain of attraction of a stable fixed point.
Next we consider equation (14) for a fixed value of J 0 but in a varying field h. Starting from a sufficiently negative field h = h min where the stable configuration has most spins pointing along the negative x-axis, the field is increased in small steps ∆h to h = h max where most spins point along the positive x-axis. Figure 4 shows the fixed point magnetization m as the field is increased from h min = −1.5 to h max = 1.5 and back to h min = −1.5 in steps of size ∆h = .01. Data for three representative values of J are shown: J 0 = 2 (red), J 0 = 1.25 (green), and J 0 = .25 (blue). J = 2 shows a familiar looking hysteresis loop but J = 1.25 shows a somewhat unfamiliar behavior. In this case, there are two symmetrically placed windows of positive and negative applied fields where the system shows hysteresis but there is no hysteresis in the intermediate region near zero applied field. For J = .25 there is no discernible hysteresis on the scale of figure (4).
The variety of behavior seen in figure (4) may be understood as follows. We saw in figure (3) that spontaneous magnetization is possible only if J 0 > J c1 ≈ 1.498. Spontaneous magnetization in zero applied field gives rise to the possibility of hysteresis as the applied field is cycled up and down across the value h = 0. Therefore the hysteresis loop for J 0 = 2 (red curve) in figure (4) centered around h = 0 is to be expected. We do not expect the green curve (J 0 = 1.25), or the blue curve(J 0 = .25) to show a hysteresis at h = 0. This is born out by figure (4). What is surprising at first sight is that the green curve in figure (4) shows two small hysteresis loops in applied fields centered around h ≈ ±0.2. We can understand this with the help of figure (5)

B. Heisenberg model
Making a transformation of variables µ i = cos α i , and u t = J 0 m(t) + h, equation (11) may be rewritten as: The integral in equation (16) is easily evaluated and yields This is rather unusual in a mean field theory. Normally, the spontaneous magnetization in a mean field theory is determined by the intersection of a straight line with an S-shaped curve. In the present case the S-shaped curve is itself a straight line in the interval −1 ≤ J 0 m ≤ 1. This means that in the absence of an applied field, the zero temperature magnetization in the random field Heisenberg model is zero if J 0 < 3 2 , and can have an arbitrary value in the range − 2 3 ≤ m ≤ 2 3 if J 0 = 3 2 . For J 0 > 3 2 , |m| > 2 3 and increases with applied field h. Figure (7) shows the magnetization curves in a cyclic field (varying infinitely slowly in the sense explained earlier) for J 0 = 1(pink), J 0 = 1.5(blue), J 0 = 1.75(green), and J 0 = 2(red). As expected from the above analysis, there is no hysteresis in the case J 0 = 1 and J 0 = 1.5, although in the case J 0 = 1.5 the magnetization shows a finite jump at h = 0. There is hysteresis for J 0 = 1.75 and J 0 = 2 with the area of the hysteresis loop increasing with J 0 .

C. Peculiar criticality
In this section we attempt to place our calculations in the context of extant work on critical hysteresis in ncomponent vector spin systems with quenched disorder. The extant work employs soft continuous spins and a Gaussian distribution of quenched field, while we have used hard continuous spins and random fields in the form of randomly oriented unit vectors. Vector spins with n ≥ 2 are called continuous spins. These can be hard or soft. Hard continuous spins have a fixed length but can make any angle from a reference axis. Computer simulations commonly use hard spins on a lattice. We have used hard spins for numerical as well as analytic work. Soft spins are continuous in angle as well as magnitude. Momentum space renormalization group uses soft spins [13,18]. It usually starts out with hard spins on a lattice but transforms them into soft spins that can take any real value but are constrained to remain close to a fixed length. This is done by introducing an effective on-site potential. Similarly continuum limit of the lattice is taken but a cutoff on the maximum wave-vector is introduced. There is some evidence that the critical behavior of models in the renormalization group theory is independent of the additional parameters introduced by the effective on-site potential and momentum cutoff [23]. However, it is not independent of the form of the random field distribution.
Hartmann et al [25] have shown numerically that the critical exponents of three dimensional random field Ising model with Gaussian distribution of random fields are significantly different from those of the same model with a bimodal distribution of random fields. Analytic results in three dimensions are not available. What is available is a perturbation series for critical exponents in 6 − ǫ dimensions for a Gaussian distribution of the quenched field [13]. It shows that similar to the random field Ising (n=1) model [18], random field XY (n=2) and Heisenberg (n=3) models have a critical point on the hysteresis loop. The critical exponents depend on n below 6 dimensions, but are independent of n in 6 and higher dimensions. The significance of 6 dimensions is that it is equal to the upper critical dimension of the model with a Gaussian distribution of the random field [24]. In 6 and higher dimensions the action is adequately described by a quadratic term, and higher order terms become irrelevant in the renormalization group sense. The quadratic action can be solved exactly, and the solution is often called the mean field solution (presumably because it gives the same critical exponents as a mean field solution based on infinite range interactions). In this variant of the mean field theory, the hysteresis loops for n=1, 2, and 3 exhibit a critical point as the width of the Gaussian disorder increases, but the critical exponents do not depend upon n.
In our variant of the mean field theory based on infinite range interactions, the critical behavior for n = 2 is different from that of n = 3. This should not raise a serious concern because we use a different distribution of the random field than used in reference [13]. Our random field distribution for n= 2 and 3 is a continuous analog of bimodal distribution in the case n=1. For n = 1, we know that Gaussian and bimodal distributions of random field give two different sets of critical exponents in 3 dimensions. This difference may persist even above the upper critical dimension although to our knowledge the upper critical dimension for a bimodal distribution is not known precisely. Nevertheless the striking difference between the nature of criticality for n = 2 and n = 3 in our mean field theory is interesting and could not have been anticipated beforehand. We may therefore take a closer look at the algebraic mechanism producing this difference and also the difference from the mean field theory of the random field Ising model based on infinite range interactions and a Gaussian distribution of the random field. In these exactly solved cases the equations of motion have a similar form but the signs of various terms depend on the details of the model. This produces distinct critical behavior in each case. We have, where u t = J 0 m(t) + h. The first equation is for the random field Ising model with a Gaussian distribution centered at zero and having unit variance. It can be easily derived and at its fixed point it reduces to the mean field equation studied by Sethna et al [18]. The equations for the XY and Heisenberg models were derived in sections III A and III B where the functions f (u) and g(u) are also defined. For simplicity, let us set h = 0 and confine to u ≥ 0. Now suppose we start with a small value of u t and iterate the above equations of motion till we reach a fixed point u * . We focus on the behavior of u * as a function of J. In each case, there is a threshold J c such that u * = 0 for J < J c . We get J c = π 2 , 2, and 3 2 for n=1, 2, and 3 respectively. At J = J c , there is a transition to a non-zero value of u * . This transition is continuous for n = 1, discontinuous for n = 2, and peculiarly discontinuous for n = 3 in the sense that u * can have any value in the range 0 ≤ u * ≤ 1. Thus the transitions for n = 1, 2, and 3 are distinct from each other.
It is not difficult to understand the above results analytically. The functions Erf(u), f (u), and g(u) are all zero at u = 0 and their first derivatives with respect to u are positive. Erf(u) is concave down for u ≥ 0; f (u) is concave up for 0 ≤ u < 1 and concave down for u > 1; g(u) has zero curvature for 0 ≤ u < 1 and concave down for u > 1. It is also instructive to write the leading terms in the series expansion of the right hand side. We get the following expressions for the Ising, XY, and Heisenberg spins respectively, n = 1 : The leading terms of the recursion relations in the limit u t → 0 show that the fixed point u * → 0 if 2 π J 0 < 1 for n = 1, J0 2 < 1 for n = 2, and 2J0 3 < 1 for n = 3. This yields the critical values J c mentioned above. For n = 1 and J > J c , u * ≈ (J − J c ) 1/2 i.e. it has a square root singularity characteristic of mean field critical behavior. In the case n = 2, the cubic term has a positive sign and therefore a physically acceptable solution does not grow continuously from u * = 0 at J = J c . The solution of the full equation shows that u * has a first order jump in this case. The recursion relation for n = 3 is peculiar because it does not have any non-linear terms. At J c = 3/2 any value of u * (0 ≤ u * ≤ 1) satisfies the fixed point equation. For J > J c , u * increases with J but remains bounded below J.  (9) show magnetization curves for the random field XY and Heisenberg models respectively in a slowly varying cyclic field. The data is obtained from simulation of the model on a simple cubic (sc) lattice with nearest neighbor (nn) interactions. In order to keep the computer time within reasonable limits, the XY model is simulated on a lattice of size 100 3 , and the Heisenberg model on a lattice of size 50 3 with periodic boundary conditions. Graphs are presented for various values of J as indicated in the captions for the figures. For each value of J, the applied field is cycled in small steps between two large values that saturate the magnetization along negative and positive x-axis respectively. For clarity, the figures depict only a part of the simulation data in a small range of the applied field where variation in magnetization is most pronounced. At each step of the applied field the system is allowed to relax till it reaches a fixed point. We assume that the system has reached a fixed point if the projection of each spin along x-axis remains invariant within an error of 10 −5 .
The mean field theory predicts the absence of hysteresis in the XY model if J 0 < 1.498. The energy scale in our model is set by the disorder term. Thus the behavior of the mean field model at J 0 may be compared with the behavior of the nn model on a sc lattice at 6J. As an order of magnitude estimate, we expect the absence of hysteresis on a sc lattice if 6J < 1.498 or J < .25 approximately. This is qualitatively in accordance with the result of simulations shown in figure (8). The magnetization curves for J = .1 show no discernible hysteresis on the scale of the figure. At J = .2, we find two isolated hysteresis loops separated by a region of zero hysteresis near h = 0. This is qualitatively similar to the prediction of the mean field theory. With increasing J the two isolated loops widen, gradually merge with each other, and the overall shape of the hysteresis loop evolves as indicated in figure (8). For much larger values of J the hysteresis loop becomes narrower and more vertical. Within numerical errors, magnetization curves in increasing and decreasing fields approach a step function at h = 0, and hysteresis appears to vanish for J ≥ 1. The large J regime marks a qualitative difference between the prediction of the mean field theory and the simulations. The mean field theory predicts hysteresis but simulations on cubic lattices with nearest neighbor interactions show no hysteresis. This discrepancy may be attributed to the use of infinite range interactions in the mean field theory. The energy barrier for rotation of a strategically placed spin may be significantly smaller if its nearest neighbors alone are taken into account rather than all spins in the system. The dynamics based on nn interactions initiates a rotation at the least stable site and gradually spreads it on adjacent sites in the neighborhood. Large J simulations take an enormously long time to reach a fixed point in the neighborhood of h = 0, but the end result appears to be simply a reversal of saturation magnetization when the sign of h is reversed. In the limit J → ∞, the system effectively acts as a single spin having the total magnetization of the system. Just as an isolated spin in the limit J = 0 does not show any hysteresis so also the entire system in the limit J → ∞. The main difference between the magnetization curves in the limits J → 0 and J → ∞ lies in their shape, but this is understandable if we rescale the applied field appropriately with the total magnetization of the system.
Simulations of the Heisenberg model are also in reasonable agreement with the predictions of the mean field theory except for large values of J. The mean field theory predicts hysteresis if J 0 > 3/2. This corresponds to J > .25 approximately. Simulations do not show any significant hysteresis if J ≤ . 25. Figure (9) shows a magnetization curve for J = .25 that reverses itself when the field is reversed. The magnetization is linear in the applied field over a wide range around h = 0. This is in qualitative agreement with the prediction of the mean field theory. Simulations for J = .4 and J = .5 show typical hysteresis loops although the range of applied field over which perceptible hysteresis is observed is an order of magnitude smaller than the range predicted by the mean field theory. The main difference between the simulations and the mean field theory lies at large values of J. The magnetization curves shown in figure (9) for J = 1 appear to form a narrow nearly vertical hysteresis loop. However size of the steps used to increase and decrease the applied field in figure (9) is of the order of the width of the apparent hysteresis loop. Simulations based on smaller steps and higher accuracy in determining the fixed points suggest that the hysteresis loop vanishes for J ≥ 1 and the magnetization has a first-order jump at h = 0.

IV. CONCLUDING REMARKS
We have analyzed a simple model to study the effect of quenched disorder on hysteresis in magnetic systems of continuous symmetry. The model is obtained by adding quenched disorder and zero temperature dynamics to the well established XY and Heisenberg models of ferromagnetism. The quenched disorder is in the form of randomly oriented fields of unit magnitude. Is this model applicable to experiments? We have argued that thermal fluctuations are of secondary importance in disorder-driven hysteresis. Therefore the use of zero temperature dynamics may not be serious. It has the virtue of being deterministic and therefore easier to analyze theoretically. A large number of studies on disordered systems employ zero temperature dynamics for these reasons. Randomly oriented crystal fields are also not uncommon in amorphous materials. These are dipolar or quadrupolar but if the activation barriers are large, may act like quenched random fields as a spin or domain pointing one way gets hard to dislodge. Thus the basic ingredients of our model are chosen to make a minimal model for understanding experiments. The parameters of the resulting model are: n components of vector spins, exchange interaction J, and the applied field h. Effectively, there are just two parameters; integer n and real J. This is because the middle term in Hamiltonian (1) does not have a tunable value, and the field h is cycled between −∞ and ∞. A two parameter model may not capture details of hysteresis in various materials but it provides a caricature of experimental observations. The variety of shapes of hysteresis loops are particularly striking for the XY model (n=2). As J is varied, we get familiar as well as rather unusual shapes of loops. The unusual shapes have been noted earlier in magnetic and other materials. These are known as wasp-waisted [26] or double-flag shaped loops [27]. These shapes have a kind of weak universality in the sense that they are seen in the mean field theory, simulations on three dimensional lattices, and experiments in diverse systems. Similar shapes are also seen in the random field Blume-Emery-Griffiths model and other models of plastic depinning of driven disordered systems.
Soft continuous spins with Gaussian random fields have been used earlier to study critical hysteresis in 6 − ǫ dimensions in the renormalization group theory. Where does our mean field calculation sit in this context? We note that our mean field results do not match the renormalization group results in any limit. There are two possible reasons for this. Firstly, we have used a different distribution of random fields than the one used in 6 − ǫ expansion. The form of random field distribution appears to be important in determining critical hysteresis. Secondly, we do not have the benefit of an appropriate renormalization group study of our variant of the model, nor do we know the upper critical dimension of our model precisely. Before the renormalization group theory, mean field theory was viewed as an approximate but self consistent theory of critical behavior in 3 dimensions because it neglected fluctuations. This variant of mean field theory was based on infinitely weak but long ranged interactions in the system. It represented the effect of the entire system on an individual spin by an effective field while keeping the energy of the system extensive. The renormalization group has given another connotation to mean field theory. In its framework, the mean field theory becomes a reduced theory based on a quadratic action that is exact at and above an upper critical dimension where fluctuations are negligible. The upper critical dimension for pure (non-disordered) magnetic systems is 4, and in this case the two variants of the mean field theory predict the same critical behavior. This is understandable because the effective field is proportional to the order parameter and the effective action is therefore quadratic. The upper critical dimension for a disordered spin model with a Gaussian random field is 6. In this case also, an explicit calculation for the random field Ising model with Gaussian field shows that the two variants of mean field theory yield the same critical behavior. However, when the randomness is of the form of randomly oriented unit vectors, we have analyzed only one variant of mean field theory that is based on infinitely weak but long range forces. It predicts strikingly different critical behavior in XY and Heisenberg models respectively. The case n = 2 has a first order transition, and n = 3 an unusual transition as discussed in section III C. A somewhat similar case of first as well as second order depinning transition in the mean field theory occurs in a viscoelasic model of driven disordered systems [28,29]. Thus we have a number of model specific results. Evidently more work is required to make any general connection between the random field distribution and the nature of criticality in the model, and to connect a conventional mean field theory to a limiting form of renormalization group theory above an upper critical dimension.
The new framework for understanding critical behavior also uses the idea of a lower critical dimension. Below the lower critical dimension the fluctuations are so great that the system does not order at all and therefore there is no question of a phase transition. The lower critical dimension for an equilibrium transition in an n-component spin system in a Gaussian random field is 2 for n = 1, and 4 for n ≥ 2. To our knowledge, the lower critical dimension for the case of randomly oriented unit vectors is not known. However, we mention a few issues that may bear on experiments in three dimensions irrespective of the form of randomness characterizing the system. It has been argued that critical hysteresis in a Gaussian random field Ising model is in the same universality class as the corresponding equilibrium critical point [13,18]. There is a reasonable experimental evidence for this [18]. In the case of Gaussian random field XY and Heisenberg models, the lower critical dimension lies above 3. Does it necessarily mean the absence of critical hysteresis in these models in 3 dimensions? The situation is not entirely convincing either theoretically or experimentally. The argument for a lower critical dimension is based on spontaneous symmetry breaking in the absence of an applied field. If the critical point on the hysteresis loop were to occur at a non-zero value of magnetization or the applied field then a unique direction is already chosen by the corresponding magnetization or the applied field. In this case we may observe critical hysteresis in 3 dimensions with critical exponents appropriate for the Gaussian random field Ising model. Much of critical hysteresis seen in experiments may belong to this case but there is discrepancy between some experiments and theory [14]. It may be that quenched disorder in experimental systems is not characterized adequately by Gaussian random fields or randomly oriented unit vectors. The presence of demagnetizing fields and dipolar forces in materials used for experiments are likely to change the simple theoretical picture based on on-site random field disorder. This applies equally to mean field theory and the renormalization group approach based on expansions around a quadratic action.
In the absence of exact solutions in three dimensions, simulations of models may be more relevant to experiments. Our simulations produce rather smooth hysteresis loops for small and moderate values of J suggesting the absence of jumps in the magnetization. Phase transitions in complex systems are difficult to decide on the basis of numerical work alone, and therefore we have focused on the shape of hysteresis loops. The shapes are not universal but this does not diminish their importance in the application of magnetic materials. The relationship between the shape of hysteresis loops and the defect mediated process of magnetization reversal is also interesting. This has been studied at zero temperature numerically in a two-dimensional XY-model with weak random anisotropy [12]. We hope future studies on these lines will clarify the relationship of the hysteresis loops to the underlying energy landscape as well as the patterns of spin configurations such as vortex loops in three dimensional XY model.  The random field has a fixed magnitude equal to unity. The hysteresis loop for J0 = 2 has a familiar shape but J0 = 1.25 shows rather unusual hysteresis in two small windows of applied field situated at h ≈ −.2, and h ≈ .2 respectively but no hysteresis outside these windows. In particular, there is no hysteresis at or near zero applied field. For J0 = .25, there is no discernible hysteresis in any region of the applied field on the scale of the above figure.