First-Order Phase Transition in Perovskites Pr$_{0.67}$Sr$_{0.33}$MnO$_{3}$ - Magneto-Caloric Properties -- Effect of Multi-Spin Interaction

We show by extensive Monte Carlo simulations that we need a multi-spin interaction in addition to pairwise interactions in order to reproduce the temperature dependence of the experimental magnetization observed in the perovskite compound Pr$_{0.67}$Sr$_{0.33}$MnO$_{3}$. The multi-spin interaction is introduced in the Hamiltonian as follows: each spin interacts simultaneously with its four nearest-neighbors. It does not have the reversal invariance as in a pairwise interaction where reversing the directions of two spins leaves the interaction energy invariant. As a consequence, it competes with the pairwise interactions between magnetic ions. The multi-spin interaction allows the sample magnetization $M$ to increase, to decrease or to have a plateau with increasing $T$. In this paper we show that $M$ increases with increasing $T$ before making a vertical fall at the transition temperature $T_C$, in contrast to the usual decrease of $M$ with increasing $T$ in most of magnetic systems. This result is in an excellent agreement with the experimental data observed in Pr$_{0.67}$Sr$_{0.33}$MnO$_{3}$. Furthermore, we show by the energy histogram taken at $T_C$ that the transition is clearly of first order. We also calculate the magnetic entropy change $|\Delta S_m|$ and the Relative Cooling Power (RCP) by using the set of curves of $M$ obtained under an applied magnetic field $H$ varying from 0 to 5 Tesla across the transition temperature region. We obtain a good agreement with experiments on $|\Delta S_m|$ and the values of RCP. This perovskite compound has a good potential in refrigeration application due to its high RCP.


I. INTRODUCTION
There is always a challenge for theorists to elaborate models which yield quantitative agreements with experiments in materials science.On the one hand, this is because experimental samples are not always "clean" due to various causes such as the presence of impurities, dislocations, domains etc. in addition to the method of sample preparation.On the other hand, theoretical models are often too simple to describe experimental samples.Of course, if models contain some main ingredients, they can at best get qualitative agreements with experiments.In magnetism, theories often start with a Hamiltonian containing pairwise interactions with a spin model (Ising, XY, Heisenberg, Potts, ... spins).The interaction can be of short range, long range or can be of various kinds causing frustration [1] or topological spin structures (skyrmions) [2].The nature of the phase transition depends on these elements: it can be of second or first order or unknown critical properties.Most of transitions due to short-range pairwise interactions in two or three dimensions (2D or 3D) are known by the modern theory of phase transitions such as the renormalization group [3] and highly performant numerical simulations such as the histogram techniques and the Wang-Landau method of simulations [4,5].So far, in most cases, the order parameter such as the magnetization M decreases as the temperature T increases and vanishes at the transition temperature T C , except in some frustrated systems where there is no ordering at any T [1], or there is a co-existrence of order and disorder at non-zero T [6,7].
In this paper we are interested in a family of perovskite compounds.Perovskite compounds have recently used in many applications ranging from solar cells with high-efficiency photovoltaics [8,9], lasers [10], light-emitting diodes [11][12][13], ... to colossal magnetoresistance [14,15] and magnetic refrigeration devices [ 16].In what which concerns their magnetic properties, there is a book edited by E. Dagotto [15] where experimenttal and and theoretical progress have been presented.The reader is referred to this book for a complete review up to 2003.More recently, it has been shown experimentally [17][18][19][20][21][22][23] and theoretically [24,25] that perovskites of manganite family with various doping atoms possess very high Relative Cooling Powers (RCP) which can be used as clean sources for refrigeration [16].The so-called "magnetic refrigeration" is based on the magneto-caloric effects (MCE) discovered by Warburg in 1881 [26] and widely investigated since then (see the review given in [27]).
In the particular case of the perovskite compound Pr x Sr 1x MnO 3 , we have however seen that experimentally M can have a plateau from T = 0 up to the transition temperature T C as seen in Pr 0.9 Sr 0.1 MnO 3 [24 ], or M decreases with increasing T before making a sharp fall at T C in the strongly bond-disordered Pr 0.55 Sr 0.45 MnO 3 [25].Very recently, it has been experimentally observed that M increases with increasing T before making a sharp fall at T C in Pr 0.67 Sr 0.33 MnO 3 which is the subject studied in this paper.Note that in the theory of phase transitions and critical phenomena, the temperature dependence of the order parameter and the order of the phase phase transition depends only on a few parameters: the nature of the interaction (short range, long range, competing interactions,...), the spin model (Ising, XY, Hesenberg, Potts models) and the space dimension [1,3].In view of the unusual behavior of M mentioned above, we have introduced a multi-spin interaction in addition to the standard pairwise ones in order to explain these experimental cases [24,25].This multi-spin interaction consists in taking into account the simultaneous interaction of 5 spins at each lattice site.In the case Pr 0.9 Sr 0.1 MnO 3 we succeeded to reproduce the magnetization plateau between T = 0 and T C experimentally observed.We have also succeeded to reproduce the experimental magnetization M as a function of T in the case of Pr 0.55 Sr 0.45 MnO 3 .Note that in this latter case, the sample is strongly disordered due to the mixing of Mn 3+ of spin S = 2 with Mn 4+ of spin S = 3/2, in addition to the random dilute positions of Pr 3+ .In such a strongly disordered system, it is rare that the transition is very sharp as experimentally seen in spin glasses.We have shown that with the multi-spin interaction, we could reproduce this result [25].
In this paper, we use again the multi-spin interaction to study the case of the compound Pr 0.67 Sr 0.33 MnO 3 .This compound has been experimentally studied.It shows an unusual behavior of M different from the previous cases: M increases as T increases up to T C .before making a very sharp transition [23].With a fine tuning of parameters, we have reproduced the experimental M versus T as seen below.We have also used the energy histogram method to show evidence that this transition is of first order in agreement with the experimental sharp transition.We have calculated the magnetic entropy change and deduced the RCP, both are in agreement with experiments.Note that the present compound has a high RCP which is interesting for magnetic refrigeration applications.
In section II, we describe our model Hamiltonian and discuss about the multi-spin interaction.The numerical simulation method is also explained.In section III, we show our results and compare to the experimental data on the temperature dependence of M .In section IV, we calculate the magnetic entropy change |∆S m | under an applied magnetic field µ 0 H up to 5 Tesla.We deduce the RCP for µ 0 H=1, 2, 3, 4 and 5 Tesla.These results are in good agreement with experimental data.Concluding remarks are given in section V.

