Magnetoelastic coupling associated with vacancy ordering and ferrimagnetism in natural pyrrhotite, Fe7S8

Magnetoelastic coupling associated with the hexagonal—monoclinic transition in a natural sample of the mineral pyrrhotite, Fe7S8, has been analysed in terms of separate coupling of spontaneous strains with two discrete order parameters, qv for Fe/vacancy ordering and qm for magnetic ordering. Coupling of the two order parameters separately with strain gives rise to two terms for coupling between them, λqm2qv2 and λqm2qv8, and a pattern of evolution in which qv varies continuously and qm discontinuously through a single transition point. The transition is ferrimagnetic and ferroelastic but the relatively slow relaxation rate for Fe/vacancy ordering, in comparison with magnetic ordering, results in elastic and anelastic properties which are quite different from those observed in other ferroic or multiferroic materials with two instabilities. Instead of classical elastic softening, there is stiffening of the elastic constants which scales with qm2 and qv2. Instead of the normal pattern of acoustic loss associated with the mobility and subsequent freezing for ferroelastic twin walls, the loss is consistently low throughout the temperature range 300 K–875 K.


Introduction
A notable feature of the mineral pyrrhotite, Fe 1−x S, is that it displays wide variations in magnetic and thermodynamic properties over a narrow range of compositions corresponding to 0 > x > 0.125 [1][2][3][4][5][6]. These have resulted in it having signicance in contexts as disparate as palaeomagnetism on Earth [7] and Mars [8,9], mineral exploration [10] and phase change magnetic memories [11]. Much of the diversity observed in natural and synthetic samples arises from the fact that small changes in iron content lead to a variety of superstructures Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. associated with ordering of vacancies on the cation site. The Néel temperature of ∼590 K [2] is nearly constant across the solid solution, but the symmetry of the magnetic ordering schemes which develop is controlled by the particular pattern of commensurate or incommensurate vacancy ordering at each composition. The expectation is that crystals will be ferrimagnetic if the vacancies are con ned to alternate iron layers of the NiAs-type parent structure and antiferromagnetic if they occur equally in all layers. As a consequence, it is generally held that 'monoclinic' pyrrhotite will have a magnetic moment whereas 'hexagonal' pyrrhotite will not (but see Horng and Roberts 2018 [12]). Changes in the hard and soft directions for magnetic ordering as functions of composition and temperature lead to a further set of magnetic transitions below the primary Néel point, as summarised, for example, by Schwarz and Vaughan (1972) [2]. Table 1. Order parameters with symmetries of different irreducible representations for the changes in space group P6 3 /mmc → C2 ′ /c ′ → P1, from Haines et al [19]. The basis vector is (2,−2,0), (2,2,0), (−1,1,2) and the origin is (0,1/2,0). a and b signify non-zero components of the multi-component order parameters. P1, P2, P4, C1 and C5 indicate the relevant order parameter directions, as given in the group theory program ISOTROPY [22]. The present study was initiated as part of a wider effort to characterise the strength and mechanisms of magnetoelastic phenomena in ferroic, multiferroic and superconducting phases through the in uence of strain coupling on elastic and anelastic properties [13][14][15][16][17]. Pyrrhotites provide an interesting example of magnetoelastic behaviour in which magnetic ordering is coupled with vacancy ordering either directly or indirectly via common strains. The magnetic transitions are accompanied by small but measurable changes in lattice parameters, but the dynamics of strain coupling and relaxational processes are constrained by the timescale required for changes in Fe/vacancy order to occur, in comparison with the much shorter time scale of changes in local spin con gurations. Experimental data of Herbert et al [18] for a sample with the 11C superstructure suggest that measurable changes in Fe/vacancy ordering can occur on a timescale of 10's to 1000's of seconds at ∼450-500 K and, by extrapolation, on a timescale of hours to 10's of hours at ∼350 K or years at ∼300 K. These results suggest that equilibrium might be maintained on a laboratory timescale between 600 K and ∼400 K, say, and on a geological time scale perhaps down to room temperature.
Symmetry relationships provide a framework of theory by which paramagnetic-antiferromagnetic/ferromagnetic transitions, Fe/vacancy ordering and Morin-type transitions across the whole range of pyrrhotite composition can be rationalised in terms of a well-de ned set of macroscopic order parameters [19]. Here the focus is on coupling of ferrimagnetism and vacancy order with strain in the 4C monoclinic structure of Fe 7 S 8 below ∼590 K. Because the symmetry changes are hexagonal → monoclinic → triclinic, this includes consideration of the development and properties of ferroelastic twin walls. Relevant aspects of symmetry changes associated with the hexagonal-monoclinic-triclinic sequence of phase transitions in Fe 7 S 8 are summarised in section 2, below. A formal analysis of coupling between strain and the relevant driving order parameters, using lattice parameter data from the neutron powder diffraction study of Powell et al [21], is presented in section 3. Section 4 contains a description of the natural single crystal used for measurements of elastic and anelastic properties by resonant ultrasound spectroscopy (RUS). Implications of the strain and elastic properties are considered for the hexagonal (P6 3 /mmc)-monoclinic (C2 ′ /c ′ ) transition in section 5. Magnetoelastic coupling associated with the Besnus transition at ∼35 K in the same natural sample is described in an accompanying paper (see reference [20]).

