Magnetoelastic properties and behaviour of 4C pyrrhotite, Fe7S8, through the Besnus transition

Pyrrhotite, Fe7S8, is a commonly occurring carrier of magnetic remanence and has a low temperature transition, the Besnus transition, involving a change in spin state. Variations of the thermodynamic, magnetic and elastic properties through this transition at ∼33 K in a natural sample of 4C pyrrhotite have been tested against a group theoretical model for coupling between order parameters relating to Fe/vacancy ordering (irrep U1(1/2,0,1/4)) and magnetic ordering (irreps mΓ4+ and mΓ5+). Magnetoelastic coupling is weak but the pre-existing microstructure of ferroelastic and magnetic domains, that develop as a consequence of Fe/vacancy and ferrimagnetic ordering during slow cooling in nature (P63/mmc → C2′/c′), causes subtle changes in the low temperature transition (C2′/c′ → P1¯). The Besnus transition involves a rotation of magnetic moments out of the a–c plane of the monoclinic structure, but it appears that the transition temperature might vary locally according to whether it is taking place within the pre-existing domain walls or in the domains that they separate. Evidence of metamagnetic transitions suggests that the magnetic field–temperature phase diagram will display some interesting diversity. Low temperature magnetic transitions in minerals of importance to the palaeomagnetism community have been used to identify the presence of magnetite and haematite in rocks and the Besnus transition is diagnostic of the existence of pyrrhotite, Fe7S8.