II. MODEL AND SIMULATION METHOD
A. Experimental results of Pr0.67Sr0.33MnO3 Let us describe the structure of Pr 0.67 Sr 0.33 MnO 3 .It is crystallized in an orthorhombic structure with Pnma space group comprising lattice parameters on three crystalline axes a = 5.4734 Å, b = 7.7284 Å and c = 5.4880 Å.This structure can be modelized as a body-centered tetragonal lattice with the longest axis is the b crystal axis.The Mn 3+ and Mn 4+ occupy the basal planes (a, c), the centered sites are occupied at 67% by Pr 3+ and 33% by Sr. Since Sr has no spin, the centered sites are magnetically dilute.The outermost spins of Mn 3+ occupy the orbital 3d 4 , so its spin is S = 4 × 1/2 = 2 by the Hund's rule.The outermost spins of Mn 4+ occupy 3d 3 , so its spin is S = 3 × 1/2 = 3/2.As for Pr 3+ its outermost electrrons occupy the orbital 4f2, so its spin is S = 1 by the Hund's rule.Note that the Mn 3+ and Mn 4+ ions occupy randomly the basal sites while Pr 3+ occupy random centered sites .We show in Fig. 1 an example of ion distribution in a bct cell.
Experiments have been performed on this compound [23].The data show that M in zero-field cooling slighly increases with T and makes a sharp fall at T C = 291 K.As said in the Introduction, this behavior is strikingly unusual.Note that the compounds with the percentage of Pr close to 0.67 show different temperature dependence of M : Pr 0.63 Sr 0.37 MnO 3 [19 ], Pr 0.60 Sr 0.40 MnO 3 [17 ], Pr 0.7 Sr 0.3 MnO 3 [17 ].Ref. [18] shows curves of M (T ) in various applied fields, so there is strictly speaking no phase transition for ferromagnets in an applied field.Refs.[17,22] show a tendency of M increasing with increasing T as in Ref. [23] but the transition is not sharp due to a small applied field.We think that the work in [23] is more recent with probably higher performant experimental techniques.We will focus on the results of this work to elaborate our theoretical model.
The Pr x Sr 1−x MnO3 compounds have complicated structures due to the random mixing of spins S = 2 of Mn 3+ , S = 3/2 of Mn 4+ and S = 1 of Pr 3+ in addition to the zero-spin Sr sites .These are strongly disordered systems of mixed spins, specially when x is between 0.4 and 0.8.To our knowledge, there are so far no theoretical works taking into account the above-mentioned factors, except our two previous [24,25].We have tried various pairwise interactions and various spin models but none gives sharp first-order-like fall of M at T C .In statistical physics, few pairwise interactions give rise to a first-order transition: we can mention what is known: frustrated systems in 3D, q-state Potts models in 2D with q > 4 or in 3D with q ≥ 3 in 3D.The perowskite compounds we study are not frustrated because, in spite of the mixed spins, the ferromagnetic interaction between Mn 3+ and Mn 4+ due to double exchange interaction via intercalated oxygen ions [28][29][30][31] is dominant, leading to a ferromagnetic ordering.So the sharp fall of M at T C in these compounds may be due to other types of interaction.We have proposed in previous works [24,25] a multi-spin interaction which reproduce not only the sharp fall of M at T C but also its temperature dependence from 0 to T C .As will be shown below, the multi-spin interaction creates spin fluctuations under control by tuning its strength with respect to the pairwise interactions.
There are 6 kinds of pairwise interactions in addition to the multi-spin interaction: J 1 : Interaction coupling of a Mn 3+ ion with a NN Mn 3+ ion, J 2 : Interaction coupling of a Mn 3+ ion with a NN Mn 4+ ion, J 3 : Interaction coupling of a Mn 4+ ion with a NN Mn 4+ ion, J 4 : Interaction coupling of a Pr ion with a Mn 3+ ion, J 5 : Interaction coupling of a Pr ion with a Mn 4+ ion J 6 : Interaction coupling between two Pr ions on the adjacent bct units.
The pairwise Hamiltonian is written as where S i is the spin at the lattice site i, <i,j> is made over spin pairs coupled through the exchange interaction J ij .H is a magnetic field applied along the z axis.Depending on the kinds of the ions at lattice sites i and j, we have J ij equal to one of the above six kinds of interaction.Using only these interactions and with spin models Ising and Heisenberg, we did not succeed to reproduce the curve M versus T without a multi-spin interaction, as shown in [24,25].At this stage, let us note that the pairwise interaction used in magnetism comes from model Hamiltonians of statistical physics.All of them, except the Heisenberg case, are introduced by "hand", but their validities are verified by many experiments: we can mention the Ising model [32] and various Potts models [33] among others.The Heisenberg case is exceptional: starting from the overlap of two neighboring orbitals, and using their spin-dependent wave functions in the Hartree-Fock approximation it was shown that the second-order perturbation gives rise to the coupling between the spins of neighbor atoms (the reader is referred to the detailed demonstration given on pp.53-55 of Ref. [34]).It is however intuitively, the interaction between a spin with all of its neighbors should be simultaneous.Strictly speaking it is not the sum of interacting pairs.However, mathematically it is not possible to demonstrate a multi-spin interaction from the first principles.It is nevertheless, starting from vertex model one can exactly show that Ising models with m-spin interactions (m > 2) can be derived from the square-lattice eight-vertex model [35] which can be mapped onto an Ising model with two-and four-spin interactions by Wu [36] and by Kadanoff and Wegner [37].Let us mention also the Ashkin-Teller model [38] which was also reformulated as an Ising model on the square lattice with two-and four-spin interactions by Fan [39].Note in passing that the Ising model with three-spin interactions on the triangular lattice was solved by Baxter and Wu [40,41].A recent work on the m-spin interaction in 1D has been studied by Turban [42].The reader is referred to this work for mathematical details.So, the formulation of various multi-spin interactions is not new, but their effects are not widely studied, although there is a renewed interest seen in some very recent works [42][43][44].On this aspect, our multi-spin interaction used to explain various unusual behaviors of M in real materials [24,25] is the first attempt to enrich the family of magnetic interactions in real materials.
The multi-spin interaction is written as where K is the interaction strength and the sum runs over all Mn sites i.The spins S i interacts simultaneously with its four nearest neighbors (NN) S i1 , S i2 , S i3 and S i4 on the xy plane (see Fig. 2).
FIG. 2: Spin of the Mn ion at the site i interacts simultaneously with four Mn spins at nearest sites i1, i2, i3 and i4 in the multi-spin interaction given by Eq. (2).See text for comments.
Note that the Ising pairwise interaction obeys the reversal symmetry, i. e. the properties of the system are invariant by the operation (S i → −S i ,S j → −S j ).Also, the properties of the system do not change under the Mattis transformation (S i → −S i ,J ij → −J ij ).This transformation just changes the feromagnetic ordering to antiferromagnetic ordering, but the transition temperature and the critical exponents remain the same (see Ref. [34], p. 143).In contrast to pairwise interactions, the multi-spin interaction given above does not have the "overall" reversal symmetry (reversing all spins) because the energy changes its sign.However, if an even number (2 or 4) of spins changes their sign, the energy remains invariant.This property is interesting because it creates a large degeneracy which gives rise to first-order transitions as seen in q-state Potts models mentioned above.We return to this point below.
Let us note that in our model the spin of amplitude S is a multi-state Ising spin, taking its values among −S, −S + 1, .., S − 1, S (2S + 1 values).This is contrast to the standard Ising model ±S.The multi-state Ising spin allows for a smooth variation of energy in the spin updating.

