Surfaces , Interfaces and Crystal Growth of Marcasite FeS

Marcasite mineral is a metastable iron sulphide, α-FeS2, less known and less abundant than pyrite β-FeS2. Their chemical compositions are similar, formed by S2 2dimers and Fe2+ species, but the marcasite crystal structure is orthorhombic while pyrite is cubic. The thermodynamic stability of a range of selected surfaces and their kinetic growth are analysed with molecular mechanics simulations implemented in METADISE code. The most stable surface of marcasite corresponds to the (101) surface having the lowest surface energy (marcasite 101 γ = 1.06 Jm-2) and the most favoured kinetic growth rate is expected for the [010] direction (the attachment energy of (010) surface is 0.20 |eV|/at). Mirror twinning is another characteristic of the mineral contributing to a distinctive crystal growth forming different shapes such as striated needles or pseudo-pentagonal flowers. Similarities, differences and overgrowth between the two dimorphs are highlighted.


Introduction
Marcasite-FeS 2 , called formerly 'white iron pyrite', is a metastable phase with an orthorhombic crystal structure. This sulphide grows under acidic conditions and low temperature [1][2][3][4][5][6] . The mineral has a chemical composition similar to pyrite (known as 'fool's gold'), more pale-yellow silvery colours, and its predominant morphologies are well distinguishable. It is usually found in layers of twinned surfaces forming a flower with a pseudo-pentagonal arrangement or acicular while pyrite adopts preferentially euhedral and platonic shapes 7 . Some controversies are concerning crystals with radiating nodules, reported as pyrite or marcasite 8 and mainly X-ray diffraction (XRD) should be able to differentiate them. Marcasite is less known and experimental data are scarce. This mineral is anisotropic 9,10 and decomposes more easily than pyrite by surface attack 11 . Some potential applications have been proposed for lithium and sodium batteries 12,13 .
Most of the properties have been investigated with DFT (Density Functional Theory) methods (such as calculating band gap, elastic properties, thermoelectric properties [14][15][16] ). Marcasite has been frequently disregarded as it was believed that its band gap was low (E g =0.34eV 17 ), incompatible for photovoltaic applications which has been first dismantled by first-principle simulations 14 and then by diffuse reflectance.
The computed band gap is slightly higher than pyrite and it was later confirmed experimentally 18 .
This article aims to bring new insights on surfaces and interfaces using force field methods in order to unpuzzle crystal growth and twinning. A comparison of structural properties between marcasite and pyrite will also enable to better understand overgrowth with preferred orientations.

Methodology
Pyrite-FeS 2 has been one of the first minerals studied by W.H and W.L Bragg (father and son) in 1913 using X-ray diffraction, giving at that time some difficulties to elucidate as it belongs to the imaginary Fedorov system. Fe 2+ /S 2 2ions occupy a face centred cubic sites (space group Pa3 , number 205) but the centre of the cubic cell is not occupied by an atom. There are separated atoms at the centre of the cell, a dimer of sulphurs orientated along the <111> direction. The Wyckoff positions for S are 8c(uuu) and 4a(000) for Fe. Iron atoms are in an octahedral environment. The octahedra are only edge sharing.
The crystal structure of marcasite-FeS 2 (see Figure 1) is orthorhombic with Pnnm space group (space group number 58). The Fe 2+ /S 2 2ions occupy body centered positions. Their Wyckoff positions are 2a(000) and 4g(uv0) using lattice parameters with c<a<b (n.b: some articles might report another lattice orientation 7 ). The dimer of sulphurs is then aligned within the (001) plane. The octahedra are corner sharing and also edge sharing.
The methodology to simulate marcasite is the same than the one described in the previous study by Arrouvel and Eon 19 , using identical interatomic potentials that have been previously validated on pyrite by de Leeuw et al. 20 [2] where r 0 is the separation between the two sulphurs (dimer S1-S2) at equilibrium.
A 3-body S-S-Fe angular interaction (the sulphur atoms do not belong to the same dimer) is represented by the potential ϕ ijk as follows: where k ijk, is the constant force and θ 0 the angle between the three atoms at equilibrium. The force field parameters are listed in Table 1.
The extraction of elastic constants is useful in geology as it provides some data about the composition in Earth's interior from seismic studies. The bulk modulus Voigt average B indicates the structural stability of the material by the resistance under a uniform pressure and it is deduced following the relation: The shear modulus G measures the material's response to the shear stress (formula in Equation 5 The Young's modulus (E) and Poisson's ratio (n) depend on B and G as follows: ( ) All the simulated (hkl) surfaces are generated automatically using the Tasker 23 approach which enables to identify the energetics and the structure of the most stable terminations. The choice of the surfaces is done in prioritising low Miller indices, observed surfaces/interfaces and XRD peaks. The kinetic growth rate is related to the attachment energy, following the theory of Hartman-Perdok 24 . E att is the energy released per mol by the growth of a layer (or slice) of thickness d hkl (interplanar distance) according to the relation: With E bulk the lattice energy of the crystal, E slice the energy of the slice (hkl) per mol. The surface energy in vacuum, hkl γ , for the (hkl) surface follows the expression: With E hkl the energy of the (hkl) surface, E bulk the energy of the bulk normalized to the number N of FeS 2 units in the supercell and hkl A the surface area. The morphology of an ideal crystal at the thermodynamic equilibrium is constructed using the law of Gibbs-Curie-Wulff 25 .
h hkl corresponds to the distance from the centre of mass of the solid to the (hkl) surface and α is a characteristic constant of the material.  Some mirror twin grain boundaries (GB) are constructed. To select the most stable GB, the unrelaxed slab is reflected and scanned in moving its position to the respect of the other slab by steps of 0.2 Å along the 2 dimensions of the (hkl) interface plane. At each point of the grid, the system is relaxed. The lowest interface potential energies E GB give the most stable GB per surface area. The formation energy (E f ) and the cleavage energy (E c ) of GB are calculated using the relations: To stabilise a grain boundary, E f is minimized and E c is maximized.
The present paper is based on force field methods to simulate surfaces and interfaces of marcasite; however, supplemental ab initio simulations have been undertaken on the pyrite (001) surface using DFT (Density Functional Theory) method to support the cogrowth between the two dimorphs. The bulk of pyrite has been optimized with VASP software 26 at the gamma-centred 12×12×12 k-points grid with RPBE+D3 functional and PAW pseudopotential (the cell parameters and ions have been fully relaxed). The cutoff on the kinetic energy is 500 eV. The k-points grid of the surface is 6×6×1, with a slab thickness of 14.8 Å and a vacuum thickness of 13.2 Å (the ions have been fully relaxed, the cell kept fixed). The scanning tunneling microscopy (STM) simulation is done using the Tersoff and Hamann 27 approximation in calculating the energy-resolved charge density, which is proportional to the tunneling current, from -2.5 eV to the Fermi level. The STM image is obtained with HIVE-STM program 28 .

Bulk
The force field simulations have been applied successfully to pyrite 19,20 and prove to be transferable to marcasite having the same chemical nature. The optimized structural parameters, elastic and dielectric properties of the dimorphs are listed in Table 2. Experimental data are rarer on marcasite and no measurement on dielectric constants has been found. The deviation of its simulated cell volume from the experimental values 29 (a=4.436Å, b=5.414Å, c=3.381Å, volume of 81.20 Å 3 ) is 5%. The orthorhombic system has 9 independent values for elastic constants and 3 for dielectric constants which reflects anisotropic properties. Elastic properties from classical simulations are in very good agreement with DFT/PBE+U calculations 12  The mineral is then easiest to distort upon a shear stress rather than compression. Some properties extracted from force field methods on pyrite have been already discussed by Arrouvel and Eon 19 , are extended in Table 2 and compared to experimental data [31][32][33][34] . This dimorph is structurally isotropic, with 3 independent elastic constants. The elastic constants C 11 = 31.99 × 10 11 dyn/cm 2 (i.e. 320 GPa), C 12 = 5.54 × 10 11 dyn/cm 2 (i.e. 55 GPa), and C 44 = 10.80 × 10 11 dyn/cm 2 (i.e. 108 GPa) are also in good agreement with experiments and theoretical results (see Liu et al. 34 and references therein). The calculated bulk and shear modulus are smaller on pyrite compared to marcasite indicating that the compressibility if higher on pyrite.
Stability and elastic properties between the two dimorphs are pressure and temperature dependent; the order of stability and elastic modulus is inverted at high pressure 34 . FeS 2 is not an interesting material for the construction due to fracture and erosion. In jewellery, pyrite is preferred to marcasite because marcasite erodes more quickly, tarnishes and is more fragile.
The calculated dielectric constants in marcasite are slightly lower than those in pyrite, in other words, pyrite should be slightly more polarizable. However, another investigation using ab initio simulations (not considered in the present study) should be more reliable to compare those properties and to improve force field reliability if necessary.

Calculated surfaces and morphologies
The marcasite (hkl) facets have been simulated using METADISE to evaluate the structure and stability of each surface. Attachment energies and surface energies are listed in Table 3. Most of the surfaces are Tasker type II with a sulphur termination. Results are compared to available experimental and theoretical data however some discrepancies are expected in particular for cases with strong relaxations and reconstructions. Indeed, reductions and oxidations of the surface species are not considered by the force field method. A surface Sspecies, especially exposed from a cut sulphur dimer, is expected to be an unstable species, to be strongly stabilised by a reduction into S 2-. The reduction of the sulphur species should also be linked to an oxidation of Fe 2+ into Fe 3+ in considering stoichiometric slabs. Another possibility would be to create sulphur vacancies, with S 2surface species, to keep the oxidation state of Fe 2+ . Such non-stoichiometric termination becomes similar to a FeS termination. Possible reconstructions stabilising surface energies with steps/kinked exposing other facets are not either considered in the present study. In marcasite, the dimer crossing the a axis is approximately oriented along the [110] direction and the other set of dimers is oriented along the [ ] 110 direction. The highest discrepancy on surface energies using DFT and force field methods corresponds to the (110) surface which is the plane perpendicular to S-S dimers. This is supporting the fact that force field methods underestimate the relaxation from Sspecies.
Strong stabilisations have also been observed on pyrite. Similarly, the surface energy obtained with force field for the (111) surface is higher compared to the surface energy obtained with DFT. The discrepancy corresponds to the family of planes perpendicular to S-S dimers which are orientated along the <111> directions. The {110} facet can also be stabilised through steps made of stable {001} facets (see Arrouvel and Eon 19 and references therein). However, the main exposed surfaces, thermodynamically stable and with higher growth rates, are in excellent agreement with ab initio and experimental studies. Force field simulations indicate that the most favourable kinetic growth on marcasite is along the [010] direction followed by the [101] direction. This result explains that acicular minerals are likely to be elongated along the kinetic direction.
The energy required for (hkl) growth rate is in the following order: (010) < (101) < (001) < (100) < (102) < (111) < (103) < (011) < (110) ~ (121) < (021) < (211) < (201) < (031) < (120) < (210) < (130) < (305) The most favoured kinetic growths are for lower Miller indices. The second surface with higher rate is the (101) surface, which is also the most stable surface, the most common twinned surface and it has an epitaxial relationship with pyrite (001). This surface will be then discussed in more detail. Richards et al. 36 Table 3 for comparative purpose. The order of surface energies after relaxation, calculated with force field, is as follows: < (010) < (001) < (111) < (121) < (102) < (103) < (031) < (305) < (100) < (130) < (211) < (011) < (201) < (120) < (021) < (210) < (110) As expected, we obtain the same results from Ngoepe et al. who have calculated 4 surfaces with force field 37 . The order on the surface energies is also consistent with DFT + U simulations 15 . All theoretical studies agree with the fact that (101) is the most stable surface followed by (010) surface. At the thermodynamic equilibrium under vacuum, only those two surfaces are exposed in the constructed morphology ( Figure 2). The most stable surface, the (101), has structural similarity with the pyrite (001) surface and energetically equivalent with surface energies of 1.06 Jm -2 and 1.04 Jm -2 19 respectively. Five-fold iron ions are exposed with the base of each polyhedron linked by sulphur dimers. STM is a powerful tool that enables to evidence some chemical resemblance between the two dimorphs. An STM image has been simulated with DFT on marcasite (101) 15 and due to lack of STM data on pyrite (001) surface, the latter has been simulated applying the same bias voltage on the S-terminated surface for comparison (n.b. DFT simulations are not the focus of the present paper. A work comparing the influence of the functionals on structural and physical properties on MS 2 pyrite structures, M=Transition Metal, is under investigation). In Figure 3 is reported the simulated STM image with a tip at 2.1 Å from the surface. The white spots forming columns in zig-zag are the most external sulphurs, upper to the iron plane. The same topology is noticeable on the STM image of marcasite (101) simulated by Dzade and de Leeuw 15 . The distances between sulphurs are also matching with those measured on marcasite (the distance between 4 columns is 16.2 Å, e.g. 3a of pyrite, similar to the distance reported by the authors 15 ). The epitaxial overgrowth of marcasite {101} and pyrite {001} has already been observed on minerals 38,39 and calculations emphasize that the overgrowth is likely to be unidirectional. The good matching, called 'leaf epitaxy', is illustrated in Richards et al. 36 in which the pyrite, on the top of a rib, has the (h00) facet orientated parallel to the (101) facet.
The (010) surface is with higher symmetry. The exposed   classical methods are for (110), (011) and (130) surfaces which appear to be more stable with DFT 15,35 . The {011} facet can play a role in marcasite crystal growth as it is reported as a possible twin. Twinning on the {101} facet is the most common phenomenon in this mineral. Mirror twins of those two facets are then simulated and their interfaces are discussed.

Interfaces
The {101} facets of marcasite have a lot of similarity with the {001} facets of pyrite. Rotational twinning is well known in pyrite, giving the iron-cross shape from pyritohedrons for example (pictures and data of minerals are accessible in https://www.mindat.org/article.php/3872/ PYRITE+-+Iron+Cross+Twins, and a classification of twins is displayed by Tanako 40 ). It is wildly admitted from macroscopic observations that iron-cross twins are formed by a rotation of 90º around one of the [001] axis. Other authors 41,42 proposed another interpretation involving mirror twinning of {110} facets observable by high resolution transmission electron microscopy. A forthcoming theoretical work with classical mechanics and ab initio methods will detail rotation and mirror glide twinning in pyrite. The apparent C 4 rotation axis is not existing in marcasite and only mirror twinning is considered in the present investigation. No experimental data is available at the atomic scale on marcasite. The twinning of the {101} plane is known as Sperkise twinning 8 . In marcasite, alternated bands of polarization colours come from twinning of anisotropic structures 9 . From simulations, the energy profile of the stable mirror S-terminated (101) slab scanned on the top of the original slab gives a global minimum and some possible local minima. The lower formation energy and the higher cleavage energy would be the most favoured grain boundary. Low Miller (hkl) indices, being with higher symmetry, have lesser distortions at the interfaces and than less possibilities of displacement. The lattice parameters parallel to the (101) interface are 5.55 Å (corresponding the b axis of the bulk) and 5.63 Å. The displacement from the energetic minimum is at 0.2 Å and 3.6 Å. The slabs are relaxed and the optimized interface is in Figure 4. Each sulphur bonds to the opposite Fe; the local environment geometry at the interface is similar to the bulk. The Fe-S distance stretches by 0.04 Å (dist(Fe-S) = 2.33 Å). The formation energy E f is 0.13 Jm -2 and the cleavage energy E c is 2.0 Jm -2 . The good continuity between twins has been mentioned by Nespolo and Souvignier 43 using crystallographic orbit analysis. The (305) plane, being quasi-perpendicular to the (101) surface, has been simulated and it appears to be relatively stable, at 2.54 Jm -2 . This twinning is the most common in marcasite exposing the flat stable (010) surface. The angle between (101) and ( ) 101 planes is 74.6º (from experimental lattice parameters), slightly larger than π/5 (72º). The mineral, made of successive {101} twins, is enrolling in 3D forming the pseudo-pentagonal flower-like with possible other facets involved in the crystal growth. The mirror (011) twin has also been calculated as it can be a possible rarer twinning. The energy formation is 0.07 Jm -2 and the cleavage energy of 5.89 Jm -2 which places this surface as a good candidate for twinning. The (101) and the (011) surfaces are then two favourable surfaces for twinning, with higher symmetry and keeping the bulk characteristic with dimers of sulphur without formation of vacancies, then without changes of the oxidation state of the species. Force field methods are well adapted to simulate such twins. For other types of twins involving defective interfaces with possible chemical changes, ab initio methods will be required.
Contact boundaries are frequent in this mineral and might be responsible for being more fragile than pyrite and breaking apart. Different morphologies are mentioned in mineral data bases (such as in mindat.org), acicular, prismatic, cockscomb, platy, however none of them proposes a clear identification of the exposed surfaces and striations. Striations have been characterised onto pyrite and SEM (Scanning Electron Microscopy) images were key to certify their directions. A forthcoming article combining an experimental and theoretical study on pyrite will detail the role of the sulphur network on the directional kinetic growth of the striations. The identification of striations becomes also useful to point out twins. Parr and Chang 10 did notify various shapes of marcasite and their difficulties to determine the surfaces are due to no defined angle and the coexistence of (010) and (hk0) surfaces. Based on some features commonly reported, Figure 4 illustrates a mirror {101} twin and striations typically observed on crystals extracted from mines in Belgium (Limites quarry, Ave-et-Auffe, Rochefort, Namur, Wallonia). In this graphical example, the twinned (101) and its pair ( ) 101 have a higher surface area than the pair ( )

101
and ( ) 101 . The striations are along the [001] direction visible on (010)/(hk0) surfaces, which is also a favourable kinetic growth direction. The same orientation of the striations is seen in Richards et al. 36 . Since the lack of reliable mineral data, other cases are not excluded. For example, it has not been possible to verify the existence of striations on other surfaces.

Conclusion
Marcasite FeS 2 is a mineral that can adopt various shapes, depending on the geological conditions. Its anisotropy, twinning and oxidation rate enable to distinguish it to pyrite. Radiating nodules are however more ambiguous and XRD is the method to characterize the phase unquestionably. Striations are also common on the two minerals, observable at naked eye onto specific facets. While their directions are evident in pyrite (along <001> directions), no experimental data had clarified the crystal growth and its striations on marcasite faces. The present study has identified structural and energetic characteristics of a wide range of marcasite surfaces using classical mechanics methods. The (101), (010) and (001) surfaces are playing an important role in the crystal growth as they are stable and with higher growth rates. The (101) surface has the lowest surface energy; it is structurally and energetically similar to the (001) surface of pyrite. The sulphur terminations, highlighted by STM simulations, are aligned on both surfaces along a specific direction which explains the orientated overgrowth of the dimorphs. Mirror twinning is also favourable at the marcasite {101} interface. The distortion at the boundary is minimal in keeping the specificity of the disulphide (i.e. dimers of sulphurs). The {011} facet is as well likely to twin, having favourable formation and cleavage energies. A typical twinned specimen from Belgium has been schematized and the striations are estimated to be aligned along the [001] direction. A forthcoming DFT study will enable to characterize electronic properties of strongly reconstructed surfaces and to consider other types of twinning (e.g. rotational twinning).

Acknowledgements
The work has been done without any financial support. The author thanks CENAPAD-SP for computational resources. The author is grateful to the two anonymous reviewers for their constructive comments.