Pyrrhotites, characterized by the chemical formula Fe1–δS (0 < δ ≤ 1/8), represent an extended group of minerals that are derived from the NiAs-type FeS aristotype. They contain layered arrangements of ordered Fe vacancies, which are at the origin of the various magnetic signals registered from certain natural rocks and can act as efficient electrocatalysts in oxygen evolution reactions in ultrathin form. Despite extensive studies over the past century, the local structural details of pyrrhotite superstructures formed by different arrangements of Fe vacancies remain unclear, in particular at the atomic scale. Here, atomic-resolution high-angle annular dark-field imaging and nanobeam electron diffraction in the scanning transmission electron microscope are used to study natural pyrrhotite samples that contain commensurate 4C and incommensurate 4.91 ± 0.02C constituents. Local measurements of both the intensities and the picometer-scale shifts of individual Fe atomic columns are shown to be consistent with a model for the structure of 4C pyrrhotite, which was derived using X-ray diffraction by Tokonami et al. (1972). In 4.91 ± 0.02C pyrrhotite, 5C-like unequally sized nano-regions are found to join at anti-phase-like boundaries, leading to the incommensurability observed in the present pyrrhotite sample. This conclusion is supported by computer simulations. The local magnetic properties of each phase are inferred from the measurements. A discussion of perspectives for the quantitative counting of Fe vacancies at the atomic scale is presented.


Pyrrhotite (Fe1–δS; 0 < δ ≤ 1/8) is one of the most common metal sulfide minerals in the Earth's ore deposits, as well as in a range of meteorites (Rochette et al. 2001, 2005; Weiss et al. 2002; Lorand et al. 2005; Yu and Gee 2005; Louzada et al. 2007). It has attracted extensive research interest over the past century and has been considered as “geologically interesting, industrially important, and crystallographically complex” (Pierce and Buseck 1974).

The complexity of pyrrhotite results from the spatial arrangement of Fe vacancies, which are present in the corresponding aristotypic structure of stoichiometric FeS. At high temperature (=483 K), pure FeS has a NiAs-type hexagonal structure with a = b ≈ 0.344 nm and c ≈ 0.588 nm (Töpel-Schadt and Müller 1982; Kruse and Ericsson 1988), consisting of alternating layers in the c plane of pure Fe or pure S. Upon cooling below ~413 K it transforms to the room-temperature-stable structure of troilite with a (3a,3a,2c)unit cell (Töpel-Schadt and Müller 1982), which exhibits macroscopic antiferromagnetism due to the anti-parallel alignment of the ferromagnetically coupled spins of the Fe2+ atoms that are located in successive basal planes perpendicular to the c axis (Ricci and Bousquet 2016). In pyrrhotite, the slight change in composition from stoichiometric FeS due to the incorporation of vacancies and the resulting vacancy ordering lead to the formation of different superstructures, which are referred to as NC pyrrhotites and are typically characterized by unit cells that are multiples of the parent NiAs-type subcell. In this notation, N is the coefficient of the c lattice parameter [see also Pósfai et al. (2000) for a description of the structural relationships].

Depending on the superstructure, pyrrhotite can exhibit ferrimagnetism [e.g., in the 4C end-member Fe7S8 (δ = 1/8) (Powell et al. 2004)], making it an important carrier of magnetic remanence in both the Earth's crust (Kontny et al. 2000; Pósfai et al. 2000) and extraterrestrial materials (Rochette et al. 2001, 2005; Keller et al. 2002; Weiss et al. 2002; Lorand et al. 2005; Yu and Gee 2005; Louzada et al. 2007). Pyrrhotite also possesses peculiar magnetic properties, including for instance the Besnus transition (Besnus and Meyer 1964), which was only recently explained by out-of-plane shifts of the magnetic moments of Fe atoms close to vacancies (Koulialias et al. 2018c, 2019). A recent study has also shown that ultrathin nanosheets of 4C pyrrhotite can act as efficient electrocatalysts in oxygen evolution reactions because of the presence of active Fe vacancy sites of mixed valence (Chen et al. 2017).

Pyrrhotite can be grouped into two categories, depending on the constituent superstructures: commensurate pyrrhotites with integer values of N and incommensurate (or intermediate) pyrrhotites with non-integer values. The chemical formulas of commensurate NC superstructures take the form Fen–1Sn (n ≥ 8), where N = n/2 for n = even and N = n for n = odd (Elliot 2010). In incommensurate pyrrhotites, the Fe concentration is intermediate between that in the commensurate counterparts. Since differences in Fe concentration between different pyrrhotites are small (the atomic percentages are between 46.7 and 50 at% across the entire group), it is difficult to distinguish between superstructures using compositional measurements alone (Morimoto et al. 1975). Based on studies of natural samples, four stable phases have been proposed: 4C (Fe7S8; δ = 1/8), 5C (Fe9S10; δ = 1/10), 11C (Fe10S11; δ = 1/11), and 6C (Fe11S12; δ = 1/12) (Carpenter and Desborough 1964; Corlett 1968; Morimoto et al. 1970; Vorma 1970; Koto et al. 1975). The incommensurate pyrrhotites appear to be metastable and can take the form of intergrowths with stable pyrrhotites (e.g., Nakazawa and Morimoto 1971; Harries et al. 2011). Hereafter, NC denotes uniquely incommensurate pyrrhotites.

