PHYSICAL REVIEW B

VOLUME 58, NUMBER 9

1 SEPTEMBER 1998-I

Magnetic hysteresis in an Ising-like dipole-dipole model

Gyorgy Szabo and Gyorgy Kadar
Research Institute for Technical Physics and Materials Science, POB 49, H-1525 Budapest, Hungary Received 6 April 1998 Using zero-temperature Monte Carlo simulations we have studied the magnetic hysteresis in a threedimensional Ising model with nearest-neighbor exchange and dipolar interaction. The average magnetization of spins located inside a sphere on a cubic lattice is determined as a function of magnetic eld varied periodically. The simulations have justied the appearance of hysteresis and allowed us to have a deeper insight into the series of metastable states developed during this process. S0163-1829 98 02133-X

I. INTRODUCTION

H The description of hysteresispredominantly for magnetization curveshas been the aim of numerous papers for more than a century now. The early phenomenological model of Lord Rayleigh was the fundamental idea of the well-known Preisach model of ferromagnetic hysteresis, which has been further developed and widely discussed together with other descriptive phenomenological models in a long series of works1 appearing up to the present day. The physical explanation for the lag of the magnetization component behind the external magnetic eld varying along a given line was elaborated by Stoner and Wohlfarth2 in a simple micromagnetic model for the case of a single-domain particle of uniform magnetization. With the advance of the techniques of statistical physics recently great interest seems to be oriented towards the investigation of hysteretic phenomena in realistic multiparticle and/or multidomain systems. The aim after all would be narrowing the gap between the phenomenological top-down and the physical bottomup approaches to the description of hysteresis in general. In this paper Monte Carlo simulation of a multiparticle Ising system of pointlike elementary magnetic dipoles will be presented. The dipoles contained inside a sphere are arranged in a three-dimensional simple cubic lattice, they are parallel or antiparallel to the z axis of the Cartesian system, and interact by the nearest-neighbor exchange as well as the long-range dipole-dipole couplings. It will be shown that this model exhibits magnetic hysteresis if the external magnetic eld is varied periodically. The present approach allows us to visualize and analyze the evolution of spin congurations.
II. MODEL

J
r,r

h
r

1 2

V r r
r,r r r

r ,

where (r) 1 and 1 for up and down spins. In the rst term the sum is over the nearest-neighbor pairs and J denotes the strength of the usual exchange interaction measured in the above-mentioned energy unit. The second term describes the effect of the external magnetic eld h. The dipole-dipole interaction between two spins having only a z component is dened as V r r 2 3z 2 . r5 2

Notice that this dipole-dipole interaction provides ferromagnetic coupling along the z axis and the coupling becomes antiferromagnetic in the x-y plane. Furthermore, the average value of the eld of a given dipole over the points located at a prescribed distance vanishes due to the cubic symmetry. Conversely, up or down spins on a spherical shell produces zero magnetic eld at the central lattice point. This is the reason why the resultant magnetic eld is zero at the center of a spherical sample if all the spins point to up or down .36 At a given site r the magnetic eld produced by the remaining dipoles is given as hd r
r r r

V r r

r .

