Modeling Premartensitic effects in Ni2MnGa: A Mean Field and Monte Carlo Simulation Study

The degenerate Blume-Emery-Griffiths model for martensitic transformations is extended by including both structural and magnetic degrees of freedom in order to elucidate premartensitic effects. Special attention is paied to the effect of the magnetoelastic coupling in Ni2MnGa. The microscopic model is constructed and justified based on the analysis of the experimentally observed strain variables and precursor phenomena. The description includes the (local) tetragonal distortion, the amplitude of the plane-modulating strain, and the magnetization. The model is solved by means of mean-field theory and Monte Carlo simulations. The results show that a variety of premartensitic effects may appear due to the magnetoelastic coupling. For large values of the magnetoelastic coupling parameter we find a premartensitic first-order transition line ending in a critical point. This critical point is responsible for the existence of large premartensitic fluctuations which manifest as broad peaks in the specific heat, not always associated with a true phase transition. The main conclusion is that premartensitic effects result from the interplay between the softness of the anomalous phonon driving the modulation and the magnetoelastic coupling. In particular, the premartensitic transition occurs when such coupling is strong enough to freeze the involved mode phonon.


I. INTRODUCTION
Many metals and alloys undergo a so called martensitic transition (MT) from an open cubic phase at high temperatures to a more closed packed phase at lower temperatures [1]. It is a displacive, diffusionless, first-order phase transition, accompanied by incomplete softening of certain transverse phonon modes. For Zr, which belongs to an ideally simple class of martensitic materials, the pure group IV-metals [2,3], it was demonstrated that the first-order character can be understood as an effect of a coupling between two simultaneous strains: an internal two-plane shuffle strain and a uniform strain [4].
A rich variety of precursor phenomena [5] have been observed in (weakly) first-order MT. Some of them, as the intermediate tweed structures [6,7], are not common to all materials but others, intimately related to the transition mechanism, are present in almost all bcc systems studied so far. The most significant is the anomalously low T A 2 [110] phonon branch ([110] propagation, [110] polarization]), accompanied by a low value of the elastic constant C ′ = (C 11 − C 12 )/2. Moreover, both the phonon branch and the corresponding elastic stiffness soften with temperature. Recently, a lot of interest has been focused on the intermetallic Ni-Mn-Ga alloy close to the stoichiometric composition Ni 2 MnGa. It is the only known ferromagnetic fcc Heusler alloy undergoing a MT on cooling. Besides its theoretical interest, it may be of technological importance too since it opens the possibility of controlling its shape memory properties (intimately related to the MT) by applying an external magnetic field [8]. It has a transition from a bcc (neglecting the atomic order) to a low temperature tetragonal phase bct, which is modulated by a five-plane shuffle strain [9,10]. Particularly intriguing is that for nearly stoichiometric composition the full MT is preceded by an intermediate phase in which apparently only the shuffle strain is activated, but not the tetragonal strain [11,12,13,14]. This intermediate phase consists in a micromodulated domain structure, without resulting macroscopic tetragonal deformation so that the cubic symmetry is preserved [14]. This is accompanied by a significant, although not complete, softening of the T A 2 phonon branch at a wave vector ξ 0 = 0.33. Only at lower temperatures, at the martensitic transition point, the homogeneous tetragonal strain is activated (and the modulation changes slightly). This particular behavior observed in Ni 2 MnGa seems to be related to the influence of magnetism. The point is whether or not this intermediate modulated structure is a genuine new phase. Different behaviors have been reported in the literature. For some samples [11,13,14] there exists evidence for a true phase transition of very weak first-order which is driven by a magnetoelastic coupling. The main proof of that is the fact that the intermediate transition (IT) shifts with the external applied field [11] while no (significant) magnetic field dependence has been found for the MT temperature [15,16]. In other studies [17], the authors could not find any indication for an IT although precursors, clearly related to the magnetization of the sample, have been observed. Apparently, the only relevant difference in the samples used by Planes et al [11], Zheludev et al [13,14] and Stuhr et al [17] is the content in Mn. Very recently, it has been suggested [8] that the tetragonal phase can be suppressed by increasing the concentration of Ni at expenses of the content of Mn. Moreover, the magnetoelastic effects have been confirmed from the experimentally observed dependence of the elastic constants on an external magnetic field [18]. In spite of the experimental evidence for the magnetoelastic coupling in Ni 2 MnGa, its microscopic origin has not been established yet.
In the present paper we have theoretically investigated the nature of the bcc to bct transition and constructed a model in order to solve the puzzling behavior and to elucidate the role of the magnetic coupling, which counterintuitively seems to couple the ferromagnetic order stronger to the modulating strain (η) than to the uniform tetragonal strain (ǫ). We shall demonstrate that the situation can be described by the degenerate Blume-Emery-Griffiths model (DBEG) [19] extended to include coupling to magnetic degrees of freedom, and with a interpretation of the variables, appropriate to the present case.
The plan of the paper is the following: First, in section II, we provide an analysis of the experimental facts and a theoretical explanation of the observed phenomena, which is used to justify a microscopic model presented in section III. The model is studied using first mean-field theory (section IV) and next Monte Carlo simulation (section V). The discussion from the comparative study is presented in section VI where we provide also our conclusions taking into account the available experimental data.

