Random Field Potts model with dipolar-like interactions: hysteresis, avalanches and microstructure

A model for the study of hysteresis and avalanches in a first-order phase transition from a single variant phase to a multivariant phase is presented. The model is based on a modification of the Random Field Potts model with metastable dynamics by adding a dipolar interaction term truncated at nearest neighbors. We focus our study on hysteresis loop properties, on the three-dimensional (3D) microstructure formation and on avalanche statistics.


I. INTRODUCTION
Microstructure formation in first-order phase transitions is a phenomenon that has been studied by physicists, mathematicians and engineers [1,2,3]. It is important not only from a fundamental point of view, but also for applications, due to its relationship with material properties. Microstructures occur in ferroic systems (ferromagnetic, ferroelectric and/or ferroelastic) which are driven through a first-order phase transition (FOPT) in which some symmetry operations of the parent phase are lost. The product phase (usually less symmetric) may appear in energetically equivalent variants which are related by the symmetry operations that have been lost in the transition. The obtained microstructures correspond to the arrangement of such equivalent variants and are decided by the interplay of many energetic terms: interface energy, surface energy and long-range interactions.
Until now many of the studies on microstructures have focused on the determination of the optimal variant configuration minimizing a certain thermodynamic potential that takes into account the above factors and external conditions. Nevertheless, in many cases, when real materials are studied, such optimal microstructures are not observed. This is mainly due to two important factors: (i) the existence of disorder sources of very different natures both in the bulk and on the surfaces and also (ii) the athermal character of the phase transition dynamics: when the temperature is not very high, the energetic barriers that separate the optimal solutions from the parent phase cannot be overcome. Thus, the system evolves following metastable paths which locally optimize the system energy, but are far from the trajectories obtained from a global minimization principle. An interesting suggestion on the behavior of microstructure formation comes from the glass-jamming transition framework (see Refs. 4, 5 and references therein), which is associated with the so-called kinetically constrained model. * Electronic address: benedetta@ecm.ub.es These models are stochastic lattice gases with hard core exclusion, with the addition of some local constraints, which mimic the geometric constraints on the possible rearrangements in physical systems. Similar behavior is quite likely to arise in microstructure formation.
The use of continuum models derived from elasticity theory [6,7,8,9,10] has been proposed as another approach to microstructure formation. Some of these models have been successful in explaining microstructures, hysteresis and avalanches. Nevertheless, they are very time consuming from a computational point of view. For this reason, in many cases only 2D problems have been addressed, and even this item presents difficulties connected with large statistics and with the scanning of the model parameters in order to study their influence. Our aim here is to find a statistical mechanical lattice model, easy to simulate and which allows for the study of statistics of microstructure sequences dynamically generated in athermal systems, under the influence of disorder.
The Random-Field Ising model (RFIM) [11] with metastable dynamics is one of the simplest models for the study of the combined effects of disorder and athermal evolution. It is formulated in a magnetic language for a spin reversal transition, driven by an external magnetic field H. The only degrees of freedom are spin variables defined on a lattice (i = 1, ..., N ), which take values S i = ±1 on the i-th lattice site. The RFIM enables computation of hysteresis loops M (H) corresponding to the behavior of the order parameter M = i S i as a function of H as well as the analysis of the intermediate states between the negatively saturated initial state {S i = −1} and the final positively saturated state {S i = +1}, and vice versa. In particular, the RFIM has been successful in understanding the avalanche dynamics (Barkhausen noise [12]) that joins the intermediate metastable states and shows absence of characteristic scales for a critical amount of disorder. Moreover, the RFIM displays several other interesting properties [11,13]: it exhibits a well-defined rate-independent trajectory, it shows return point memory, it satisfies the abelian property and, from a computational point of view, it is fast to simulate trajectories in relatively large systems [14].
Nevertheless, the usefulness of the RFIM for the study of microstructures is almost null. This happens because the parent and product phases in the RFIM are the positively magnetized phase and the negatively magnetized phase. These two phases are single variant and totally equivalent from a symmetry point of view. Consequently, the obtained hysteresis loops are symmetric under the exchange H → −H and M → −M , there is absence of latent heat associated with the FOPT and the domains are spherically symmetric (except for some short-range correlations due to lattice symmetries).
Within this framework, the aim of the present work is to explore some modifications that should be introduced in the RFIM in order to obtain 3D microstructures without losing, as far as possible, some of the useful RFIM properties that we have mentioned above. A first step in this direction was done several years ago by defining the Random Field Blume-Emery-Griffiths model, in which the 'spin' variables take three different values (S i = −1, 0, 1) with metastable dynamics [15]. In this case, the FOPT takes place from a single variant parent phase, represented by S i = 0, and a product phase with two variants S i = +1 and S i = −1. In the cited work, the hysteresis loops, phase diagram and avalanche distribution were studied for this type of simple case. In the present paper we would like to go one step forward. This will be done by starting from the Random Field Potts model [16] with metastable dynamics. In the Potts model the spin variables can take an arbitrary number of values. This model allows phase transitions to be studied from a non-degenerate phase S i = 0 to a multivariant product phase. We will explore the effect of an extra interaction term (of a dipolar nature, but truncated to the nearest-neighbor approximation), which will be necessary in order to produce microstructures (for the introduction of dipolar interactions in RFIM, see [17]). It is not our aim to focus on the detailed analysis of any particular transition. We will study a model that, from the point of view of symmetry, would correspond to a transition from a cubic phase (single variant) to a tetragonal phase (with three equivalent variants), although neither the detailed interactions nor the external constraints will be tuned for the particular modeling of such transitions in ferroelastic systems. This will be the aim of a future work [26].
The paper is organized as follows: in section II we introduce the hamiltonian of the model. In section III we detail the metastable dynamics that has been used for the simulations. Section IV is devoted to the discussion of the obtained results: in section IV A we present our analysis on the shape of hysteresis loop cycles as a function of the model parameter values. In this section it is shown that the loops happen to be unsymmetric, in contrast to the RFIM results. The microstructures are analyzed in section IV B, where we discuss three different regimes corresponding to different ranges of the parameters values. Moreover, in section IV C we present the statistical analysis of avalanche behavior. Finally, we summarize and discuss the future perspectives in section V.

