Theory and application of the vector pair correlation function for real-space crystallographic analysis of order/disorder correlations from STEM images

The manner in which local structure and chemistry deviate from the average crystalline unit cell has a great influence on the properties of many functional materials. This is particularly true for technologically important material categories, such as classical and relaxor ferroelectrics. These materials find applications in energy storage capacitors, computer memory, transducers, and actuators, among others. Tetragonal tungsten bronze (TTB) is a metal-oxide structure compatible with a wide range of compositions, similar to the widely studied perovskite structure, including many lead-free compositions with ferroelectric or relaxor behavior.^1,2 Additionally, various TTB compositions are known to exhibit other unique properties, such as very high Curie temperatures relative to most perovskite-based ferroelectrics,³ pyroelectricity,⁴ photorefractivity,⁵ and second harmonic generation.⁶ TTBs are, however, much more structurally complex than perovskites and commonly exhibit modulated superlattices.^1,7 The prototypical tetragonal unit cell has four main types of cation sites, as shown in Fig. 1. These are often fully occupied, though many “unfilled” TTBs exist, which have a large fraction of vacancies, primarily on the A1 sites. Octahedral connectivity is a complex corner-shared network in the a–b plane with direct corner-sharing along the c axis. In modulated TTBs, the true c lattice parameter is doubled from the prototype cell, primarily as a result of alternating octahedral tilts along that axis.^2,7 Other observed deviations from the high symmetry tetragonal structure result in an orthorhombic unit cell due to symmetry breaking from octahedral tilts, displacements of the A2 sites along the [110] type directions, and a small (4 pm) lattice parameter difference between the a and b axes.² 

As shown by Zhu et al., the A-site cation sizes are predictive of either commensurate or incommensurate modulations.¹ They demonstrated that, with increasing average A-site cation size and A1 site tolerance factor, the onset of incommensurately modulated structures coincides with relaxor behavior. Relaxors are known to have small polar regions with length scales on the order of nanometers. It is thought that chemical or structural disorder in the lattice disrupts polar ordering, preventing the formation of long-range ferroelectric domains.⁸ In addition to chemical disorder on the A or B sites, in TTBs the structural disorder inherent to incommensurate modulations may also contribute to the disruption of long-range ferroelectric correlations and favor relaxor behavior.¹ 

Woike et al. refined the structure of the unfilled TTB, Sr_1−xBa_xNb₂O₆ (SBN) with x = 0.39, from primary x-ray reflections and first-order satellite peaks. They used the super-space approach to quantify the modulation, refining the structure as tetragonal, and found it has two distinct wave vectors, which are related to each other by the lattice symmetry: $q1=α(a*+b*)+12c*$ and $q2=α(a*−b*)+12c*$⁠. Furthermore, in agreement with previous studies of TTB materials, they concluded that the structural modulation is primarily on the oxygen and the A2 sub-lattices. Since they were not able to integrate higher-order satellite spots in their refinement, it remained an open question whether the material was truly tetragonal, having two continuous modulation vectors, or if it had an orthorhombic structure with two possible modulative orientations.⁷ For SBN with x = 0.5, Bursill and Lin determined the presence of two orthorhombic variants using a combination of electron diffraction and dark-field imaging.⁹ 

Due to their structural complexity, TTBs have primarily been studied experimentally by low spatial resolution methods, such as x-ray diffraction,^2,7 parallel-beam electron diffraction, and conventional dark-field transmission electron microscopy (TEM).^1,9 While these techniques are valuable, they tend to average the structure over large length scales and cannot directly observe local correlations, which are key to understanding the origin of relaxor ferroelectricity. Here, we apply aberration-corrected scanning transmission electron microscopy (STEM) to study Ba₅SmSn₃Nb₇O₃₀ (BSSN) directly at the atomic scale. We use BSSN as a model material to develop analysis methods for quantifying and visualizing short- to long-range structural correlations, such as modulations, from atomic resolution STEM images.