A. Experimental Facts
The structural properties of Ni 2 MnGa have been investigated in a series of papers [9,10,11,12,13,14,17,20]. At high temperatures the alloy has the fcc (L2 1 ) Heusler structure which, neglecting the atomic order, can be regarded as a bcc lattice. It is paramagnetic at high temperatures with the magnetic moments mainly on the Mn sites [21]. At temperatures below T m , it orders ferromagnetically with no particular easy direction of the moment. At the temperature T M (< T m ) there is a (first-order) structural phase transition of the martensitic type to an average tetragonal structure, which additionally is modulated by a transverse five-layer shuffling strain. Prior to this transition precursor structures of that phase as well as of the fcc having a six-plane modulation have been observed in neutron scattering experiments [14,17]. This may happen also as a transition (first-order) at a temperature T I (T M < T I < T m ) giving rise to a genuine intermediate phase [22] without any macroscopic tetragonal deformation [14].
The above temperatures T m , T M and T I are extremely sensitive to the composition and atomic ordering of the sample. Thus T m may vary between 360K to 395K, whereas T M may vary from 175K to 450K. In the sample studied by Stuhr et al. [17], T m ≈ 364K and T M ≈ 284K with precursor signs for T > T M , but no intermediate phase was found. In the sample studied by Zheludev et al. [13,14] T m ≈ 380K and T M ≈ 220K and an intermediate phase is observed below T I ≈ 260K. Similarly, Planes et al [11] studied a sample with T m ≈ 381K and T M ≈ 175K and found an intermediate transition at T I ≈ 230K. In spite of the large variation in the temperatures, we assume the basic physics is the same, and hence we shall use the information obtained by observations on different samples. In particular, we shall attempt to explain the precursor and the intermediate phase phenomena by constructing an effective microscopic Hamiltonian, which allows analysis beyond a mean-field or Landau expansion treatments [8,11,23].

B. The two-strain model
Although Ni 2 MnGa is a metallic alloy for which both the structure and the magnetism is determined by the conduction electrons, it is instructive to consider a model system for the mismatch between cubic crystals, similarly as it is found in KCl grown on NaCl (001) [24,25]. In this case, the very large (∼ 17%) lattice constant mismatch makes the interface to buckle simultaneously in the [110] and [110] directions. This produces superstructures consisting of seven NaCl and six KBr layers, or multipla thereof, perpendicular to the interface, while preserving the square symmetry. Both the NaCl and the KBr crystals are modulated since they have similar elastic properties. This gives rise to superstructure peaks in the X-ray and Helium scattering spectra, precisely as those observed in the intermediate phase of Ni 2 MnGa, where one also would expect more higher order satellites in off-symmetry directions to be found, but not so far looked for experimentally. The buckling gives rise to a variation of the local [001] direction and a change in lattice constant perpendicular to the interface (corresponding to a local tetragonal distortion). We stress that for the ionic crystals the forces by no means triggers the local modulating strains, which are caused by the forced contact between unequal crystals. For Ni 2 MnGa the situation is rather different. Here, a nesting feature of the Fermi surface causes a strong electron-phonon coupling and an incipient soft phonon mode at q = ξ, ξ, 0 positions in the fcc phase, where ξ ∼ 1 3 [14,17]. This is presumably the driving mechanism in Ni 2 MnGa, and giving, as a consequence, the tetragonal distortion. As precursor phenomena, quasi elastic peaks are observed at ξ = 1 3 and ξ = 1 6 [17] corresponding to a similar six-plane modulation as discussed for the alkali salt interface. To match this to the fcc (001) lattice plane it is advantageous to make a lattice mismatch and expand the lattice and to rumple the interface, or in other words making the [001] direction fluctuate in epitaxial-like angles. Stuhr et al [17] have shown that the six-plane fcc modulation corresponds precisely to a five-plane modulation of the tetragonal phase, and the latter is found as a precursor phenomenon. An 5:6 expansion would be very drastic. However it suffices to create a superstructure of a common divisor, for example a 30 or 60 plane repeat distance. The latter would correspond to a lattice mismatch of 1.67%, which is very close to the mismatch observed between the fcc and tetragonal phase of Ni 2 MnGa: ∼ 1.6% [21]. Hence we argue that the electronically driven six-plane modulation in turn also causes the tetragonal distortion. The theory for why it is advantageous for mismatched crystals at epitaxial interfaces to develop mutual superstructure peaks was discussed in more detail by Vives and Lindgård [26]. The modulation occurs simultaneously in the q = [ 1 3 , 1 3 , 0] and q = [ 1 3 , − 1 3 , 0] directions, thus preserving the square symmetry of the (001) plane and yielding a modulated (001) plane, rather than a direction. In the fcc structure there are three such equivalent planes (100), (010) and (001). In order to (essentially) preserve the volume, V , a change of which is expensive for the electrons, the c-axis [001] must shrink (from T = 4.2 K to 295 K it is indeed observed that V 1 3 only increases by 0.5% [21]). Therefore a model must include a coupling between the plane modulation η and the tetragonal strain ǫ. Then, the question is: how can there be apparent separate temperatures for the onset of ordering of the two kind of strains? This is possible if in the intermediate phase the tetragonal strain is only local, varying in direction and also in the allowed directions in the fcc crystal. Moreover, only in the martensitic phase microcrystals with a resulting tetragonal strain should be formed. The corresponding tetragonal structure is observed as highly mosaic with a large variation of the [001] directions [17], in agreement with the above picture.

