Misﬁt-generated structural and optical anisotropies of the natural MoS 2 -PbS van der Waals heterostructure merelaniite

.


I. INTRODUCTION
The ability to control the structure of materials and thereby to tailor their properties has long been the goal of materials scientists and engineers. Following the discovery of graphene, vertically stacked van der Waals (vdW) heterostructures have provided a vast space for exploration and utilization of new chemistries and structures for materials design [1]. While tremendous advances have been made in fabricating such heterostructures amid challenges presented by issues of stability, naturally occurring minerals that are layered vdW heterostructures provide another route to studying such complex materials [2,3]. In particular, the structural, electronic, optical, optoelectronic, and electrocatalytic properties detectivity for infrared photosensor applications [13,[21][22][23][24][25][26]. For solar-cell applications, MoS 2 layers have been shown to be an effective hole-transport material in combination with light-absorbing PbS quantum dots [27]. Optoelectronic devices fabricated with nanoplate PbS grown on few-layer MoS 2 have demonstrated nonvolatile optical memory characteristics with infrared pulses [28]. High surface area MoS 2 clusters modified with PbS quantum dots have shown enhanced performance as a NO 2 detector over MoS 2 clusters alone [29].
Among naturally occurring vdW materials [2], merelaniite (Mo 4 Pb 4 VSbS 15 ) is the newest member of the family of misfit-layer sulfosalts of the cylindrite structure type, but is composed primarily of alternating MoS 2 and PbS layers [30]. The new mineral is found as fine cylindrical whiskers of circumferentially wound layers [ Fig. 1(a)] that can be associated with graphite, alabandite, zoisite (variety tanzanite), prehnite, quartz, calcite, and other minerals from the tanzanite gem mines near Merelani, Lelatema Mountains, Manyara Region, Tanzania. Rarely, it also occurs as nonwound thin sheets [Figs. 1(b)-1(c)], which are sometimes attached to cylindrical whiskers [ Fig. 1(d)]. To date, the Merelani tanzanite mines are the only confirmed locality for merelaniite.
Merelaniite is unique among related misfit-layer compounds in regard to its semicommensurability along the modulation direction. The H and Q lattice parameters along that direction (b H and b Q , respectively) have been observed [31,32] to be in the ratio of two successive integers (Vernier principle) such that b H /b Q ≈ (n + 1)/n ≡ n Q /n H , with corresponding ratios, for example, n Q /n H = 13/12 for cylindrite and lévyclaudite, 12/11 for Sn 2+ -rich franckeite, 16 In this study we present the crystal structure of merelaniite, as determined from single-crystal x-ray diffraction, demonstrating large modulations in bond distances and atomic positions along the modulation direction q, with consistent results from high-resolution scanning transmission electron microscopy (STEM) data. X-ray photoelectron spectroscopy (XPS) was used to study the oxidation states of the primary elements. The large in-plane structural anisotropy of merelaniite was probed by angle-resolved polarized Raman spectroscopy. We also demonstrate that merelaniite crystals can be mechanically exfoliated down to a thickness of 8 nm, corresponding to only around seven layer pairs of the heterostructure. By performing polarization-dependent thirdharmonic generation (THG) measurements, the anisotropic nonlinear optical response from the exfoliated merelaniite thin flakes is demonstrated, and the third-order nonlinear susceptibility of merelaniite crystal is estimated.

A. Crystal structure
The structure of merelaniite has been determined by single-crystal x-ray diffraction on the basis of the (3 + 2)dimensional superspace approach (see Ref. [3] for a review). This misfit-layer compound, of the cylindrite type, results from the combination of two heavily modulated C-centered triclinic H and Q subsystems (Fig. 2). The Q pseudotetragonal layer can be described as derived from the NaCl archetype, while the H pseudohexagonal layer was derived from the CdI 2 archetype. The modulations along q 11 are readily apparent, as illustrated in Fig. 2(a). Refined crystallographic parameters [33] are summarized in Table I. A thorough analysis of the reflection distribution revealed the two misfit subsystem cell orientations: the overall pattern can be interpreted as the combination of two heavily (i.e., with satellites up to third order) modulated triclinic subsystems with a common incommensurate q wave vector (q 11 = q 21 ) and only one shared reciprocal axis (a * 13 = a * 23 ). Thus, all the reflections were indexed with six integer indices, and the structure was refined starting from the atomic coordinates reported for lévyclaudite [32] using the program JANA2006 [34]. Details of the embedding of reflections in superspace can be found in Ref. [32]. Unfortunately, the number of independent reflections did not allow us to refine the occupancy of the site positions of the two layers. We therefore constrained the metal atom positions in the H and the Q subsystems as fully occupied by Mo and Pb, respectively, neglecting the minor but significant amounts [30] of V and Sb. In the final stage, the refinement converged to R = 0.0419 (wR = 0.0498) for 33 794 (6728 main reflections and 27 066 satellite reflections) independent reflections (3σ level) and 157 parameters.
Within the H layer, Mo atoms occupy octahedral sites [ Fig. 2 [35], and even closer to that observed in bulk 1T -MoS 2 [2.389(5) Å] by Fang et al. [19]. The Mo-S distances vary from 1.94 to 3.07 Å in the structure, which could be symptomatic of different valence states for Mo (see XPS results below)-and likely different coordinations-in different portions of the modulated structure, as is made evident by the calculation of the bond valence through the modulated structure [ Fig. 3(a)] as a function of the phase coordinate t along the fourth reciprocal superspace vector (whose real-space projection is q 11 ) [32,36]. Alternatively, the strong variation of the Mo-S distances could be related to a Mo-for-V substitution, as is indicated by the STEM data (see below).
In the Q layer, resembling a two-atom-thick slice of a NaCl-type structure, the Pb atoms are on the outside of the corrugated layer [ Fig. 2(c)] and may therefore make possible bonding interactions with the S atoms of the adjacent H layer. Each metal atom is then coordinated to five S atoms within the Q layer and to one or two S atoms in the neighboring H layer. The Pb-S bond distance is 2.869(5) Å on average, with variations of ±0.9 Å along the internal t superspace coordinate. It is remarkable to observe such a large dispersion of the bond distances as a function of t [ Fig. 3(b)], in accord with the huge displacement of the atoms perpendicular to the layers ( z of the Pb position ranging from −0.64 to 0.65 Å). It is noteworthy that a similar displacement of atoms was observed in lévyclaudite [32], with z of the Sn position ranging from ∼ −0.75 to 0.75 Å; however, this produced only a small dispersion of the Sn-S bond distances. Figure 4(a) shows a high-angle annular dark-field (HAADF) STEM image of a merelaniite crystal projected along the [100] zone axis. The contrast in HAADF STEM images is proportional to Z 2−δ , where Z is the atomic number and δ is an instrument factor that depends on the HAADF detector geometry [37,38]. In this image, the brightest contrast comes from Pb atoms followed by Mo atoms. Figure 4(b) shows a magnified section of Fig. 4(a), with a portion of the refined modulated structure model superimposed onto the STEM image projected along the same zone axis (with Pb, Mo, and S atoms shown in green, gray, and yellow spheres, respectively). Considering that there are strains induced by the bend that is visible at the top right of Fig. 4(a), and the fact that the structure is incommensurate, which leads to a lack of periodicity in real space, there is reasonably good agreement between the STEM image and the refined structure model. The differences visible in Fig. 4(b) may also be related to the fact that the two isolated single-crystal domains used for the x-ray and electron diffraction experiments, which were extracted from two different whiskers, may not be precisely representative of the same bulk. Indeed, natural materials (i.e., minerals) can show subtle chemical variations responsible for minor structural deviations and disorder, even within a single-crystal domain. Additionally, some shearing may possibly have taken place in the ultramicrotome preparation of the TEM foil, which could also account for some of the observed differences with the refined structure model, such as differences in some of the Pb/S atom positions normal to the stacking direction, within the Q layers shown in Fig. 4(b). A solution to the problem could be the application of a total scattering technique such as the atomic pair-distribution function analysis of powder-diffraction data (e.g., Ref. [39]). Unfortunately, the paucity of available material and the chemical inhomogeneity precluded such investigation.
To confirm the atom positions within the crystal structure, we further have performed energy-dispersive x-ray spec-trometry mapping at the atomic scale. Overall composition in atomic fraction (at. %) obtained from nine areas are as follows: Mo 17.5(±3.5), Pb 14.8(±1.3), V 4.0(±0.4), Sb 1.7(±0.3), and S 62.0(±4.9). Figure 5(a) shows an energy-dispersive spectroscopy (EDS) x-ray map where the HAADF image is superimposed with x-ray maps of Pb, Mo, and S, respectively, as well as the maps of these individual atoms. The corresponding line-scan plot is shown in Fig. 5 [33]). There, along with the similar trends with Pb, Mo, and S atoms, we have found that the V-atom positions anticorrelate with Pbatom positions, and partly coincide with the Mo-and S-atom positions.
Next, XPS was performed on merelaniite flakes to determine the surface chemical composition and the oxidation states of the elements. According to the XPS spectra, the main elements of Mo, Pb, and S are present, whereas V and Sb are not apparently observed. Figure 6 shows the high-resolution XPS spectra at the binding-energy regions of Mo 3d, Pb 4 f , and S 2p, where the peak fittings are used to determine the oxidation states of individual elements. Results indicate that lead is present in the material in Pb 2+ (88%) and Pb 4+ (12%) oxidation states, which is the same as that found for Pb in cylindrite [10] and very similar for lengenbachite (91% Pb 2+ and 9% Pb 4+ ) [11]. Molybdenum is present in the material in Mo 4+ (91%) and Mo 6+ (9%) oxidation states, while sulfur was found to be exclusively present as S 2− .

B. Studying structural anisotropy with Raman spectroscopy
As illustrated in the inset of Fig. 7(a), the Raman spectra from the top surface area of a bulk merelaniite cylinder (spot 1) and the connected unrolled flat area (spot 2) are collected and analyzed. It is found that the Raman spectra from these two spots are almost the same ( Fig. S2 [33]), indicating that the crystal structure remains the same as the merelaniite cylinder is unrolled into the flat flake.   [16,40], galena PbS [41], stibnite Sb 2 S 3 [42,43], and vanadium disulfide VS 2 [44,45]. The peaks at 122 and 142 cm −1 are assigned as the two-phonon process modes of VS 2 (121 and 150 cm −1 ). The 163-cm −1 peak represents a combination of the A g mode of Sb 2 S 3 (156 cm −1 ) and the two-phonon process mode of VS 2 (168 cm −1 ). The 189-cm −1 peak is assigned as the A g mode of Sb 2 S 3 (191 cm −1 ). The peak at 215 cm −1 corresponds to the longitudinal optical phonon modes of PbS (204 cm −1 ). The 230-cm −1 peak is assigned as the B 1g /B 3g mode of Sb 2 S 3 (238 cm −1 ). The 286-cm −1 peak is attributed to the A g mode of Sb 2 S 3 (283 cm −1 ), as well as the E 1g mode arising from the octahedral coordination of Mo atoms in 1T -MoS 2 [17,19,20,46]. We note that 1T -MoS 2 with distorted octahedral coordination also has three additional Raman bands (J 1 at 150-156 cm −1 , J 2 at 212-226 cm −1 , and J 3 at 324-333 cm −1 ) that are thought to arise from superlattice zone folding, and may contribute to respective Raman bands in merelaniite [17,20,46]. The 326-cm −1 peak represents the A 1g out-of-plane vibration mode of VS 2 (325 cm −1 ) and the peak at 354 cm −1 is attributed to one Raman mode of VS 2 (360 cm −1 ). The peaks at 383 and 392 cm −1 are attributed to the E 2g in-plane vibration mode of MoS 2 (384 cm −1 ). The prominent Raman peak at 404 cm −1 is assigned as the A 1g out-of-plane vibration mode of MoS 2 (409 cm −1 ). The 459-cm −1 peak is related to the first overtone of the longitudinal optical phonons of PbS (454 cm −1 ) and the second-order 2LA(M ) mode due to longitudinal acoustic phonons of MoS 2 (450 cm −1 ). Finally, the broad Raman shifts about 611 and 777 cm −1 are related to the first-order and second-order Raman modes due to the resonant Raman scattering in MoS 2 (572, 600, 643 cm −1 and 756, 782, 820 cm −1 , respectively).
Furthermore, angle-resolved polarized Raman spectroscopy is conducted to study the structural anisotropy and identify the crystal axes of merelaniite. Figure 7 As shown in Figs. 7(c)-7(n), the measured data agree well with the theoretical fits to Eq. (1). Similar to results for other naturally occurring incommensurate vdW heterostructures franckeite [9], cylindrite [10], and lengenbachite [11], the in-plane structural anisotropy leads to anisotropic twofold patterns in the polar plots of Raman intensity in parallel polarization configuration as a function of incident polarization angle for all of the examined modes (Fig. 7). As observed for cylindrite and franckeite, the intensities for merelaniite reach primary maxima at 0 • and 180 • along the ±x direction for all modes. Secondary maxima are pronounced along the ±y direction for all four of the examined Raman modes for cylindrite, whereas secondary maxima along the ±y direction are observed for only 2 of 12 studied modes in merelaniite, and are much less prominent than the primary maxima [Figs. 7(f) and 7(h)]. Secondary maxima are also absent or very weak along the ±y direction for franckeite. On the other hand, depending on the mode studied, primary maxima for lengenbachite can occur along either the ±xor ±y-directions, with secondary maxima evident for most of the studied modes.
Polarization-resolved Raman scattering in parallel polarization configuration has been used to infer the directions of the in-plane crystal axes and modulation directions in franckeite [9], cylindrite [10], and lengenbachite [11]. Since, in the present study, the x direction (θ = 0 • ) was experimentally set up to be perpendicular to the axial direction of the merelaniite cylinder, the above results indicate that the modulation direction q 11 of the merelaniite crystal is either along the cylinder axis (±y direction) or normal to it (±x direction) in the plane of the layers. However, since in the case of lengenbachite it was shown by complimentary TEM studies that the modulation direction can be along the direction of the primary maxima or along the direction of the secondary maxima, depending on the Raman mode examined [11], the orientation of the modulation direction relative to the cylinder axis for merelaniite is not determined unambiguously. Nevertheless, this is consistent with prior studies that have shown that in cylindrite, the modulation vector is parallel to the cylinder axis [48,49].
According to Eq. (1), the Raman intensity ratio Raman tensor |u| 2 /|v| 2 in the merelaniite crystal. It is found that the ratios of |u| 2 /|v| 2 are 1.265 for the 230-cm −1 mode and 1.424 for the 404-cm −1 mode from the top surface area of the cylinder (spot 1), whereas almost the same ratios of 1.297 and 1.487 are obtained for these two modes from the unrolled flat area (spot 2). This indicates that the structural anisotropy of the merelaniite crystal remains almost unchanged between the cylinder and the unrolled flake.

C. Mechanical exfoliation
The vdW-heterostructure nature of merelaniite crystal leads to weak interlayer interactions, which indicates the feasibility of obtaining thin flakes by mechanical exfoliation. Here, the bulk merelaniite crystal is mechanically exfoliated using Nitto tape (SPV 224) to get thin flakes with various thicknesses on a quartz substrate. Figure 8(a) shows the reflection optical microscope image of an isolated ultrathin merelaniite crystal on a quartz substrate. Figure 8(b) presents an atomic force microscope (AFM) image of the scanned region marked by the yellow dashed box in Fig. 8(a)

D. Anisotropic nonlinear optical response
The effect of the in-plane structural anisotropy on the nonlinear optical response of merelaniite crystal was further investigated by measuring the polarization-dependent anisotropic THG emission. Figure 9(a) plots the recorded THG spectrum from the 8-nm-thick merelaniite flake showing an expected peak emission at the wavelength of 520 nm, which is exactly at one-third of the fundamental wavelength of 1560 nm. The inset of Fig. 9(a) shows the transmission optical microscope image with the green color emission. The log-scale plot of the THG power dependence on the fundamental laser power in Fig. 9(b) shows a cubic power law with the slope of 2.97 ± 0.04, which further confirms the THG process. To determine the crystal axes of the probed exfoliated flakes, polarization-resolved Raman spectra are recorded. would be the cylinder axis, consistent with the orientation of the sample in Fig. 7.
Next, anisotropic THG responses from the exfoliated merelaniite flakes are characterized. Considering that merelaniite belongs to the triclinic crystal system, the contracted form of the third-order nonlinear susceptibility tensor χ (3) can be approximately expressed as [50] where the first subscript 1, 2, 3 represents x, y, z respectively, and the second subscript refers to the following combination of the three components as = E 0 (cos θx + sin θŷ), wherex and y are the unit vectors along the x-and y-axes. Since the excitation electric-field polarization always remains in the x-y plane in our experimental configuration, no contribution from the components relating to the z terms of the χ (3) tensor is expected in the THG signal. Therefore, the x-and y-polarized components of the resulting THG signal can be expressed as x ∝ (χ 11 cos 3 θ + χ 12 sin 3 θ + 3χ 18 cos θ sin 2 θ y ∝ (χ 21 cos 3 θ + χ 22 sin 3 θ + 3χ 28 cos θ sin 2 θ The effect of the incident linear polarization angle on the THG emission power from the merelaniite thin flakes is then characterized. The desired linear polarization of the pump beam is obtained by placing a linear polarizer oriented along the x axis and a rotating half-wave plate. Figure 9(d) plots the angular dependence of the THG emission power as a function of the excitation linear polarization angle with respect to the x axis for the 8-nm-thick flake. The black and blue data points are the recorded x-and y-polarized components of the THG power, while the red data points are the measured total THG power. The solid curves are the theoretical fits of the experimental data using Eqs.  Table II lists the extracted relative magnitudes of the elements of the χ (3) tensor for merelaniite, cylindrite, franckeite, and lengenbachite. It shows that χ 11 and χ 22 dominate for all these naturally occurring vdW materials. It is noteworthy that merelaniite has relatively larger values of χ 18 and χ 29 compared to cylindrite and lengenbachite, while for franckeite the values of χ 18 and χ 29 are much lower than all the others.
Furthermore, the dependence of THG emission as a function of merelaniite flake thickness is explored. Figure 9(i) plots the evolution of the measured average THG power P (3ω) as a function of the flake thickness d at a pump power of P (ω) = 3.6 mW. The black data points are the measured THG power from five different flakes of thicknesses, ranging from 8 to 121 nm, when the incident linear polarization is along θ = 0 • . The THG power increases quadratically up to around 2.11 pW as the flake thickness is increased up to 49 nm, giving the maximum conversion efficiency of 5.86 × 10 −10 , and then decays exponentially. The significant attenuation of THG power is attributed to the strong absorption through thick merelaniite flakes. The measured thickness dependence of the THG response can be fitted with an exponentially decaying function P (3ω) (d ) = Cd 2 exp(−2d/δ), where δ ≡ λ 3 2πk 3 , C is a constant, and k 3 is the imaginary part of the refractive index at λ 3 = 520 nm, which gives k 3 = 1.93. Finally, the effective scalar value of χ (3) for merelaniite crystal is estimated by using the following relationship [11]: where n 1 and n 3 are the real parts of the refractive index of merelaniite at the fundamental wavelength λ 1 = 1560 nm and the third-harmonic wavelength λ 3 , respectively. The parameters of the fundamental pump laser are given by spot size W = 1.5 μm, repetition rate f rep = 80 MHz, and pulse width τ = 90 fs. By assuming [11] the real part of the refractive index to be n 1 = n 3 = 3.5, the χ (3) eff value of merelaniite crystal is estimated to be 4.82 × 10 −20 m 2 /V 2 , which is just slightly smaller than those of the other naturally occurring vdW materials, as shown in Table II.

III. CONCLUSION
We have determined the crystal structure of merelaniite using single-crystal x-ray diffraction. Using the (3 + 2)-dimensional superspace approach, we find that merelaniite's structure is of the cylindrite type, being composed of two C-centered triclinic subsystems: pseudotetragonal (Q layer) of primarily PbS, and pseudohexagonal (H layer) that is primarily composed of MoS 2 . The layers show large misfit-induced bond-distance modulations and atomic displacements along a single common incommensurate q wave vector. Within the H layer, Mo atoms occupy octahedral sites. High-resolution scanning transmission electron microscopy shows reasonable agreement with the x-ray diffraction studies, Minor yet visible differences between x-ray and electron diffraction results might be related  to associated differences in sample preparation and the fact that the two isolated single-crystal domains used for the two experiments came from different whiskers and may not be precisely representative of the same bulk. High-resolution scanning transmission electron microscopy further indicates V atoms are concentrated in the H layers.
The in-plane structural anisotropy of merelaniite crystal is revealed through angle-resolved polarized Raman spectroscopy. We have also shown that merelaniite flakes can be mechanically exfoliated down to 7 layer-pair thickness. Anisotropic THG response from the exfoliated merelaniite thin flakes is further demonstrated with the extracted anisotropic χ (3) tensor elements, together with the estimated χ (3) eff value of merelaniite crystal, indicating that merelaniite exhibits strong nonlinear optical response similar to other anisotropic nonlinear vdW materials.
Yang et al. [51] have described nanoscale layered PbS inclusions 2-, 4-, and 6-atomic layers thick, epitactically aligned and parallel to MoS 2 layers within 2H 1 -polytype molybdenite from the Huanglongpu Mo-Pb ore deposit, Qinling orogenic belt, China. They suggest that the PbS layers most likely formed by the diffusion and subsequent exsolution of Pb in Pb-bearing molybdenite due to the different ionic radii of Pb ions as compared to S and Mo ions. While it is unlikely that merelaniite formed by such a mechanism, the observations of Yang et al., together with the cylindrite-group structure of merelaniite, suggest that it may be possible to fabricate a homologous series with different stacking sequences of MoS 2 and PbS layers, and further demonstrate an expanded structure-chemistry space for engineering stable layered materials for potential device applications.

A. Materials
Cylindrical merelaniite whiskers in this study were taken from the same specimen as holotype material [30]. They were associated with graphite and occurred on a large alabandite crystal from an unspecified mine in the tanzanite gem mines near Mererani (Merelani), Lelatema Mountains, Simanjiro District, Manyara Region, Tanzania. Separate whiskers were used for each of the different measurements, using the methods described below.

B. Single-crystal x-ray diffraction
A small merelaniite fragment (28 × 25 × 11 μm in size) from the holotype was mounted on a 5-μm-diameter carbon fiber, which was, in turn, attached to a glass rod. Singlecrystal x-ray diffraction intensity data were collected with a Bruker D8 Venture Photon 100 CMOS equipped with graphite-monochromatized Mo Kα radiation. The detectorto-crystal distance was 70 mm. Data were collected using ωand ϕ-scan modes, in 0.5 • slices, with an exposure time of 30 s per frame. The data were corrected for Lorentz and polarization factors and absorption using the software package APEX3 [52].

C. Scanning transmission electron microscopy
TEM foils were prepared by embedding a merelaniite whisker in Epofix embedding resin, curing at 60 • for 12 h, sectioning with a DiATOME Ultra 35 • diamond knife (feed 40 nm, sectioning speed 2 mm/s), and mounting on EMS C-flat TM holey carbon TEM grids. The atomic resolution imaging and x-ray spectroscopy maps were obtained using an FEI Titan Themis aberration-corrected microscope at Michigan Technological University. The imaging and spectroscopic measurements of the sample were performed in scanning mode of the TEM (STEM). The point resolution in STEM mode is 0.8 Å. The microscope was fitted with superX quad detectors for fast mapping with increasing capability to count x rays. For x-ray mapping, typically 25-μs dwell time and 100 picoamps of beam current were used. For both STEM and EDS mapping drift-correction modes were used. STEM and EDS data acquisition and image analysis were performed using Thermo Fisher VELOX.

D. XPS
The XPS spectra were obtained using the Thermo Scientific Nexsa XPS spectroscopy system with a monochromatic Al x-ray source.

E. Optical spectroscopy
The Raman spectra were collected by exciting the sample with a 632.8-nm He-Ne laser using an objective lens and recording the back-reflected light via the same objective lens in a spectrometer (Horiba, iHR 550). The linear polarization of the incident laser beam was adjusted by a linear polarizer and a half-wave plate. An edge filter (Semrock, LP02-633RE-25) was used to reject the laser light. Another linear polarizer in the collection path was used to resolve the parallel component of Raman spectra. The THG signal was collected from the sample excited with a femtosecond laser (wavelength 1560 nm, pulse width 90 fs, repetition rate 80 MHz) using an objective lens. The transmitted THG signal was collected by another objective lens. A short-pass filter was placed in the collection path to block the laser light. The THG signal was then recorded in a spectrometer and a charge-coupled device camera for the analysis.
samples were selected for this study. We thank Owen P. Mills (Michigan Technological University) and Stacie Kirsch (Electron Microscopy Sciences) for preparation of TEM foils.
L.B., X.Y., and J.A.J. conceptualized this study. X-ray diffraction experiments and analysis were performed by L.B.
Mechanical exfoliation, XPS, AFM analysis, and anisotropic Raman and third-harmonic-generation studies were carried out by A.D., J.G., and X.Y. STEM and EDS studies were conducted by P.M. All authors contributed to writing the manuscript.