C. Method of simulation
We first use the Metropolis algorithm [45] to perform simulations.The sample sizes are 30 3 bct cells, namely 2 × 30 3 lattice sites.The total number of spins is N .We generate a random spin distribution on the lattice of Pr x Sr 1−x Mn 3+ x Mn 4+ 1−x with x = 0.67 We use periodic boundary conditions in all directions to reduce surface effects.The first 10 5 MC steps are used to equilibrate the system, and the averaging of physical quantities is taken during the following 10 5 MC steps.The calculated physical quantities are: the average internal energy E per spin, the specific heat C V per spin, the magnetization M and the susceptibility χ per spin.They are defined by where M ℓ (ℓ = 1, 2, 3) is the average of spins of kind ℓ: ℓ = 1 for Mn 3+ , ℓ = 2 for Mn 4+ , and ℓ = 3 for Pr 3+ , defined as where N ℓ is the number of spins of kind ℓ.Note that the Landé factor g in ( 5) should be understood as an "effective" g since it concerns three types de spin.
Let us emphasize that in the updating of spins, we use the over-relaxation method by Creutz [46,47] which is known to accelerate the convergency to equilibrium [48].The reader is referred to [25] for more details of the implementation.In short, at a given spin S i , we calculate its interaction energy using Eqs.( 1) and ( 2).We take a new value of S i among its (2S + 1) values, we calculate its new energy.This new state is accepted or not using the Metropolis algorithm.Now, instead of taking another spin to proceed, we stay with the spin S i and try to update it several times: this is called over-relaxation updating by Creutz [46,47] which optimizes the equilibrium state.The reproductivity of the results is checked by using independent runs with different random spin distributions.
Before showing our results, let us emphasize that for the purpose of the present paper, we do not need to use highperformance MC techniques such as the multi-histogram method [4] (see an example in [49]) or the Wang-Landau algorithm [5] (see an example in [50]).We shall show below however a simpler energy histogram technique which indicates evidence of the first-order character of the transition observed the the present compound.