We consider a three-dimensional Ising model with spins located on the points of a simple cubic lattice [r (x,y,z); x, y, and z are integers within a sphere of radius R (r 2 x 2 y 2 z 2 R 2 ). Dimensionless quantities will be used throughout the paper. Namely, the length will be measured in terms of lattice constant a, the dipole moments in terms of the unique dipole moment , the external magnetic eld h will be expressed in units of /a 3 , and the energy per spins in terms of 2 /a 3 . In this case the Hamiltonian is given by
0163-1829/98/58 9 /5584 4 /$15.00 PRB 58

Our analysis is restricted to spherical systems because here the early theories predict zero eld inside the ferromagnetic sphere.36 More precisely, Cohen and Keffer6 have shown that the contribution of pointlike dipoles to the local eld may differ from zero in the close vicinity of the surface. We have numerically studied the local eld because this phenomenon plays a crucial role in the magnetic hysteresis as will be described in the next section. In the ferromagnetic state the average h d (r) energy per sites due to dipolar interactions is zero5 independently of the system size. In this case the local eld satises the con5584 1998 The American Physical Society

PRB 58

MAGNETIC HYSTERESIS IN AN ISING-LIKE . . .

5585

FIG. 1. Radial dependence of the local magnetic eld if m 1 and R 15.33 for choosing different directions: 1,0,0 , diamonds; 0,0,1 , open squares; 1,1,1 , open circles; 5,0,1 , ; and 1,0,5 , .

FIG. 2. Magnetic hysteresis for J 0 and R 15.33.

ditions of reection h d (x,y,z) h d ( x, y, z) and rotation h d (x,y,z) h d ( y,x,z) symmetries. The numerical results see Fig. 1 demonstrate clearly that inside our sphere h d (r) 0 in agreement with the predictions of analytical calculations.36 At the same time large variations are indicated on the outer shells. For (r) 1 the highest values of h d appear along the periphery of the top and bottom layers. Along the 1,1,1 directions the local eld vanishes see the open circles in Fig. 1 . The lowest eld values are found at the eight sites symmetrically equivalent to 15,0,3 if R 15.33. The probability distribution of the lattice points as a function of h d exhibits a sharp peak around h d 0, and its maximum value increases with R. In a cube-shaped sample such calculation yields a signicantly different wide probability distribution that causes drastic changes in the hysteresis too. The dipolar energy of ordered spin congurations was determined previously in innite systems of cubic symmetry. Using the Ewald method Luttinger and Tisza5 evaluated the energy per sites for all the periodic ordered structures characterized by a spin conguration within the 2 2 2 unit cell. In the energetically favored states the up or down spins form vertical columns as expected. For h 0 and J 0 the spin conguration of lowest energy is a twofold degenerated chessboardlike antiferromagnetic arrangement in the x-y plane of ferromagnetic columns along the z direction, that is, (r) 1 ( 1) if x y is odd even . In this columnar antiferromagnetic CAF structure the energy per sites is given as 5 E CAF 2.676 J. 4

15.33. For both the above-mentioned spin congurations it is found that nite-size corrections are proportional to 1/R. The numerical technique has allowed us to study periodic antiferromagnetic structures different from those studied by Luttinger and Tisza.5 These calculations indicate that double or multilayer structures similar to LAF become stable for sufciently strong ferromagnetic coupling and the preferred layer width increases with J. These indications, however, should be considered as preliminary results because more systematic analyses are required to investigate the effect of anti ferromagnetic coupling on the spin congurations in the ground state, the size effects, etc. Henceforth, we will concentrate on the model with J 0. In this case the slow cooling Monte Carlo technique has conrmed that the twofold-degenerated ground state is equivalent to the CAF spin conguration in zero external eld (h 0).
III. SIMULATION OF HYSTERESIS

A series of zero-temperature Monte Carlo simulations has been performed to study hyteresis phenomena in the model described above. In these simulations the system is started from a random spin conguration. For an elementary process a randomly chosen spin is ipped if h d (r) h (r) 0, otherwise the spin value remains unchanged. This process is repeated until all the spin signs become equivalent to the sign of the local magnetic eld. Then the external magnetic eld (h) is increased decreased by h and by repeating the above-mentioned spin-ip processes the system is allowed to relax toward a new local energy minimum. In such a way the external eld is varied periodically with an amplitude of 10. During this procedure we monitored the magnetization dened as m 1 N r ,
r

Furthermore, for the fourfold-degenerated layered antiferromagnetic LAF spin conguration, where (r) 1 ( 1) if x or y) is odd even , the energy per sites is E LAF 2.422 J. 5

It is obvious from Eqs. 4 and 5 that the CAF state is preferred to LAF if the ferromagnetic coupling J 0.127. In our spherical model numerical calculations are performed to check the nite-size effects choosing R 10.25, 12.25, and