BSSN is a TTB exhibiting an incommensurate superstructure and relaxor behavior prototypical of its class. The relaxor behavior is evidenced by a broad and frequency-dispersed maximum of the dielectric permittivity below 250 K. At room temperature, however, BSSN displays linear PE loops, indicating a paraelectric phase, though the incommensurate modulation is maintained.¹⁰ The incommensurate nature of BSSN is apparent from the superlattice reflections in the [110] selected area diffraction pattern (SADP) in Fig. 2. Atomic resolution STEM offers unique and complementary insights for understanding local structure when compared to other methods that may achieve higher accuracy but tend to average the structure over large length scales. For thin samples, STEM high-angle annular dark-field (HAADF) images are directly quantifiable and atomic columns can be located with sub-pixel accuracy by fitting with a two-dimensional Gaussian function.¹¹ Scanning techniques, however, are subject to a number of image distortions, such as drift (from thermal or charging effects), scan-line origin noise, and scan vector misalignment.^12,13 Therefore, care must be taken to correct these to enable accurate measurements beyond a range of 1–2 nm, particularly important for modulated structures.¹⁴ The processes used for drift and distortion correction in this study, as well as atomic column position determination, are detailed in Sec. V. The question we seek to answer in this work is whether STEM imaging is able to quantify structural information accurately from short- to long-range length scales, and on a sublattice basis, while retaining vectorized orientation information. We are motivated to perform this analysis in real space so that it can be applied to heterogeneous materials for which ensemble information is not sufficient.

Our approach to answer this question is to develop and implement the vector pair correlation function (vPCF) method to quantify the structural modulation and other displacive correlations in BSSN from STEM images. To place this in context, we will first describe the well-known radial PCF (rPCF). The rPCF is the ratio of the local to the average density [ρ(r) and ρ₀, respectively] as a function of distance (r) from a reference atom,¹⁵ 

where N is the number of atoms in the system and δ is the Dirac delta function. The sums being over all the atoms in the system, the PCF yields an ensemble average representation of the local “neighborhood” at distance r from a reference atom. A rPCF derived from a perfect lattice exhibits delta functions at the exact lattice spacing distances (indicating perfect ordering), while a completely random system would manifest as g(r) = 1 (indicating no correlation between positions, regardless of distance). An amorphous atomic structure will exhibit some order in the rPCF, as indicated by an initial intense peak (due to the nearest-neighbor bond length) and a few subsequent, increasingly broad peaks. Longer correlations, however, are absent in amorphous structures, and g(r) tends to unity with increasing r. In an imperfect crystalline structure, peaks in the PCF will broaden. In the case of a modulated crystal, or for structures with discrete displacements, peak splitting may occur. Overall, PCFs offer a simple and concise way to summarize local structural deviations from the global structure. Historically, PCFs have been primarily developed and used to analyze the structure of amorphous materials, such as liquids and glasses from x-ray or neutron scattering data. PCFs are increasingly applied to crystalline scattering data as the influence of local structure on many materials’ properties becomes ever more apparent.¹⁵ 

The one-dimensional rPCF described above is helpful in understanding the degree of crystallinity, quantifying short-range order in a system, or for structure determination from total scattering data using Monte Carlo methods.¹⁶ The concept has also been applied to real-space STEM data, by directly computing the projected rPCFs from specific zone-axis lattice images. Radial PCFs, however, do not explicitly preserve orientational information.^17,18 Additionally, interpretation of an rPCF can be complicated by overlapping peaks, especially at increasing distances. rPCFs, therefore, are of limited utility for studying complex or long range disorder phenomena. An example of radial partial PCFs from [110] zone axis STEM images of BSSN is shown in Fig. 3. Note that the peaks have a significant width and that some are wider than others, indicating more variability in displacements along the corresponding crystallographic directions. Additionally, note that the apparent single peak at 12 Å is actually two peaks separated by only 4 pm in the average BSSN structure.