III. RESULTS -COMPARISON WITH EXPERIMENT
The temperature dependence of M is the main experimental data which is used for our modeling.The energy is not experimentally accessible.The specific heat and the susceptibily are not experimentally available.
To fit with the experimental M , the best set of parameters we found are where J is the energy unit for the simulations.Its real value is obtained when we fit the simulation M (T ) curve with the experimental M (T ).C is the reduction coefficient applied for interactions along the z axis whose lattice constant is much longer than those in the xy plane.
Using the above interaction values we obtained the MC transition temperature T C (M C) ≃ 3.684J.To find the value of J, we use the equality T C = 291K = 3.684J since the transition temperature is proportional to the energy unit.We have J ≃ 79 Kelvin.This J is however an effective interaction since there are several kinds of interaction.Let us show the MC curve and the experimental magnetization M in Fig. 3a.By fitting the value of MC magnetization at T = 0 with experimental value of Ref. [23], one obtains gµ B ≃ 0.74 emu/g.Figure 3b shows the MC result for the susceptibility χ vs T (two colors are two independent runs in order to have more data near T C ).
Several remarks are in order: • The MC data coincide with the experimental magnetizatio for the whole temperature range from 0 to T C , • The magnetization increase slighly with increasing T up to T C , this unusual behavior was also experimentally observed in early works [17,22] but the transition was not sharp as in the recent work [23], • The fall of M at T C is vertical, suggesting a first-order transition.We will show evidence of this below.
• The peak of χ is very sharp.
Knowing J, we can calculate the interactions J 1 , J 2 , ..., J 6 in real unit.For example the dominant interaction J Now, if we view the system as an "effective" ferromagnet with an effective exchange interaction J ef f , we can use the following expression for the ground-state (GS) energy where Z = 6Mn+8 × 0.67Pr=11.36 is the effective coordination number at a Mn lattice spin, S ef f is the effective spin length given by S ef f = (0.67 × 2 + 0.33 × 1.5 + 0.67 × 1)(0.67 + 0.33 + 0.67) ≃ 1.503.
Note that this value of the effective exchange interaction, very close to J (=79 K), is in the order of magnitude of interactions in magnetic materials with T C at the room temperature [51,52].
Let us show the internal energy per spin E and the specific heat C V in Fig. 4. Two colors correspnd to two runs to have more data around T C .Some remarks are in order: • The GS energy is 14.1 in unit of k B J in the figure which is -96 meV [see (11)], • E has a very stiff slope at T C .We will show below by energy histograms that in fact E is discontinuous at T C , • The peak of C V is extremely high as a consequence of strong fluctuations due to the muli-spin interaction.We will give a discussion below.
We show now in Fig. 5 the z components of three sublattice magnetic ions ⟨S i ⟩ (i = 1, 2, 3).We see here that the Mn 3+ (S 1 ) and Mn 4+ (S 2 ) have the same sign, namely they order ferromagnetically, while Pr 3+ ions (S 3 ) are ordered antiferromagnetically with the Mn ions.These data are MC results, they were not available by experiments.Note that all curves fall sharply at T C .This results in the vertical fall of M shown in Fig. 3a.