Order parameters and symmetry
For Fe 7 S 8 , the space groups of the sequence of structures with falling temperatures which emerges from the analysis of Haines et al [19] is P6 3 /mmc-C2 ′ /c ′ -P1, in agreement with the possibility of P1 or P1 proposed by Wolfers et al [21] for the low temperature structure. The key order parameters and their non-zero components are listed in table 1. The P6 3 /mmc → C2 ′ /c ′ (4C structure) transition is driven by a combination of vacancy ordering, for which the order parameter has the symmetry of irreducible representation (irrep) U1(1/2,0,1/4), and magnetic ordering, for which the order parameter belongs to magnetic irrep mΓ + 5 . This combination gives secondary order parameters belonging to irreps mΓ + 2 , mΓ + 4 and mΓ + 6 . The direction P2 of mΓ + 5 has non-zero components (−a,0.577a) for the setting given in table 1. By itself, the equivalent mΓ + 5 (0,a) would give a magnetic structure with space group Cm ′ c ′ m, which would have moments aligned parallel to [100] * H in ferromagnetic layers and antiferromagnetic alignments between layers. Subscript H is used here and below to signify indices given with respect to the parent hexagonal structure. In combination, mΓ + 4 and mΓ + 5 allow individual moments (and, hence, the net ferrimagnetic moment) to rotate within the a-c plane of the monoclinic structure, as reported by Powell et al (2004) [23]. It follows that the change in order parameter which would drive a transition from the C2 ′ /c ′ structure to P1 is mΓ + 5 (−a,0.577a) → mΓ + 5 (a,b). This allows rotation of the moments out of the a-c plane [19].
On its own, mΓ + 2 would give a ferromagnetic structure with moments aligned parallel to [001] * H while mΓ + 6 with direction P1 (a,0) would give a ferromagnetic structure with moments aligned parallel to [100] * H . Neither of these ferromagnetic ordering schemes has been observed in the monoclinic 4C (C2 ′ /c ′ ) structure, however, implying that the components of mΓ + 2 and mΓ + 6 are zero for reasons of thermodynamic stability.