II. MODEL
The model can be defined on any regular lattice. We will consider a simple cubic lattice of size N = L × L × L with periodic boundary conditions. At each lattice site we define a variable S i (i = 1, . . . N ), which can take four different values that we will call 0, x, y and z. We can choose different representations for our variables but it is convenient to consider a vector S i having three components: we will indicate the four possible values as 0 = (0, 0, 0), x = (1, 0, 0), y = (0, 1, 0) and z = (0, 0, 1).
The order parameter for the phase transition under study here is M = i ( S i ) 2 , where the sum spans over the whole lattice. M represents the amount of the system that has transformed from the cubic to the tetragonal phase. By following the analogy with the magnetic case, we will refer to M as the total magnetization of the system. Moreover, we define the normalized magnetization m, as m = M/N . We will drive the system by an external field H coupled to M , since we are interested in the transition from the 0 phase to the multivariant phase which will be composed of regions (variants) in the states x, y and z. The field H would correspond to the driving effect of the temperature in athermal structural transitions. We will start by decreasing H, from the M = 0 state. We consider the following hamiltonian: The first term is a Potts exchange term extending to nearest-neighbor (n.n.) pairs. The parameter k will be always positive, favoring the ferromagnetic interaction between spins in the same state. In the following, without loss of generality, we will consider k = 1. The second term is a dipolar interaction truncated to n.n. pairs. As we have already said, the aim of adding this term is to generate some simple microstructures. To include higher order terms would lead to more realistic ones. The vector r ij is the lattice vector joining 'spins' S i and S j . We will study the cases with λ < 0 and λ > 0 separately. As can be easily seen from Eq. 1, in fact, the λ < 0 case corresponds to favoring the growth of prolate (needle-like) domains parallel to the spin direction. On the other hand, for λ > 0 such a growth is not favored, but as it is partially compensated by the exchange term it basically corresponds to the formation of oblate (disklike) domains, perpendicular to the spin direction. We will illustrate these features in Sec. IV B.
The third term of the hamiltonian accounts for the interaction between the system and the external field H. This field will be driven from very high positive values to very negative values (and vice versa) in steps ∆H, mimicking an adiabatic triangular driving force (i.e. field frequency ω → 0). By a deliberate abuse of language we will refer to the step ∆H as the driving rate. One can notice that it is possible to add a second driving term G S i which would be convenient for the study of the transitions from one variant to another, mimicking, for instance, the effect of an applied external stress.
The last term H dis accounts for the quenched disorder of the system. We will restrict ourselves to a random field type with zero averages. However, there are still several possibilities for such a hamiltonian term because the random fields can couple either to S i or to the order parameter ( S i ) 2 . We will thus consider: where g i is a three-component vector random field, whose components are extracted from a gaussian distribution N (0, 1) with zero mean and unitary standard deviation, and h i is a scalar field, again extracted from a gaussian distribution N (0, 1). The parameters σ and ρ control the amount of quenched disorder in the system. In order to compare our model with the standard RFIM, we define the total amount of disorder σ 2 0 = σ 2 + ρ 2 . In practice, the two disorder terms can be understood as arising from a local random field f i , whose components are correlated, being