A. Energy Histograms: Evidence of First-Order Transition
In the theory of phase transitions and critical phenomena, a transition of first order is characterized by a discontinuity of the first derivatives of the free energy (see chapter 7 "Phase Transitions" in Ref. [34]).In our problem, it is M and E which should be discontinuous at T C .Take E for example.At T C , theoretically there is a co-existence of the ordered phase of energy E 1 and the disordered phase of energy E 2 .The distance E 2 − E 1 is called "latent heat" in first-order transitions.Let us expand a little bit how in MC simulations a first-order character is observed: due to the updating dynamics, the co-existence of two phases can take place in two ways: i) at a time t, the two phases spatially coexist (half space for each phase for example), ii) at time t one of the phases occupy the whole sample and this lasts for a lapse of time, it changes to the other phase during a next lapse of time, and repeats this cycle.For the first scenario, if one averages the energy on the sample at t, we find the average value Ē = (E 1 + E 2 )/2 which is located at the middle of the vertical slope at T C , one can have a wrong feeling that the energy is not discontinuous.In fact, it is.For the second scenario, if we again average the energy of the system during a long time, as we do in MC simulations, the time-averaged energy is also Ē = (E 1 + E 2 )/2 which is located in the vertical slope of E at T C (we suppose the two lapses of time are equal).So, how to see the discontinuity of E at T C ?The answer is we have to record the energy at t for a very long time and make a histogram h(E(t)).The second scenario will be seen by a double-peak structure of h(E(t)), one peak is centered at E 1 and the other peak centered at E 2 .
We have realized the energy histogram at and near T C .We show in Fig. 6  Fig. 6a and Fig. 6b correspond to the second scenario in which the system goes between the ordered (left peak) and disordered (right peak) phases.As said above, if we take the energy average of Fig. 6a, we have < E >= −6.42686 which is located in the middle of the energy gap.For Fig. 6b, we have a slightly different < E >= −5.89695 located also in the energy gap.The points in the vertical slope of Fig. 4a are results of the average: in fact the system does not stay at the averaged points, but at E 1 and E 2 as shown by the double peaks in these figures.Below and above T C , we have only one peak, the averaged energies are near the positions of the peaks, taking into account the tiny peaks.
To conclude this section, we emphasize on two points: i) the multi-spin interaction allows for an excellent agreement on M (T ) between our model and experiments, ii) the energy histogram technique clearly shows evidence that the transition in the present perovskite compound is of first order.