Strain analysis
It is not possible to measure the values of order parameters directly and by far the most straightforward way of extracting their interdependence and variations with temperature are through measured variations of lattice parameters. If the order parameters are strictly zero in the high symmetry structure, it is necessary to choose reference states for calculation of the symmetry-adapted strains which also must be zero above the transition point. The reference system used here for analysis of spontaneous strains which couple with order parameters for vacancy ordering and magnetic ordering is shown in gure 1. In this setting, the basis vectors are (−2,−4,0), (2,0,0), (1,2,2) with origin (1/2,1/2,1/2) and the non-zero components of the U1(1/2,0,1/4) order parameter are (−0.707a,−0.707a,0,0,0,0.707a,−0.707a). For the purposes of producing a Landau expansion to describe the P6 3 /mmc-C2 ′ /c ′ transition, the U1(1/2,0,1/4) order parameter can be represented as q v and the non-zero component of the mΓ + 5 order parameter as q m . Coupling between the two order parameters and strain and between themselves then has the form (after Haines et al [19]) T cm and T cv are critical temperatures for the magnetic and vacancy ordering transitions, respectively, a v , a m , b v , etc, are standard Landau coef cients, e i (i = 1-6) are components of the spontaneous strain tensor, λ's are coupling coef cients and C o ik (i, k = 1-6) are elastic constants of the hexagonal parent structure. The lowest order coupling between the two order parameters, either directly or indirectly via these strains, is biquadratic, λq 2 v q 2 m . Coupling between the monoclinic shear strain e 4 and q v has the unusual form, λe 4 q 4 v , whereas it couples with q m as λe 2 4 q 2 m . These two terms would lead to indirect coupling between q v and q m with the form λq 2 m2 q 8 v . The equilibrium condition, ∂e/∂G = 0, gives strain/order parameter relationships as Individual strains are determined from lattice parameters according to (eg see Carpenter and Salje [24])  [25], shown as lled triangles, in comparison with those for the pseudohexagonal cell of Fe 7 S 8 from Powell et al [23], shown as red lines. The FeS data have been scaled by factors given in the legends so that they overlap with the Fe 7 S 8 data at temperatures above 600 K. Dashed lines represent reference parameter, a o , for the FeS structure, as obtained by tting the expression a o = a 1 + a 2 Θ s coth(Θ s /T), and similarly for c o , V o . The saturation temperature, Θ s , was xed at 150 K. In this sample, a second magnetic transition occurred at ∼390 K. (d) Spontaneous strains for FeS, calculated from the data in (a)-(c), compared with the intensity of a superlattice re ection, I (0001) , from neutron powder diffraction of Andresen et al [31]. The second transition occurred at ∼440 K in the sample of Andresen et al [31]. The right axis is scaled to demonstrate I (1000) (∝ q 2 m ) ∝ e 3 . (e) Comparison of I (1000) for the 4C structure of Fe 7 S 8 from Powell et al [23] with data for FeS from Andresen et al [31], showing that magnetic ordering has essentially the same temperature dependence in each phase. (f) Comparison of V and c parameters for Fe 7 S 8 , scaled so that they overlap at T 600 K. Open and lled symbols represent data from two different instruments in the original work of Powell et al [23]. Data for c show the same form of steep change at ∼600 K as seen in the data for V, but the trend is then immediately reversed.
Here a, b, c, β * and V refer to the F2/d cell of the monoclinic structure, rather than the C2/c cell. a o , c o and V o are reference parameters of the parent hexagonal structure, extrapolated into the stability eld of the monoclinic structure. Separation into strain contributions due to coupling with q v and q m is made possible by comparison of the lattice parameter data of Powell et al [23] for Fe 7 S 8 with data of Tenailleau et al [25] for a natural sample of FeS, troilite, which has antiferromagnetic ordering without vacancy ordering. The change in symmetry expected on the basis of group theory for FeS is P6 3 /mmc → Pnma (moments parallel and antiparallel to [010] H in alternate layers) or Pnm ′ a ′ (moments parallel and antiparallel to [100] * H in alternate layers) [19]. Tenailleau et al [25] were able to t their neutron powder diffraction patterns with hexagonal lattice parameters down to at least ∼400 K, . (e 1 − e 2 ) values remain low at all temperatures, while cos β * shows a steep increase at ∼600 K, consistent with rst order character for the hexagonal-monoclinic transition. The single negative value of (e 1 − e 2 ) at 593 K is presumed to be an artefact. (c) Comparison of I (0001) with the variations of 'e 3 scaled', showing that part of the e 3 strain has a similar variation with q 2 m to that of e 3 in FeS. (d) Variations with temperature of (I (0001) ) 1/2 , representing q m , and cos β * 1/4 , representing q v . cos β * 1/4 has a steep discontinuity between 603 and 593 K, while (I (0001) ) 1/2 varies continuously and has a tail due to short range ordering above 600 K.
indicating that the distortion from hexagonal geometry was small, ie that coupling of q 2 m with the symmetry-breaking shear strain (e 1 − e 2 ) was weak. The two sets of lattice parameters are reproduced in gure 2, with the data for FeS scaled so that they overlap with data for the hexagonal cell of Fe 7 S 8 above 600 K. It is worth noting that this scaling produces close overlap for the two samples when q v = q m = 0, implying that their evolution ahead of the phase transitions is almost identical, i.e. including precursor effects due to uctuations.
Also shown in gure 2 are baselines of the form a o = a 1 + a 2 Θ s coth(Θ s /T) [26][27][28][29][30] which were obtained by tting to values for a, c and V of FeS in the interval 693-873 K. The saturation temperature, Θ s , was xed at 150 K in order to provide a semi-quantitative representation of the physically correct form of variation as T → 0 K. The a parameter of FeS decreases through T N while c increases, giving the continuous variations of linear strains e 1 and e 3 shown in gure 2(d).
Although not required by symmetry, |e 3 | is almost equal to |e 1 |. Overall, there is a negative volume strain, V s . The intensity of the (0001) magnetic superlattice re ection, I (0001) , in powder neutron diffraction is expected to scale with the square of the magnetic order parameter, and data of Andresen et al [31] show that it also scales with e 3 ( gure 2(e)). This is consistent with e 3 (∝ e 1 ∝ V s ) ∝ q 2 m . Fe 7 S 8 differs from FeS by having additional negative and apparently discontinuous reductions in a, c and V at the transition point ( gures 1(a)-(c)). Figure 2(f) highlights this aspect of variations of V and c which have been scaled to show that the increase in c below the transition point is preceded by the same form of reduction, though it is smaller in magnitude.
Variations of non-symmetry-breaking strains, e 1 , e 2 , e 3 , V s , and symmetry-breaking shear strains (e 1 − e 2 ) and cos β * for Fe 7 S 8 , using the lattice parameters of Powell et al [23] and baselines from FeS, are shown in gures 3(a) and (b), respectively. Each of e 1 , e 2 , e 3 and V s (≈e 1 + e 2 + e 3 ) display the same steep and apparently discontinuous reduction between 603 and 583 K, followed by smooth variations with further lowering of temperature. (e 1 − e 2 ) remains small at all temperatures, consistent with the weak coupling already inferred  [34]. (c) Expanded view of a section of (b), showing that the 5C structure does not t the tell-tale peaks at low diffraction angles.  [23] for the 4C pyrrhotite structure and the model of Liles and de Villiers [34] for the 5C structure.
Recovery of c following the initial steep reduction can be separated out as a strain parameter if the overall volume change is assumed to contain a component in the c-direction which continues to scale with a and b in the same way down to lower temperatures. A strain parameter 'e 3 scaled' has therefore been obtained by rst scaling the average values of e 1 and e 2 so that they overlap with e 3 from 773 down to 593 K, and then subtracting this scaled average from the observed values of e 3 . Figure 3(a) shows that 'e 3 scaled' is indistinguishable from e 3 for FeS, while gure 3(c) shows that it also varies in proportion to I (0001) from gure 4 of Powell et al [23]. These relationships indicate, rstly, that magnetic ordering in Fe 7 S 8 follows the same temperature dependence as in FeS and, secondly, that the strength of coupling of the magnetic order parameter with strain in the crystallographic c-direction is independent of composition.
Powell et al [23] reported a hysteresis of ∼5 K between the P6 3 /mmc-C2 ′ /c ′ transition temperature measured during heating and cooling of their synthetic sample, while Herbert et al [18] have reported a hysteresis of 2 K for a natural sample. This is consistent with the rst order character for the transition they inferred from the steep changes in lattice parameters between 603 and 593 K. However, there is a marked contrast between the evolution of the two order parameters, as highlighted in gure 3(d), which shows (I (0001) ) 1/2 and, hence q m , varying continuously through the transition point, while cos β * 1/4 , representing q v , shows the rst order discontinuity at ∼600 K. Direct measurement of magnetisation also appears to show a continuous variation of the saturation remanence [32].