Pyrrhotite superstructures have been investigated intensively using several reciprocal-space-based techniques, including X-ray diffraction (Corlett 1968; Vorma 1970; Tokonami et al. 1972; Morimoto et al. 1975; de Villiers and Liles 2010; Liles and de Villiers 2012), electron diffraction (Dódony and Pósfai 1990; Pósfai et al. 2000), and neutron diffraction (Powell et al. 2004), as well as using real space methods, including high-resolution (HR) transmission electron microscopy (TEM) (Pierce and Buseck 1974; Morimoto et al. 1974; Nakazawa et al. 1975; Harries et al. 2011). The most thoroughly studied superstructure to date is monoclinic 4C pyrrhotite (a = 1.1902 nm, b = 0.6859 nm, c = 2.2787 nm, and β = 90.4° based on the non-conventional F2/d space group (Tokonami et al. 1972); see also ICSD 42491). In this phase, Fe vacancies are arranged in layers parallel to the cm plane [i.e., the (001)m plane, where “m” refers to monoclinic] and alternate with full Fe layers along the cm axis (Bertaut 1953; Tokonami et al. 1972). Powell et al. (2004) performed refinements of the 4C pyrrhotite structure and described it as belonging to the conventional space group C2/c. The stacking sequence in the direction normal to the basal planes can be written …FAFDFBFC…, where A, B, C, and D are layers that have spatially different configurations of Fe vacancies, while F describes a fully occupied Fe layer [see Figs. 1a and 1c in (Elliot 2010) for schematic diagrams].

In contrast to 4C pyrrhotite, studies of other pyrrhotite phases are less conclusive. For example, in an early study of 5C pyrrhotite it was proposed that its structure could be created by adding two fully occupied Fe layers to the 4C structure if these layers are arranged as far apart as possible (Corlett 1968), which would correspond to stacking sequences such as …FAFBFFCFDF… (Vaughan et al. 1971). Although this sequence was also reported by Dódony and Pósfai (1990), the fact that the spacings in the diffraction pattern were uneven suggests that the structure was probably incommensurate. In further investigations by Elliot (2010), de Villiers et al. (2009), and Liles and de Villiers (2012), 5C structures were described by models that contained partially occupied vacant sites, i.e., …FA0.5A0.5C0.5C0.5FD0.5D0.5B0.5B0.5… for Fe9S10 (Elliot 2010), …A0.914A0.285A0.891B0.484B0.468C0.914C0.285C0.891D0.484D0.468… for Fe9.021S10 (see also ICSD 190011), and …B0.938C0.302(C0.803B0.947)A0.598(A0.475C0.950)A0.938D0.302(D0.803A0.947) B0.598(B0.475D0.950)… for Fe9.007S10 (see also ICSD 190012) (Liles and de Villiers 2012). In this notation, subscripts denote the occupancies of the Fe positions described by Elliot (2010), while symbols in brackets indicate that they are located in the same layer. Similar results have been reported for 6C pyrrhotite (Nakazawa and Morimoto 1971; Koto et al. 1975; de Villiers and Liles 2010).

Pierce and Buseck (1974) studied incommensurate pyrrhotites using HRTEM dark-field imaging and described the superstructures in terms of disordered sequences of anti-phase domains. This description was extended by Harries et al. (2011), who proposed a translation interface modulation (TIM) model on the basis of electron diffraction and HRTEM dark-field imaging results. According to the work of van Smaalen (1995), incommensurate phases can be grouped into three categories: incommensurately modulated structures, incommensurate intergrowth compounds and quasicrystals. The studies by Pierce and Buseck (1974) and Harries et al. (2011) suggested that NC pyrrhotites are incommensurate intergrowth compounds. However, since they did not consider structural features, i.e., displacive shifts (e.g., Ge et al. 2015) and occupancies of individual atomic columns (e.g., van der Lee et al. 1991) that are decisive for identifying incommensurately modulated structures (van Smaalen 1995), the possibility that NC pyrrhotites are incommensurately modulated structures cannot be excluded. A four-dimensional (4D) superspace formalism has been developed for the description of all incommensurate pyrrhotites (Yamamoto and Nakazawa 1982; Izaola et al. 2007). However, it takes the form of graphical models, which are difficult to compare with real structures. Its validation also requires a large number of crosschecks with high-quality experimental data, which are still rare. Thus, a detailed understanding of pyrrhotite superstructures at the atomic level is still far from complete.

The primary difficulty in characterizing pyrrhotite superstructures results from the intrinsic complexity of the Fe vacancy ordering, in particular in incommensurate superstructures. There are also fundamental limitations to obtaining large single-phase crystals for studies using reciprocal-space-based techniques. Since results that are obtained using the latter techniques are typically averaged spatially, their interpretation is also not straightforward. The structural complexity of pyrrhotite can then lead to multiple interpretations of the same observation (Koulialias et al. 2018b), or to oversimplified explanations (Dódony and Pósfai 1990). Earlier HRTEM results (Pierce and Buseck 1974; Morimoto et al. 1974; Nakazawa et al. 1975; Harries et al. 2011) were usually not directly interpretable, due to the effect of lens aberrations. To the best of our knowledge, the only image that has shown Fe vacancy ordering with atomic spatial resolution was obtained from 4C pyrrhotite using high-angle annular dark-field (HAADF) scanning TEM (STEM) (Xu et al. 2015). Further quantitative analysis, however, was not carried out in the latter study.

It is important to note that relevant structural details in pyrrhotite include not only the positions of Fe vacancies within the FeS framework, but also the Fe atom shifts that are associated with the presence of these vacancies. Both parameters are important for understanding incommensurability (van Smaalen 1995). However, the importance of the latter measurements has been underestimated, despite the fact that a limited number of studies have been carried out using X-ray and neutron diffraction (e.g., Tokonami et al. 1972; Powell et al. 2004; Liles and de Villiers 2012; Koulialias et al. 2018a). For example, according to the model proposed by Tokonami et al. (1972), Fe atoms in 4C pyrrhotite are not expected to line up coincidently when they are viewed along major zone axis directions perpendicular to the cm axis (see, e.g., atomic columns 1–4 in Figs. 1a1c and Supplemental1 Fig. S1), which leads to the observed column splitting in the projected plane (i.e., the atomic columns are projected as atom pairs, as can be clearly seen, e.g., for column 4). The Fe sublattice is then no longer rectangular, as it is in NiAs-type FeS. However, direct local experimental observations are lacking.

Here, aberration-corrected STEM is used together with image simulations to perform an atomic-scale study of a natural pyrrhotite specimen that contains both commensurate 4C and incommensurate 4.91 ± 0.02C superstructures. Quantitative analyses of image intensities and measurements of shifts of the positions of individual Fe atomic columns are used to characterize the superstructures locally with atomic spatial resolution. The present study of the 4.91 ± 0.02C superstructure provides a high-quality model that can be used to assess the use of the 4D superspace formalism for pyrrhotite, as well for other omission/defect structures.