III. DYNAMICS
There are several possibilities for the choice of metastable dynamics. In Fig. 1 we show examples of hysteresis loops obtained with three possible choices of dynamics. At first sight, the three loops look very similar. In all the cases we start from a metastable state, we increase or decrease the field by a ∆H step and then, at constant field, we recursively minimize the system energy by using a local rule based on single-spin changes. Only after a new metastable state is reached we proceed with a new field change ∆H.
In the extremal selection + extremal update case, we scan the whole system and check which variable S i can change to a new value with the minimum energy difference ∆H. The proposed change is accepted if this minimum ∆H is negative. In the random selection + random update case, we randomly choose a spin on the lattice and propose a random change to a new value. If the proposed change implies ∆H < 0 the change is accepted. Finally, in the random selection + extremal update dynamics, we randomly choose a spin on the lattice and check among the three possible new values which one represents a minimum ∆H. If this minimum value is negative we accept the change. From a computational point of view the first choice is much more time consuming than the other two, since the effort scales with L 3 . Although the loops are very similar, detailed analysis reveals that the obtained hysteresis loops, as well as the sequence of metastable states, are not identical. This tells us that the proposed model is not abelian and that the final state will depend on the order in which unstable spins will be changed. In order to ensure some robustness of the results, we are thus forced to choose extremal selection + extremal update dynamics, that is: to propose the optimal spin change among the whole lattice and among all the three possible final values at each time step. This kind of dynamics is deterministic and thus, by definition, independent of the updating order. We will keep to this dynamics for the rest of the paper.
We now study the effect of changing the value of ∆H. In Fig. 2 we show three hysteresis cycles obtained for different values of the driving rate ∆H using the extremal selection+extremal update dynamics. A detailed analysis reveals that the differences between the three loops can be attributed to the fact that driving with a smaller ∆H allows more metastable intermediate states to be found, but for the same applied field values, in the three realizations not only the magnetization, but also the final microscopic configurations reached are the same. The independence from the field rate is an important prop- erty from the point of view of the simulations, since it allows us to use a relatively large ∆H for the study of the properties of hysteresis loops.