C. The Landau models
Very recently, Landau models for the MT in cubic ferromagnetic materials have been proposed [8,23]. In these models, the magnetoelastic coupling between the uniform strain tensor and the (vector) magnetization is fully considered. Nevertheless, they seem more appropriated to the study of the magnetic properties of the martensitic phase rather than to the analysis of the IT itself. Here, accordingly with the discussion given above, we shall adopt a different strategy and study the structural transition in Ni 2 MnGa in terms of a Landau expansion of the only most relevant strain and magnetization variables, including non-uniform strains or modulations. Then, similarly to the case for the bcc to hcp transition in Zr [4], we include the following terms, which are allowed by symmetry: where η is the discussed plane modulation strain and ǫ the local uniform contraction perpendicular to that plane, but we do not consider higher order uniform strain terms. Here ω 2 s = a(T − T s ) is the squared frequency [27] for the incipient soft mode phonon with q = 1 3 , 1 3 , 0 , other constants are positive, and C ′ = 1 2 (C 11 − C 12 ), which is small and temperature dependent. By eliminating the local tetragonal strain ǫ we can [4], write the free-energy along the optimum energy path involving both ǫ = −∆ 2 / C ′ η 2 and η as The coupling between the two strains therefore makes it possible forB to become negative and hence to cause a first-order transition at T M before the soft mode transition at T = T s , even without the coupling with the magnetism. Next let us consider the influence of the magnetism. A ferromagnetic moment m can influence an itinerant magnet as N i 2 M nGa in two ways. Firstly, giving rise to a splitting of the electronic energy bands, which is proportional to the amplitude |m|, but which is not sensitive to the moment direction. This will cause a change in the soft phonon frequency below T m proportional to |m| 2 : where a, b and u are positive constants, and we have assumed a mean-field behavior for |m| 2 . A measurement of the soft mode frequency-squared should therefore exhibit a kink at the magnetic ordering temperature. This is exactly the behavior observed by Stuhr et al [17], who found a = 0.018 meV 2 /K and b = 0.020 meV 2 /K. Unfortunately for the samples showing the IT [11,13,14] measurements in the paramagnetic phase are not available. Nevertheless, in what follows, we shall assume the behavior discussed above for all the samples. Hence the magnetic free-energy part can be written as: where α = A(T c − T ) and T c is the magnetic transition in absence of magnetoelastic coupling. A and β are positive parameters and γ = −µu is yielding the coupling between the amplitude-squared of the magnetization and the effective modulating strain. By eliminating |m| 2 one can write an effective free-energy in the form of the eq.(2). This yields a further temperature dependence of ω 2 s andB, as discussed by Planes et al [11]. The other possible coupling between the magnetization and the structure is by magnetostriction, which deforms the crystal in the direction of the magnetization [8]. Experimentally, the easy direction of magnetization is not known with certainty, even in the tetragonal phase. Webster et al [21] propose it might be in the 111 directions of the L2 1 phase [28], but that other directions are almost as likely. Analysing their data perhaps allows the conclusion that it is at least not in the tetragonal [001] direction, which would normally have been the obvious choice. If so, it is confined to the (001) plane, say, which has four-fold symmetry and therefore not yielding a strong easy axis. The lowest order coupling would then be of the form + m 2 ǫ 2 . An interesting possibility is if the easy direction is along the [100] direction, because this may belong to two different modulation planes, (001) and (010), therefore yielding a minimum magnetostriction because it cannot distinguish between a [001] and a [010] tetragonal strain. The effect is further reduced because of the equivalence between moments in the [100] and [010] directions in the modulated plane (001).
Inclusion of the coupling can be done in the eq.(4), with no qualitative change (except for an induced dependence on the magnetization of C ′ [11]), and hence we shall neglect it in the following. If the coupling is sufficiently strong, it would no doubt prevent the existence of an intermediate phase in Ni 2 MnGa.
The conclusion of the above analysis is that the phases may be characterized by the minimum path modulation strainη in eq.(2) which includes a finite, but local tetragonal strain which is on average zero because of fluctuating directions. This is consistent with the observation of a significant broadening of the fcc (002) peak in the neutron scattering data [17]. This strain is coupled to the magnetization such that a ferromagnetic moment on neighboring sites will favor a modulation also at these sites, due to a change in the local band structure, but without differentiating between the three possible planar modulations and accompanying tetragonal local strains. In order to give meaning to the notion of a local band structure, we consider sites not as the Mn atom positions, but at a more coarse grained level, e.g. at the unit cell level, say. At the MT a resulting tetragonal strain, ǫ, may arise as the interfaces between the modulated fcc and the modulated and locally tetragonal crystallites grow larger and thicker. This strain is on the average in the 001 directions but with considerable variations in the epitaxial angles relative to those. To construct a simplified model, we shall consider one with only two variables such that each has only two degrees of freedom. One variable with only two values is representing the plane modulationsη, and the other also with only two values representing the resulting tetragonal strain ǫ. This does not correspond to considering the behavior of a subspace of variants, but constitutes a simplified statistical analogue to the physics of Ni 2 MnGa.