Materials and methods

Sample origin and preparation

The natural pyrrhotite crystal examined in this study was obtained from the ETH Zurich mineral collection (No. 5573). It was obtained from the locality Auerbach in Hesse, Germany (Charilaou et al. 2015; Koulialias et al. 2016) and had a crystal morphology that can be described schematically as a hexagonal prism with (pseudo-) hexagonal basal planes (Supplemental1 Fig. S2). Lamellar TEM specimens were prepared from its polished top surface using focused ion beam (FIB) milling with Ga in an FEI Helios Nanolab 400s dual beam system (ER-C 2016a). The lamellae were cut in the following two directions for systematic investigations: cross-sectional cuts with a normal of [100]m and 30° pre-tilted cuts with a normal parallel to [110]m (or [110]m), as shown in Supplemental1 Figure S2 in red and green, respectively.

The as-prepared FIB lamellae were thinned using 2.5 kV Ar ion milling in a Bal-Tec Res-120 system, followed by final cleaning using a focused 500 eV Ar ion beam (Fischione Nanomill, Model 1040), to remove damaged layers that may have been introduced during the previous milling steps. Final sample thicknesses of between 30 and 50 nm could be achieved without detectable ion beam damage. However, the possibility of redeposition of Cu from the Omniprobe grid during ion milling could not be excluded. The lamellae were stored in a vacuum environment and cleaned for 2 × 5 s (with an interval of 1 min) using a shielded specimen holder port and an Ar and O2 gas mixture in a Fischione plasma cleaner (model 1020), to minimize surface oxidization during plasma cleaning and the build-up of carbon contamination during atomic-resolution STEM.

Measurement techniques, imaging conditions, and analysis methods

Selected-area electron diffraction and energy-dispersive X-ray spectroscopy (EDXS) were carried out at 200 kV in an FEI Tecnai F20 microscope equipped with a field emission gun (FEG) and an EDXS attachment (ER-C 2016b). The areas of the specimen that were analyzed had diameters of ~200 nm. Kinematical simulations of electron diffraction patterns were generated using SingleCrystal software (CrystalMaker Software Limited, U.K.).

Distributions of 4C twin variants and other superstructures in the samples were studied using scanning nanobeam electron diffraction (Gammer et al. 2016; Jin et al. 2016) at 300 kV in an FEI Titan 80-300 microscope equipped with a probe spherical aberration (CS) corrector (ER-C 2016c). The semi-convergence angle and the diameter of the electron probe were 0.5 mrad and ~3 nm, respectively. The probe was scanned across the sample and used to generate an annular dark-field STEM image with a size of 1024 × 1024 pixels. At each position, a nanobeam electron diffraction pattern was recorded using a dwell time of 0.1 s, to obtain a 4D data set, which was processed using principal component analysis (PCA) (Hotelling 1933a, 1933b). Spectral unmixing (SU) (de la Peña et al. 2011) was used to identify and separate the contributions of different twin variants and superstructures to the diffraction patterns. A non-negative matrix factorization algorithm for SU was used, since the resulting non-negative principal components are correlated with the most prominent intensity features in the diffraction patterns. As the present analysis focused on the positions of the diffraction spots, their intensities (i.e., count rates) were equalized prior to SU. This approach allows the principal components to be related to changes in the positions and not the intensities of the diffraction spots. Furthermore, the shift of the Laue center, which was caused by unavoidable electron beam tilt when scanning the probe across a large field of view, as well as sample tilt, was compensated by aligning the diffraction patterns relative to each other using common alignment procedures (Schaffer et al. 2004). A brief workflow is given in the Supplemental1 Information. Based on the principal components, a diffraction image could be constructed, with color coding used to show the twin and phase distributions.

Atomic-resolution HAADF STEM images were recorded at 200 kV in an FEI Titan G2 80-200 ChemiSTEM microscope equipped with a high-brightness FEG (XFEG), a probe CS corrector, and a super-X EDXS system (ER-C 2016d). The incident electron beam convergence semi-angle for HAADF STEM imaging was ~25 mrad, while the collection semi-angle was 70–200 mrad. During HAADF STEM imaging, “view” mode was used at high magnification (e.g., 20.5 M on an FEI instrument; image window: 512 × 512 pixels; pixel size: 8.6 pm; dwell time: 0.4 μs/pixel) instead of “focus” mode, which would result in a longer dwell time, while tuning the defocus of the image as fast as possible. The microscope was then switched quickly to low magnification (e.g., 1.8 M) to record the image (dwell time: 0.5 μs/pixel). The image size was 4096 × 4096 pixels to achieve a large field of view with sufficient image sampling (e.g., 12.1 pm at 1.8 M). In this way, an image with dimensions of approximately 50 × 50 nm could be obtained, in which the central area (approximately 4.5 × 4.5 nm in size) may be influenced by more prolonged illumination during focusing, while the rest (the majority) can be regarded as a fresh area. For quantitative analysis, the use of the central (approximately 10 × 10 nm) area was generally avoided. Selected HAADF STEM images were denoised using a nonlinear filtering algorithm (Du 2015).

Multislice HAADF STEM image simulations were carried out using Dr. Probe software (Barthel 2018) for the experimental imaging parameters. A series of sample thicknesses up to ~100 nm was considered in the simulations. The sample tilt angle and all aberrations of the microscope were set to zero. Rigid supercell models based on the structure of 4C pyrrhotite (Tokonami et al. 1972; ICSD 42491) were built using home-made scripts. Structural models were visualized using VESTA software (Momma and Izumi 2011).

The mean intensities of individual atomic columns in experimental and simulated images were measured using iMtools software (Houben 2009), following a procedure similar to that described elsewhere (Barthel 2018). Specifically, atomic column positions were measured by fitting two-dimensional Gaussian functions to the image intensity with sub-pixel precision. The intensities of pixel positions in the circular area assigned to each atomic column were then summed and divided by the measured numbers of pixels in this circular area.

