Observation of a network of structural modulations

FeSe presents an excellent playground to explore the interplay of electronic nematicity and symmetry breaking strain due to its structural simplicity and the absence of magnetic ordering that is present in many other Fe-based superconductors^29,30. In principle, various experimental methods can be used to apply strain to a material, such as voltage-controlled piezoelectric setups^{11,13,25,31,32}, mechanical actuators²⁷, differential thermal contraction^24,33, and heteroepitaxial film growth^34–37. In this work, we use molecular beam epitaxy (MBE) to grow FeSe thin films (a = 3.8 Å) on SrTiO₃(001), a substrate with a ~2% lattice mismatch (a = 3.9 Å) (Fig. 1a, b, see “Methods” section). We find a 2D network of modulations emerging at the surface of FeSe, propagating approximately along the Fe–Fe lattice directions (Fig. 1c, d). As we will subsequently show, this in turn leads to a spatially varying strain at the surface. The spatial distribution of modulation lines is qualitatively similar to those in heteroepitaxially-grown heterostructures of other chalcogenides^35–37 and arsenides³⁴. In our FeSe films, ranging from 3 to 6 monolayers in thickness, this distance between neighboring modulation lines is approximately 15–20 nm, consistent with the spacing determined from cross-sectional transmission electron microscopy³⁸ and roughly consistent with the expected value based on the lattice constant mismatch between FeSe and SrTiO₃(001) (Supplementary Note 1).

Fig. 1

Structural characterization of FeSe thin films.

a Reflection high-energy electron diffraction (RHEED) images (15 keV) of SrTiO₃(001) (top) and FeSe (bottom). b Top view schematic of FeSe and SrTiO₃ crystal structures, in which a-axis and b-axis denote the nearest-neighbor Fe–Fe directions. Black squares in b outline the Se (left) and Ti (right) unit cells. Brown, dark green, and yellow spheres in the FeSe crystal structure denote Fe, top Se, and bottom Se atoms, respectively. Light green, blue and red spheres in the SrTiO₃ crystal structure denote Sr, Ti, and O atoms, respectively. c 3D rendered large-scale STM topograph that shows SrTiO₃(001) substrate, 1 monolayer (ML) FeSe, 3 ML FeSe and 4 ML FeSe terraces. The periodic structural modulations can be more easily distinguished on the 3 ML and 4 ML FeSe layers. d Magnification of the region outlined by the black square in c, showing the structural modulations in the STM topograph T(r). Box smoothing over ~0.36 nm is applied to the topograph in d. STM setup conditions: c I_set = 10 pA, V_sample = 1 V; d I_set = 10 pA, V_sample = 600 mV.

Visualizing electronic nematic domains

To characterize electronic properties of our thin films, we use low-temperature spectroscopic-imaging scanning tunneling microscopy. Although a single-unit-cell thick FeSe grown on SrTiO₃(001) can exhibit superconductivity with ~10–15 meV pairing gap (Supplementary Fig. 7), the surface of thicker FeSe films (~a few to ~50 nm thick²¹) grown on the same substrate typically does not show superconducting behavior (Supplementary Note 8). By acquiring dI/dV(r,V) or I(r,V) spectra (where I is the tunneling current, V is the voltage applied to the sample, and r is the relative xy-position of the tip) on a densely-spaced pixel grid, we are able to visualize spatial variations in electronic density of states as a function of energy and position. We focus on an area shown in Fig. 1d, where we observe two striking features not immediately obvious from STM topographs (Fig. 2b). First, we can discern dark irregularly-shaped contours enclosing parts of the sample (denoted by white dashed lines). Second, we observe horizontal (an example denoted by green arrows) or vertical stripes (purple arrows) oriented along Fe–Fe lattice directions, with ~1.8 nm nearest-neighbor distance, which do not disperse as a function of energy in dI/dV maps (Supplementary Fig. 1). This has been interpreted as the formation of charge-stripes in the electronic nematic state^21,28, with the direction of electronic nematicity rotating by 90° across the domain boundaries. This interpretation is further supported by dispersive C₂-symmetric modulations pinned to individual dumbbell-shaped impurities, which also rotate by 90° across the same boundaries (Fig. 2d). Putting this information together, we conclude that the sample consists of two types of electronic nematic domains. We note that the smallest domains observed here are only ~100 nm², significantly smaller than those in bulk single crystals^14,23. As we will show, the reduced domain size can be attributed to the rapidly varying strain landscape (Fig. 3).

Fig. 2

Visualizing electronic nematic domains.