III. THE MODEL
The desired degrees of freedom can be represented by the p-degenerate Blume-Emery-Griffiths Hamiltonian (DBEG) [19]. It was first introduced with the aim to account for the entropy stabilization of the high temperature phase in martensitic transitions. Recently, it has been shown [29] that it is equivalent, with respect to universal properties, to that of the ordinary BEG model with the crystal field shifted by a term k B T ln p.
Although the DBEG model is very simple it includes most of the relevant physical ingredients to understand MT, namely a multivariant low temperature deformed phase and a high temperature average cubic phase with enhanced entropy. The model is an extension of the ordinary three-state Blume-Emery-Griffiths Hamiltonian defined on a lattice, which we shall take as simple cubic (or square), as motivated above. On each lattice site i = 1, .., N a variable σ i = 1, 0, −1 represents the deformation state near each site on the lattice. The state σ i = 0 represents the undistorted phase, and it is chosen to be p-fold degenerate (p ≥ 1), in order to approximately account for the high entropy of vibration of the cubic phase. The states σ i = ±1 represent the distorted phase. The Hamiltonian accounting for the energy gain in having the same structure on neighboring sites was written as [19] where the sums are performed over all nearest-neighbor pairs. In what follows we will take J > 0 as the unit of energy and work in terms of reduced magnitudes defined as H * = H/J, K * = K/J . . .. Two order parameters can be defined, φ 1 = σ i /N and φ 2 = σ 2 i /N . The model in (5) was solved, for K * ≥ 0, by mean-field and Monte Carlo simulation techniques [19]. We found a phase transition from a cubic (disordered) phase with φ 1 = 0 to a tetragonal (ordered) phase with φ 1 = 0. The behavior of the secondary order parameter φ 2 in the disordered phase depends on K * and T * , but exhibits non-analytical behavior at the same temperature as φ 1 . Large values of K * stabilize the cubic phase. Moreover, the K * parameter and the degeneration p of the 0-state control the order of the transition, which changes from being of second order (for low values of K * ) to first-order.
The DBEG model [19] was introduced as a simplification and a further generalization of the Lindgård-Mouritsen model [4], initially designed to mimic the bcc to hcp transition in Zr. Consequently, it naturally inherits the identification of the two order parameters suitable for Zr: φ 1 is the two-plane shuffle strain and φ 2 the homogeneous strain. In the previous work [19] this was not emphasized since both φ 1 and φ 2 exhibit a phase transition at the same point.
For Ni 2 MnGa it is more natural to identify the order parameters in the opposite way. Therefore, in the present work, we define: The breaking symmetry order parameter ǫ corresponds to the tetragonal distorsion in the sense that if ǫ = 0 (equal population of σ i = +1 and −1) it corresponds to having all the variants equally populated, and hence an average cubic phase. For q the relation is more involved because of the complicated physics of the modulating strainη. In the high temperature cubic phase, when the σ i variables distribute at random q = q 0 = 2/(p + 2). Let us assign the difference q 0 − q with the amplitudeη of the plane modulating strain without distinguishing between the three possible modulating planes. We now include the magnetic degrees of freedom by means of spin variables S i = ±1 (defined on the lattice site i = 1, .., N ) having a ferromagnetic Ising interaction. Thus, the purely magnetic contribution is: where J * m > 0. The total Hamiltonian should further include a coupling term between the structural and magnetic variables. We have argued in section II that the magnetic influence of electronic properties gives rise to a coupling between the magnetic moment and the plane modulation. On a microscopic level let us assume that the presence of a moment of neighboring sites gives rise to a modulation on neighboring sites as motivated previously. To describe this, let us consider the following symmetry allowed [30] interaction contributions: that may be rewritten as: In the pure ferromagnetic phase (S i = 1) this Hamiltonian can be viewed as an Ising Model for the variables 1/2 − σ 2 i . Thus a phase transition exists for (U * 11 + U * 00 − 2U * 10 ) > 0 and U * 11 = U * 00 . For simplicity, in what follows we shall take U * 11 = U * 00 = 0 and denote U * = U * 10 < 0. Then, the coupling Hamiltonian becomes: As we shall see in the next section, in the mean-field approximation, the first term becomes of the form of the coupling term in eq.(4) whereas the last term gives a simple modification of the purely magnetic interaction J * m defined in eq.(8). Furthermore, (11) shows that it may be particularly convenient to choose p = 2, which gives q 0 = 1/2. The total Hamiltonian model for Ni 2 MnGa can then be written as with H * M , H * m and H * int respectively given by expressions (5), (8) and (11). We shall demonstrate that it is possible to split up the structural transition into one determined by the order parameter q, which we will associate with the IT and another one, determined by ǫ, to be associated with the tetragonal deformation occurring at the MT.