To capitalize on the usefulness of STEM image-derived rPCFs while avoiding the noted drawbacks, we have extended the concept to multiple dimensions as the vPCF. This is similar in concept to 3D pair distribution functions (PDFs), which may be derived from 3D single crystal scattering data acquired with synchrotron or neutron sources.¹⁶ In this paper, we present a two-dimensional vPCF applied to atomic-column locations in a STEM image, in principle, however, the method could be applied to 3D atomic position data, resulting in a 3D vPCF. Adapting from Eq. (1), a vPCF is mathematically described by

Here, $r̄−r̄vu$ is the inter-atomic column vector relative to the vu atom column. To compute the vPCF from a finite-sized image, two-dimensional column-to-column vectors are binned into pixels of dimensions Δx, Δy. $g(r̄)$ is, then, normalized by the chosen pixel size.

From distortion-corrected STEM images, it is relatively straightforward to calculate the vPCF once atomic column locations have been determined. Depending on the types and ordering of displacements, different methods to quantify the vPCF peak shape may be appropriate. In general, it is likely that PCF peaks do not represent Gaussian distributions, so determining shape by Gaussian function fitting could lead to analytical errors. vPCF peak quantification can be accurately and efficiently accomplished using central moments.¹⁹ First-order moments are used to calculate the mean of a distribution, while second-order moments give information about the width of the distribution and orientation of anisotropy. From the second-order moments, the covariance matrix can be constructed and eigenvalue decomposition can be used to find the axes of the equivalent ellipse. The orientation and eccentricity of this ellipse will be primarily used in this study. Figure 4 illustrates these concepts.

To illustrate how various structures will affect the resulting vPCF, a model lattice system was generated. A random, Gaussian distribution of displacements was added to the lattice to mimic experimental measurement error and the presence of defects. For the purposes of explanation, we use large displacement distributions (relative to those observed experimentally in BSSN) for these models. All vPCFs are presented as a axis projections of an orthonormal lattice so that b and c are in the plane of the image. The resulting vPCFs are shown in Fig. 5.

Figure 5(a) presents a model vPCF for a lattice with only random displacements. The geometry of the peaks is the same as the lattice geometry, and isotropic peak spreading is present due to the random displacements in both the x and y directions. In Fig. 5(b), a vPCF from a lattice with a structural modulation vector of $q=12c*$ with alternating displacements in the b direction is illustrated. Peak splitting is observed for vPCF peaks at odd multiples of c from the origin due to the anti-correlated relationship to the displacement of the origin atom column. At even multiples of c, which are highlighted by the horizontal phase lines, displacements are correlated, so no splitting or additional broadening occurs. Figure 5(c) presents a vPCF for the case of a commensurate modulation vector having components in both the x and y directions, $q=13b*+12c*$⁠. Here, the modulation phase lines are shown and intersect lattice points at multiples of 2b and 3c from the origin. In this vPCF, peaks midway between the constant phase lines show the broadest spreading or splitting, while those progressively closer have narrower, more isotropic distributions. The moments described previously were used to calculate the eccentricity of each peak, with the coloring of the peaks indicating this value.

It should be noted that atomic pair correlation vector information is included in the electron diffraction pattern (Fig. 2) and that it is possible to calculate a vPCF directly from such data. However, a number of challenges make this difficult experimentally, particularly the presence of multiple scattering of electrons and the limited dynamic range of the detectors available in most instruments. Moreover, it is not possible to produce partial vPCFs by this method, and we will show that the ability to quantify sublattice-specific information is a substantial advantage of vPCFs calculated from direct lattice images.