a Atomically-resolved STM topograph T(r), and b tunneling current map I(r,V = 60 mV) over an identical region, showing nematic domain boundaries denoted by white dashed lines. Horizontal charge-ordered stripes (denoted by green arrows) rotate into vertical stripes (denoted by purple arrows) across the boundaries between A and B. c Average differential conductance dI/dV spectra over nematic domains A and B, vertically offset for clarity. d Magnification of the regions near three different impurities in a outlined by the light gray, dark gray and black squares, which are located in the electronic nematic domain A, B and on the boundary, respectively. The last four figures in each row are dI/dV(r,V) maps encompassing each impurity, showing unidirectional dispersion from −20 to −50 mV. The blue and white arrows serve as guides to the eye for the dispersions. We note that the direction of this signal is not dependent on the impurity shape, as all three impurities shown are located at the equivalent Fe site. Symmetry of the electronic signal around impurities at the boundary of the two regions is broken down further, with one peak along each a-axis and b-axis, further supporting the intrinsic symmetry-broken electronic state of the system. STM setup conditions: a I_set = 110 pA, V_sample = −100 mV; b I_set = 110 pA, V_sample = −100 mV; d I_set = 110 pA, V_sample = −100 mV, V_exc = 5 mV.

Fig. 3

Strain analysis and comparison between experimental and theoretically predicted strain maps.

a Atomically resolved STM topograph, where a-axis and b-axis denote the nearest-neighbor Fe–Fe directions. Dashed lines are guides to the eye for the structural modulations. Strain tensor components b $$ u_{\mathrm{aa}}\left(\mathbf{r}\right)$$ and c $$ u_{\mathrm{bb}}\left(\mathbf{r}\right)$$ derived from a. d Symmetric strain map $$ S\left(\mathbf{r}\right)\equiv u_{\mathrm{aa}}\left(\mathbf{r}\right)+u_{\mathrm{bb}}\left(\mathbf{r}\right)$$. e Antisymmetric strain map $$ U\left(\mathbf{r}\right)\equiv u_{\mathrm{aa}}\left(\mathbf{r}\right)-u_{\mathrm{bb}}\left(\mathbf{r}\right)$$. f Schematic of an edge dislocation. Vector b denotes the Burgers vector. Blue (red) circles schematically represent the FeSe (SrTiO₃) lattice. Theoretical maps of g $$ u_{\mathrm{aa}}\left(\mathbf{r}\right)$$, h $$ u_{\mathrm{bb}}\left(\mathbf{r}\right)$$, i $$ S\left(\mathbf{r}\right)$$ and j $$ U\left(\mathbf{r}\right)$$. The inset in g shows the strain created by a single dislocation line as a function of the distance from it. In our model described in detail in Supplementary Note 7, Burgers vector b is set to be 0.53 nm, which is the Se–Se distance along [110] direction on the top or the bottom of an FeSe monolayer (ML), and the thickness d is set to be 1.6 nm, or 3 ML of FeSe. STM setup condition: a I_set = 110 pA, V_sample =−100 mV.

Strain analysis

To measure local structural distortions, we start with an atomically-resolved STM topograph (Fig. 3a and Supplementary Fig. 3), and apply a geometric phase analysis method^33,37 based on the Lawler–Fujita drift-correction algorithm⁴. This method allows us to determine the displacement of atoms $$ \mathbf{u}\left(\mathbf{r}\right)=u_{\mathrm{a}}\left(\mathbf{r}\right)\widehat{\mathbf{a}}+u_{\mathrm{b}}\left(\mathbf{r}\right)\widehat{\mathbf{b}}$$ with picoscale resolution. The four-component strain tensor $$ u_{\mathrm{ij}}\left(\mathbf{r}\right)\equiv \mathrm{d}{u}_{\mathrm{i}}\left(\mathbf{r}\right)/\mathrm{d}{r}_{\mathrm{j}}$$ (where i, j = a, b) can be used to extract different types of strain deformations. For example, $$ u_{\mathrm{aa}}\left(\mathbf{r}\right)$$ represents the change of the lattice constant along the a-axis (relative to the average lattice constant in the field-of-view), with positive (negative) values denoting local tensile (compressive) strain. Taking into account strain along both lattice directions, it is convenient to define symmetric strain component: $$ S\left(\mathbf{r}\right)\equiv u_{\mathrm{aa}}\left(\mathbf{r}\right)+u_{\mathrm{bb}}\left(\mathbf{r}\right)$$ and antisymmetric strain component: $$ U\left(\mathbf{r}\right)\equiv u_{\mathrm{aa}}\left(\mathbf{r}\right)-u_{\mathrm{bb}}\left(\mathbf{r}\right)$$. The latter is particularly useful as a quantitative measure of structural anisotropy between the two lattice directions. We apply the strain analysis algorithm to the same area as in Fig. 2 to obtain strain maps (Fig. 3b–e). Tensile strain is observed along the modulation lines in both $$ u_{\mathrm{aa}}\left(\mathbf{r}\right)$$ and $$ u_{\mathrm{bb}}\left(\mathbf{r}\right)$$ maps, which is sandwiched by two ribbons of compressive strain. Further away from the modulation lines, there is tensile strain again in broader areas. To support the robustness of the strain algorithm, we note that strain maps calculated from STM topographs acquired in a range of different biases look qualitatively indistinguishable (Supplementary Figs. 2 and 4). Theoretically calculated^39,40 strain maps based on a network of edge dislocations (Fig. 3g–j), Supplementary Note 7) also show a close resemblance to our experimental data (Fig. 3b–e). Therefore, we can conclude that the observed strain can be modeled well by a misfit dislocation network, with small differences that could be attributed to intrinsic orthorhombic distortion accompanying each electronic nematic domain.