Sample characterisation
The natural pyrrhotite sample selected for detailed investigation was a shiny, bronze coloured, cm-sized cluster of single crystals with well-developed crystal faces. It came from the mineral collection of the South Australia Museum and had originated from a mine in Mexico.
One large single crystal had faces consistent with a hexagonal growth habit. Indexing of the faces as (001) and (100) was con rmed by x-ray diffraction. Samples with faces cut perpendicular to [001], [100] and [120], as de ned with respect to hexagonal axes, were prepared from this using a ne annular saw. The sample used for high temperature RUS measurements was nearly in the shape of a rectangular parallelepiped with two pairs of parallel faces 1.7 mm and 1.5 mm apart and a third pair of non-parallel faces separated by between 1 mm and 1.5 mm. A small piece of the original single crystal with oriented faces was mounted and polished for chemical analysis using a Cameca SX100 electron microprobe in the Department of Earth Sciences, University of Cambridge (instrumental conditions: 20 keV, 10 nA, 1 µm beam diameter). With respect to eight sulphur atoms, the average of eight analyses gave the number of Fe atoms as 7.00 ± 0.06. Co, Cu, Mn, Ni and Zn were measured for, but all were below the detection limit on all of the eight analysis spots.

X-ray diffraction
An x-ray powder diffraction pattern was obtained from a ground up piece of the original pyrrhotite sample, using Bragg-Brentano geometry on a D8 Bruker diffractometer equipped with a primary Ge monochromator, Cu Kα 1 radiation and a Sol-X solid state detector. Conditions for the data collection were: 5-70 • 2θ, 0.02 • step size, 12 s/step, 0.2 mm divergence slit, 0.2 mm receiving slit, rotating sample. Rietveld re nements were performed with software Topas v6. [33], starting with 4C [23] and 5C [34] structures of Fe 7 S 8 retrieved from the Inorganic Crystal Structure Database [35]. NIST standard LaB6 660b was used to model the instrumental parameters, following a fundamental parameters approach [36]. A shifted Chebyshev function with six parameters and pseudo-Voigt TCHZ function were used to model the background and the peak shape respectively. Figure 4(a) shows the match of the re ned 4C structure, based on the model of Powell et al [23], with the observed diffraction pattern (C2/c, r wp = 26.98, r exp = 21.2, χ 2 = 1.27). Most natural samples contain a mix of 4C and 5C phases. In order to check whether this was the case for the sample used in this study we also performed a re nement on the basis of a 5C structure. Figure 4(b) shows the same comparison between observations and a re nement starting from the 5C structure given by Liles and de Villiers [34]. (P2 1 , r wp = 38.60, r exp = 21.2, χ 2 = 1.82). Lattice parameters from both re nements are given in table 2. 4C and 5C superstructures share the same NiAs substructure so that most of the diffraction peaks, especially at higher angles, will be closely similar. The low angle peaks correspond to large d-spacings and provide the most important means of distinguishing between different superstructures. Re nement on the basis of the 5C structure produces a worse t to the data, statistically, and does not match the low angle peaks ( gure 4(c)). Re nement on the basis that both structures were present increased the number of tting parameters but did not lead to a signi cant improvement in the t (χ 2 = 1.26) in comparison with the result for the 4C structure alone (χ 2 = 1.27). It has therefore been concluded that any 5C pyrrhotite in the sample was below the level of detection.