Figure 2 shows the [110] SADP from BSSN exhibiting incommensurate satellite spots. Note that the pairs of satellite diffraction spots visible in this pattern are of the same modulation vector type (i.e., either q₁ or q₂). The opposite modulation vector type lies in a perpendicular plane of the reciprocal lattice, as shown in Fig. 2(d), and so corresponding spots do not appear in the diffraction pattern. In Fig. 2(b), a conventional dark-field TEM image formed using a pair of satellite spots reveals a mottled intensity distribution in the grain. This indicates that only the bright regions have the specific modulation vector corresponding to the diffraction spots used to form the image. Figure 2(c) shows a section of the first-order Laue zone ring from an [001] zone axis SADP. In this case, sets of four satellite spots are visible, indicating that the diffracting region of the sample possesses both types of modulation vectors. Thus, it is most likely that a modulation is present throughout the entire grain, but that it fluctuates between q₁ and q₂ locally corresponding to two variants of orthorhombic symmetry. The ideal orientation to observe the effects of the modulation at the atomic scale is the [110] zone axis, where one of the modulation vectors will be parallel to the image plane. The other modulation vector, however, is out-of-plane, and corresponding regions will appear unmodulated.

HAADF STEM images of BSSN were acquired along the [110] zone axis. At this orientation, two projections of A2 columns are present, and we refer to these as the A2₁ and A2₂ sublattices in this work. (See Fig. 1 for an illustration of this lattice projection.) The A2₁ columns are those that have the long axis of their distorted pentagonal coordination polyhedron (largest corner-to-corner distance) perpendicular to the zone axis. Displacements are known to occur primarily along this axis and therefore are visible in the STEM images.² A2₂ columns’ long polyhedral axes are parallel to the viewing direction, and, as a result, the imaging mode is not sensitive to the expected displacements on these sites.

Atomic column displacement vectors are a common way to display position data in STEM images. Often, this vector is measured relative to the average position of near-neighbor columns of a different sublattice type. For example, in perovskites, B-site displacements may be measured from the average position of the four nearest-neighbor A-site columns. Some of the atomic columns in [110]-oriented BSSN images, however, have no obvious sublattice from which to measure the relative displacement. Therefore, in this study, we define displacements as the vector from a rigid global reference lattice point to the atomic column position found by Gaussian fitting. Rough registration of the reference lattice was made using the FFT and an atom column selected from the image as the origin. To achieve an accurately reference lattice, its basis vectors and origin position were subsequently refined by least-squares minimization of the displacements on the A1 sites. This sublattice was chosen because it appears as prominent peaks in the image and it is expected that the site has minimal modulative displacements, as evident in Fig. 3. Displacement vector plots shown here are based on the final, refined reference lattice. Figure 6 shows representative displacement vector plots for the A1, A2₁, and B2 sublattices. A1 exhibits random and relatively small displacements. A2₁ shows a large degree of displacement ordering, which we will show to be consistent with the modulation. Finally, on the B2 sublattice, displacements are small overall, but we measure short-range correlated displacements primarily along [001].

Figure 7 presents the quantified local structure mapped onto the STEM image where colored lines show deviations of the A1 and A2₁ near-neighbor distances from the average structure. Comparing the two plots, A1 distances vary much less than those of A2₁. No clear pattern of distance variations is seen for A1. On the other hand, it is immediately obvious that in the leftmost portion of the image, the A2₁ sublattice is modulated, while weak or no modulation is observed in other regions of the image. Note that the alternating long and short NN distances (indicated by red and blue lines) along the [001] direction agree with the $12c*$ component of the modulation. Columns of unit cells with NN distances close to the ideal separation appear approximately every three unit cells along the $[11̄0]$ direction in the modulated region, agreeing with the ≈0.3(a* + b*) component. These observations are consistent with our hypothesis of modulation domains corresponding to two orthorhombic variants. Also note the much longer correlation length of the modulation along [001] as compared to $[11̄0]$⁠. While the expected modulation can be observed in this plot, it is possible to quantify it more precisely using vPCFs.