STEM image simulations of the detectability of vacancy ordering in 4C pyrrhotite

The present study relies on an a priori knowledge of the structure of 4C pyrrhotite. In this section, the structure model for 4C pyrrhotite proposed by Tokonami et al. (1972) is used to perform simulations that assess the detectability of Fe vacancy ordering in HAADF STEM images. In comparison with the model proposed by Powell et al. (2004), in which the cell is non-orthogonal, the present model, which is based on a pseudo-orthorhombic cell, provides better visibility of vacancy ordering. Atomic-resolution HAADF STEM images are considered because they have superior tolerance for the examination of the thicker samples that are studied here. The use of this imaging mode is crucial to the present research, as discussed in “Atomic-resolution imaging of 4C pyrrhotite and Fe column shifts.”

Figures 1a1c show the atomic model projected along the [010]m, [110]m, and [110]m directions of the 4C monoclinic structure, respectively (corresponding to the hexagonal <2110=-type axis in the NiAs-type subcell). The atomic columns that contain vacant Fe sites are denoted Fe(V) and are marked by orange circles. In these columns, every second atomic site is vacant (corresponding to 50% Fe occupancy on average over the whole column), in contrast to columns that are fully occupied and denoted Fe(F), as shown on the right of Figure 1c. The stacking sequence along the cm axis can be described using a similar notation to that used by Pierce and Buseck (1974) as …FVAFVBFVAFVB… in Figure 1a, …FVAFVBFVBFVA… in Figure 1b, and …FVAFVAFVBFVB… in Figure 1c, where F are (001)m atomic planes in which all of the Fe positions are filled (as in the “Introduction”), while VA (VB) are atomic planes that include Fe(V) columns occupying A-type (B-type) positions2. The (001)m interplanar spacing d between neighboring Fe layers is ~0.29 nm (i.e., 1/8 of the c lattice parameter of 4C pyrrhotite labeled in Fig. 1b). Closer examination reveals that the projected structures in Figures 1b and 1c are essentially equivalent and yield identical diffraction patterns, as shown in Figures 1e and 1f. As a result of this equivalence between the ±[110]m and ±[110]m directions, unless otherwise specified, they are referred to as [110]m-type or type II below. In contrast, the diffraction pattern along [010]m shown in Figure 1d is different, containing only one reflection in addition to the primary NiAs-type reciprocal lattice, which is outlined by a red rectangle. It is referred to as type I below.

Simulations of HAADF STEM images were carried out along both the [010]m (type I) and the [110]m-type (type II) direction. Figures 2a2d show simulated images for sample thicknesses of 11 nm (Figs. 2a and 2c) and 88 nm (Figs. 2b and 2d). It is clear from the simulations that the visualization of Fe vacancy ordering depends strongly on sample thickness (Figs. 2a2d).

Figure 2e shows mean intensities I measured from the simulated images for atomic columns of Fe(F) (black: [010]m and red: [110]m-type projection), Fe(V) (blue) and S (yellow) plotted as a function of sample thickness. For the Fe(V) and S columns, the values are essentially the same for both projections. The inset to Figure 2e shows the ratio IFe(V)/IFe(F) for sample thicknesses up to ~90 nm. It can be seen that when the sample is thinner than ~10 nm this ratio is only slightly higher than the nominal value of 0.5 used in the simulations (as a result of electron dechannelling, which is discussed below). With increasing sample thickness, this ratio displays a monotonic but nonlinear increase in both projection directions, reaching a value of ~0.9 at a thickness of ~90 nm, resulting in an intrinsic difficulty to resolve Fe vacancy ordering in thicker samples.

It should be noted that, even for fully occupied Fe columns such as 1–4 in Figure 2c (for which the same notation is used as in Fig. 1), the simulated intensities differ from each other (i.e., they are successively lowered from column 1 to column 4) as a result of electron dechannelling caused by differences in the splitting of atomic columns projected in the electron beam direction (by 7.2, ~18, 30.9, and 46.7 pm for columns 1–4, respectively), as shown in the atomic models in Figures 1a1c. This situation results in reduced mean intensities for Fe(F) in the [110]m-type orientation (red line in Fig. 2e), when compared with intensities measured along [010]m (black dotted line in Fig. 2e), as the splitting in the latter case is ~18 pm for all Fe(F) columns (see Fig. 1a). For the same reason, the differences between the Fe(F) intensities are smaller in the simulated image projected along [010]m (Fig. 2a), as evidenced by the smaller error bar (i.e., standard deviation) than in the [110]m-type orientation (Fig. 2e).

The same analysis is shown for hexagonal <1100=-type axes in Supplemental1 Figure S1. Although vacancy ordering can still be discerned, lamellae fabricated in this orientation are not recommended for analysis, as a result of the more compact arrangement of Fe and S atomic columns in this projection direction. In particular, the Fe(V) atomic columns, which are marked by orange arrows in Supplemental1 Figure S1a (identical to our cross-sectional sample; see “Materials and methods”), show a strong shift toward neighboring Fe(F) columns, thereby hindering quantitative analysis.

After taking the above considerations about the sample geometry into account (see “Material and methods” and Supplemental1 Fig. S2), the analysis presented below was performed on lamellae that were prepared in [110]m-type orientations.

Experimental STEM results

Twinning and phase distribution

The occurrence of twinning perpendicular to the cm axis in 4C pyrrhotite is well known (e.g., Bennett et al. 1972; Putnis 1975; Pósfai et al. 2000; Xu et al. 2015) and is found in the present sample. Figure 3a shows an annular dark-field image, in which bands of dark contrast (~300 nm in width) are visible in a brighter matrix. Figure 3b shows a “diffraction” image [analogous to an electron energy-loss spectrum image (Jarausch et al. 2009)], which was obtained from the same area by summing the weight maps of the main diffraction components, as shown in Figures 3c and 3d, based on the results of PCA. The constituent diffraction patterns shown in Figures 3c and 3d are both [110]m-type at first glance. However, by using least-square fitting of the diffraction spots, it was determined that they are associated with 60°-rotated or 180°-rotated 4C twins, between which there is a small angular difference of ~1.1° (Figs. 3c and 3d; see also Figs. 1e and 1f for reference). The different angles can be recognized in the composite diffraction pattern shown in Figure 3e.