Microstructure
The con guration of ferroelastic twins due to a change in point group symmetry 6/mmm → 2/m is accompanied by shear strain which, in tensor form for one orientation variant, can be expressed as [37] The transition will produce six twin orientations with W and W ′ twin walls in three pairs, parallel to (1, √ 3,0) and (−3p/2, √ 3p/2,−q), parallel to (100) and (001), and parallel to (−1, √ 3,0) and (3p/2, √ 3p/2,q), using Miller indices de ned with respect to orthogonal reference axes [37]. From the lattice parameter data shown in gure 2, values of |p| and |q| are ∼1%. This means that most of the twin walls arising from vacancy ordering will be aligned nearly parallel to the [001] * , though there will also be one set aligned parallel to (001). Because of the increase in the repeat distance parallel to c by a factor of four, there will also be antiphase domains with boundaries that involve translational errors of the ordered structure [38]. Segments of the boundaries between the antiphase domains which are oriented parallel to (001) can be non-conservative in the sense of having a local change in Fe:S ratio [38].
Electron back-scattered diffraction (EBSD) observations on a freshly cleaved surface of a pieces of the original sample used in the present study con rmed the presence of ferroelastic twin walls, as shown in gure 5(a). In this image, twins with three distinct orientations can be observed. The orientations are as expected for a hexagonal system, i.e. they are separated by 60 degrees in a pole plot in the reference axes of the low symmetry monoclinic phase ( gure 5(b)). The twin individuals are separated by straight walls and have widths of a few micrometres up to several tens of micrometres. The EBSD data have a resolution of 0.4 µm and were collected on an FEI Quanta 650FEG SEM equipped with a Bruker e-Flash HR EBSD detector, operating at 20 kV and beam spot size 5.8 (beam current ∼10 nA). The detector resolution was 320 × 240 pixels, while working distance and sample to detector distance were 26.4 mm and 13.0 mm respectively. The data collection and indexing was perfomed with Bruker QUANTAX CrystaAlign software. Data were analysed using MTEX V5.0.3 [39], a freeware toolset for the commercial software package MATLAB. The misorientation angle precision was estimated to be ∼0.2 degrees using a linear t of a misorientation pro le from the dataset [40] . There must also be six possible directions for the ferrimagnetic moment in the monoclinic structure. In the accompanying paper [20], the orientation of the net magnetic moment is used to con rm the presence of multiple magnetic twins, with a net moment that lies in the (001) plane. In order to visualise the domain structure, we used the Bitter technique, in which a thin layer of a magnetic colloidal suspension, or ferro uid, is applied to the surface and the magnetic particles are attracted to the high eld gradients present at the domain walls [41]. Bitter patterns from the polished face of the sample used in the electron microprobe analysis are shown in gures 6(a) and (b). Previous reports of magnetic 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 [42][43][44][45][46]. 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.