IV. MEAN-FIELD TREATMENT
In this section we solve the presented model (12) by using standard mean-field techniques. The state of order of the system depends on the occupation numbers N S σ . This stands for the number of points in the structural σ = −1, +1, 0 and magnetic S = +, − state. There are six different occupation numbers which should fulfil the following normalization condition: where N is the total number of points in the lattice. We define the following order parameters: The corresponding entropy can be written as: where p ≥ 1 is the degeneracy factor of the 0-state. The mean-field expression of the free-energy per particle is: where T * = k B T /zJ, and z is the coordination number of a given site. Notice that m 1 appears only in the entropic contribution to the free-energy. Standard minimization with respect m 1 renders the following relationship m 1 = ǫ(m−m0) q , which has to be fulfilled at all temperatures. Then, after substitution in eq.(20) we obtain the following expression for the free-energy as function of ǫ, q, m, m 0 : Further minimization with respect the other four order parameters yields to the next set of coupled equations: . Then, the values of K * for which the IT exists are determined by eq. (23). Indeed, by setting ǫ = 0, it follows that for K * (T * ) = −2T * ln(p/2) the order parameter q has a continuous phase transition. The results will be presented for two values of the degeneracy factor, p = 2 and p = 4. In Table I we give the identification of the different phases in relation to the problem of interest here and their corresponding abbreviated notation. Figures  1 and 2 show the temperature behavior of the order parameters for a given value of U * = −3.50 along the path determined by the condition K * (T * ) = −2T * ln(p/2) (coexistence line). In particular, for p = 2 one has K * (T * ) = 0. For the sake of clarity, m 0 (T * ) is not shown. At high temperature, a magnetic transition appears at T * m = J * m + U * /2 from a paramagnetic cubic phase (PC) (ǫ = 0, q = 1/2, m = m 0 = 0) to a ferromagnetic pure cubic phase (FC) (ǫ = 0, q = 1/2, m = 2m 0 = 0). From equation (23) it is easy to see that, for q = 1/2 (T * < T * I ), m and m 0 are not independent but m = 2m 0 . At lower temperatures, the order parameter q separates into two branches q + (q > 1/2) and q − (q < 1/2). This occurs at the critical point T * Ic and separates two different ferromagnetic phases both with ǫ = 0 but with different values of q: a pure cubic phase (ǫ = 0, q = 1/2, m = 2m 0 = 0) and an average cubic phase (ǫ = 0, q = 1/2, m = m 0 = 0). The two branches of q are identified with the tetragonal-like modulation (q + ) and the fcc-like modulation (q − ) in the sense that when all are equally populated the cubic symmetry is preserved. Finally, at a temperature T * M , a martensitic-like transition to a ferromagnetic martensitic (tetragonal) phase (FMT) (ǫ = 0, q = q + , m = 0, m 0 = 0) occurs.
In figures 3 and 4 we show the phase diagram as function of K * . The intermediate transition (IT) is of first-order and it is represented by the dashed line (T * I ) separating, below the critical point (black diamond), the regions with q + and q − . For p = 2 (Fig. 3), this boundary is a vertical straight line and therefore cannot be crossed by sweeping T * at constant K * . We advance here that this is a mean-field artifact. The exact treatment by Monte Carlo simulation will show that this transition line is always bent, even for p = 2. The mean-field solution renders a real first-order IT only for p > 2, as it can be seen in Fig. 4 for p = 4. As an example, in figure 5 we show the temperature behavior of the order parameters for p = 4 and three values of K * around the critical value K * c ; K * < K * c (a), K * = K * c (b) and K * > K * c (c). The insets shows an enlarged view of the magnetization behavior around T * I . In figure 6 we show the location of the critical point (T * Ic ) for p = 2 (K * c = 0) with respect to the magnetic (T * m ) and the martensitic (T * M ) transitions as function of the coupling parameter U * . One observes that, as the strength of the coupling decreases, the critical point approaches the martensitic line in such a way that the IT disappears well before U * becomes zero. This is consistent with some experimental observations [18] indicating that strong coupling is required for the IT to appear. The same behavior is found for p = 4, as shown in Fig. 7. We note that, in this case, the critical points correspond to different values of K * (K * c = −2T * ln p 2 ). Figure 8 shows the location of the critical points (thick line) and the region of first-order IT in the U * -K * plane. For very weak coupling, the IT does not exists. For p = 4, as the strength of the magnetoelastic coupling increases, the IT exists for a larger interval of values of K * above K * c . The magnetic transition is of second order for all values of model parameters studied. The MT is found to be first-order only for low values of the coupling strength for which ǫ and q order simultaneously. Notice that this does not contradict the results obtained previously for the DBEG with U * = 0 and K * ≥ 0. Actually, for p = 4, it was found that the MT transition is discontinuous only for K * > 0 [19]. Thus, the mean-field solution of the present model does not reproduce the evident first-order character of the MT in Ni 2 MnGa. At this point, this could be attributed to the insufficient accuracy of the mean-field treatment but, in the next section, we will see that Monte Carlo studies render the same feature. Thus, it is more plausible that it is due to the incompleteness of the model, addressed preferently to the study of the IT and related effects. Actually, our model description is done in terms of a single modulation strain whereas it is known [31] that the MT and the IT phases have different modulations.
In conclusion, the main effect of the magnetoelastic coupling parameter U * is to generate a critical point (T * Ic , K * c = −2T * ln p 2 ), between the magnetic (T * m ) and the martensitic (T * M ) transitions. It is the end point of a first-order transition line T * I (U * ) that, emerging from the martensitic phase boundary, separates two average cubic (ǫ = 0) ferromagnetic (m = 0) phases with different modulation amplitude (η = q 0 − q): q + (η < 0) and q − (η > 0). Provided the coupling strength is large enough, the IT exists for a limited range of K * (below the critical value) which, in turn, depends on U * .
Before ending this section, we would like to show that the magnetoelastic coupling behind the model under discussion is consistent with the one considered in the Landau expansion (4), which in turn has been inspired by the experiments. Let us consider the simplest case of p = 2 (and K * = 0). Furthermore, we shall assume that m 0 may be approximated by m 0 ≃ (1 − q)m. Although this decoupling (see eq. (17)) is exact only for T * ≥ T * I , it provides the first coupling term between m and q. Indeed, from eq. (21) we see that the magnetoelastic contribution to the internal energy becomes of the form: whereη = q 0 − q = 1 2 − q. As it was already discussed at the end of section (III), the first contribution represents a simple correction to the Curie temperature, and the second is the magnetoelastic coupling.