Nanoscale 4C [010]m twin variants of limited size (<30 nm), which are marked in yellow in Figure 4, were also observed. They are usually surrounded by a [110]m-type matrix (marked in red in Fig. 4b), resulting in a composite nanobeam diffraction pattern (Fig. 4d) that includes contributions from both the [010]m (Fig. 1d) and the [110]m-type axes (Figs. 1e and 4c). These observations are consistent with the work of Pósfai et al. (2000). Further analysis is presented in Supplemental1 Figure S3.

The abundance of twinning in this sample suggests that the unit cell doubling reported by Koulialias et al. (2018b) may result simply from the superposition of two twin-related zone axes [150]m and [120]m.

In addition to twinning, an extra superstructure (marked in green in Fig. 4b) associated with a 5C-like nanobeam diffraction pattern (Fig. 4e) was found. This superstructure was denoted 5C* by Charilaou et al. (2015) and Koulialias et al. (2016), based on powder X-ray diffraction results. Closer inspection of a selected-area electron diffraction pattern (Fig. 5a) reveals that this 5C-like superstructure is, in fact, an incommensurate 4.91 ± 0.02C phase (4.91C for short; see also Supplemental1 Fig. S4 for indexing of the spots). Its fraction is found to be significant in this sample. It is discussed in detail below.

Statistical EDXS data yielded an Fe concentration of 45.2 ± 1.8 at% in the 4C matrix. These measurements, which are consistent with results reported by Koulialias et al. (2016), do not allow the commensurate 4C and incommensurate 4.91C phases to be distinguished from each other chemically.

Atomic-resolution imaging of 4C pyrrhotite and Fe column shifts

Vacancy ordering in pyrrhotite was found to be sensitive to electron beam irradiation at 200 kV. Figure 6a shows an atomic-resolution HAADF STEM image recorded from a region of 4C pyrrhotite along a [110]m-type axis (i.e., from the red region in Figs. 3 and 4). This image was selected because it contains not only a region of 4C pyrrhotite that had not been illuminated by the electron beam prior to image acquisition, but also a modified region (denoted MA in Fig. 6a), in which vacancy ordering was no longer visible after prolonged electron beam illumination. This single image therefore allows us to study not only the atomic structure of 4C pyrrhotite, but also the consequence of irradiation by a 200 kV electron beam. The boundary between the two regions is marked by a dashed line.

Figures 6b and 6c show magnified versions of the modified and unirradiated 4C regions after image averaging. The right part of Figure 6c is color-coded, to better visualize Fe vacancy ordering. Two features can be seen directly. First, the Fe vacancy ordering, which is indicated by bold orange circles in the unirradiated 4C pyrrhotite region (Fig. 6c), is no longer visible in the modified region. Second, in addition to the disappearance of Fe vacancy ordering, the atomic column shifts in the MA region (marked using thin yellow circles in Fig. 6b) are much less visible than in the 4C region, resulting in a quasi-square arrangement of the Fe sublattice.

Quantitative measurements of the Fe atomic column shifts were performed on a unit-cell-by-unit-cell basis directly from the HAADF STEM image, first for the simulated image shown in Figure 7a (reproduced from Fig. 2c; sample thickness: 11 nm) and then for the experimental image shown in Figure 6a. For clarity, two global axes, z and x, which are always parallel to the 4C cm axis and the axis nearly perpendicular to cm (e.g., [130]m in Figs. 6a and 7a), respectively, are defined. The measured shifts are then given by their components along these two axes, as described below.

The Fe atomic column positions were determined by locally fitting two-dimensional Gaussian functions to the image intensity, as shown schematically in Figure 7b, in which the centers of the marked circles (i.e., the positions of intensity maxima) denote the column positions. Two methods were used to determine the atomic column shifts from the measured positions. In the first method, the deviation (ΔX, ΔZ) of each atomic column (denoted as individual shifts) from the nearest (calculated) site of the average Fe lattice, which is marked by a white grid in Figure 7b, was calculated. The laterally averaged shifts (in the global x direction) were then plotted as a function of Fe layer stacking along the cm axis (in the global z direction). Depending on the periodicity in the x direction, results from the cyan and yellow vertical lines are presented separately in Figures 7c and 7f and Figures 7d and 7g, respectively. In the second method, only the relative shift (Δx, Δz) of each yellow column was measured with respect to the calculated center (marked by cyan dots) of its two closest cyan neighbors marked in Figure 7b. The averaged results are shown in Figures 7e and 7h.

For a given Fe layer, the relative shifts shown in Figures 7e and 7h (i.e., Δx and Δz, respectively) are equal to the difference of the individual shifts between the yellow (Figs. 7d and 7g) and cyan (Figs. 7c and 7f) columns. Hence, the relative method can be regarded as a simplification of the individual approach. It is evident from Figures 7c7h that the tendency of the shift variations along the cm axis, which is the most important factor for understanding incommensurability, is the same for both methods. However, an advantage of using relative shift analysis is that it can reduce the influence of the unavoidable sample drift that is present in experimental STEM images. Therefore, only relative shifts are used hereafter.

Figures 8a8d show laterally averaged (in the global x direction) atomic column shifts Δx (Figs. 8a and 8c) and Δz (Figs. 8b and 8d) plotted as a function of Fe layer stacking sequence along the cm axis (in the global z direction) for the unirradiated 4C and modified regions denoted by white rectangles in Figure 6a. As a reference, the shifts that are derived from the simulated HAADF STEM image of 4C pyrrhotite (see Figs. 7e and 7h) are shown as asterisks in Figures 8a and 8b. The expected Fe atomic column shifts are present in the 4C pyrrhotite matrix. However, they are close to zero in the modified area. It should also be noted that in Figure 8a large deviations in Δx are present (see, e.g., layer 8), suggesting that the investigated area is still affected by the prolonged electron beam illumination nearby. Δx is more sensitive to electron beam damage than Δz (Fig. 8b). Such measurements could therefore be used to assess the influence of electron beam damage on vacancy counting at the atomic scale.