Resonant ultrasound spectroscopy
RUS involves excitation of acoustic resonances of a mmsized sample and measurement of the resulting spectra [47][48][49][50][51]. Individual resonance modes are dominated by shearing motions so that the results for their frequency, f, and peak width at half height, ∆f , provide information predominantly relating to shear elastic constants. The elastic constants scale with f 2 , and acoustic loss is expressed in terms of the inverse mechanical quality factor, Q −1 , which is usually taken to be ∆f /f . In the experimental set up used in the present study [54], the parallelepiped sample was held lightly across its corners between alumina buffer rods inserted into a horizontal tube furnace. Temperature is measured using a thermocouple located within a few mm of the sample and a nal calibration is made based on the transition temperature of quartz.
Spectra were collected in an automated heating sequence from room temperature up to 875 K, with a settle time of 20 min to allow thermal equilibration at each set point. Variations of f 2 and Q −1 extracted from the primary spectra are shown in gure 7. Some oxidation of the sample occurred at the highest temperatures, so that it was not possible to follow the reverse sequence during cooling. f 2 of a resonance peak with frequency near 625 kHz at room temperature decreased smoothly (elastic softening) with increasing temperature. The peak broadened and disappeared altogether in the interval 611-684 K, signifying stronger attenuation. Relatively sharp peaks (Q −1 ∼ 5 × 10 −3 ) reappeared in the primary spectra at 684 K and above, with a net elastic softening of ∼8% across the transition point. Data for a resonance peak with frequency near 1420 kHz at room temperature allowed the pattern of f 2 variations through the transition temperature to be more clearly de ned. There is a minimum in f 2 at 614 K, with softening by ∼5% between temperatures below and above this. As de ned with respect to the parent P6 3 /mmc structure, represented by extrapolation of a straight line to data above 684 K (dashed line in gure 7), the C2 ′ /c ′ structure stiffens abruptly at the transition point and then more gradually with falling temperature. The overall pattern for the change in f 2 is closely similar to the form of variation shown by the strains V s ( gure 3(a)) and cos β * ( gure 3(b)). The most striking thing about the evolution of Q −1 is that it is rather low across the whole temperature range. There may be a small increase in Q −1 in the vicinity of the phase transition, but it is never larger than fractions of one percent.