V. MONTE CARLO SIMULATIONS
In this section we solve numerically the model (12) by using Monte Carlo simulation techniques [32]. Our objective is to find system configurations ({σ i } , {S i }) distributed according to the canonical ensemble probability. The corresponding equilibrium simulations have been carried out using the standard Metropolis algorithm. The changes in the σ i and S i variables are proposed independently and accepted or rejected according to the single-site transition probability W = min 1, e −∆H * /T * . We have used a 2d-square lattice with N (= L 2 ) sites subjected to periodic boundary conditions. Different lattice sizes ranging from L = 20 to L = 100 have been studied. The unit of time is the Monte Carlo step (MCS) and consists in N attempts of changing the σ i and S i variables. The simulations have been carried up to ∼ 30 · 10 3 MCS per site. Runs have been performed starting from two initial conditions: (i) a perfect FMT phase (σ i = 1,S i = 1, i = 1, 2, · · · N ) and (ii) a perfect FC phase (σ i = 0,S i = 1, i = 1, 2, · · · N ). This is very convenient in order to detect metastability and hysteresis when crossing first-order transition lines. Notice that, from the mean-field solution presented in the previous section, we already have an idea of the range of the space of parameters we have to explore. Accordingly, we shall fix J * m = 4.0 and use different values of U * < 0. Concerning the degeneracy factor we restrict the Monte Carlo simulations to p = 2. Nevertheless, we have verified that other values of p > 2 render qualitatively similar results. Most of the simulations have been performed at fixed values of the model parameters (U * and K * ) and sweeping the temperature T * , but few have been performed at fixed T * and sweeping the parameter K * .
The different quantities measured after each MCS are: the internal energy H * , and the order parameters m, ǫ and q. These quantities have been averaged over ∼ 200 configurations taken every 100 MCS and discarding the first 10 4 MCS for equilibration. Such averages will be denoted by · · · . We have computed H * , |m| , q , |ǫ| . Moreover, the specific heat and the susceptibilities associated with the fluctuations of the order parameters have also been measured: All these definitions correspond to intensive quantities. In many cases the specific heat c * has also been obtained from the numerical derivative (1/N )d H * /dT . The agreement between this and the estimation obtained from eq. (27) gives confidence that the equilibration times used are appropriate. From the behavior of the order parameters, the specific heat and the corresponding susceptibilities the phase diagram can be obtained. The phase transitions associated with ǫ, q and m have been determined from the location of the peaks in either the specific heat or in the corresponding susceptibilities [33]. This method is more accurate than to look for singularities directly on the behavior of the order parameters. Moreover, we have checked whether of not the peaks correspond to a true phase transition by studying their dependence with increasing the system size L. As an example, in Fig. 9 we show the temperature dependence of the specific heat for U * = −3.5, K * = 0.15 and four different system sizes (L = 10, 25, 50, 100). The smooth peak at T * ∼ 4.3 does not correspond to a true phase transition since it does not exhibit scaling behavior. The second order phase transitions (at T * M ∼ 2.25 and T * m ∼ 5.2) exhibit peaks which shift and become narrower and higher as one increases the system size L [34]. Besides, first-order phase transitions (T * I ∼ 3.9 ) exhibit sharp discontinuities which, although they can also increase, neither shift nor become smoother with increasing L. In this last case if L is not very large or averages are not taken for long enough MCS, hysteresis may appear. Figure 10 shows three sections of the phase diagram as a function of K * corresponding to three different values of the coupling parameter: U * = 0.0(a), U * = −2.5(b) and U * = −3.5(c). We notice that, for each value of U * , the magnetic (T * m ) and martensitic (T * M ) transitions are almost independent of K * . The overall conclusion emerging from Fig. 10 is that the premartensitic effects are more important the larger the strength of the magnetoelastic coupling U * is. As it was anticipated by the mean-field calculations, its main effect is the showing up of a critical point (black diamond), and a first-order transition line (dashed line) separating two FC phases, q + and q − , with different values of the modulation amplitudeη. Contrarily to the mean-field solution, for p = 2 , the T * I (K * ) is now bent due to σ i σ j fluctuations, as will be discussed in the next section.
The second interesting point manifested by the Monte Carlo results is the existence of large fluctuations close to the critical point which, we stress, do not correspond to true phase transitions. These are revealed by anomalies in the response functions defined in equations (27)(28)(29)(30). In the case of the specific heat, such anomalies appear in the form of smooth peaks that become difficult to resolve as we move away from the critical point. In figure 11, which is an enlarged view of figure 10(c), we denote the position of such smooth peaks by two dotted lines (with white squares) that, from the critical point, extend towards both sides of the FC phase. We have also indicated the metastability limits associated with the IT (points inside triangles). Actually, the position of the IT (denoted by black points along the dashed line representing T * I (K * )) is determined as the middle point of these limiting lines. The metastability limit points have been obtained by performing some runs at constant T * (< T * Ic ) and sweeping K * (either increasing K * from the q + phase or decreasing K * starting from the q − phase) and some others at constant K * and increasing the temperature from the MT phase. The first-order IT, as will be discussed below, is not found by decreasing temperature.
Keeping Fig. 11 in mind, now we shall study the temperature behavior of the specific heat at constant K * and U * = −3.5, above and below the critical point. The corresponding results are shown in Fig. 12. Different phenomenology may be observed when, increasing the temperature, we move from the bottom to the top of the figure.
a) For K * = −0.1, a martensitic transition and a magnetic transition with no sign of an intermediate transition.
b) For K * = 0.12, a martensitic transition, an anomaly (indicated by ↑) due to the proximity of the critical point and a magnetic transition. As we mentioned before, this anomaly is due to fluctuations and appears when crossing the dotted line (with squares) in Fig. 11. c) For K * = 0.14, 0.15 and 0.18, an additional peak (to the left, also indicated by ↑) due to the intermediate transition shows up. As K * increases one observes that the fluctuation peak gets smoother (since we move away from the critical point) while the IT shifts towards lower temperatures. The entire temperature behavior for K * = 0.15 has been previously discussed in Fig. 9. It is difficult to be exactly at the critical point (that we have estimated K * c ≃ 0.130 ± 0.005) but it would correspond to the value of K * at which the anomaly peak begins to split up.
It is interesting to note that such double peak behavior found around the IT has been observed experimentally in Ni 2 MnGa. Figure 13a shows an example corresponding to MC simulations with p = 2, U * = −3.5 and K * = 0.14 The two peaks are found for heating runs only. Due to the metastability of the first-order intermediate phase transition the transition to the q + phase it is not found when cooling. A similar behaviour is found when performing calorimetric measurements as it is illustrated in figure 13b. The corresponding experimental details can be found in Ref. [11]. The lines correspond to thermograms obtained by heating and cooling (as indicated by the leanning arrows). Actually, Fig. 13b corresponds to an enlargement of figure 1 in Ref. [11] that has been reproduced with permission of the authors. Here, additional calorimetric runs are shown in order to reveal the systematic character of such double peak obtained when heating. Indeed,in the original published figure the two peaks are almost unobservable and were not considered by the authors who treat both peaks as a single one. When comparing Figures 13a and 13b one observes that the highest peak occurs in different order. We do not have an explanation for this yet, but it could be related to the calorimeter inertia.
Besides the behavior of the specific heat, it is also instructive to look at the behavior of the susceptibilities. Figure  14 shows the evolution of χ m , χ q and χ ǫ with temperature for K * = 0.15 and U * = −3.5 (this is one of the cases shown in Fig. 12 and discussed in point c)). The smooth peak in the specific heat (marked by ↑) is associated only with fluctuations of q, while at the IT (large peak) the discontinuity in χ q (modulating strain) is accompanied by an increase of the fluctuations χ ǫ (tetragonal homogeneous distortion).
Very recently, measurements of magnetic and transport properties of N iM nGa alloys have been reported [35]. Magnetization measurements as a function of temperature, for small values of the applied field, reveal the existence of premartensitic anomalies at two different temperatures (separated ∼ 20 o ). Only the low temperature anomaly is found when performing resistivity measurements. Nevertheless, it has been suggested [35] that both correspond to different true pre-martensitic transitions. In the light of the present results we suggest that the high temperature anomaly could be a signature of the critical fluctuations. We agree with the authors that more careful experimental studies are needed in order to clarify the results.