The Δx shift of 4C pyrrhotite in Figure 7e (as well as in Fig. 8a), which shows a shift sequence with alternating N-shaped (red) and inverse-N-shaped (green) patterns, is now considered. In Figure 7e, each 4C pyrrhotite unit cell is associated with an N-shaped shift sequence followed by an inverse-N-shaped shift sequence along the cm axis, corresponding to a periodicity of 8d (where d ≈ 0.29 nm; see Figs. 1b and 7b), which is identical to that of the image intensity modulation indicated by red and green rectangles in Figure 7a. In the relative shift measurements, the central Fe column is located in every second vertical line of Fe columns (yellow circles in Fig. 7b). In the shift sequence shown in Figure 7e, the lower and upper halves of the 4C pyrrhotite unit cell (i.e., the red and green frames in Fig. 7a) are linked to each other by an anti-phase operation, with a displacement between them of R ~ [130]m/8 = 0.297 nm. There is an additional displacement of ~[110]m/8 = 0.172 nm in the viewing direction to match the S sublattice.

The fully occupied Fe layers corresponding to column 4 (marked by red arrows in Fig. 7a; see also Figs. 1b and 2c) have zero shift, as shown in Figure 7e. These fully occupied layers are marked using the symbol F (slim) and those with maximum |Δx| shifts using the symbol F (bold) in Figures 7a and 7e (and Fig. 1c). They are used as references for the determination of incommensurate structures below.

Planar defects in 4C pyrrhotite

Planar defects were also found in the 4C pyrrhotite sample. The HAADF STEM images shown in Figure 9a (along the [010]m direction in the defect-free matrix, i.e., type I) and Figure 9c (along the [110]m-type direction, i.e., type II) illustrate two examples of such faults in the same crystal. Figures 9b and 9d show, respectively, color-coded versions of the same images after noise subtraction (Du 2015). The faulted areas were identified by tracing the stacking sequence of the positions of individual Fe(V) atomic columns (orange circles) along the cm direction. The centers of the faults are indicated by open arrows. According to the stacking sequence (see also Fig. 2 as a reference), the faulted area in Figures 9a and 9b can be regarded as a half-unit-cell intergrowth (yellow rectangle) of the [110]m-type variant (type II; see the full unit cell marked by white rectangles in Fig. 9d) in an [010]m matrix (type I matrix; see the full unit cell marked by white rectangles in Fig. 9b), whereas in Figures 9c and 9d the intergrowth and matrix structures are reversed.

From a crystallographic perspective, the faults shown in Figures 9a and 9b and Figures 9c and 9d are intrinsically the same. This is also apparent from the measurement of the Δx displacement as a function of layer stacking sequence shown in Figure 9e, in which the half-unit-cell twin variant (i.e., the defective region) is marked in yellow. In the type I matrix, the twin variant shows an inverse-N-shaped shift sequence (green), as expected for the type II structure, whereas in the type II matrix the defective band shows a shift sequence that is identical to that in the type I structure (blue). This analysis demonstrates that quantitative measurements of Fe atomic column displacements can be used to provide additional insight into the local structure of pyrrhotite.

Atomic imaging of incommensurate structures

Figure 10a shows a representative raw HAADF STEM image taken from an undamaged region that includes both commensurate 4C and incommensurate 4.91C pyrrhotite superstructures, between which the boundary is marked by a white dotted line. Figure 10b shows the same image after moderate noise subtraction (Du 2015). The difference between the raw and denoised images is shown in Figure 10c. The denoised image was further processed using Gatan DigitalMicrograph software by color coding and γ tuning (value = 0.32), to better identify the Fe columns visually. Based on the atomic column intensities in the processed image, all possible Fe-rich columns could be marked. Unit-cell-by-unit cell shift measurements were also performed. The interface between 4.91C and 4C was found to be atomically sharp, and the 4.91C pyrrhotite structure was found to be almost periodic along the [130]m direction of 4C pyrrhotite (i.e., along the global x direction). The observed incommensurate modulation vector q is oriented only in the 4C cm* (|| cm) direction (see Fig. 5a; i.e., along the global z direction). Therefore, it is reasonable to average the image along the global x direction, to obtain more statistical information, i.e., to determine the average of each atomic column intensity in every unit cell. The averaged image is shown in the upper panels of Figures 10d and 10e (note the 90° clockwise rotation), to provide a direct correspondence with the analysis shown in the lower panels.

The lower panel in Figure 10d shows the Δx shift plotted as a function of Fe layer stacking in the cm direction of 4C pyrrhotite. In Figure 10e, the lower panel shows a simplified stacking sequence of the Fe layers determined on the basis of the image intensity. (For clarity, the VA and VB layers are indicated only using orange dashes). The following features are apparent:

  1. In Figure 10d, the shift sequence in the 4C region (highlighted by a blue shadow) is characterized by an alternating stacking of red N-shaped and green inverse-N-shaped patterns with a periodicity of 8d, where d ≈ 0.29 nm, as indicated in Figure 1b and discussed in “Atomic-resolution imaging of 4C pyrrhotite and Fe column shifts.”

  2. In the incommensurate 4.91C pyrrhotite region, the 8d shift sequence is no longer present. Instead, N-shaped (or inverse-N-shaped) patterns of width 4d are separated by 1d-wide intervals, resulting in the formation of a 5d displacive sequence, as indicated by green and red vertical lines. In combination with the apparent 10d image intensity modulation (see, e.g., the areas marked by green and red rectangles in Fig. 10e), the 5d displacive sequence is associated with the presence of nanoscale regions of 5C periodicity (highlighted by green and red shadows in Fig. 10d). These nanoscale regions have lengths of up to 10 nm in the z direction.

