Atomic scale chemical ordering in franckeite—a natural van der Waals superlattice

The mineral franckeite is a naturally occurring van der Waals superlattice which has recently attracted attention for future applications in optoelectronics, biosensors and beyond. Furthermore, its stacking of incommensurately modulated 2D layers, the pseudo tetragonal Q-layer and the pseudo hexagonal H-layer, is an experimentally accessible prototype for the development of synthetic van der Waals materials and of advanced characterization methods to reveal new insights in their structure and chemistry at the atomic scale that is crucial for deep understanding of its properties. While some experimental studies have been undertaken in the past, much is still unknown on the correlation between local atomic structure and chemical composition within the layers. Here we present an investigation of the atomic structure of franckeite using state-of-the-art high-angle annular dark-field (HAADF) scanning transmission electron microscopy (STEM) and atom probe tomography (APT). With atomic-number image contrast in HAADF STEM direct information about both the geometric structure and its chemistry is provided. By imaging samples under different zone axes within the van der Waals plane, we propose refinements to the structure of the Q-layer and H-layer, including several chemical ordering effects that are expected to impact electronic structure calculations. Additionally, we observe and characterize stacking faults which are possible sources of differences between experimentally determined properties and calculations. Furthermore, we demonstrate advantages and discuss current limitations and perspectives of combining TEM and APT for the atomic scale characterization of incommensurately modulated von der Waals materials.