Figure 8 presents partial vPCFs for the A1, A2₁, B1, and B2 sublattices. Peaks are represented by ellipses proportional to the size and shape of each peak’s calculated equivalent ellipse. The color of the markers indicates the ellipse eccentricity. In the A1 and B1 partial vPCFs, peaks are generally small and no strong peak shape pattern emerges, as expected for a sublattice without correlated displacements. In the A2₁ partial vPCF, the orientation of the ellipses is along the $[11̄0]$ direction, as expected from the known displacement direction on the site. The $12c*$ component of the modulation is readily obvious, as indicated by alternating round and elliptical peaks along the c axis of the partial vPCF. Example peaks are shown in the insets. Note the multimodal distribution of the highly elliptical peaks, which is consistent with strongly modulated displacements. The peak eccentricity pattern along the $[11̄0]$ axis exhibits minima at multiples of three unit cells, consistent with α ≈ 0.3. Furthermore, we note that along the [001] direction, the eccentricities retain a strongly alternating pattern for at least 20 unit cells (note: Fig. 8 only displays the vPCF for eight lattice vectors in this direction). Along the $[11̄0]$ axis, on the other hand, the pattern rapidly decays, showing an anisotropy in the correlation length of the modulation. This anisotropy is likely a result of long-range coordinated tilts along [001] oxygen octahedral chains, while the complex octahedral connectivity in the a–b plane frustrates the corresponding tilt coordination. To measure the modulation vector components from the vPCF, we performed Bayesian inference, modeling the periodic eccentricity pattern as

where $A$ is the amplitude, b is the bias, and $ϕ(r̄,q̄)$ is the phase angle found by the projection of vPCF lattice vectors onto the modulation vector according to

Only peaks in the region of the vPCF 0 ≤ x ≤ 4 lattice cells (along the $[11̄0]$ direction) were used for inference because of the anisotropic and limited correlation length observed. The inference revealed α = 0.314 with 95% predictive interval [0.310, 0.318], in close agreement with the SADP. The eccentricity phase angle data are plotted in Fig. 8(e) along with the model (using the mean inferred parameter values). This plot highlights the strong periodicity of the vPCF peak shapes. Note the loss of coherence between the model (red line) and experimental data for the set of peaks at x > 4 lattice cells along $[11̄0]$ [gray data points in Fig. 8(e)]. For x ≤ 4 [blue data points in Fig. 8(e)], the model and data are highly coherent to more than nine multiples of the modulation vector. This further illustrates the correlation length anisotropy of the modulation.

Finally, considering the B2 partial vPCF, we observe a short-range correlation in displacements along the [001] axis over a distance of approximately two unit cells. This is indicated by the small eccentricities of the [001] and [002] peaks, but with predominantly eccentric peaks elsewhere with the major axis along the [001] direction. Despite the near-zero remnant polarization for BSSN at room temperature,²⁰ short-range order on the B2 sublattice indicates the presence of residual local polar distortions persisting from the lower-temperature relaxor behavior. This is consistent with other relaxors, which have been shown to retain polar regions several hundred degrees above their permittivity maximum.²¹ 

Here, we have developed the vPCF as a real-space method for visualizing and quantifying the local to long-range disorder of complex structures from STEM lattice images. Furthermore, we have demonstrated its utility when applied on a sublattice basis as a partial vPCF. We have shown that the globally observed tetragonal structure of BSSN is locally broken by two variants of orthorhombic symmetry, which are related to the incommensurate modulation. From quantified STEM images and electron diffraction, we deduce that the incommensurate modulation is discontinuous, with both orthorhombic symmetries coexisting in the same grain. By analysis of the vPCFs, we observe that, among the cation sublattices, only the A2 site is modulated. Furthermore, the vPCF reveals the correlation length of the modulation to be anisotropic, and we infer that this is a result of the complex connectivity of the oxygen octahedral network in the a–b plane. In the B2 partial vPCF, we observe short-range displacement correlations in the [001] direction, indicating the presence of residual polar distortions.