At this stage, the boundary regions are considered. Three boundaries, which are present in Figure 10d (and Fig. 10e), can be identified on the basis of both the measured Δx shifts (Fig. 10d) and the atomic column intensities (Fig. 10e). First, the shift pattern (Fig. 10d) is discussed. When the N-shaped (red) and inverse-N-shaped (green) nano-regions meet, the structure corresponds locally to two 4C-like shift sequences that are overlapped with each other (green and red dashed lines), leading to a new M-shaped or W-shaped shift pattern (black lines). The width of the boundary, which is characterized merely by the Δx shift measurements, is 5d (yellow shadows in Fig. 10d). The boundaries are also characterized by the presence of two fully occupied Fe layers, which are marked FF in Figure 10e.

The stacking sequence that includes the FF layers is less localized, effectively widening the boundary regions, as shown in Figure 10e. In the 4.91C region, two types of stacking sequences can be observed: …FVAFVAVAFVBFVBVB… (green and red shadows) and …FFVAFVAFFVBFVB… (yellow shadows). By comparing Figures 10d and 10e, the former sequence is found to occur only in parts of the green and red nano-regions in Figure 10d, leading to much wider boundary regions, which are characterized by the latter sequence, in Figure 10e. In the following discussion, regions that include both the 5d shift sequence (Fig. 10d) and the …FVAFVAVAFVBFVBVB… stacking sequence (Fig. 10e) are referred to as 5C-like nano-regions. The presence of two fully occupied (or at least Fe-rich; see below) FF layers suggests that in the boundary region the Fe concentration is higher than in the 5C-like nano-regions, thereby resulting in variations in composition within 4.91C pyrrhotite.

Both the shift measurements and the intensity analysis are consistent with the formation of a 1d (or 4d) anti-phase-boundary-like junction in the z direction between two nano-regions, as illustrated by shifted green and red wavy lines in Figures 10d and 10e. It should be noted that: (1) The spatial overlap of two 5C-like nano-regions cannot reproduce the observed shift and intensity features simultaneously. (2) The structural boundary features are consistent between all of the observed boundary regions, including pure 4.91C areas (i.e., where no 4.91C/4C interface is present), suggesting that this boundary structure is intrinsic to 4.91C pyrrhotite. (3) Although the boundary areas behave like anti-phase boundaries consisting of 4C cells connected by FF layers (i.e., …FFVAFVAFFVBFVB…), similar to the TIM model reported by Harries et al. (2011), the Δx shift sequence in such areas suggests that the boundary structures are more complex than that produced by a simple 4C anti-phase operation.

A structural model for 4.91C pyrrhotite

Two representative segments of 4.91C pyrrhotite, comprising an individual 5C-like nano-region and a boundary region (red and yellow rectangles in Fig. 10e), are enlarged in Figures 11b and 11c, respectively. The 4C matrix (blue rectangle in Fig. 10e) is also shown in Figure 11a as reference for the following analysis. The cm axis is aligned vertically for visual reasons.

The following observations are based on an analysis of the Fe column intensities in the 5C-like and boundary regions:

  1. The Fe(V) columns, which are marked by orange pentagons in Figures 11b and 11c, have higher intensities than those in 4C pyrrhotite (Fig. 11a), suggesting that the Fe vacancy sites that are present in 4C pyrrhotite are partially occupied in 4.91C pyrrhotite.

  2. In the FF layers (Fig. 11c), an intensity reduction is observed in the columns that are marked by red dashed circles in the 4.91C structure, when compared to the 4C structure, suggesting that the fully occupied Fe columns in 4C pyrrhotite are no longer fully occupied in 4.91C pyrrhotite, but contain vacancies.

To understand the modified Fe occupancies/vacancies in the 4.91C superstructure, HAADF STEM image simulations were performed. Starting from the 4C structure (Tokonami et al. 1972), the occupancies of the Fe vacancy sites in some of the Fe(V) columns were changed, as illustrated in the form of partially filled orange circles in Figure 11d. The resulting modified columns are referred to as Fe(M). Five occupancies of the vacancy sites, i.e., 0.1, 0.3, 0.5, 0.7, and 0.9, were used in the simulations, leading to Fe concentrations of 55, 65, 75, 85, and 95% in the whole Fe(M) columns.

The situation with no additional occupancy, i.e., pure 4C pyrrhotite, is first considered. A simulated image for a sample thickness of 33 ± 1.4 nm provides the best match to the atomic column intensities in the experimental results obtained from the 4C structure, as shown in the form of cyan and blue bars in Figure 11d. All of the intensities (including Fe(M); right cyan bars in Fig. 11d) are normalized with respect to those of fully occupied S columns. The intensities of the Fe(1–3) and Fe(4) columns are shown separately because the latter intensities are reduced, as they are associated with the largest Fe column splittings in projection (46.7 pm; see Fig. 1).

The Fe occupancies in 4.91C pyrrhotite were estimated based on the fit between the experimental image and the simulation for the 4C pyrrhotite matrix. The observed enhancement in Fe(V) column intensity (orange pentagons) in the 5C-like (red) and boundary (yellow) areas was determined to correspond to an Fe occupancy of approximately 0.3 at the partially occupied Fe sites (resulting in ~65% average column occupancy). Similarly, for the Fe(1–3) and Fe(4) columns that are fully occupied in 4C pyrrhotite, the reduction in column intensity is determined to correspond to an Fe occupancy of ~0.7 at every second Fe site (resulting in a column average of 85%).

It should be noted that the simulations assume that the atomic sites that contain partial Fe occupancies are ordered in every second site along the column. The presence of an ordered arrangement of partially occupied Fe atoms in the electron beam direction is very difficult to verify using electron microscopy. However, it has been reported to be energetically favorable in natural 5C pyrrhotite (Elliot 2010; Liles and de Villiers 2012).