IV. RESULTS
We have performed numerical simulations of systems with sizes L = 8, 16, 32, 40 and 60, averaging over 10 2 − 10 3 realizations of the quenched random fields. We have focused our analysis on hysteresis loops behavior, on microstructure formation and on the statistical properties of the avalanches. We consider the cases with λ < 0 ( Fig. 3(a)) and λ > 0 separately ( Fig. 3(b)), because they show a clearly different behavior. In the first case, which corresponds to the formation of prolate domains, the more negative λ is, the larger the width of the loop. In the case of very negative values, the loops start to exhibit a plateau in the increasing field branch: the retransformation to the 0 phase is done in two separate steps. (This effect will be discussed below). For the second case, increasing lambda towards positive values increases the tilt of the hysteresis loop so that saturation in the transformed phase can only be obtained when the field is very negative. This effect is In Figs. 4 and 5 we show the effect of the two disorder parameters σ (Fig. 4, in the low (a) and high (b) ρ regimes) and ρ (Fig. 5, in the low (a) and high (b) σ regimes). In all the cases it can be seen that increasing the amount of disorder increases the tilt and decreases the width of the loop. Moreover, as expected, for low values of the amount of disorder (σ or ρ) the loops exhibit sharp discontinuous (ferromagnetic-like) behavior. This feature is in agreement with recent results on the standard RFIM, concerning the observation that the transition from sharp to smooth loops can be induced by different kinds of disorder parameters: not only the random field variance σ but also random anisotropy [18], the vacancy concentration [19], etc.. Our model shows that the correlation with the random fields of intensity ρ can also act in a similar way.
Let us now discuss the plateau observed in Fig. 3(a) in the increasing field branch. As shown in the examples in Fig. 6, this plateau occurs at smaller magnetizations when the system size is increased. This suggest that it may be due to the stabilization of "slab" domains that cross the whole system from one face to the other and that due to the periodic boundary conditions behave as infinitely large. Such slabs become less and less frequent by increasing the system size. This suggestion has been confirmed by analyzing sequences of configurations. An example will be discussed in Sec. IV B.
In order to perform a quantitative analysis of the hysteresis loops it is important to measure some of their properties. One of the most studied hysteretic features is loop area. In fact, it represents the amount of energy dissipated during a cycle and thus it is an important quantity to be controlled both from the theoretical and materials application point of view. In Fig. 7 we show the loop area, averaged over several disorder configurations, as a function of the parameter λ, for various values of σ.
As can be seen, the area shows a much more important dependence on λ for negative than for positive λ. This behavior can also be seen by studying the coercivity (amplitude) of the hysteresis cycles at m = 0.5, which displays a very similar dependence on λ.
In Figs. 3, 4 and 5 we can see that the hysteresis cycles obtained with our model are asymmetric, i.e. the decreasing field branch (transformation) cannot be related by an inversion operation to the increasing field branch (retransformation). This is an interesting property since, experimentally, materials displaying a transi- tion to a multivariant phase show such behavior, which cannot be reproduced with the RFIM (see section I). In our model, this feature descends from the intrinsic difference of the physical processes occurring in the two branches (transition from the 0 state to the three variant phase in the first branch, and the opposite process in the second branch). In order to study this feature more quantitatively, we define an asymmetry factor A as: where (dM/dH) 1 and (dM/dH) 2 are the derivatives of the hysteresis cycle at the coercive fields (defined as the two fields at a height m = 0.5) in the two branches. As can be seen in Fig. 8, A is greater than zero for negative values of λ, while for λ > 0 the effect of the 'dipolar' term is screened by the disorder and the Potts term and A is essentially 0, irrespective of the value of σ. We have cut the curve with σ = 3 in Fig. 8 at the value λ = −7. In fact, for the considered range of parameter, lower λ hysteresis curves begin to show the plateaux explained above, and thus our definition of asymmetry loses sense.
The effect of disorder on hysteresis loop properties is illustrated in Fig. 9, where we show the width W of the loops at m = 0.5 (coercivity) as a function of increasing correlation ρ, for different values of σ. As mentioned in