Sample Preparation: Details of the ceramic processing methods have been described previously.²⁰ TEM samples were cut from bulk ceramic pellets using a low speed diamond saw and mechanically wedge polished with diamond lapping films down to a grit of 0.1 µm until distinct optical thickness fringes were apparent. A Gatan PIPS II was used to further thin the samples and remove damage from mechanical polishing. Samples were liquid nitrogen cooled during ion milling; accelerating voltages of 1 kV, 500, and 100 V were used sequentially to achieve the thinnest possible sample with minimal residual damage.

STEM imaging was performed using an aberration-corrected ThermoFisher Titan S/TEM operated at 200 kV with the camera length set to give a HAADF collection angle of 70–200 mrad. Drift and scan distortion corrections were done using a non-linear drift correction code.¹² Sets of five images were collected in a serial acquisition with the scan orientation being rotated 90° between subsequent images in the set. Strong charging-related drift was present in the first acquisition of each set but thereafter reached equilibrium. While the code was capable of correcting this drift, it slightly reduced the useful area of the final image. Since including the first image was not necessary to correct distortion or improve the signal to noise ratio, only the remaining four images were used to produce the final drift-corrected image. Drift correction alone does not remove scan coil calibration errors so, to enable long-range correlations to be measured in BSSN, an image series of a silicon standard sample was also collected to measure the misalignment of the scan vectors.¹⁴ This was done during the same instrument session and with identical imaging parameters to those used for the BSSN images. A peak finding algorithm was applied to the drift-corrected silicon images by and the average crystallographic basis vectors were determined from the peak positions. The required transformation to recover the known silicon crystallographic basis was calculated as that needed to transform the basis vectors to be orthogonal and of equal magnitudes. This correction was then applied to the BSSN images. The final error in the scan vector distortion after correction of the silicon images was <10⁻³ for both the scaling and shear components.

Atomic column location and indexing was performed using a custom Python^TM code. First, the image was filtered using a Laplacian of Gaussian filter to differentiate closely spaced columns. Then, the watershed algorithm²³ was used to segment the image along minima between the atomic column peaks. This produced a segmentation map, labeling the pixels that belong to each intensity peak in the image. In order to simultaneously index the atomic columns and avoid fitting peaks that represent noise in the background intensity, a structural reference lattice was registered to the image. The reference lattice was oriented and scaled using the FFT and then fixed laterally by selecting an appropriate atom column in the image as the origin. Next, for each projected column in the reference lattice, the segmentation region with the closest center of mass was found and used to mask the original, unfiltered, distortion-corrected image. A 2D Gaussian function was then fit to the isolated peak by least squares minimization. The center of the Gaussian fit was taken to indicate the atom column location. The fitting parameters along with the reference lattice point coordinates and sublattice indexing label were all stored in a data frame format using the pandas²⁴ module. This data structure enabled many filtering, sorting, and plotting options for further analysis of the data, as well as the storage of additional data derivatives.

This material is based on the work supported by the National Science Foundation, as part of the center for Dielectrics and Piezoelectrics under Grant Nos. IIP-1841453 and IIP-1841466. This work was performed, in part, at the Analytical Instrumentation Facility (AIF) at North Carolina State University, which is supported by the State Of North Carolina and the National Science Foundation (Award No. ECCS-2025064). The AIF is a member of the North Carolina Research Triangle Nanotechnology Network (RTNN), a site in the National Nanotechnology Coordinated Infrastructure (NNCI). Work at the Molecular Foundry was supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The present work was partially supported by the National Natural Science Foundation of China under Grant No. 51790493.

There are no conflicts to declare.

S.D.F. and E.C.D. conceived the idea and developed the methodology. Z.J.Y. and X.M.C. synthesized the materials. M.J.C. prepared the TEM samples. S.D.F. collected the data with supporting contributions from M.J.C. S.D.F. designed, developed, and implemented the image analysis code, interpreted the experimental results, performed data visualization, and drafted the manuscript. C.O. advised on data collection and processing. E.C.D. and S.D.F. reviewed and edited the manuscript with supporting contributions from all other authors.

The data and code that support the findings of this study are openly available in the KiltHub Repository at http://doi.org/10.1184/R1/14318765.²²