Discussion
The Néel temperature for FeS, troilite, is ∼590 K (e.g. reference [52] and references therein), and the temperature of the onset of magnetic ordering is essentially the same across the Fe 1−x S solid solution to at least Fe 7 S 8 , irrespective of whether or not it is accompanied by vacancy ordering [5,12,53]. Based on the neutron diffraction data of Powell et al [23], the onset temperatures for magnetic ordering and vacancy ordering are indistinguishable, though the two order parameters have markedly different trajectories. The magnetic ordering transition is essentially the same as occurs in FeS with respect to transition temperature, the temperature dependence of q m and the strength of coupling between q m and strain. Strains that re ect the variation of q v have a discontinuity at the transition point while q m varies continuously through it ( gure 3(d)). The overall behaviour is of a long-range magnetic ordering transition accompanied by signi cant secondary vacancy ordering that gives the rst order character. Superposition of the two contributions in this way is consistent also with the form of the excess heat capacity shown by Grønvold and Stølen [53] which, for Fe 7 S 8 (their gure 6), looks like the excess heat capacity for FeS (their gure 2) combined with an additional latent heat.
Strictly speaking, conventional Landau potentials refer to phase transitions at the displacive limit and will not reproduce the temperature dependence of the two order parameters in detail for order/disorder processes. However, an expansion of the form of equation (1) is expected to at least give the correct overall form for the observed variation of q m with respect to Figure 8. Variation of (I (0001) ) 1/2 , representing the variation of q m , with respect to cos β * 1/4 , representing the variation of q v . The initial increase in (I (0001) ) 1/2 at cos β * 1/4 = 0 is due to short range ordering above the transition point. This is followed by an apparently discontinuous increase in cos β * 1/4 , followed by a smooth variation of both parameters down to 30 K. Between 593 and ∼300 K the variation is linear.
q v , as represented by the variations of (I (0001) ) 1/2 and cos β * 1/4 in gure 8. The lowest order term for coupling between q m and q v is biquadratic and is most likely to be the result of a common strain mechanism. Coupling would be favourable via (e 1 + e 2 ) since both order parameters induce negative values for this strain, but unfavourable via e 3 because q m causes the cdimension to increase and q v causes it to decrease. The volume strain is negative in both cases, however, so that the coupling is likely to be favourable overall.
Generic models of biquadratic order parameter coupling based on two 246 potentials were treated in some detail by Salje and Devarajan [54]. If the coupling is weak, two successive transitions are expected unless the transition temperatures are coincident. Here we propose that there is one transition temperature set by the magnetic ordering. The vacancy ordering then takes place due to coupling with the magnetism through common strain. It is certainly not the case that q m varies independently of q v in Fe 9 S 10 because abrupt changes in temperature produce relatively slow changes in moment as the degree of Fe/vacancy order re-equilibrates [55]. Analogous kinetic evidence of coupling, also in Fe 9 S 10 , has also been reported by Marusak and Mulay [56] and Herbert et al [18].
The models of Salje and Devarajan [54] for strong coupling included only situations where the coef cient of the fourth order term was positive for both order parameters and did not produce any pattern of variations resembling that shown in gure 8 for the C2 ′ /c ′ structure of Fe 7 S 8 . Bilinear coupling between two order parameters, i.e. with the form λq 1 q 2 , and linear-quadratic coupling, i.e. of the form λq 1 q 2 2 , can result in two separate transition temperatures being renormalised to one [54,57], but neither are allowed in this case. Instead, there is an unusual higher order coupling with the form λq 2 m q 8 v that arises as a consequence of the coupling terms λe 4 q 4 v and λe 2 4 q 2 m in equation (1). The implications of this unusual term have not been explored but, given that the magnetic ordering is the same as in FeS (without vacancy ordering), the outcome represented in gure 8 is presumably that symmetry breaking is driven by the magnetic instability and that coupling of the strain with the magnetic order parameter generates a eld which induces vacancy ordering. The normal expectation for changes in elastic constants associated with phase transitions is that softening will occur in the low symmetry phase due to coupling of strains with the driving order parameter [24,58,59]. The paramagnetic-antiferromagnetic transition in CoF 2 provides an example of this classical mechanism arising from coupling of the form λeq 2 [60]. In Fe 7 S 8 , elastic stiffening is observed instead of softening ( gure 7), but this is readily understandable in terms of the slow timescale of the relaxation of the order parameter with respect to the timescale of the applied stress. Changes in vacancy order give rise to a large shear strain, e 4 , but applying a stress at a frequency of ∼1 MHz, does not allow time for q v to relax. Instead, the elastic constants depend primarily on the second derivative of the excess free energy with respect to strain, δ 2 G/δe 2 , and therefore scale with combinations of q 2 m and q 2 v due to terms of the form λe 2 q 2 . For example, C 44 will soften or stiffen according to C o 44 + 6λ em3 q 2 m + 2(λ ev5 + 3λ ev6 )q 2 v , and similarly for the other elastic constants. The step in f 2 at ∼600 K and further non-linear increases with falling temperature are thus expected to follow the variation of V s , which scales with q 2 m and q 2 v , as appears to be the case. Stiffening rather than softening re ects positive values of λ in each of the λe 2 q 2 terms.
When measured at RUS frequencies, the mobility of ferroelastic twin walls in perovskites through some interval below a ferroelastic transition temperature and above the pinning/freezing temperature of the walls gives rise to a characteristic pattern of high attenuation [61]. The P6 3 /mmc-C2 ′ /c ′ transition in Fe 7 S 8 is, in principle, ferroelastic but switching of the resulting ferroelastic domains in response to an applied stress requires local reordering of the vacancies at the twin walls and will be slower by many orders of magnitude than switching of domains in materials where the ferroelasticity arises from a displacive phase transition. Even close to the transition temperature, it seems unlikely that this mechanism for acoustic loss will operate and the observed high attenuation is more likely to be intrinsic. Given that there is signi cant coupling between strain and the magnetic order parameter at a macroscopic scale, one source of the high loss would be changes in local spin state in response to the dynamic stress within a mechanical resonance of the sample. The degree of loss is highest close to and immediately below the transition point and, instead of there being evidence Debye peak in Q −1 to mark a discrete freezing event, the amplitude of the loss process then decreases smoothly with falling temperature. Q −1 values remain high at temperatures above the transition point, signifying the presence of strain coupling to some local dynamical changes, perhaps of vacancy distribution in the disordered structure. There is no strain contrast across 180 • magnetic domain walls, so they cannot contribute directly to the observed variations in Q −1 .