B. Microstructures
As we have already mentioned (see section II), when the dipolar term is large enough compared to the Potts term, the transformed domains lose their spherical symmetry and start to show a non-trivial microstructure.
The microstructure of a system is defined as the arrangement of the variants of the product phase. In Fig. 10 we can see an example of these threedimensional microstructures. We represent the views of the yz, xz and xy surfaces, when the sample has reached saturation (fully transformed state). In this case (λ < 0), the domains tend to be prolate. For instance green domains (corresponding to S i = (0, 0, 1)) tend to be elongated along the z direction both in the xz and zy planes. This effect can be quantitatively measured as will be explained below.
For λ > 0, we observe the formation of oblate domains, as shown in Fig. 11, developing in the plane perpendicular to the spin direction. This effect generates a sort of "chessboard" correlation as can be seen, for instance in the yz plane by observing the red and green domains. In order to quantify the shape of the domains in such microstructures, we have calculated the average linear size D x , D y and D z , of the domains of the three variants x, y and z, along the three spatial directionsx, y andẑ, at the saturation configuration. For instance, the average size matrix corresponding to Fig. 10 and Fig. 11, are shown in Tables I and II. In the first case (corresponding to the prolate domains), the diagonal elements of the matrix are sensibly larger than the others, confirming the growth tendency of domains along the orientation of each variant. As is quite obvious, for symmetry reasons, the diagonal elements can be averaged giving what we will call the average linear size in the parallel direction D and the off-diagonal elements can also be averaged giving the average linear size in the perpendicular direction D ⊥ . For the case of Fig. 10 (λ < 0) we get D = 4.44 ± 0.20 and D ⊥ = 2.19 ± 0.05. For the case of Fig. 11 (λ > 0), the tendency of the system to form oblate domains is confirmed by the values D = 1.175 ± 0.005 and D ⊥ = 1.89 ± 0.03.
For the sake of completeness, as a final microstructure example, in Fig. 12 we show a system configuration with high disorder and λ = 0. As expected, no domain asymmetry arises and every spin just aligns with its local random field. The values of D = 1.51 ± 0.14 and D ⊥ = 1.50 ± 0.20 are equal to within statistical errors.
The formation of microstructures such as those in Fig. 10 and 11 is quite clearly affected by the dynamics due to the effect of kinetic constraints. In fact, when a domain of one variant starts to grow, it necessarily blocks the growth of neighboring domains of other variants and vice versa. Thus the first variant to locally break symmetry will facilitate the nucleation of large domains, thus restricting the dimensions of the other variants. This effect could be seen, for instance, by analyzing the decreasing field branch in Fig. 13: the formation of the domains of type x and y that cross the system block the growth of domains of the z variant.
Moreover, with the help of the microstructure representation, we can analyze in more detail the bump formation due to finite size effects, discussed in section IV A (see Fig. 13): in fact, if the system size is finite there is the formation of domains spanning a whole system side, such as the y domains in microstructures labeled with 4 and 5 in Fig. 13. As already pointed out, these kind of slab domains are actually infinite due to the periodic boundary conditions and are thus very stable. As we can see in the figure, they keep existing up to high driving field, giving rise to the existence of the plateau and disappear when the field overcomes a certain threshold.
As we already mentioned in section II, the truncation of the dipolar term does not allow elastic effects to be reproduced, which would lead to more realistic microstructure. In particular, it would be interesting to control the tendency for the different variants to exhibit a preferred habit plane, as would be the case of a real cubic to tetragonal transition. We expect that including a nextnearest-neighbor dipolar interaction will allow for such an interesting property to occur. The model presented here is, therefore, a promising starting point for modeling materials with a phase transition from a single variant to a multivariant phase.