Based on the measured atomic column shifts and statistical intensity measurements, a rigid supercell model containing a total of 3712 atomic sites was constructed. (Details of the supercell construction are given in Supplemental1 Fig. S5.) The supercell model shown in Figure 12 includes green and red 5C-like nano-regions, as well as their boundary (magenta dashed rectangle in Fig. 10d). The occupancy values determined for the partially occupied Fe atomic sites are also given (with a total of 1640 Fe atoms in the model).

A comparison of a simulated HAADF STEM image calculated on the basis of this supercell model for a specimen thickness of 33 nm (inset superimposed on Fig. 10b) with the experimental image shows that the features are reproduced well, including the presence of uneven spacings between spots in simulated electron diffraction patterns (Figs. 5b and 5c), suggesting the validity of the model.

Based on this model, the Fe concentration can be calculated to be 1640/(1640+1856) ≈ 46.9%, which is slightly higher than that for ideal 4C pyrrhotite (≈46.7%). Assuming that the Fe atoms in this model have the same spin configurations as in all Fe1–δS omission structures (i.e., ferromagnetic coupling between spins within Fe layers vs. antiferromagnetic coupling between adjacent Fe layers), the net magnetic moment will be proportional to the number of uncompensated Fe atoms between adjacent layers. The difference is <2 Fe atoms per supercell (2.377 × 1.374 × 16.63 nm in size), suggesting that the incommensurate 4.91C pyrrhotite is macroscopically antiferromagnetic.

In the present supercell, the occupancy values of the partially occupied Fe sites were not modulated incommensurately. Nevertheless, almost the same diffraction pattern could be obtained when removing the local chemical variations from the supercell model (i.e., with the Fe occupancies set to 1). The contribution from structural (i.e., displacive) anti-phase-like stacking between unequally sized nano-regions is therefore thought to be the dominant origin of incommensurability in 4.91C pyrrhotite.

By comparing the 5C-like segments in the present supercell model with those for natural 5C pyrrhotite proposed by Elliot (2010) and Liles and de Villiers (2012), one may conclude that 5C-like nano-regions in natural 4.91C pyrrhotite differ from 5C Fe9S10 described by Elliot (2010), whereas they have structural and chemical similarity to 5C Fe9.021S10 described by Liles and de Villiers (2012).

Perspectives for Fe vacancy counting

By using the sample preparation method described above, a specimen thickness of between 30 and 50 nm could be achieved without noticeable ion beam damage, thereby allowing superstructures in both commensurate and incommensurate pyrrhotites to be explored in detail. However, according to the image simulations shown in Figure 2, such specimen thicknesses are still far from being optimal (<10 nm at least) for Fe vacancy counting. Further improvements in sample preparation could include the use of focused He or Ne ion beam systems (Winston et al. 2011; Pekin et al. 2016), together with sample cooling during subsequent ion milling. In addition, Mo grids could be used to eliminate the re-deposition of Cu from the grids that are currently used.

A further issue that has to be considered is electron beam irradiation damage (within several seconds at 200 kV), which, according to Figure 6, can alter the Fe vacancy distribution. Similar observations have been reported by Nakazawa et al. (1975) and Pósfai et al. (2000) due to ion milling and electron irradiation and are expected to be stronger when the sample is sufficiently thin for Fe vacancy counting. Further investigations of the influence of electron dose and/or dose rate on irradiation damage at different accelerating voltages are required. The use of a newly developed low-dose imaging method (Yu et al. 2016; Zhu et al. 2017) based on direct electron detection would also be beneficial.

Last but not least is the issue of the precision of Fe vacancy counting, which, according to the image simulations shown in Figure 2, is affected by electron dechannelling caused by atomic column splitting, in particular along the [110]m-type axis. If the sample geometry is suitable, then the use of [010]m variants could be a possible solution. Another approach could involve support from first-principles calculations, to predict values of atomic column splitting in thin samples.


Significantly, in addition to determining the local Fe vacancy concentration, this paper provides an indirect method for inferring the local magnetic properties of pyrrhotite at close to the atomic scale. In agreement with the model for the 4C structure (Tokonami et al. 1972), the results are consistent with the fact that 4C pyrrhotite is ferrimagnetic. The inferred stacking sequence and changes to the occupancies of atomic columns in the 4.91C super-structure indicate that it is antiferromagnetic, both in nanoscale 5C-like regions and at the boundaries. This indirect information is important for studying pyrrhotites and other minerals because neither electron magnetic circular dichroism (Schattschneider et al. 2006) nor phase contrast techniques such as off-axis electron holography (Lichte and Lehmann 2008) or 4D-STEM (Ophus 2019) are currently able to measure magnetic properties at this spatial resolution, in part because the acquisition time to record such a weak magnetic signal would be so long that the material would be affected by electron irradiation damage.

Since quantitative atomic-resolution transmission electron microscopy is a well-established technique that has been applied to a wide variety of materials, such as perovskite oxides, it will be possible to apply the present approach to a vast group of incommensurate pyrrhotites and other minerals. For example, it may be possible to use atomic-resolution transmission electron microscopy combined with quantitative analysis to determine whether incommensurate 5.5C and commensurate 11C superstructures are structurally identical.

Acknowledgments and Funding

The authors are grateful to Doris Meertens for FIB preparation of specimens for HAADF STEM, Dongsheng Song and Xiankui Wei for valuable discussions, and Juri Barthel, Lothar Houben, and Hongchu Du for providing software programs or scripts. Part of this project was supported by the Swiss National Science Foundation (Grant No. 200021-153173). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreements 856538 (3D MAGiC) and 823717 (ESTEEM3), as well as from the DARPA TEE program through grant MIPR no. HR0011831554 and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project ID 405553726-TRR 270.


Deposit item AM-21-17479, Online Material. Deposit items are free to all readers and found on the MSA website, via the specific issue's Table of Contents (go to http://www.minsocam.org/MSA/AmMin/TOC/2021/Jan2021_data/Jan2021_data.html).
In projection, only two of the four types of positions can be distinguished (Elliot 2010).
This is an open-access article distributed under the terms of the Creative Commons Attribution CC-BY 4.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.Open access: Article available to all readers online. This article is CC-BY.