IV. MAGNETOCALORIC EFFECT -MAGNETIC ENTROPY CHANGE
In order to calculate the magnetic entropy change, we apply a magnetic field H i ∈ (0, ..., H) on the sample at a given T .We calculate the magnetization M (T, H).This is shown in Fig. 7.These results are to be compared with the experimental data taken from [23] of sample P1200 shown in Fig. 8.We see that the agreement between experiment and simulation is good except at very low H where simulations give higher magnetizations at temperatures around T C : experimental curves have stronger slopes, while numerical curves attain horizontal slopes at lower fields in this temperature range.This is expected since simulation samples are exempt of defects, domains,... which are more or less present in experimental samples.These defects prevent an increase of M at low H in experiments, but they do not resist to H at stronger values.
Let us calculate the magnetic entropy change |∆S m | defined by the Maxwell's formula FIG.8: Effect of magnetic field (Tesla) on the magnetization experimentally obtained for T = 260 K, 262 K, ..., 320 K (from top to bottom line).This figure is made using data of Fig. 8 for sample P1200 of Ref. [23].
where δM is the variation of the magnetization at H when T varies from T to T + δT .This formula is discretized and used in experiments and simulations as The experimental magnetic entropy change |∆S m | using this formula is shown in Fig. 9a (sample P1200 of [23]): the peak temperature does not change significantly with increasing H but the peak height is higher with larger H.For comparison |∆S m | obtained by MC simulations is shown in Fig. 9b.We note that the peak heigth is slightly higher than the experimental one, for all H.We believe that this is a consequence of higher M at very low H mentioned above.
. We show in Table I the experimental Relative Cooling Power (RCP) taken from Ref. [23] and our MC results of RCP.Note that Ref. [23] gives only the RCP for 1 and 2 Tesla for the sample P1200.But using their curves shown in Fig. 9, we calculated the experimental RCP for 3, 4 and 5 Tesla shown in the Table .Our RCP values are sommewhat higher than the experimental values but they are within the same order of magnitude.We consider this as a good agreement.

V. CONCLUSION
In this paper, we proposed a model Hamiltonian which includes a multi-spin interaction in addition to the pairwise interactions to study the magnetic properties of Pr 0.67 Sr 0.33 Mn 3+ 0.67 Mn 4+ 0.33 .We compared our results with experimental measurements of the magnetization M versus temperature T which shows that M (T ) unusually increases with increasing T before making a vertical fall at the transition temperature T C .With our model Hamiltonian, we performed numerical simulations and obtained results of M (T ) in good agreement with experiments for the whole range of temperature from 0 to T C .In addition, we showed a clear evidence of the first-order character of the transition, using the energy histogram technique.
We also calculated the magnetic entropy change and the Relative Cooling Power (RCP) for applied fields from 1 to 5 Tesla.These results are in agrrement with experimental data.The high coefficients RCP make the present compound a potential candidate for magnetic refrigeration.
We note that with the multi-spin interaction, we previously succeeded to obtain results for Pr 0.9 Sr 0.1 Mn 3+ 0.9 Mn 4+ 0.1 [24] and Pr 0.55 Sr 0.45 Mn 3+ 0.55 Mn 4+ 0.45 [25] which, in both cases, are in agreement with experiments.Rarely, modeling using a microscopic interaction Hamiltonian leads to good agreement on macroscopic properties experimentally observed.
To conclude, we believe that the multi-spin interaction is necessary for studying materials with complex structures such as in the above-mentioned compounds whose behaviors are unusual.

FIG. 3 :
FIG. 3: (a) MC result (violet squares) and experimental magnetization (green void circles) versus temperature T in Kelvin, are shown for comparison, (b) Magnetic susceptibility versus temperature T .See text for comments.

FIG. 4 :
FIG. 4: (a) Internal energy per spin versus temperature T in Kelvin, (b) Specific heat versus T .Line is guide to the eye.