Introduction
Two-dimensional (2D) materials fascinate with their unique and exciting properties originating from quantum effects. 2D materials such as graphene, phosphorene and hexagonalboron-nitride have been successfully introduced in device applications [1][2][3][4][5]. However, it remains desirable to find means to tailor material properties.
Designing vertical stacks of 2D materials has been shown as a possibility to tune material properties [6][7][8][9][10]. The stacked layers are mainly connected by weak van der Waals forces and the class of materials is therefore called van der Waals heterostructures. In particular, the periodic layering (vertical stacking) of crystal lattices connected by a van der Waals interlayer is described as a class of van der Waals superlattice [11].
Besides the progress in designing new materials, nature already formed van der Waals superlattices long before. Those naturally occurring van der Waals superlattices provide templates for growths and synthesis in the laboratory [12], inspire first-principle studies on new van der Waals heterostructure combinations to predict electronic structures and properties of newly designed materials [13,14] or can also be applied directly in applications according to their extraordinary properties [15,16].
One of those naturally occurring van der Waals superlattices is franckeite, a sulfosalt mineral composed of periodically stacked layers. Franckeite has gained substantial attention because of its promising properties and the relative ease with which it can be mechanically exfoliated. Density functional theory suggests that the HOMO and LUMO bands in franckeite are confined to alternating layers [17,18], analogous to type-II staggered gap alignment in multiple quantum well structures. Charge-transfer and energy-transfer in van der Waals heterostructures with carrier confinement due to staggered gap alignment is itself the subject of intense interest [19][20][21][22]. In contrast with drawbacks of designing van der Waals heterostructures artificially (limited control of relative crystal orientation of layers, risk of damage, air or impurities between stacked layers), the alignment of layers in the superlattice of franckeite has shown to be preserved across nanoflakes and appears to be free from impurities at the surface/interfaces of the stacked layers [17,18]. Its stability in air together with the reported possibility to exfoliate fewlayer flakes [17,18] raises its intriguing potential in versatile energy conversion and energy storage applications [18], optoelectronic applications [16,17,[23][24][25][26] as well as in the field of chemical and biosensing applications [27][28][29]. Moreover, the franckeite constituent PbS is a material with electronic structure of interest for thermo-electric materials [30]. Of especial relevance to the work reported here, franckeite is imbued with non-trivial structural features, including incommensurate constituent layer lattices and a spontaneous, anisotropic, out-ofplane, lattice modulation [31]. Our observations of franckeite reveal the complexity in structure and chemical order that can be present within a van der Waals superlattice, and which may require refinement of electronic structure calculations.
A description of the basic structures of the franckeite layers is given in the x-ray diffraction (XRD) refinement by Makovicky et al [32]. The pseudo tetragonal Q-layer has been described as being four atomic planes thick (two tightly bonded bilayers) MS layer with M = Pb 2+ , Sn 2+ or Sb 3+ 6 The stacking direction is perpendicular to the ( a × b) plane. Neighbouring layers are connected by weak van der Waals bonds while atoms have strong chemical bonds within the layers. This stacking implies that the direction of the van der Waals gap is defined to be parallel to the stacking direction [ a × b] and the van der Waals plane is parallel to the ( a × b) plane.
Furthermore, Makovicky et al [32] determined the modulation vectors of the layers to q Q = −0.001 29 (8) The displacement is dominated by a first order transversal sinusoidal (mainly parallel c, compare figure 1, or figures 3-7 within [32]). The occupational modulation, which occurs on the Pb/Sb sites in the Q-layer and on the Fe/Sn sites in the Hlayer, is more complicated due to important contributions up to the third order (compare insert in figure 1(a), or figures 8 and 9 within [32]). For a confirmation of these XRD results, local measurements at atomic resolution are necessary.
A method that has been successfully applied for the study of modulated structures is transmission electron microscopy (TEM) [33][34][35]. Williams et al [33] studied both franckeite and cylindrite (variation of franckeite with only two atomic planes thick Q-layers) with very fine details using selected area electron diffraction and high resolution transmission electron microscopy (HRTEM). They were able to support the structure determination from previous XRD measurements by local analysis, for example the difference in Q-layer thickness between cylindrite and franckeite as well as the sinusoidal displacive modulation. Furthermore, Q-layer stacking faults were identified showing displacive modulations in images from [010] zone axis while this is only expected in images from [100] zone axis (Q = A, H = B, C = Q-layer stacking fault: '. . . ABABCBABA. . . '). These faults were characterized as 90 • rotations about the stacking direction.
While insightful, the imaging techniques used previously have not allowed to draw any direct conclusion from the image contrast on the composition or structure of the samples. The contrast in HRTEM images depends on the relative phase of diffracted beams and is thus sensitive to sample thickness, objective lens defocus and atomic structure of the sample [36]. Contrast interpretation is only possible through comparison with image simulations, such as those presented in the cited studies [33,34]. An imaging mode that is more directly interpretable is high-angle annular dark-field (HAADF) scanning transmission electron microscopy (STEM) [37]. The signal collected on an annular dark-field detector with large inner angle (high-angle scattering) depends strongly on the atomic number (Z to the power of 1.6-1.9) and linearly on the thickness of the locally scanned region [38,39].
Wang et al [35] used HAADF-STEM imaging for their studies on franckeite. They presented images along the zone axes within the van der Waals plane. In their work, they confirmed the microstructure of the two stacked layers from the direct interpretation of the HAADF-STEM image contrast. Furthermore, they showed that the sinusoidal modulation might have a large compositional component. However, the spatial resolution enabled by their instruments was not sufficient to distinguish single atomic columns within the layers at the time.
The spatial resolution of STEM has significantly improved since the work was published. Ongoing microscope developments and major technological breakthroughs such as aberration correctors have enabled sub-Å spatial resolution for HAADF-STEM imaging. Therefore, characterization studies of complex structures such as franckeite can be further enhanced. Recent works [17,18] have demonstrated the possibility of more extensive (S)TEM studies providing atomically resolved cross-section imaging of the stacked layers and chemical analysis using energy dispersive spectroscopy (EDS) analysis in STEM.
Another chemical analysis technique, complementary to EDS-STEM, which has gone through a remarkable development in recent years is atom probe tomography (APT). APT can provide atomically resolved analysis in reconstructed 3D volumes [40] and thereby, in principle, investigate displacive or occupational modulation in 3D at the atomic scale. In particular, investigations at the atomic scale of such modulations can become accessible even though they are not perpendicular to a main crystallographic direction that can become a requirement for atomically resolved measurements using (S)TEM, while generally APT reconstructions profit from structural information as input given by (S)TEM. In the structure of franckeite, the displacive modulation vector is approximately perpendicular to a of both layers enabling the possibility for atomically resolved analysis and comparison in both STEM and APT. It should be noted, however, that APT analysis of franckeite and other van der Waals heterostructures is challenging due to a strong anisotropic atomic density and anisotropic electrical properties. Only very limited literature exists about the application of APT on van der Waals heterostructures with similar structure properties to franckeite (e.g. [41,42]).
The demonstration of exfoliation of franckeite flakes with single unit cell thickness opens the path to incorporating such van der Waals heterostructures in future devices [31]. Knowledge of the structure and chemical order of franckeite at the atomic scale is essential to the development of accurate models of electronic structure and physical properties. Therefore, our work focuses on the characterization of franckeite using aberration corrected HAADF-STEM imaging in various low index zone axes within the van der Waals plane and APT analysis. Through this work, we provide new insights into the local atomic structure and chemistry within the layers and demonstrate that the stacking is influenced by the modulation which constrains α while β varies for both layer types from stack to stack.

Materials and methods
The franckeite crystals used in our experiments originate from a single ore specimen from the San José mine in Oruro, Bolivia [43]. Franckeite flakes were prepared by extraction of franckeite crystals from the ore sample. The crystals were mechanically exfoliated and transferred with Nitto Denko thermal release tape onto a degenerately doped Si substrate with a 300 nm oxide layer, pre-patterned with Au alignment marks to facilitate flake indexing (supplemental figure S.2 (https://stacks.iop.org/JPCM/34/055403/mmedia)). Both exfoliated flakes and bulk franckeite crystals were used for the experiments.
TEM and APT sample preparation was carried out in a focused ion beam (FIB) scanning electron microscope Thermo Scientific Helios G4 UXe DualBeam Plasma-FIB equipped with a Xe + source and a Zeiss NVision 40 equipped with a Ga + source. Common crosssection lift-out procedures were applied followed by parallel milling for the preparation of TEM lamella [44] or annular milling for the preparation of needle shaped APT tips [45]. Selection of the prepared viewing direction is facilitated by a characteristic surface morphology (supplemental figures S.1 and S.2).
STEM work was performed on a double spherical aberration (C S ) corrected FEI Titan 80-300 Cubed TEM operated at 200 kV with a semi-convergence angle of 19 mrad and an inner semi-collection angle of the dark-field detector of 64 mrad. To minimize beam damage, the beam dose was chosen as the lowest possible beam current that provided a still acceptable signal-to-noise ratio (screen current <50 pA in vacuum). In order to reduce the noise of the scanning unit, several frames were acquired with a pixel time of 0.1 μs and a rotation of the scanning directions by 90 • between frames followed by post-processing (rigid and non-rigid alignment of the frames and averaging) using Smartalign [46]. Further image analysis was performed using Gatan's GMS 3.0 and Matlab. Structure models were visualized using VESTA [47].
APT work was performed on a Local Electrode Atom Probe (LEAP) 4000X HR [48]. Experiments were run under a vacuum of approximately 3 × 10 −11 Torr and at an approximate temperature of 24 K. Field evaporation was promoted by illuminating the tip with a pulsed UV laser (λ = 355 nm) at a rate of 50 kHz and a pulse energy of 5 pJ. The rate of evaporation was set to a target value of 0.002 ions/pulse. Reconstruction of the APT data was done using CAMECA's Integrated Visualization and Analysis Software (IVAS) 3.8.0.
The reported structure of franckeite is composed of two triclinic systems (basic units of Q-layer and H-layer). Due to the presence of the two layers there can be some ambiguities in the indexation of zone axes (hk0) with h, k = 0, because of different unit lengths and angles between crystallographic axes. For clarity, we use in our work Cartesian coordinates as reference for indices of viewing directions. The crystallographic direction a of both basic units is defined to be parallel to [100]. Since  [33] and Wang et al [34]) with crystallographic data given by Makovicky et al [32].

Results
The HAADF-STEM images in figure    Translation symmetry can be found between stacks along c while intensity patterns repeat within layers along b approximately after one modulation. (b) Q-layer visualized in temperature contrast scale (dark blue low intensity, yellow high intensity). Highest intensity is detected in the convex regions (arrows). Additionally, contrast alteration is present between neighbouring atomic columns in the two inner atomic planes. (c) H-layer in temperature contrast scale as (b) and in original HAADF-STEM contrast. Only slight intensity variation related to the modulation appears. Signal of lighter atoms is measured close to the heavy atomic columns in agreement to the S atoms in the model (see scheme at the right).
( figure 2(b)). This averaging leads to a qualitative blurring of the projected atomic columns in particular when the sample thickness is in the order of, or larger than, the modulation wavelength.

[100] viewing direction
A more thorough analysis of the [100] viewing direction is visualized in figure 3. Figure 3(a) shows an HAADF-STEM image similar to figure 2(a), but obtained at higher magnification and shown in an inverse contrast for better visibility of the projected atomic columns which therefore appear dark. The significant contrast variation in the Q-layer can be used as a fingerprint for the identification of translation symmetries. For example, shifts of a basic unit length along c repeat the contrast pattern (compare schematically drawn blue boxes in figure 3(a)). On the other hand, basic unit shifts along b are not applicable due to the modulation. An approximate repeat of the contrast pattern is estimated to take place after one modulation wave (compare violet boxes in figure 3(a)) which indicates the adjustment of the modulation to the basic units of the structure.
Moreover, the intensity signal in the Q-layer is observed to adjust to the sinusoidal modulation. This effect is highlighted in figure 3(b) for a representative part of a Q-layer drawn in a temperature contrast scale for which blue represents low intensity and yellow high intensity. The intensity of atomic columns in the two outer atomic planes of the Q-layer appears to be highest at convex (expanded, black arrows) regions of the Qlayer. This is suggestive of the fact that there are elements with higher Z-number on average in the expanded regions than in concave (compressed, opposite side of layers at black arrows) regions where the intensity is weaker or that there are more subtle effects in the contrast due to further displacements of the atoms within the depth of the TEM lamella [49]. The transition in intensity between the convex and concave regions is smooth for the two outer atomic planes. In the two inner atomic planes, there is another superimposed contrast oscillation visible between neighbouring atomic columns. This contrast suggests that these atomic columns consist of elements with an alternating higher and lower Z-number on average in addition to the effects of the modulation.  Stacking of layers appears to be not aligned (compare green and orange lines connecting equivalent atomic sites in adjacent stacks). Variations in the projected inner structure of the Q-layers are observed (compare A and C and yellow line inside the Q-layers and less clearly resolved D).

[010] viewing direction
A detailed analysis of the atomic structure in the [010] viewing direction is visualized in figure 4. The Q-layers can be separated in different forms, here defined as A, C and D as follows. A and C differ in the stacking of the two tightly bonded bilayers. Yellow lines in the Q-layer are used to schematically present, for clarity, the direction of the shortest connection between atomic columns in the two inner atomic planes. The alternating direction between A and C blocks can be explained by a rotation of 180 • about the stacking direction or about a. The D block shows a stacking of the two tightly bonded bilayers similar to A. The difference between A and D is an out of optimal zone axis condition for D. This observation is suggestive of a small rotation relative to the projection obtained in A. Overall, the A block is the most representative structure in the field of view where the stacking sequence can be written (B = H-layer), from top to bottom as ABCBABABDBA. C and D can thus be seen as stacking variants which are further investigated by additional zone axes imaging presented in the next section.
The projected c vector of the Q-layer and the H-layer of adjacent stacks can be drawn by the connection of closest equivalent atomic sites. This is schematically visualized in figure 4 by orange lines for the Q-layer and by green lines for the H-layer. From these observations, it becomes clear that neither c Q nor c H is aligned over all stacks in this projection. This implies that β Q and β H can both vary for each stack. In contrast, no variation of β within a single stack is observed. Furthermore, no correlation between the variation of β and the identified Q blocks (A, C, D) is found. Figure 5(a) presents the selection of an H-layer region of figure 4. It is qualitatively visible that there is a specific pattern of two projected columns with high intensity followed by one with low intensity. This pattern is consistent with chemical ordering in the H-layer. The projected atomic column with weaker intensity could contain, on average, elements with lower Z-number than those with higher intensity. The line profile in figure 5(b) confirms this qualitative observation. The minimum intensity of the line scan was subtracted prior normalization on the maximum.
A similar investigation is given in figure 6 for the Q-layer. Figure 6(a) shows a representative Q-layer region of figure 4 and the line profiles for each atomic plane in figures 6(b)-(e) quantify observed intensity variations by relative normalization of the intensities on the maximum of all profiles after subtracting the offset. The intensity of the two outer atomic planes appears uniform (figures 6(b) and (e)) while the intensity alternates between neighbouring projected atomic columns in the two inner atomic planes (figures 6(c) and (d)). This is similar to the qualitative observations in the [100] viewing direction ( figure 3(b)). This pattern suggests a chemical ordering in the two inner atomic planes of the Q-layer such as atomic columns of elements with lower Z-number alternating with atomic columns of elements with higher Z-number on average.
Furthermore, the minima of the linear profile between atomic columns are not identical in the two inner atomic planes (figures 6(c) and (d)). This effect suggests that every two projected atomic columns in the two inner atomic layers of the   Since γ varies minimally from 90 • , such a stacking fault can explain the structural blurring of D observed in figure 4. In contrast, the stacking fault labeled as C in figure 4 is caused by a rotation of 180 • about the stacking direction. Otherwise, no optimal zone axis condition would be present in contradiction to our observation. In the viewing direction of figure 7(a), C cannot be distinguished from A.

Atom probe tomography
In the previous sections, information about local chemistry are given qualitatively without identifying the present elements. Such local chemistry information with identification of the present elements is revealed from APT analysis on franckeite as presented in figure 8. Figure 8(a) shows the reconstructed volume of a tip. Samples were prepared from flakes with a surface parallel to the van der Waals gap plane using standard lift-out procedures, such that the final tip axis would be normal to this plane in order to ensure optimal spatial resolution conditions perpendicular to the layers [50]. Precisely, such an optimal spatial resolution perpendicular to the planes is not achieved over the entire tip volume, but only in a so called pole region [40,48,50]. Therefore, a volume including only the pole region of the H-and Q-layers was reconstructed by selecting the according region on the 2D detector hit map (lower density) with a diameter of 8 nm. The reconstruction parameters were adjusted so that the average spatial distribution of Fe atoms perpendicular to the layers fits the stacking distance of H-layers, in agreement to our STEM observations. The global atom distribution in the so achieved volume is similar to the model from literature and compatible with our STEM observations (see volume projection along planes in figure 8(b) with atomic sites related to Q-layer in orange and atomic sites related to H-layer in green).
The diameter of 8 nm of this cylindrical reconstruction volume is longer than one displacive modulation wavelength (q ≈ 4.5 nm). Thus, more than one displacive modulation wave is expected along each layer in the reconstructed volume. Even though the absolute amplitude of the modulation is only 1 Å, it is of interest to identify the displacive modulation within the reconstruction in order to investigate relations to elemental distribution in 3D. For a direct visibility, the viewing direction has to be set perpendicular to the modulation vector, similar to the HAADF STEM investigation. In the structure of franckeite, the modulation vector coincides with the crystallographic direction b for both Q-and H-layer. However, crystallographic directions in addition to the pole of H-layers and Q-layers are not accessible unambiguously in the present franckeite APT data. That is why, in order to find a viewing projection perpendicular to the modulation vector, rotations by 2 • of the reconstructed volume around the stacking direction (identified pole, z-axis of reconstructed volume with pole region as center) was performed. Variable clipping (down to less than 1 nm) of the volume along the viewing direction was applied. Despite the effort of this manual data treatment, the displacive modulation did not become uniquely visible. Reasons could be found in (i) the quality of the acquired raw data, (ii) the reconstruction algorithm, (iii) the optimization of the available reconstruction parameters or (iv) the method of analyzing the reconstruction.
The quality of the raw data of a van der Waals heterostructure can be limited by an inhomogeneous evaporation in comparison to homogeneous materials due to varying bonding types. In particular, the structure of franckeite consists of varying bonding types perpendicular to the van der Waals gap. Weak van der Waals bonds connect the H-and Q-layers while strong chemical bonds act between atoms within the single layers. The resulting varying evaporation barrier can cause a change in the order of evaporation in comparison to materials with homogeneous bondings [40,48] which cannot be considered in the available reconstruction software, yet. Additionally, homogeneous and isotropic atomic volumes and densities are generally assumed in APT reconstruction algorithms. These assumptions are not fulfilled in van der Waals heterostructures. Local volumes of higher atomic densities within the layers alternate with lower atomic densities in the van der Waals gaps. Thus, the in-depth parameter Δz along the reconstructed tip apex does not only depend on the atomic species, but is a function of the in-depth position within the complex structure itself.  In the presented reconstruction, the choice of parameters was set manually. This comprises the pole center, the diameter of the selected pole region as well as the adjustment of the reconstruction parameters to achieve spatial distributions of elements in agreement to STEM observations. Also the qualitative search for the viewing direction perpendicular to the modulation vector was done manually. Principally, all steps can be automatized enabling a powerful quantitative analysis which enhances the possibility to identify the modulation within the 3D reconstruction. However, its development is not part of the here presented work. Figure 8(c) visualizes 1D atomic fraction plots of the elements Sb and Pb (left) and Fe and Sn (right) perpendicular to the qualitatively resolved layers of the franckeite structure of figure 8(b). The concentrations of (Pb, Sb) and (Fe, Sn) show an anticorrelation which identifies the (Sb, Pb)-rich Q-and the (Fe, Sn)-rich H-layers, as described qualitatively in figure 8(b).
In figure 8(c), the measured local atomic fractions are quantified. In the Q-layers, Pb and Sb are present in concentrations of approximately 30 at% and 12 at%, respectively. In the H-layer, Sn has an average concentration of approximately 22 at% and Fe of 10 at%. Additionally, S is present in both layers (not shown). The values can be compared to the EDS-STEM quantification by Velicky et al [18]. There, the Pb concentration is in average slightly less than 30 at% within the identified Q-layers while the average concentration of Sb is higher than 12 at%. The concentration of Sn and Fe within the identified H-layer is approximately the same as measured by APT.
Although the nature of the signal that leads to the quantification is fundamentally different for both techniques EDS and APT, the results are approximately the same giving more confidence in the local concentrations. The EDS-STEM quantification depends mainly, besides spectrum treatment, on the proportionality factors, determined by standard samples, which convert the ratio of acquired and integrated signal of the elemental ionization edges into compositions. In contrast, single atoms are counted in APT so that no proportionality constants are needed to be determined. However, the determined composition depends mainly, here again besides spectrum treatment, on the evaporation probability per pulse for each element, which is assumed to be identical for all elements. Such an evaporation probability ratio between elements can vary for different experimental parameters (e.g. temperature, tip shape, voltage pulse fraction or laser power) and in particular, local variations can appear at pole regions due to influences of the local atomic structure on the effective electric field. Nevertheless, the fact that the two compared complementary techniques reveal similar local concentrations ensures the reliability in the APT data and motivates for improvements in order to reveal displacive modulations within the reconstruction, as discussed above.
In the concentration plots of the APT measurements, finer details than the anticorrelation of compositions in Q-layer and H-layer can be deduced. In the H-layer, the components Fe and Sn are expected to be distributed on only one atomic plane (average positions highlighted by horizontal dotted lines). Therefore, differences in the broadening of the concentration peaks could be expected for local compositions relative to the displacive modulation: a sharper concentration peak indicates a higher concentration of an element at the inflection point while a broader peak can be caused by a higher concentration at regions of the displacive modulation where its amplitude is reached. Comparing the concentration plots of Sn and Fe, there is no systematic difference that means that no preferential site occupation relative to the displacive modulation can be derived.
In the Q-layers, the shape of the concentration peaks of the components Pb and Sb, however, varies significantly. The Pb concentration shows a double peak (sometimes shoulders beside a main peak) while the Sb concentration has one broader peak with its maxima in between the double peak of the Pb concentration. A higher Pb concentration in the outer atomic planes of the Q-layers than in the inner atomic planes can be one explanation. Additionally, a relation to the displacive modulation of the Q-layers can be suggested from the observation when the highest concentration of Pb correlates with the convex regions of the Q-layer, even though a clear resolution of the displacive modulation was not visible within the reconstruction.

Layer stacking
Our results show that the angles β Q and β H change from stack to stack (figure 4) and that reported values may refer to average values. Contrariwise, α seems to be constrained by the modulation (figure 3(a)). Although not directly reported, these observations can be estimated from the comparison of both angles for reported structures (table 1 in [32]) in which α varies only within 1 • while β covers a range of up to 7 • . This result of constrained α and varying β shows a possible effect of the modulation on the stacking of incommensurate layers in such a van der Waals heterostructure. Parallel to the modulation wave vector, the stacking is constrained when the modulation adjusts the structure of both incommensurate layers. Perpendicular to the modulation wave vector within the layer, no constraint of the incommensurate structures occurs since there is no adjustment of the structures.
In addition to the stacking angles variability, we characterize stacking faults in the Q-layer such as rotations of 180 • about either the stacking direction or a (figure 4). Twins (180 • rotational stacking faults) have been previously reported by Makovicky et al [32], but characterized with b as their rotation axis. If those twins were characterized by observations along a, they could be also described as twins about the stacking direction ( a × b) enabling the possibility to describe all stacking faults in franckeite by rotations about the stacking direction. A confirmation of the rotation axis would be necessary by imaging single twins along various crystallographic directions.
Furthermore, we observed Q-layer stacking faults with a rotation of ±γ about the stacking directions (figure 7). 90 • stacking faults similar to our observations have been characterized before by Williams and Hyde [33]. In contrast to their observations, however, we do not find a resolved modulation in the [010] viewing direction for the Q-layer stacking faults. From our observation, we expect that the modulation in the Q-layer is generally independent of rotational stacking faults and only resolved in the [100] viewing direction. Since a Q and b Q are very close to each other, mismatch between a Q and b H could cause similar modulations as observed between b Q and b H while this is not the case for mismatch between a H and a Q or b Q .
The alignment of symmetrically equivalent direction between Q-layer and H-layer, but in different directions of the sample, would also be maintained by 60 • stacking faults in the Q-layer or 90 • stacking faults in the H-layer. However, they are not observed. A possible reason can be the modulation which breaks the in-plane symmetry of the basic structures. Therefore, studying the influence of crystallographic parameters on the formation of stacking faults and modulations in synthetically produced incommensurate van der Waals heterostructures, such as synthetic tin-selenium [51], could shed new light on the formation of both stacking faults and modulations.

Occupational modulation
The modulation has not only a relation to the stacking, but also to the measured characteristics within the layers. The intensity between projected atomic columns in the H-layer changes slightly over one modulation (in the HAADF-STEM image figure 3(c)). Interpreting the measured intensity as pure Z-contrast, Sn (Z Sn = 50) is in slight excess at regions of largest modulation displacement along c while Fe (Z Fe ) is in excess at the inflection points. This is in contrast to the suggestion made by Makovicky et al [32] that Sn is in excess at the inflection points while Fe is segregated at regions of largest modulation displacement along c in order to compensate varying bonding distances (compare figure 1(a)). Nevertheless, the variation is small. No direct elemental analysis, for example our APT investigation, identified significant occupation modulation in the H-layer, making clear statements difficult.
Significant intensity variations, however, are measured in STEM between projected atomic columns along one modulation within the Q-layer ( figure 3(b)). In the two outer atomic planes, a smooth transition between the high intensity at convex regions of the Q-layer and low intensity at concave regions is observed. The measured intensity modulation of one outer atomic plane is shifted by π relative to the other outer atomic plane of one Q-layer. The intensity pattern in the two inner atomic planes is more complex. Highest intensities are measured at the inflection points of the single planes. Additionally, there is an intensity oscillation between neighbouring atomic columns superimposed on the intensity variations related to the modulation. This effect will be discussed in the next section while we focus here on occupational modulation correlated to the displacement. The intensity variations of both inner atomic planes are approximately in phase (maxima/minima at the same location relative to the displacement wave). In contrast to our results, the XRD refinement of Makovicky et al [32] shows that the occupation modulations within all the atomic planes of the Q-layer follow approximately the same behavior with respect to the modulation displacement (no shift between outer and inner atomic planes with respect to their maxima/minima, figure 8 in [32]). In their refinement, the two outer atomic planes have not been considered separately and thus an occupational modulation phase shift between both outer atomic planes could not be accounted for. Our result shows a symmetrical correlation of the occupational modulation to the sinusoidal displacement, but with phase shifts between the atomic planes of the Q-layer.

Chemical ordering in Q-layer
We showed in figure 6 that the intensity of neighbouring atomic columns alternates in the two inner atomic planes of the Q-layer for the [010] projection indicating a chemical ordering. For example, the high intensity is explained by an higher composition of Pb (Z Pb = 82) in comparison to the weaker intensity which is explained by a higher composition of Sb (Z Sb = 51). A similar chemical ordering can be deduced for the [100] viewing direction in figure 3(b), even though it is more concealed due to the modulation and the possibly related elemental segregation discussed above and the projection effects within the thin lamella. The observed chemical ordering in both viewing directions implies a symmetry for rotations of γ (or γ − 180 • ).
On the base of these results, we propose in figure 9 an example of a Q-layer model which better describes our observations by chemical ordering in the two inner atomic planes of the Q-layer. Figure 9(a) is a plan view on a single inner atomic plane of the Q-layer. The composition of the sites within such a plane is chosen as follows. Every second heavy atomic site (not sulphur) has the same composition as described by Makovicky et al [32] for the two inner atomic planes. The other sites are chosen to have the same occupation probability as the outer atomic planes for Pb and Sb. Those occupation probabilities are adjusted in order to keep a composition consistent with that reported for the Q-layer [32]. Occupation modulations are not considered in this scheme.
The cross-section views of such a model, shown in figures 9(b) and (c), are in qualitative agreement with our experimental observations for STEM in both viewing directions [100] and [010] and for APT concentration plots. In particular, we see that the proposed ordering in the two inner atomic planes is more consistent with the contrast change such as the intensity in every second atomic column is lower than in the two outer atomic planes ( figure 6). Furthermore, highest Pb concentration is present in the outer atomic planes of the Qlayers while the concentration of Sb is highest in the inner part of the four atomic planes thick Q-layer that is consistent with our APT observations ( figure 8(c)). As such a chemical ordering in the structure, the basic unit cell of the Q-layer needs to be increased along a and b in order to keep the symmetry of the superspace group of franckeite characterized by Makovicky et al [32] or another superspace group needs to be assigned. Note that the varying distance between the projected atomic columns in the inner atomic planes is also respected in our proposed model.

Chemical ordering in H-layer
We observed intensity variations in the heavy atomic columns of the H-layer for the [010] viewing direction (figure 5) as an alternation between one atomic column of lower intensity and two with higher intensity. Similar contrast variations in the H-layer have been imaged by Velický et al [18]. This is anomalous with homogeneous (Fe, Sn) site occupation of 0.57 Sn and 0.43 Fe for the [010] viewing direction as reported previously [32]. We propose, based on our Z-contrast observations, a structural model example of the H-layer in figure 10. Figure 10 , respectively, in agreement with the hexagonal structure. The proposed ordering within the H-layer is one of the possibilities to explain the contrast observations. Mixed atomic occupation probability with preferential Fe occupancy instead of pure Fe sites would also explain the observed contrast variations.

Conclusion
We have investigated the atomic structure of the mineral franckeite using HAADF-STEM which allows more direct interpretation of image contrast. In particular, our work focused on structure projections of the Q-layer and the H-layer within the van der Waals gap. Complementary to the HAADF-STEM investigation, we applied APT as a chemical analysis method revealing local chemistry down to atomic scale. APT offers additionally a high potential for new insights into such modulated van der Waals heterostructures due to 3D analysis at atomic scale. To realize its full potential, further improvements are required such as the consideration of the influence of van der Waals gap structures on the evaporation and, as a consequence, on reconstruction algorithms. We compared our observations to a structure model of a recent XRD refinement and former (S)TEM studies at lower spatial resolution. We conclude that the overall microstructure is confirmed including the modulation of the incommensurate layers. Beyond this, however, we show essential details about the stacking of the Q-layer and the H-layer, the effects of the modulation and chemical ordering within each of both layers not previously observed.
Stacking faults in the Q-layer have been characterized as rotations of ±γ or 180 • about the stacking direction. Furthermore, the basic units of Q-layer and H-layer stacks are shown to be precisely fixed in α while β varies for both layers.
In both layers, Z-contrast variations related to the modulation are observed and interpreted as a response on a segregation of differently sized elements to expanded and compressed regions which originate from the displacements along the modulation.
Chemical ordering in the two inner atomic planes of the Q-layers is interpreted from HAADF-STEM image contrast and a structure model is proposed. Similarly to the pseudo tetragonal Q-layer, chemical ordering is observed in the pseudo hexagonal H-layer and a structure model is presented. Both orderings cause an increase of the basic unit cells, or a change in the superspace group needs to be defined, and is, as such, suggested for further work.
Our results can be used as a base for an extensive study of the electronic properties of franckeite comparing experimental investigations with simulations. Specifically, we provide new insights into the internal structure of the Q-layer and H-layer and their stacking which are expected to significantly influence the electronic structure. This information is essential to enable further work to provide a better understanding of the formation and orientation alignment of the incommensurate layers and the origins of the modulation.