Correlation of nematic domains and antisymmetric strain

To investigate strain inhomogeneity further, we superimpose the outlines of electronic nematic domain boundaries on top of the antisymmetric strain map U(r) (Fig. 4c). For the electronic nematic domain A, where the charge-stripe wave vector is oriented along the a-axis, it is expected that the lattice constant along the a-axis (a₀) is larger than that along the b-axis (b₀)²¹. Indeed, the average antisymmetric strain within this region is consistent with this expectation (orange color in Fig. 4c). Similarly, within the electronic nematic domains B, where the charge-stripe wave vector propagates along the b-axis, we find that on average, b₀ is greater than a₀ (purple color in Fig. 4c). This is consistent with the global picture revealed in elasto-resistance experiments of bulk single crystals, where electronic nematic response followed the direction of externally applied anisotropic strain^11,12. However, our ability to probe both local anisotropy and electronic nematicity at the nanoscale enables us to explore their correlation at previously inaccessible atomic length scales.

Fig. 4

Interplay of electronic nematicity and local strain.

a STM topograph T(r), and b tunneling current map I(r, V = 60 mV) acquired over an identical region of the sample. Dashed lines in a, c denote the network of structural modulations. Solid lines in b, c outline the electronic nematic domain boundaries. c The antisymmetric strain map U(r) calculated from a. d A schematic of electronic nematic domains determined from b, where A and B refer to orthogonal electronic nematic domains. e Histogram of the antisymmetric strain magnitude in electronic nematic domains A and B in c. f Magnification of a single nematic domain B and the antisymmetric strain distribution within it. White dots and double lines depict the Fe impurities and electronic stripes associated with the direction of electronic nematicity. Dashed lines in f denote the borders between regions of positive and negative antisymmetric strain. STM setup condition: a I_set = 110 pA, V_sample = −100 mV; b I_set = 110 pA, V_sample = −100 mV.

Interestingly, we find that the distribution of antisymmetric strain values within each electronic nematic domain is highly inhomogeneous (Fig. 4c, e). Moreover, not only do we find variations in magnitude, but also in the sign of the anisotropic strain. In other words, within electronic nematic domain A where the direction of electronic nematicity would dictate that a₀ should be larger than b₀, we observe sizeable regions (31% of the area) where this trend is opposite (an example denoted by blue arrow in Fig. 4c). The same observation is also apparent in the orthogonally oriented electronic nematic domain B (yellow arrow in Fig. 4c), where regions with a₀ > b₀ comprise 21% of the whole area. Thus, our experiments reveal a local decoupling of structural and electronic anisotropy in an electronic nematic system. We emphasize that this result does not rely on the theoretical strain model in Fig. 3 or the nature of strain modulation lines. The local strain on the surface is directly determined from atomically-resolved STM topographs, and then correlated with simultaneously acquired dI/dV maps where we can see electronic nematic domains.

MBE growth

FeSe films were grown on Nb-doped (0.05 wt%) SrTiO₃ (001) (Shinkosha). The substrates were sonicated in acetone and 2-propanol, followed by annealing in O₂ supplied tube furnace at 1000 °C for 3 h. This step created √13 × √13 R33.7° surface reconstruction, which has been observed in both RHEED images and STM topographs (Supplementary Fig. 6). The substrates were then introduced into our MBE system (Fermion Instruments) with a base pressure of ~4 × 10⁻¹⁰ Torr. Continuously monitored by a pyrometer, the substrates were slowly heated up to ~400 °C for growth. Fe (99%) and Se (99.999%) were co-evaporated from two Knudsen cells held at 1100 and 145 °C, respectively, corresponding to flux rates of 9.87 × 10⁻⁵ atoms/(sec*Ų) for Fe and 3.37 × 10⁻³ atoms/(sec*Ų) for Se measured by the quartz crystal microbalance. At these relatively low flux rates, it takes about 28 min to form each monolayer, followed by post-growth annealing at ~450 °C for 2–3 h. After growth, the samples were either quickly transferred to the STM using a vacuum suitcase chamber held at ~1 × 10⁻⁹ Torr or capped with ~50 nm thick amorphous Se layer, and de-capped in the STM chamber at ~500 °C for 2 h. We note that FeSe films with modulations were observed in both films transferred by suitcase and de-capped thin films (Supplementary Note 4). We hypothesize that the modulations in multilayer FeSe may be related to the √13 × √13 R33.7° surface reconstruction of SrTiO₃ (001)³⁸.

STM measurements

STM data was acquired using a Unisoku USM1300 STM at the base temperature of ~4.5 K. Spectroscopic measurements were made using a standard lock-in technique with 915 Hz frequency and bias excitation as detailed in figure captions. STM tips used were home-made chemically-etched tungsten tips, annealed in UHV to bright orange color prior to STM imaging.