Implications for palaeomagnetism
A key feature of pyrrhotite in the context of preserving palaeomagnetic signals over long periods of geological time is that the ferrimagnetism is coupled with Fe/vacancy ordering. The coercivity of the C2 ′ /c ′ structure might be weak but the orientation of the net moment is likely to return to its original orientation in a crystal even after the application of a large external eld because the underlying pattern of vacancy ordering will remain xed unless the temperature is raised sufciently to allow Fe/vacancy diffusion to occur. In other words, the orientation of a remanent signal is xed by the underlying pattern of ferroelastic domains. However, although this should give long term stability, it raises the question of what initial signal is recorded. In the limiting case, some arbitrary con guration of domains with Fe/vacancy order could determine the net magnetisation independently of the magnetic eld in which the crystal grew. To address this issue it will be necessary to develop a formal model of the P6 3 /mmc-C2 ′ /c ′ transition, to include the form and strength of coupling between q v and q m , together with the kinetics of vacancy ordering.

Conclusions
Analysis of the hexagonal-monoclinic transition in a natural sample of pyrrhotite with composition close to Fe 7 S 8 from the perspectives of strain and elasticity has shown: (a) There is signi cant coupling of strains with the order parameters for vacancy ordering and magnetic ordering. Monoclinic lattice geometry occurs because of strong coupling between the order parameter for vacancy ordering and the shear strain e 4 . By itself, magnetoelastic coupling would cause only small deviations from hexagonal geometry. (b) Magnetic ordering appears to have the same continuous temperature dependence as in FeS, and the strength of coupling with strain also appears to be the same. First order character for the transition is due to the vacancy ordering. An additional and unusual coupling term, λq 2 m q 8 v has emerged from considerations of symmetry and perhaps has a role in determining that magnetic ordering and vacancy ordering occur simultaneously but, individually, have markedly different trajectories with changing temperature. (c) Although the transition is, in principle, ferroelastic, the dynamical response to applied stress depends on the relatively long time scales required for Fe/vacancy diffusion in comparison with spin reorientation. As a consequence, changes in elastic and anelastic properties measured at RUS frequencies involve stiffening rather than softening and a pattern of acoustic loss which differs from the classical Debye freezing of ferroelastic twin walls seen in perovskites.

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).