C. Avalanches
Another phenomenon that may be analyzed with our model is the avalanche dynamics. The hysteresis loops in athermal first-order phase transitions consist of a sequence of jumps between metastable states. Such discontinuities are, in general, microscopic. However, for certain values of the disorder one or more may become comparable to the system size and then correspond to the observed macroscopic discontinuities in the ferromagnetic loops. In the magnetic case, microscopic avalanches can be detected experimentally by appropriate coils. They correspond to the so-called Barkhausen noise [12]. In structural transitions, avalanches can also be detected typically by acoustic emission techniques [20]. Knowledge of the distribution of the number of avalanches along the transition, as well as their size and duration, is an important piece of information in order to characterize athermal FOPT.
Good discrimination of individual avalanches in the simulations can only be performed in the limit of ∆H → 0. This will require a true adiabatic simulation algorithm which is not available at present, as opposed to the case of the standard RFIM [14]. In our case, after a small (but finite) change ∆H some spins can be updated. We will consider all of them as being part of a single avalanche. This is an approximation and, therefore, we should keep in mind the fact that we are slightly overestimating the size of the observed avalanches due to some unavoidable overlaps. In the experiments recording avalanches the same problem occurs [21,22].
With this definition we can study, for instance, the average number of avalanches N av as a function of the external field H. As an example, in Fig. 14(a) we present the behavior of this number for the case λ = −8, σ = 4, ρ = 0, compared with the behavior of the average hysteresis loop in Fig. 14 (b). N av presents two peaks in correspondence with the two coercive fields and goes to zero (as expected) at the m = 0 and m = 1 saturations. This kind of information is very relevant for the understanding of the measurements of acoustic emission in structural transitions using the pulse-counting technique [23].
More interesting information can be obtained by measuring the avalanche sizes S and computing the integrated probability distribution P (S), by analyzing all the avalanches in a single branch of the loop (the two branches must be analyzed separately since they are not symmetric). As a naive approximation, in our case one can define in our case the size S of the avalanche as the order parameter variation ∆M associated with an avalanche (i.e. when the field is varied by ∆H.) Nevertheless, this definition imported from the standard RFIM should be carefully adapted to our multivariant FOPT. Inside an avalanche, in fact, one can distinguish between different kinds of processes taking place, depending on their effect on the order parameter variation. Let us focus on the decreasing field branch starting from the m = 0 phase up to the m = 1 saturated configuration: there are several microscopic possibilities for a spin jump. A spin could jump from the 0 state to one of the three variants x, y and z thus giving rise to a positive contribution to the magnetization change ∆M ; it could jump from the x, y or z states to 0, causing a negative contribution ∆M < 0; or finally there is a possibility of a change from one variant to another without contributing to the change in the system magnetization. These three possible updates will be called +, − and 0. Instead of only measuring the total size ∆M of an avalanche, for each of them we will measure the three quantities n + , n − and n 0 , which are the number of spin updates of each kind. Moreover, since as we have seen in section IV A, the hysteresis loops are not symmetric, we should separately analyze the avalanches during the decreasing field branch and during the increasing field branch. This gives, therefore, 6 possible distributions: P down + (n), P down − (n), P down 0 (n) for the decreasing field branch and P up + (n), P up − (n), P up 0 (n) for the increasing field branch. It should be mentioned that in the increasing field branch the total number of events of the + and 0 kind are much smaller than the number of − events, typically by 2-3 orders of magnitude. For the decreasing branch the − events are rare, but 0 events are frequent, since an important number of transitions within variants may occur during the avalanches in the later stages of the transition. Of course such an optimization between variants cannot occur in the reverse process during the increasing field branch.
In Figs. 15 and 16 we show some examples of the probability distributions for varying values of the two disorder parameters on log-log plots, respectively σ and ρ. Actually, we show only the P down + (n) and the P up − (n) distributions which account for the majority of the events. In Fig. 15(b) it is possible to notice another finite size effect: the slab domains of size corresponding to multiples of the system size L present the tendency to disappear abruptly (and thus the P − up (n) presents some peaks in correspondence of values multiple of L), as we have already pointed out in section IV B.
An interesting result for the case of ρ = 0 (see Fig. 15(a)) is that the distribution P down + (n) shows a exponentially damped behavior but seems to become a power-law for a certain critical value of σ (which will be located around σ ≃ 3). Such a tendency has been well studied in the standard RFIM. The fitted value of the power-law exponent is 1.1, clearly smaller than the expected universal value 2.03 ± 0.03 [24], but this feature could be due to finite size effect. Only a detailed finite-size scaling analysis and the study of the geometry of the avalanches [25] (out of the scope of this paper) will reveal if the observed power law conforms to universal behavior or not. Interestingly it seems that for the increasing field branch (see Fig. 15(b)), the distribution is always exponentially damped, at least for the studied range of values of σ. Therefore, the critical point for the increasing field branch, if it exists, would be located at a different (smaller) value of σ.

V. SUMMARY
The analysis of microstructure formation in ferroic systems undergoing a first-order phase transitions is an interesting issue both from a purely theoretical and an applicative point of view. Microstructures arise since the product phase, arising from the balance of many energetic terms, may show energetically equivalent variants. Despite the interest in this issue, the models that have been used up to now for the study of the interplay between disorder and athermal evolution (for example, the Random Field Ising model [11]) are not suitable for the analysis of microstructure formation, due to the equivalence of the variants of the product phase.
In the present work, we have introduced a modification of the Random Field Potts model, which consists of adding a dipolar term truncated to the nearest-neighbor approximation, which represents a promising step to- wards the analysis of athermal transitions from a degenerate to a multivariant phase.
In our simulations we have chosen extremal updating in order to preserve the independence of the trajectory from the applied field rate. We have studied the dependence of the hysteresis loop shape on the hamiltonian parameters values. From a quantitative point of view, this has been performed by measuring the the loop area, its asymmetry and width. Our loops display a large area and asymmetry regime for very negative values of the 'dipolar' term parameter λ, associated with the formation of microstructures with prolate domains, oriented along three equivalent spatial directions. On the other hand, for λ > 0 domains are planar (oblate) and the loops show low, almost constant values of the area and the asymmetry.
We have also addressed the study of the avalanches in the hysteresis loop, distinguishing between the two branches. For certain parameter values, the probability distribution of the size of the leading process in an avalanche displays power-law behavior, which is the typical result for athermal phase transitions in disordered systems.