VI. DISCUSSION AND CONCLUSIONS
The phase diagram of the model exhibits, qualitatively, the same features in both mean-field theory and Monte Carlo simulations. A given value of the ferromagnetic interaction parameter J * m > 0 determines the distance between the magnetic T * m and the martensitic T * M (< T * m ) transitions. For appropriated values of the parameters K * and U * , the following phases are found, from high to low temperature: paramagnetic cubic (PC), ferromagnetic cubic (FC), ferromagnetic intermediate or premartensitic (FPMT) and ferromagnetic martensitic (FMT). The change from the FC to the FPMT phases occurs below a critical point (T * Ic , K * c ) and it takes place trough a true phase transition only for K * > K * c . The existence of this intermediate phase depends on both K * and U * . First one needs the magnetoelastic coupling U * (−J * m < U * < 0) be strong enough. We obtain that T * Ic decreases with U * whereas T * M remains almost unaltered so that the critical point disappears, on the T * M line, well before U * reaches the value zero (Fig.6). Moreover, provided the U * is adequate, the IT exists for a limited range of values of K * (> K * c ) across the first-oder transition line T * I (K * ). The corresponding order parameter is the modulating strain amplitude and changes from high to low values as one decreases the temperature across the IT. According to the theory of the harmonic thermal vibrations in a crystal, this is consistent with the behavior observed for the phonon frequency [14]. Moreover, in our results, the IT is accompanied by a jump in the magnetization, only visible for p > 2. Since this happens in both mean-field and Monte Carlo solutions we conclude that this has to do with the coupling rather than with the fluctuations. Experimentally, a jump in the magnetization has not been detected so far [21,23].
Some differences are obtained between the Monte Carlo simulation and the mean field solutions. First, in the simulations, the T * I (K * ) line always bends towards increasing K * for p ≥ 2 whereas the mean-field solution gives, for p = 2 a perfect vertical boundary at K * = 0. This is due to σ i σ j fluctuations and may be understood as follows. The internal energy of a pure q − phase (σ i = 0, S i = 1) is given by (see eq. (12)): while the energy of a pure q + phase (σ i = ±1, at random, S i = 1) is: In the mean-field approximation the term between < .. > is neglected and therefore the condition E − = E + gives K * = 0. In Monte Carlo simulations, the fluctuations are present. We can estimate its value by taking the results from the standard Ising model and therefore − < nn σ i σ j >≃ E * Ising (T * ). This function is zero only at T * = ∞. For T * = 0 takes the value −2N and increases monotonously. In the Ising model, at T * = T * c ≃ 2.27, E * Ising = − √ 2. Therefore, the condition E − = E + gives a transition line at K * (T * ) ≃ −E * Ising (T * )/2N . It bends towards K * > 0 values when T * decreases. Remember that in the mean-field approximation a similar behavior is found only when the degeneracy factor of the 0-state is p > 2.
Another point illustrated by the Monte Carlo simulations is the important role played by the fluctuations in describing premartensitic effects. In figure 15 we show, schematically, the K * − U * section of the phase diagram as it is obtained from the numerical results presented in the previous section. The triangle defines the region with a first-order intermediate (premartensitic) transition and it is limited by the line of critical points (thick line) and the line where the IT disappears on the MT line (thin line). The shadow region denotes the zone of large (critical) fluctuations and it has been estimated from the anomalies in the specific heat. It is interesting that, apart from the region of fluctuations, the mean-field solution renders a qualitative similar phase diagram (Fig.8) although shifted to negative values of K * .
Experimentally, the interplay between U * and K * , required for the IT to occur, is determined by the composition of the sample. Indeed, premartensitic effects in the NiMnGa alloy have only been reported for compositions around the stoichiometry (Ni 2 MnGa). Among them, the most important and common to all the samples, is a significant softening of the 1 3 [110]T A 2 acoustic mode with decreasing temperature, accompanying the formation of an intermediate structure that, while preserving the cubic symmetry, is transversely modulated [10,36] along the [110] direction with wavevector 1 3 . Moreover, this intermediate structure may appear through a true phase transition [11,14] or not [17]. We point out that both behaviors are compatible with the present results. In simple words, the samples showing the IT would fall inside the triangle in Fig. 15 while the others not. More precisely, in the samples showing the IT the strength of the magnetoelastic coupling (U * ) is enough to freeze completely the anomalous phonon with degree of softness related to K * . For deeper discussions it is imperative to have a better understanding about what is behind the composition dependence giving rise to the different behaviour in the samples. In particular, experiments in order to compare the different degree of softness of the anomalous phonon are required. Furthermore, it is important to know the strength of the magnetoelastic coupling. For this, measurements in both the ferro and the paramagnetic phases of the temperature behaviour of ω 2 s are needed. An study of the characteristics of the kink around the Curie point would be very helpful. Finally, we point out that samples in the critical region (shadow region in Fig 15) should exhibit a significant increasing of diffuse scattering when decreasing the temperature below the Curie point.
The present model renders, independently of the technique used to solve it, a MT which is continuous for the range of model parameters studied. This would be a serious setback in case we where interested in the properties of the MT itself. For this, models such as discussed in references [8] and [23] could be more appropriate. Here, we have developed a model with the aim to focus on the study of the IT and related premartensitic effects in Ni 2 MnGa. It is based on the assumption, sustained by the change of the slope dω 2 s dT at the Curie point, that the magnetism causes the freezing of the incipiently unstable 1 3 [110] T A 2 phonon and the splitting from the homogeneous strain. With this point of view, the intermediate phase is a precursor of the MT. On the other hand, in view of the fact that the actual modulation of the martensitic phase is different from that of the intermediate phase, some authors have claimed that both phase transitions have to be regarded as independent [37]. If so, additional coupling terms between the modulation and the homogeneous strains are required in order to produce a change in the modulation at the martensitic transition. In this sense, it has been observed that the dip in the T A 2 branch shifts under an external unaxial stress [38]. However, more along with the point of view adopted here, Sthur et al. [17] have demonstrated that the [ 1 From the above discussion is follows that, in spite of the appealing results obtained, the present model needs to be improved in order to reproduce the whole scenario of the structural transition in Ni 2 MnGa. Finally, it is worth mentioning that, by increasing the strenght of the magnetoelastic coupling U * , it is possible to extend the present study to samples for which the MT takes place in a paramagnetic matrix [39] VII. ACKNOWLEDGEMENTS We acknowledge A. Planes, Ll.Mañosa and E. Obradó for fruitful discussions. This work has received financial support from CICyT (project number MAT98-0315).       Fig. 10. Besides, we have indicated the position of the large fluctuation peaks with squares and the metastability limits of the first-order transition lines with points inside triangles. The orientation of the triangles indicates the direction of the MC simulation runs performed in order to locate each metastability limit. Note that we have performed increasing T * runs (△) and increasing ( >) and decreasing (< ) K * runs. Dotted lines are guides to the eyes