Introduction
Minerals which hold a remanent magnetic signal for tens to hundreds of millions of years provide a material record against which models of the evolution of the Earth's magnetic eld and past tectonic plate motions have been tested. Hematite (Fe 2 O 3 )-ilmenite (FeTiO 3 ) and magnetite (Fe 3 O 4 )-ülvospinel (Fe 2 TiO 4 ) are the most important solid solutions in this context but have proved to be of intrinsic interest also for the phenomenological richness of the magnetic and structural instabilities which they display. In addition to the initial magnetisation which provides the palaeomagnetic record, hematite has a spin op transition (the Morin transition) at ∼260 K [1], magnetite has a charge order transition (the Verwey transition) at ∼122 K [2,3] which is also ferroelastic, hemo-ilmenite undergoes Fe/Ti ordering, and glassy magnetic behaviour occurs at low temperatures in Ti-rich members of the hematite-ilmenite solid solution [4]. In ülvospinel, the magnetic transition at ∼125 K is accompanied by cooperative Jahn-Teller distortions [5].
Although not quite as widely relevant from a geological perspective, the sulphide mineral pyrrhotite (Fe 1−x S) also has the potential to preserve a valuable palaeomagnetic signal [6][7][8], including for meteorites and the crust of Mars [9,10]. However, it has an even more remarkable diversity of magnetic properties over a narrow composition range (0 < x 0.125) due to the fact that small changes in iron content lead to a variety of superstructures associated with ordering of vacancies on the cation site [11][12][13][14][15][16]. There is a uniform Néel temperature of ∼590 K across the pyrrhotite solid solution, Fe 1−x S with 0 < x 0.125, [12,14] above which the structure is hexagonal with space group P6 3 /mmc (NiAs structure). When x < 0.125 the magnetic ordering is antiferromagnetic with moments con ned to the hexagonal basal plane. The magnetic moments are on layers of iron atoms and are ferromagnetically aligned within the layer. Neighbouring layers are coupled antiferromagnetically. When x = 0.125 the vacancies occupy only one of the antiferromagnetic sublattices, leading to a bulk ferrimagnetic moment. When 0.08 < x 0.125 the iron vacancies order and it is this vacancy ordering that gives rise to the range of closely related superstructures. Two features in common across all these phases and their multiple instabilities are coupling between different order parameters and the development of spontaneous strain due to spin-lattice and spin-orbital effects. The primary objective of the present study was to investigate magnetoelastic and order parameter coupling associated with the low temperature transition at ∼35 K in a natural sample of pyrrhotite with composition Fe 7 S 8 . This transition has drawn comparisons with the Morin transition in hematite and the Verwey transition in magnetite [17][18][19][20][21], with a suggestion of in uence due to the Jahn-Teller effect [22]. Rochette et al (2011) [18] suggested that the low temperature transition in crystals with composition close to Fe 7 S 8 be named after Marie-Jeanne Besnus who, in 1964 [17], reported an anomaly in the resistivity and magnetisation at 30 K. She commented also that it was coincident with a broad feature in the heat capacity reported ve years earlier by Grønvold et al [23]. The transition came to prominence in the literature after the associated drop in magnetic moment in natural crystals was reported by Fillion and Rochette [19] Dekkers et al [24] and Rochette et al [8]. It was also seen in Mössbauer spectra, which showed a splitting of one of the crystallographic sites into two equally populated subsites [22,25]. Since then, the transition has been investigated in great detail using a variety of techniques, including x-ray diffraction [26], neutron diffraction [27][28][29][30] DC magnetometry [20,[31][32][33][34][35][36][37][38], AC magnetometry [21,32,36,37], calorimetry [21,23,34], electrical resistivity [21,34], Mössbauer spectroscopy [22,25,39], small angle neutron scattering [21], and torque magnetometry [40,41]. Reported changes in temperaturedependence of the structural parameters, magnetic properties and electrical conductivity appear to be consistent with a discrete, thermodynamically continuous phase transition [21,30,32,34,36,38,40,41]. A discrete reduction in symmetry is also suggested at an atomic scale by the Mössbauer evidence of splitting of one of the cation sites [22,25].
The crystallographic space group of 4C pyrrhotite is C2/c at room temperature. Magnetic moments are aligned perpendicular to the b-axis within a plane which is a few degrees from (001) [28]. With falling temperature, the moments rotate progressively towards the direction of the crystallographic c-axis within the a-c plane [27,28,30,34]. Wolfers et al (2011) [29] proposed that the Besnus transition involves a symmetry change from F2/d (a convenient non-standard setting of C2/c) to P1 or P1. It has not yet been possible to detect any distortion from monoclinic lattice geometry [30], but the proposed change in space group is consistent with group theoretical considerations of how the magnetism couples with the vacancy ordering [42]. At a macroscopic scale there is a change to four-fold symmetry in the orientation of easy magnetic directions within the a-b plane [29,38]. Koulialias et al [41] also reported the emergence of this four-fold symmetry in torque measurements, but they had previously considered changes in properties of their sample to be an extrinsic effect due to the presence of two epitaxially intergrown superstructures [34,36].
Changes in lattice parameters accompanying the symmetry change P6 3 /mmc → C2/c at 598 ± 5 K shown by Powell et al (2004) [28] arise by strain coupling with an order parameter for the 4C vacancy ordering scheme and an order parameter for the ferrimagnetic ordering [43]. Variations in lattice parameters below T B are much smaller [30] but are sufcient to con rm that there is signi cant coupling of strain also with the second magnetic order parameter. It follows that there must be coupling between two magnetic order parameters and a vacancy/structural order parameter. Because the symmetry changes are hexagonal → monoclinic → triclinic, there will be some contribution to the overall pattern of strain relaxation from ferroelasticity and the resulting microstructures must contain some combination of interacting antiphase domains, magnetic domain walls and ferroelastic twin walls.
Here we present new heat capacity, magnetic, elastic and anelastic data collected as a function of temperature and magnetic eld to show that the Besnus transition in a natural pyrrhotite crystal with composition Fe 7 S 8 appears to occur in stages which depend on the orientation of the eld with respect to different crystallographic twins arising from the hexagonal → 4C monoclinic transition. Static magnetic measurements as a function of temperature in three directions have revealed how the ferromagnetic moment evolves both above and below the transition temperature, T B . The role of strain is considered in respect of the transition and microstructure dynamics through determinations of the temperature-and eld-dependent variations of elastic constants and acoustic loss, measured by Resonant Ultrasound Spectroscopy (RUS). Magnetoelastic and order parameter coupling behaviour associated with the hexagonal (P6 3 /mmc)-monoclinic (C2 ′ /c ′ ) transition in the same sample is described in a separate paper (Haines et al 2020a [43]).

Sample description
A natural 4C pyrrhotite sample with composition Fe 7.00(6) S 8 was obtained from the mineral collection of the South Australia Museum and has been described in detail by Haines et al (2019b) [43] It is inevitable that the parent crystal contained multiple magnetic domains from ferrimagnetic ordering and ferroelastic twins from the rhombohedral → monoclinic symmetry change associated with vacancy ordering. There are expected to be six possible directions for the ferrimagnetic moment in the monoclinic structure. If all six twin domains were present in equal proportions, the net moment would lie in the (001) plane without a unique direction. . This falls close to the direction of the magnetisation which would represent the average, predominantly, of one pair of twin domains. Crystals 1 and 3 need not have had identical domain con gurations but both will have had a net moment either within or very close to the (001) plane at room temperature.

Heat capacity
Heat capacity measurements were carried out with and without an applied magnetic eld in a quantum design PPMS. Crystal 1 (12.8 mg) was held on the puck with Apiezon N grease. Two sets of data were collected over the same temperature range for each setting of the external magnetic eld. The rst was of the puck plus grease, to give the background heat capacity. The second was of the puck, grease and crystal, from which the measured background was subtracted in order to retrieve the heat capacity of the crystal. In all cases the sample was cooled in zero eld. At low temperatures the heat capacity of the crystal was large compared with the background. At temperatures above ∼30 K the background heat capacity was larger than that of the crystal.

DC magnetism
Magnetisation was measured using a Quantum Design MPMS. Crystal 2 (15.4 mg) was held between two identical quartz cylinders inside a brass half tube with [001] * , [100] * or [120] * aligned parallel to the magnetic eld. Zero-eld cooled (ZFC) data were collected by cooling the crystal to 2 K in zero applied eld before measuring the moment in a eld of 10 mT while heating back up to 300 K. Field cooled (FC) data were collected by rst cooling the crystal from room temperature in a 5 T eld. The eld was switched off at 2 K and the magnetisation then measured in a 10 mT eld during heating. Temperature was stabilised at each setpoint during measurements with H//[001] * direction and swept continuously at 3 K min −1 during measurements with H//[100] * or [120] * directions. Hysteresis loops were subsequently measured in the three principal directions at a selection of temperatures in elds between +5 and -5 T. In between measurements of the hysteresis at different temperatures the sample was heated well above the Besnus transition temperature and cooled to the target temperature in zero eld.

AC magnetism
AC magnetic measurements were carried out on crystal 1 (12.8 mg) using the ACMS II option of a Quantum Design PPMS, with a DC eld of 2 mT and an AC eld of 0.3 mT at frequencies of 0.01, 0.1, 1 and 10 kHz. Measurements were made on heating after cooling in zero eld.

Resonant Ultrasound Spectroscopy
The RUS technique has been described in detail in the literature [44][45][46][47][48]. It involves excitation of acoustic resonances of a mm-sized sample and measurement of the frequency, f , and peak width at half height, ∆f , of peaks in the resulting spectrum. Resonant modes are dominated by shearing motions so that the results from individual peaks provide information predominantly relating to shear elastic constants. Values of the elastic constants of the sample scale with f 2 . Acoustic loss is expressed in terms of the inverse mechanical quality factor, Q −1 , which is usually taken to be ∆f /f . For the low temperature experiments described here, crystal 3 (17.4 mg) was placed across a pair of faces between the two transducers in a holder that has been described by McKnight et al (2007) [49]. The holder was lowered into an Oxford Instruments Teslatron cryostat which has a superconducting magnet capable of delivering a magnetic eld of up to 14 T [50,51]. Data were collected in automated sequences with varying temperature at constant eld or with varying eld at constant temperature. A settle time, generally of 10 min, was allowed for equilibration at each set point of temperature or 1 min at each set point of eld. The orientation of the crystal was such that the magnetic eld was applied parallel to [120] * . We note here that due to the distance between the sample and the thermometer there is sometimes an offset between the measured and true sample temperature. The offset has been measured and is up to ∼2 K in the region of 35 K where the Besnus transition is observed. We therefore rely on the heat capacity measurments to provide the most accurate estimate of the Besnus transition temperature.

Heat capacity
Results of the temperature-and eld-dependent heat capacity measurements through T B are given in gure 1(a). Features at ∼215 K and ∼280 K seen in the insert, which shows data over the full temperature interval in zero eld, are known to be artefacts. Apart from a slight shoulder at the high temperature side, the peak at 33 K in zero eld looks like the classic excess heat capacity associated with a magnetic phase transition. There is a split into three overlapping peaks separated by ∼1 K in data collected with a eld applied parallel to [100] * . The three peaks are best seen in the data collected in an applied eld of 1 T ( gure 1(a)). Each of these peaks shifts to higher temperatures with increasing eld strength and their amplitudes also change systematically.
In order to determine the excess heat capacity, ∆C p , and excess entropy, ∆S (= C p T .dT), associated with the structural/magnetic changes which occur below T B ≈ 33 K, it is necessary to generate a baseline for the high temperature reference state. Instead of following the conventional approach of tting a polynomial function to the data above the transition point and extrapolating to lower temperatures, we have tted a 5th order polynomial to the anomaly-free heat capacity between 2 K and 60 K of a sample of 5C pyrrhotite which had been produced by annealing a piece of the 4C starting crystal at 600 • C (see Haines et al 2020b [52]). The t parameters were used to construct the baseline shown in the upper insert of gure 1(b). The resulting variations of ∆S, after integration of the differences, are shown in the main eld of gure 1(b) (see appendix A1 for further details). Variations of ∆S with temperature extracted in this way yielded a total entropy change of ∼-9 J mol −1 K and a tail above the transition point. The main differences between results obtained under different strengths of eld applied parallel to [100] * occur in an interval of ∼3 K through the transition point itself (lower insert of gure 1(b)). Heat capacity data for the same sample with the eld applied perpendicular to the (001) plane show a single peak pushed to higher temperature and broadened by the applied eld (appendix A2).

DC magnetisation
Values of the magnetisation measured in a 10 mT eld along [001] * , decreases slightly at rst and then more steeply below ∼200 K ( gure 2(b)). Variations in the orientation of M with temperature, as reconstructed from the three sets of measurements, are shown in terms of spherical coordinates in gure 2(b), where θ is the angle from Z and ϕ is the angle between X and the projection of the moment onto the XY plane. Tilting out of the (001) plane towards [001] * progressed to θ ∼ 56 • at 35 K, in comparison with previously reported values of ∼61 • at 11 K (Powell et al 2004) [28] and ∼66 • at 35 K (Koulialias et al 2018c) [30] from neutron powder diffraction. Below T B the direction of M rotates out of this plane, as found previously by Koulialias et al (2018c) [30]. The rotation observed here was from ϕ ≈ 60 • at 35 K to ϕ ≈ 90 • at 2 K ( gure 2(b)). The value of the net magnetisation, |M|, rst decreased by ∼20 % between ∼35 and ∼30 K and then increased slightly down to 2 K ( gure 2(b)). The maximum in d|M|/dT occurs between ∼31 and ∼33 K ( gure 2(c)).
It is important to realise that in the temperature interval between the Besnus transition and room temperature the system is in its fully ordered magnetic ground state with respect to the 595 K transition. The crystallographic twins are set, the vacancies will not reorder and the magnetic domains will be determined (in the absence of an applied eld) by the sample geometry and microstructural details, including twinning. The applied measuring eld is signi cant when compared to the coercive eld (H c ∼ 5-10 mT [8,32,36]), and the moment measured is therefore an appreciable fraction of the saturation moment.
Cooling in a 5 T eld resulted in signi cant changes in magnetisation in all three measuring directions below T B ( gure 2(a)). At   Blue squares show the temperature evolution of θ, the angle of the magnetisation vector from the c * axis, and green triangles the evolution of ϕ, the angle within the basal plane from the X-axis (the [001] direction of the parent NiAs structure). (c) d|M|/dT −1 for net moment calculated from the ZFC data has a single peak at 31 K ± 1 K. Coercivities reported in the literature for temperatures above T B are only ∼5-10 mT [8,32,36]. It is therefore not surprising that any poling of magnetic domains might not persist in a measuring eld of 10 mT during heating to above ∼35 K in the FC sequence, though there are slight differences between FC and ZFC results up to 300 K.
Variations of the orientation of the net moment, M, from the ZFC measurements between 2 and 300 K are shown with respect to an orthogonal set of reference axes, X, Y and Z, in gure 3. In this setting [100] * of the hexagonal and monoclinic structures has been set parallel to X.
[120] * of the hexagonal structure has been set parallel to Y, and is also parallel to [010] of the monoclinic structure.
[001] * of the hexagonal structure is parallel to Z.
[001] * of the monoclinic structure is not shown but β ≈ 90.5 • at 300 K in space group setting F2/d (Powell et al 2004 [28]), and so is within ∼1 • of [001] * of the hexagonal structure. Miller indices refer only to the hexagonal structure. The initial orientation was within the XY plane at an angle of ∼56 • from X ([100] * ) and 34 • from Y ([120] * ). This is shown as a red circle in gure 3 and is close to the orientation expected for one pair of the six possible twins. If only one twin was present, the moment would rotate towards Z along a straight line with ϕ = 60 • as temperature reduces from room temperature to ∼35 K. Below T B , the orientation of M rotates out of this plane at approximately constant θ, ending up on the plane between Y and Z (ϕ = 90 • ) at 2 K. An additional anomaly is evident in the M-H data collected with H//[120] * . At 40, 50 and 60 K there is a distinct break in slope at ∼0.6-0.7 T ( gure 4(c)). A break in slope is evident at ∼0.5 T also in the measurements made at 100 K, but the effect is weaker. A similar feature has been seen before and was attributed to the presence of 4C and 5C structures [36]. The crystal used in the present study only contained the 4C phase, however, and the breaks in trend are more likely to have been due to a metamagnetic transition.
Data for H//[001] * collected at 5, 20, 40 and 60 K have the same pattern of variations as for the other orientations, ie open loops below T B and very slightly open loops above T B ( gure 4(d)). However, the moment induced at 300 K was substantially less than at lower temperatures and in comparison with the other two directions. The insert in gure 4(d) shows that saturation was not reached even at 5 T, consistent with [001] * being a magnetically hard direction at room temperature [53].

AC magnetism
The variations of real, χ ′ , and imaginary, χ ′′ , components of magnetic susceptibility shown in gure 5 from the present study are closely similar to the AC magnetic response of other single crystal [21,36] and powdered samples [32]. Variations of χ ′ and χ ′′ measured at 10 kHz through the temperature interval 20-55 K are shown for H//[120] * in gure 5(a) and for H//[001] * in gure 5(b). Data collected at 0.1, 1 and 1 kHz were qualitatively the same and, in particular, the temperature at which the peak in χ ′′ occurred did not change with frequency. The main features of the temperature dependence of χ ′′ are an increase as the transition was approached from above, a peak at ∼37 K, followed by a steep drop to more or less constant values below ∼35 K. These were accompanied by rounding of χ ′ as a precursor effect ahead of T B , followed by a steep drop to very low values at ∼35 K. There is an additional small peak at ∼33 K in the data for

RUS
RUS spectra collected in a heating sequence from 5 K to 295 K are shown as a stack in gure 6(a). There are small increases and decreases in frequency of different resonance peaks showing that the single crystal elastic constants varied in a smooth but somewhat irregular manner through the temperature interval ∼150-295 K. More signi cant anomalies occurred near 35 K and the similarity of the trends for all the resonance peaks ( gure 6(b)) indicates that the shear elastic constants followed closely similar patterns of softening/stiffening through the Besnus transition. Broadening of the peaks also indicates signi cant attenuation in the vicinity of T B .
Variations of f 2 and Q −1 from tting of representative resonance peaks collected in zero eld are shown in gure 7(a) for the full temperature range and in gure 7(b) for small steps through T B (see the supplemental material for results of analysis on multiple peaks, available online at http://stacks.iop.org/JPhysCM/32/405401/mmedia). The Besnus transition is marked by steep softening from either side, with an eventual increase in stiffness by ∼2-6 % below the transition. The transition is also accompanied by marked increases in attenuation, with the maximum value of Q −1 occurring at ∼32 K ( gure 7(b)). Above T B there is hysteresis in the values of f 2 of some resonances between heating and cooling ( gure 7(a)). These differences and the fact that the crystal ended up elastically softer by up to ∼6% are most likely due to a change in the proportions of different twin orientations following cycling through the Besnus transition. There is also an overall upward drift in the values of Q −1 with increasing temperature above T B . Repeating the RUS measurements in elds of 2 T ( gure 7(c)) and 9 T ( gure 7(d)), applied parallel to [120] * , caused the temperature of the minimum in f 2 to shift upwards in temperature. The width of the minimum also broadened ( gure 7(c)) and became two clearly resolved minima ( gure 7(d)), apparently corresponding to two discrete phase transitions at ∼34.5 and ∼37.5 K.
Variations of f 2 and Q −1 for representative resonance peaks from spectra collected while increasing and decreasing the applied eld parallel to [120] * at constant temperature are shown in gure 8 (see the supplemental material for results of analysis on multiple peaks). In the data for 60 K ( gure 8(a)), a steep minimum in f 2 and a narrow peak in Q −1 imply that a metamagnetic transition occurred at ∼0.9 T. With increasing eld, Q −1 initially decreased from its value at zero eld and f 2 values were higher after the transition than they were before. This pattern was similar at 40 K ( gure 8(b)), though the eld at which the minimum in f 2 occurred is at ∼1.2 T. The pattern is completely different in data from the eld sweep at 20 K, ie below T B . With increasing eld, there was a steep reduction in f 2 at ∼2 T, but there was then no recovery at higher elds. This was reversed when the eld was reduced, but with a shift of the anomaly to ∼1.3 T. The loci of these anomalies were not explored for the other orientations of eld but are consistent with the evidence from M-H data collected from crystal 2 that a metamagnetic transition occurred at ∼0.9 T when a eld was applied parallel to [120] * at 40 and 60 K ( gure 4(c)). On the other hand, the M-H loop measured at 20 K does not show evidence of metamagnetic transition. The anomaly in f 2 evident at 20 K might, Figure 7. Temperature dependence of f 2 and Q −1 from tting of representative peaks in RUS spectra from crystal 3. Blue = cooling, red = heating. (a) Zero eld data for a resonance peak with f ∼ 700 kHz at room temperature. The Besnus transition is marked by a steep minimum in f 2 and a peak in Q −1 at 34 K. f 2 ended up ∼6% lower on reheating above the transition than it had been during cooling. (b) Zero eld data for a resonance peak with f ∼ 1220 kHz at room temperature. The variation of Q −1 accompanying the minimum in f 2 at 34 K is resolved into two peaks, at 33 and 35 K. (c) Resonance peak with f ∼ 700 kHz at room temperature, from spectra collected in a 2 T eld applied parallel to [120] * . The decrease in f 2 is over a wider temperature interval than seen for zero eld and is shifted by up to 3 K. (d) Resonance peak with f ∼ 700 kHz at room temperature, from spectra collected in a 9 T eld applied parallel to [120] * . therefore, be related only to a change in the con guration of ferroelastic twin domains/ferrimagnetic domain walls.
RUS measurements were also undertaken on a separate sample characterised in the earlier work of Volk et al (2016) [38] and Volk et al (2018) [21]. The pattern of elastic softening was closely similar to that shown above and the minimum in f 2 occurred at 33 ± 1 K, consistent with the value of T B estimated from other measurements such as heat capacity and resistivity [21,38] (appendix A3). Small differences in the value of T B between samples is likely to re ect differences in stoichiometry. The composition of the natural pyrrhotite sample used by Volk et al (2016) [38] and Volk et al (2018) [21] was Fe 7.2 S 8 .

Order parameters and symmetry
Anomalies in heat capacity, magnetism and elasticity data through the temperature interval ∼30-35 K presented here for a natural sample of pyrrhotite with composition close to Fe 7 S 8 are consistent with a discrete phase transition. This ts with a more general symmetry analysis in which different combinations of order parameters were considered to account for the vacancy and magnetic ordering transitions that have been identi ed in the Fe 1−x S solid solution (Haines et al 2019 [42]). Space groups for the proposed sequence of structures with falling temperature in Fe 7 S 8 were P6 3 /mmc-C2 ′ /c ′ -P1, in agreement with the possibility of P1 or P1 proposed by Wolfers et al (2011) [29] for the low temperature structure.
The key order parameters and their non-zero components for this sequence were listed by Haines et al (2019a) [42]. The P6 3 /mmc → C2 ′ /c ′ transition is driven by a combination of vacancy ordering, for which the order parameter has the symmetry of U 1 (1/2,0,1/4), and magnetic ordering, for which the order parameter has symmetry mΓ + 5 . This gives mΓ + 2 , mΓ + 4 and mΓ + 6 as secondary order parameters but, as discussed in more detail by Haines et al (2019a) [42] and Haines et al (2019b) [43], moments allowed by mΓ + 2 and mΓ + 6 are not observed, implying that they are zero for reasons of thermodynamic stability. In combination, mΓ + 4 and mΓ + 5 give the observed magnetic structure in which there are ferromagnetic moments within Fe layers and antiferromagnetic moments between them, with freedom for the direction to rotate between [100] * and [001] * . The change in order parameter which drives Steep anomalies near 1 T at 60 K (and 40 K) are taken to indicate that a metamagnetic transition occurred. Applying a magnetic eld perpendicular to the (001) plane was problematic as the sample was prone to be pushed out from between the transducers. Some data were recorded and are presented in appendix A2. the C2 ′ /c ′ → P1 transition is mΓ + 5 (−a,0.577a) → mΓ + 5 (a,b), allowing rotation of the moments out of the a-c plane.

Strain coupling
For kinetic reasons, it is unlikely that changes in vacancy ordering can occur at these low temperatures so that observed variations in spontaneous strain will be due to coupling only with the magnetic order parameters. By itself the proposed order parameter sequence mΓ + 5 (0,0) → mΓ + 5 (0,a) → mΓ + 5 (a,b) would give P6 3 /mmc → Cm ′ c ′ m → P2 1 /m. From the full Landau expansion given in Haines et al (2019a) [42], the symmetry-breaking strain at the rst transition is (e 1 -e 2 ) and at the second it is e 6 . Both couple to the relevant order parameter components, q, with the form λeq 2 . The structure is already monoclinic due to shear strain e 4 (≈cos β * in the setting used by Haines et al 2020a [43]) arising from Fe/vacancy ordering, which means that the symmetry-breaking strain e 6 causes the Besnus transition to be monoclinic → triclinic rather than orthorhombic → monoclinic. High resolution lattice parameter data of Koulialias et al (2018c) [30] from powder neutron diffraction obtained using a natural 4C sample did not reveal any measurable distortion from monoclinic lattice geometry but revealed small changes in the a-and b-parameters that are clearly related to the transition. Because of the non-linear contribution of saturation effects to thermal expansion as T → 0 K, it is not possible to extract reliable values of linear strains but it is still possible to follow variations of cos β * , where β * is the reciprocal lattice angle of the F2/d cell, and (e 1 -e 2 ). The latter can be estimated without signi cant loss of precision as (a − b)/ ((ab) 1/2 ). The temperature dependences of these two shear strains are shown in gure 9 and are consistent with variations over a wider temperature interval from the lattice parameter data of Powell et al (2004) [28] for a synthetic sample. The shear strain represented by cos β * arises primarily from coupling with Fe/vacancy ordering (e 4 in the setting of Haines et al 2019 [43]). Variations below 300 K are presumably due to a weaker, additional contribution of coupling with the magnetic order parameter but there are no obvious changes associated speci cally with the Besnus transition. In comparison, values of (e 1 − e 2 ) are small at all temperatures and arise by weak coupling with the magnetic order parameter. They also show a distinct reduction below T B from ∼0.002 to ∼0.001 ( gure 9), con rming that the magnetic order parameter for the Besnus transition also couples with shear strain. An equivalent coupling of e 6 ∼0.001, would imply distortion of one of the monoclinic lattice angles from 90 • to ∼90.06 • and it is not surprising that such a small change has not yet been detected.
Weak coupling between shear strains and the magnetic order parameter is consistent with the very slight overall differences in shear elastic constants below and above T B , as represented by changes in values of f 2 amounting to only ∼2-6% ( gure 7). The fact that the net response is stiffening rather than softening is likely to be due to the time scale of relaxation of the magnetic order parameter. If the order parameter relaxes on the time scale of an induced strain, ie ∼10 −5 -10 −6 s in an RUS experiment, coupling terms of the form λeq 2 are expected to lead to softening. If the relaxation is slower than this, changes in the elastic constants are expected scale with q 2 due to higher order coupling terms with the form λe 2 q 2 . Stiffening occurs when the coupling coef cient, λ, is positive (e.g. references [54][55][56]).
If the Besnus transition occurs as a consequence of an electronic instability, such as a cooperative Jahn-Teller distortion for example, a separate driving order parameter, q E , would be needed to describe the change from monoclinic to triclinic symmetry. This would belong to a zone centre irrep and would have bilinear coupling with the symmetry breaking shear strain, e sb , as λe sb q E (Haines et al 2019 [42]). This, in turn, would give rise to non-linear softening of the related shear elastic constant as T B is approached from above and below, characteristic of pseudoproper ferroelastic behaviour. Such a pattern of softening is not observed and it follows that, if there is an electronic component in the transition, it is a secondary effect. The observed elastic softening in a temperature interval of ∼5-10 K on either side of T B is invariably accompanied by increased attenuation ( gures 6 and 7(a) and (b)), signifying that it has a dynamical origin rather than being due to relaxation accompanying classical bilinear or linear/quadratic strain coupling.

Splitting of the Besnus transition
The transition at T B is thermodynamically continuous [21] and has a physical origin that is likely to be essentially the same as for the spin-op transition in end-member FeS and in crystals with other compositions in the Fe 1−x S solid solution (Haines et al 2019 [42]). In general, for a purely magnetic transition due only to ordering of spins, the total entropy change is expected to be nR ln(2s + 1) per mole, where the number of spin states, s, is 2 for Fe 2+ , R is the gas constant and n = 7 for heat capacities expressed with respect to one mole of Fe 7 S 8 . This gives -93.6 J mol −1 K as the total magnetic entropy difference expected between fully disordered and fully ordered states. Above T B there is already a high degree of magnetic order, however, and it is no surprise that the observed entropy change is only a fraction of this gure. If, as proposed by Wolfers et al (2011) [29] and Oddou et al (1992) [22], the transition is characterised by the breaking of symmetry speci cally associated with the vacant cation site, the total entropy change might be nR ln2. Here n = 1 for the number of vacancy sites per mole and a degeneracy of 2 is attributed to ordering which relates simply to a doubling of the number of independent Fe cations in the structure. This gives ∆S = -5.76 J mol −1 K, which, if there is also some electronic or phonon contribution, could reasonably account for the observed entropy change of ∼9 J mol −1 K from integration of the excess heat capacity ( gure 1(b)).
[001] * is a magnetically hard direction of the 4C structure at room temperature and moments of the monoclinic structure only start to rotate substantially out of the basal plane below ∼250 K [28,30]. Correlation of different properties measured on different pieces of the same original crystal show that rotation of the ferrimagnetic moment out of the a-c plane starts at ∼35 K, as measured in a eld of 10 mT ( gure 2(b)). Onset of the heat capacity anomaly in zero eld is also at ∼35 K, but the maximum excess occurs at 32.9 K and there is a slight shoulder at ∼34.1 K ( gure 1(a)). 35 K is also the temperature at which the steep reductions in both χ ′ and χ ′′ from AC magnetic data reach constant, low values ( gure 5) and there is a second small peak at ∼32.5 K ( gure 5(a)). A minimum in f 2 occurs between ∼33 and ∼36 K, and there is a peak in Q −1 near 32.5 K ( gure 7(a)). The steepest variation in net magnetisation, given by the maximum value of dM dT −1 , probably occurs between ∼31 and ∼33 K ( gure 2(c)). In detail, it thus appears that the transition is spread into two parts, with the onset of symmetry breaking due to changes in magnetic order occurring at ∼35 K and the main changes in thermodynamic properties, corresponding to the steepest change in the degree of order, occurring at ∼32.9 K. This difference is likely to be due to heterogeneity in the local magnetic structure arising from the presence of pre-existing ferrimagnetic domain walls.
Within magnetic domains of the C2 ′ /c ′ structure, the rst magnetic order parameter, q C , is large and within the walls it reduces to zero, with the possibility that coupling with the magnetic order parameter of the P1 structure, q P , could give rise to small variations in the value of the transition temperature. Biquadratic coupling of the form λq 2 C q 2 P is expected for two transitions arising from different instabilities. For transition temperatures which are widely separated, the simplest solution for the Besnus transition has q C remaining effectively constant, giving rise to the temperature dependent term of a conventional Landau expansion as 1 2 a(T − T B + 2λ a q 2 C )q 2 P . In other words, the transition temperature for the second order transition is renormalised from T B to T B − 2λ a q 2 C . The peak in heat capacity would come from the domains, which form the bulk of the crystal and have q C ≈ 1, but a small positive value of λ would cause the onset of the magnetic transition to occur at a slightly higher temperature in the middle of the domain walls where q C = 0.
The transfer of excess entropy from the main peak at 32.9 K to additional peaks at ∼34 and ∼35 K with increasing eld up to 9 T ( gure 1(b)) is most likely to be due to different parts of the crystal having slightly different transition temperatures depending on the orientation of the eld with respect to the C2 ′ /c ′ structure. Volk et al (2018) [21] showed that T B , as de ned by the maximum of the temperature derivative of measured magnetisation, varies systematically between ∼32 and ∼34 K depending on the orientation of the eld within the (001) plane. Their data for χ ′ and χ ′′ also suggest a difference of ∼0.5 K between transition temperatures measured in directions parallel and perpendicular to (001). Three distinct transition points would be consistent with the presence of ferroelastic domains with three distinct crystallographic orientations, with respect to the orientation of the applied eld, in the crystal used for the present study. Splitting of the transition in this way is consistent with the differences in excess entropy between different eld strengths occurring only in the close vicinity of the transition point ( gure 1(b)). The integrated excess entropy is the same, irrespective of eld strength, once the transition has occurred in all parts of the crystal.

Ferroelastic and magnetic domains
Previous reports of domain con gurations in natural pyrrhotite grains using the Bitter technique and the Magneto-optical Kerr effect have shown a characteristic pattern of lamellar 180 • domains on a scale of ∼2-30 µm [57][58][59]. Although not speci ed, the polarity of these domains is presumably parallel and antiparallel to the orientation of the individual moments determined by neutron diffraction, with domain walls close to (001) at room temperature. The observed values of θ and ϕ at room temperature for the orientation of the net moment are consistent with the crystal used in the present study containing a preponderance of twin domains which have ϕ = 60 • , θ ≈ 84 • and ϕ = 60 • , θ ≈ 96 • (including the expected out of plane angle of 6 • from Powell et al 2004 [28]). The observed rotation from θ ≈ 90 • to ∼58 • at T B , with ϕ remaining more or less constant ( gures 2(b) and 3), also ts with the individual moments rotating within the a-c plane as temperature is reduced.
An abrupt change of trend of the net moment from decreasing θ at approximately constant ϕ to increasing ϕ at more or less constant θ, such that ϕ ≈ 60 • at 35 K and ∼90 • at 2 K, matches the analysis of Wolfers et al (2011) [29] in relation to how the easy magnetisation direction would switch below a monoclinic → triclinic transition point. In particular, six-fold symmetry for the easy direction of a twinned crystal above T B would become fourfold due to the twins being further split into two equal proportions. Volk et al (2016) [38] also found a rotation of the easy magnetisation direction with respect to the a-b plane below T B and the same change from six-fold to four-fold symmetry.
The reduction in net magnetisation below T B under ZFC measuring conditions resulted from cancellation of components of moments which have opposing orientations within new sets of magnetic twins. Poling of these by cooling in a 5 T eld parallel to [001] * led to a net magnetisation that was larger by a factor of ∼3 at 2 K than the magnetisation acquired during cooling in zero eld ( gure 4(a)). On the basis of evidence in the hysteresis loops, poling in a 5 T eld applied parallel to [100] * and [120] * should also have caused an increase in remanent moment, M r . Instead, a signi cant reduction in the net moment was observed when measured in a eld of 10 mT after switching the 5 T eld off ( gure 1(a)). Failure of the crystal to remain poled in these two directions implies that magnetic twin domains renucleated when the high eld was removed. On heating back up to temperatures above T B , the distribution of moments between the original domain orientations was restored, re ecting the fact that the underlying con guration of twins and antiphase domains due to cation/vacancy ordering has a strong in uence on the con guration of magnetic domains.
Ferroelastic twinning due to the change in point group symmetry 6/mmm → 2/m has been analysed in detail by Haines et al (2019b) [43]. Most of the twin walls arising from vacancy ordering will be aligned nearly parallel to [001] * , though there can also be one set aligned parallel to (001). A few individual twin walls were observed in the sample used in the present study by electron backscatter diffraction. The increase in cdimension by a factor of 4 must result also in the development of antiphase domains with boundaries that can be nonconservative, in the sense that segments which are oriented parallel to (001) can have a local change in Fe:S ratio [60]. Reorientation of twin walls and non-conservative antiphase boundaries (APB's) from the con gurations they acquired during slow cooling in nature requires Fe/vacancy diffusion. On the basis of the kinetic data of Herbert et al (2015) [35] discussed above, this is unlikely to occur to any signi cant extent during cycling below room temperature. Differences between ZFC and FC magnetisation in gure 2(a) show that an external eld can induce rotation of the magnetic domains below T B but that they return to their original con guration in the stability eld of the monoclinic structure. Changes in the con guration of ferroelastic domains would cause changes in the elastic constants, as signi ed by changes in resonance frequencies of the RUS sample but, while there are changes between values measured before and after cooling through the Besnus transition ( gures 7(a) and (c)), these are very small.
The orientation of new ferroelastic twin walls below the monoclinic (2/m) → triclinic (1) transition will be the same as set out for the same change in point group symmetry in the mineral albite (NaAlSi 3 O 8 ) by Salje (1993) [61] (and see, also Carpenter and Salje 1998 [56]). One set of twin walls will be parallel to (010) of the monoclinic structure and the second set perpendicular to this, with an orientation between the x-and z-axes that depends on the relative values of the triclinic shear strains (e 5 and e 6 for the setting in which e 4 is the strain associated with β = 90 • in the monoclinic structure). Evidence from conventional diffraction measurements is that any distortion from monoclinic lattice geometry below 35 K is small, however, and, hence, that strain contrast across the twin walls would be very small. It is unlikely that there will be any associated changes in Fe/vacancy distribution at these low temperatures, so that the twin walls should be free to change orientation in the same manner as the new magnetic domains when an external magnetic eld is applied and they should not have the same memory effect as seen for ferroelastic domains arising from the hexagonal → monoclinic transition. It is not clear why resistance to recon guration of magnetic domains should be greatest when poling is induced parallel to [001] * at temperatures below T B but the local strain heterogeneity of ferroelastic twin walls of the monoclinic structure and local chemical heterogeneity of non-conservative APB's could act as effective sites for re-nucleation of magnetic domains to cause depoling in the other two directions.

Metamagnetic transitions
Evidence for a metamagnetic transition is seen in the magnetic hysteresis loops measured with H//[120] * ( gure 4(c)), and there are changes in elastic properties at 60 and 40 K with a pattern of softening/stiffening and acoustic loss ( gures 8(a) and (b)) which is closely similar to the pattern associated with the Besnus transition observed with varying temperature in zero eld. The new magnetic structure at high elds has not been characterised but is most likely to involve some different combination of non-zero components of the mΓ irreps. The M-H loop measured with H//[120] * at 20 K does not show evidence of a metamagnetic transition, but there are hysteretic changes in elastic properties at ∼1 T ( gure 8(c)) which are quite different from those associated with the Besnus transition. It remains to be seen whether the elastic anomaly is due to a rst order metamagnetic transition or to a change in the con guration of ferroelastic twin domains, but it appears that the H-T phase diagram would have an interesting topology.

Dynamics: AC magnetic and acoustic loss
The principal features of the χ ′′ variations reported here ( gure 5) are closely similar to those reported previously [21,32,36]. In particular, the temperature at which the peak occurs immediately above T B is independent of measuring frequency, indicating that it is not due to a thermally activated loss process [36]. The steep reduction to more or less constant values below T B has been attributed to a reduction of the mobility of 180 • magnetic domain walls by additional pinning effects in the stability eld of the low temperature structure [32]. A small peak observed a few degrees below the steep reduction coincides with the peak in heat capacity measured in zero eld (Volk et al 2018 [21] and compare gures 1(a) and 5(a) above). This association must be due to some aspect of changing dynamics of local moments at the temperature(s) where the temperature dependence of the new magnetisation is at a maximum. The weak anomaly is present when the measuring eld is applied in a direction within the plane (001) but not when the direction is perpendicular to (001), presumably depending on the orientation of new domain walls in the low temperature structure.
Variations of Q −1 shown in gure 8 have a different overall form in comparison with those of χ ′′ . Acoustic loss is relatively high in the stability eld of the monoclinic structure but reduces with falling temperature ( gure 8). Because ferroelastic twin walls generated at the hexagonal → monoclinic transition are tied to Fe/vacancy ordering, it is unlikely that the loss mechanism is due to their displacement under the low stress conditions and short timescales of a natural acoustic resonance. The loss mechanism is presumably due to motion of unspeci ed defects which couple locally with strain, including local variations in spin state and magnetoelastic coupling.
The onset of softening ahead of T B occurs at ∼45 K and has a pattern typical of uctuations of an order parameter with coupling to strain ahead of a co-elastic or improper ferroelastic phase transitions, such as occurs in quartz and perovskites for example [62,63]. The transition itself is then marked by an interval of increased acoustic loss (high values of Q −1 ) which coincides with a small amount of elastic softening followed by stiffening ( gures 6 and 7). A minimum in f 2 occurs between ∼33 and ∼36 K, but the variations of Q −1 become hard to follow in this interval because of the broadness of the primary resonance peaks. The most reliable measurements of peak width appear to indicate that there is a maximum in Q −1 between ∼32 and ∼32.5 K, which coincides with the peak in heat capacity ( gure 1(a)) and the small, rounded peak in χ ′′ ( gure 5(a)). In principle, the mobility of the second generation of ferroelastic twin walls discussed above would be expected to be restricted at such low temperatures. However, the thickness, w, of ferroelastic twin walls varies as w ∝ (T c -T) −1 below a second order transition [64][65][66] and thick twins are less likely to be pinned by defects that thin ones. This means that there can be a narrow temperature interval immediately below the transition point where more or less well developed but thick twin walls develop and would be mobile under the in uence of an external stress eld. It is tentatively proposed, therefore, that the anomalies in Q −1 and χ ′′ re ect dynamical magnetoelastic coupling related to uctuations of the magnetic order parameter above T B and to new twin walls immediately below T B . The mobility of the ferroelastic twin walls then diminishes steeply due to pinning when they become thinner as T reduces away from T B .
Splitting of the Besnus transition into two parts in gure 7(d), with minima in f 2 at ∼34.5 and ∼37.5 K, and peaks in Q −1 at slightly lower temperatures, might be taken to indicate that there are two discrete phase transitions in a 9 T eld with H//[120] * , but it is more likely that they represent transitions in two twin domains with different orientations with respect to the applied eld.

Conclusions
(a) Variations of the thermodynamic, magnetic and elastic properties through the Besnus transition at ∼33 K in a natural sample of 4C pyrrhotite with composition close to Fe 7 S 8 have been tested against the group theoretical model of Haines et al (2019a) [42] for coupling between order parameters relating to Fe/vacancy ordering (irrep U 1 (1/2,0,1/4)) and magnetic ordering (irreps mΓ + 4 and mΓ + 5 ). (b) The physical origin of the transition can be understood as being related to the spin-op transition at ∼450 K in FeS, modi ed by the presence of vacancies and the in uence of the Fe/vacancy ordering scheme. No evidence has been found for an additional electronic instability. (c) It is unlikely, for kinetic reasons, that changes in Fe/vacancy ordering occur through the transition point but details of the transition appear to be controlled by the pre-existing microstructure due to the P6 3 /mmc → C2 ′ /c ′ transition which occurs during slow cooling in nature. In particular, an unusual splitting of the transition could be understood in terms of differences in magnetic structure between the pre-existing magnetic domain walls and the domains that they separate. (d) The proposed change in symmetry, C2 ′ /c ′ → P1, will necessarily result in the development of a new set of ferroelastic twin walls but coupling of the magnetic order parameter with the symmetry-breaking strain is weak and the only evidence for these is a peak in the acoustic loss at a temperature which correlates with a small peak in the imaginary component of the AC susceptibility. (e) Changes in the elastic properties are consistent with weak magnetoelastic coupling and relaxation times of the magnetic order parameter, in response to an applied stress, that are longer than ∼10 −5 -10 −6 s. (f) Evidence for metamagnetic transitions under the in uence of a eld applied to [120] * suggests that H-T phase diagrams should display an interesting diversity of magnetic structures. The stability of these will depend on coupling with the underlying Fe/vacancy ordering, for which a common strain coupling mechanism is likely to be important.

Acknowledgments
This work was funded by a grant from the Leverhulme Foundation (RPG-2016-298), which is gratefully acknowledged. RUS facilities have been established and maintained in Cambridge through grants from the Natural Environment Research Council and the Engineering and Physical Sciences Research Council of Great Britain to MAC (NE/B505738/1, NE/F17081/1, EP/I036079/1). Heat capacity, DC and AC magnetisation measurements were carried out using the Advanced Materials Characterisation Suite, funded by EPSRC Strategic Equipment Grant EP1M00052411.

A1. Heat capacity entropy calculation
To obtain a reliable estimate of the excess entropy associated with the Besnus transition it is necessary to de ne a baseline heat capacity. This is not always easy and often the method of de ning the baseline can have serious implications for the extracted entropy. Following the standard procedure of tting a polynomial to the data above and below temperatures semiarbitrarily deemed to be 'far away' from the transition gave approximate results. However, the order of the polynomial and the exact cut off temperatures used had a reasonably signi cant effect on the extracted excess entropy. Fortunately, in this case, an independent baseline could be de ned from the heat capacity of a very similar pyrrhotite sample in which the Besnus transition does not occur. The measurement used was of a piece of 5C pyrrhotite that had been cut from the same crystal as the 4C sample studied here. The curve from the 5C sample was Measurements of heat capacity and RUS were also made with an external magnetic eld applied in the direction normal to the (001) plane. When the eld was applied perpendicular to the (001) plane, the peak in the heat capacity associated with the Besnus transition was not split into two or three peaks, as seen in the data presented in gure 1 with the eld applied within the (001) plane, but is pushed to higher temperatures and broadened considerably (see gure A1). The high eld data were scaled to match the zero-eld data over the range 60-180 K. There are no signi cant features in this range. The scaling is necessary due to a slight discrepancy in the addenda measurements over the two separate runs. RUS data with a eld of 2 T perpendicular to the (001) plane were also measured. The resulting variations of f 2 and Q −1 are shown in gure A2.

A3. RUS Data for the Volk et al Sample
Variations of f 2 and Q −1 from RUS spectra of a crystal from the sample described by Volk et al 2016 and 2018. [21,38] are shown in gure A3. There is no measureable hysteresis between spectra measured on cooling (blue) and heating (red) ( gure A3(a)). The average temperature of the minimum in Figure A2. Results of peak tting of RUS spectra collected from crystal 1 with a magnetic eld of 2 T applied in a direction perpendicular to (001). f 2 for three peaks around the Besnus transition is 33 ± 1 K ( gure A3(b)).