where N is the number of spins inside the sphere. These simulations have clearly indicated the appearance of the usual magnetic hysteresis. A typical result obtained for J 0, R 15.33 and h 0.1 is shown in Fig. 2. To check the size effect the simulations have been carried out for different sizes (R 10.25 and 12.25 . The comparisons have

5586

GYORGY SZABO AND GYORGY KADAR

PRB 58

FIG. 3. Spin congurations when decreasing the magnetic eld from h 10 to 2 a , 1.7 b , 0 c , 1 d , 4 e , and 6 f . Black and white spheres indicate up and down spins.

indicated that the size effect is comparable to the uncertianty of data observed between subsequent cycles see Fig. 2 . In other words, the above-mentioned system size containing 15 203 spins is sufciently large to study the hysteresis adequately in this model. Unfortunately, for larger sizes the simulation becomes very time-consuming because the computational time is proportional to R 6 . The choice of h 0.1 is also motivated by the minimization of run time. Further decrease of h does not modify the plotted curves signicantly. Figure 2 indicates the presence of avalanches when varying the magnetic eld. Similar phenomena were observed in experiments by Barkhausen7 and also in the random-eld Ising models.8 Recently different approaches are used to study the avalanches as well as its relation to the hysteresis. 9,10 In the present model the appearance of avalanches is a consequence of dipole-dipole interaction. By varying the magnetic eld several spins of the actual conguration become preferred to ip into the opposite state. Reversing a randomly chosen spin the local eld h d (r) will be modied in all sites of the system. Due to the dipoledipole interaction the neighboring spins in the same column are favored to change direction too. This effect drives the avalanche of spin ips in a given column. During the simulations the spin ips are observed in several columns simultaneously. This process leads predominantly to such congurations where complete spin columns have been reversed as demonstrated in Fig. 3. Figure 3 shows six spin congurations appeared at different magnetic elds during a half-cycle when decreasing h from 10. In order to visualize a more complete threedimensional picture on the spin congurations parts of the horizontal (z 0) and vertical (x 0) cross sections are displayed by removing a quarter of spins from the sphere. In this gure the small black spheres with white border and the white spheres with black border indicate up and down spins. Decreasing the magnetic eld the initial saturated state (m 1) remains unchanged until h 2.253 whose value depends slightly on R. For the given size the spins that can be the rst to ip are positioned at one of the eight sites equivalent to (15,0,3). Consequently, the spins in the corresponding four columns can ip simultaneously, because the short

columns are too far apart to affect each other signicantly. The result is well recognizable for a possible subsequent conguration plotted in Fig. 3 a . Notice that here the mentioned symmetries are no longer valid due to a branching process described as follows. Immediately after the decrease of h some of the spins become reversible. One of them is chosen randomly to ip. This elementary process fastened by the above-mentioned columnar spin-ip avalanche can actually prevent the ips at the rival spins. Thus the order of consecutive steps during the formation of the subsequent columnar structures see Fig. 3 is occasional. This phenomenon causes the hysteresis curves to be unreproducible for nite sizes as demonstrated in Fig. 2. Conguration c in Fig. 3 represents a typical state with a remanent magnetization when reaching h 0. The magnetization vanishes at h 1. Here, the appearance of domains with CAF and LAF structures is striking see conguration d . The CAF state dominates the vicinity of the z axis. The positions of the LAF domains are initialized by the rst columnar spin ips after saturation. These domains remain recognizable in a wide range of the magnetization processes as indicated by conguration e obtained for h 2. Further decrease of h will destroy these ordered regions leaving the spins unchanged only in a few columns positioned randomly conguration f for h 6 . Notice that this conguration differs signicantly from the initial ones compare to conguration a . Finally all the spins point to downward if the magnetic eld becomes less than 6.2. When studying a more complete series of spin congurations one can easily recognize that the back spin ips appear rarely during a half-cycle. According to our numerical investigations the number of back spin ips is less than 1% of the total number of spins. During the simulations we have also recorded the number of nonuniform columns that are found to be zero in the initial and nal stages of demagnetizing processes. This number can occasionally differ from zero and its rate remains below 1%. The low number of back ips explains the phenomena observed when we have articially prevented the free spin ips. In this case the above-described dynamics is modied by assuming a hysteretic behavior for the individual spins. Namely, the spin ips may occur for the randomly chosen sites only if the driving force exceeds a threshold value ( f t 0). Evidently, for f t 0 this modication leaves the above results unchanged. In agreement with the expectation the hysteresis loop becomes wider for positive f t . More precisely, the sides of the loops are shifted outward without causing any observable changes in the slopes. This means that the main features of the subsequent spin congurations are similar to those described above see Fig. 3 and the columnar spin ips are delayed within the cycles. However, the number of back ips during the half-cycle decreases with f t . For example, only a few back ips can be observed for R 15.33 and f t 1 and this elementary process vanishes practically if f t 2. Finally, we have investigated the effect of ferromagnetic coupling (J 0) on the hysteresis. For this purpose the hysteresis curves are recorded for different J values as shown in Fig. 4. In this case the simulations are started from the ferromagnetic state (h 10). We have observed that the extension of avalanches as well as the uncertainty of the magne-

PRB 58

MAGNETIC HYSTERESIS IN AN ISING-LIKE . . . IV. CONCLUSIONS

5587

FIG. 4. Magnetic hysteresis curves for different exhange constants: J 0 a , 0.5 b , 1 c , 2 d , and 3 e for R 15.33. The curves are averaged over 10 cycles.

tization process for the subsequent cycles increases with J. Figure 4 demonstrates that these phenomena are accompanied with the increase of average slope along the side of the loops and the broadening of the hysteresis loop. These observations are related to formation of larger and larger ferromagnetic domains as well as the more and more massive avalanches when increasing the ferromagnetic coupling. One can observe that the initial steps of the demagnetization process is only slightly modied by the nearest-neighbor couplings because the rst spin ips appear on the surface at low latitudes where h d (r) has low negative values. This is the region where the ferromagnetic force is weak because of the low number three or four of neighboring spins; therefore it is not able to prevent the spin ips driven by the local eld h h d (r). Further systematic analysis is required to study what happens when the typical ferromagnetic domain sizes become larger than R.

Numerical simulations have been performed to investigate the hysteresis phenomena in a three-dimensional Ising model involving both the nearest-neighbor exchange and the longrange dipolar interactions. For simplicity the spins are located on a cubic lattice within a sphere. This shape provides zero local eld h d (r) inside the sphere for a ferromagnetic state. The large deviations of h d (r) at the surface are found to play a crucial role at the beginning of demagnetization processes started from a saturated state. Using zerotemperature Monte Carlo simulations the magnetization process is investigated under an alternating external magnetic eld with an amplitude providing complete saturations. During this process we have observed avalanches and branching points leading to a large variation of paths along which the system can evolve. The avalanches consist of spin ips constrained into a few columns. As a result, after having varied the magnetic eld the new stationary metastable states built up dominantly from columns with uniform spins. This feature decreases drastically the number of possible metastable states characterizing the hysteretic behavior in comparison to the total number of congurations. Our simulations have justied that the number of back ips is negligible within the half-cycles of the magnetization process. We think that the present model seems to be a good candidate for exploring the relationship between a microscopic description and the phenomenological e.g., Preisach model approaches of hysteresis. For this purpose, however, we need more detailed information about the internal hysteresis loops, the effects of shape and exchange constant, the local-eld distribution, etc. Fortunately, the model simulations give us a wide scale of opportunity to have a deeper understanding.
ACKNOWLEDGMENT

This research was supported by the Hungarian National Research Fund OTKA under Grant No. T-23555.

Lord Rayleigh, Philos. Mag. 23, 225 1887 ; F. Preisach, Z. Phys. 94, 277 1935 ; A. Visintin, Differential Models of Hysteresis Springer-Verlag, Berlin, 1994 . 2 E. C. Stoner and E. P. Wohlfarth, Philos. Trans. R. Soc. London, Ser. A 240, 599 1948 . 3 J. A. Sauer, Phys. Rev. 57, 142 1940 . 4 C. Kittel, Introduction to Solid State Physics, 5th ed. Wiley, New York, 1976 . 5 J. M. Luttinger and L. Tisza, Phys. Rev. 70, 954 1946 ; 72, 257 1947 .

M. H. Cohen and F. Keffer, Phys. Rev. 99, 1128 1955 . Z. Barkhausen, Z. Phys. 20, 401 1919 . 8 J. P. Sethna, K. Dahmen, S. Kartha, J. A. Krumhansl, B. W. Roberts, and J. D. Shore, Phys. Rev. Lett. 70, 3347 1993 ; K. Dahmen, S. Kartha, J. A. Krumhansl, B. W. Roberts, J. P. Sethna, and J. D. Shore, J. Appl. Phys. 75, 5946 1994 ; K. Dahmen and J. P. Sethna, Phys. Rev. Lett. 71, 3222 1993 ; O. Perkovic, K. Dahmen, and J. P. Sethna, ibid. 75, 4528 1995 . 9 O. Narayan, Phys. Rev. Lett. 77, 3855 1996 . 10 G. Bertotti, Phys. Rev. Lett. 76, 1739 1996 .
7