Natural van der Waals heterostructure cylindrite with highly anisotropic optical responses

Abstract

The mechanical exfoliation of naturally occurring layered materials has emerged as an easy and effective method for achieving ultrathin van der Waals (vdW) heterostructures with well-defined lattice orientations of the constituent two-dimensional (2D) material layers. Cylindrite is one such naturally occurring vdW heterostructure, where the superlattice is composed of alternating stacks of SnS₂-like and PbS-like layers. Although the constituent 2D lattices are isotropic, inhomogeneous strain occurring from local atomic alignment for forcing the commensuration makes the cylindrite superlattice structurally anisotropic. Here, we demonstrate the highly anisotropic optical responses of cylindrite thin flakes induced by the anisotropic crystal structure, including angle-resolved polarized Raman scattering, linear dichroism, and polarization-dependent anisotropic third-harmonic generation. Our results provide a promising approach for identifying various natural vdW heterostructure-based 2D materials with tailored optical properties and can be harnessed for realizing anisotropic optical devices for on-chip photonic circuits and optical information processing.

Introduction

Since the discovery of graphene, van der Waals (vdW) materials are the subject of intensive research due to their fascinating optical, electrical, and mechanical properties. Hexagonal boron nitride, black phosphorus (BP), and transition metal dichalcogenides are a few examples among the other vdW materials, which have received great attention not only in the context of fundamental science but also for various practical applications^1,2,3. In these materials, individual layers of chemically bonded atoms are stacked and held by vdW interactions, and the nature of the vdW stacking is critical in defining their physical properties. One of the most promising prospects is to design vdW heterostructures with distinct stacking symmetries and functionalities by reassembling individual atomic layers of different materials in a precisely chosen layer-by-layer sequence, where the electronic band structures, carrier mobility, optical and magnetic properties can be tailored^4,5,6,7,8,9. The fabricated vdW heterostructures can enable intriguing functionalities beyond the reach of the existing materials and advance various photonic and optoelectronic applications covering a broad spectral range from visible to mid-IR for optical communication and sensing¹⁰, transistors^11,12, photodetectors^13,14,15, and ultrafast lasers¹⁶.

Currently, there are two ways of fabricating vdW heterostructures. One approach is to directly grow them through layer-by-layer sequential synthesis using chemical vapor deposition (CVD)^17,18,19. However, this is a challenging task and it has only been demonstrated for a few vdW heterostructures so far. Another approach is to mechanically exfoliate the individual layers of different materials and then deterministically place one layer on top of the other^20,21,22,23. Improper crystalline alignment between the stacked lattices and the presence of interlayer contaminants are two major issues for this method, which will affect the device performance. Alternatively, the mechanical exfoliation of naturally occurring vdW superlattices facilitates an easy and hassle-free approach for fabricating ultrathin layers of vdW heterostructures with well-defined lattice orientations of the individual layers. These materials are naturally formed with alternating layers of two different 2D materials due to the phase segregation process during their formation, and thus are free of the aforementioned issues. Several exfoliated naturally occurring vdW heterostructure materials such as cylindrite²⁴, franckeite²⁵, and levyclaudite²⁶ have been reported, which belong to the sulfosalt mineral family with the vdW superlattices formed by alternate stacking of SnS₂-like pseudo-hexagonal (H) layer and PbS-like pseudo-tetragonal (Q) layer. Out of these materials, the optical and electrical properties of bulk to single unit cell thick franckeite have been explored for optoelectronic applications in the linear optical regime^27,28. The high electrical conductivity and magnetic properties of cylindrite flakes have also been demonstrated²⁴. Although the stacked alternating layers are comprised of isotropic vdW materials, the incommensurability between the constituent SnS₂-like and PbS-like lattices introduces a strain-induced deformation in the crystals and hence induces the structural anisotropy showing anisotropic electrical, vibrational, and optical properties. The anisotropic optical properties of franckeite have recently been reported²⁷. However, to date, there is no study on the anisotropic optical responses of cylindrite yet.

The recent demonstration of broadband nonlinear optical response with ultrafast carrier dynamics exhibited in franckeite has opened up the possibility of developing nonlinear optoelectronic devices based on naturally occurring vdW heterostructures²⁹. Nonetheless, the effect of the crystal structural anisotropy on the nonlinear harmonic generation in franckeite and cylindrite has not yet been explored. The optical anisotropy in vdW materials exhibits polarization-sensitive nonlinear light-matter interactions at the nanoscale, which is important in the context of many applications in optical information processing and communication³⁰. The ultrathin thicknesses of the exfoliated vdW materials make them appealing for realizing high-performance integrated nonlinear optical devices^31,32. Previously, several types of vdW materials such as BP^33,34,35, ReS₂³⁶, GeSe³⁷, and GeAs³⁸ have been explored for anisotropic third-harmonic generation (THG). The widely-used anisotropic nonlinear 2D material of BP shows strong third-order optical nonlinearity, but BP is unstable under ambient conditions which limit its use in practical applications. Although other materials of ReS₂, GeSe, and GeAs are stable, they exhibit simple crystal structures formed by stacking only one type of molecular layer. Additionally, nonlinear harmonic generation is also regarded as a rapid and reliable characterization tool for evaluating material homogeneity and identifying crystal orientation and thickness of 2D materials^39,40.

Here, we study the highly anisotropic optical responses of mechanically exfoliated cylindrite thin flakes down to a few unit cell thickness. First, we perform transmission electron microscopy (TEM), energy-dispersive X-ray spectroscopy (EDXS), and X-ray photoelectron spectroscopy (XPS) analysis to determine the anisotropic crystal structure of cylindrite and its chemical composition. We also demonstrate that the in-plane structural anisotropy of cylindrite flakes can be revealed by angle-resolved polarized Raman spectroscopy. Furthermore, the optical anisotropy induced linear dichroism in cylindrite flakes is shown by performing polarization-dependent differential reflectance spectroscopy (DRS) measurements. Finally, highly anisotropic THG from cylindrite flakes is demonstrated by characterizing the angular dependence of THG emission with respect to the pump polarization. The third-order nonlinear susceptibility of cylindrite is also estimated to be in the same order of magnitude as other anisotropic nonlinear 2D materials. The demonstrated results are promising for identifying various natural vdW heterostructures to expand the current existing 2D material library. Our results can also be harnessed further to realize ultrathin anisotropic optical devices for polarization-sensitive photodetectors, optical modulators, frequency- and polarization-multiplexed encoding prototypes for secure optical information processing, and communication.

Results

Morphology and chemical composition of cylindrite

Cylindrite was first discovered by Frenzel from a tin-mining region in Bolivia in 1893. It attracted attention due to its peculiar triclinic crystal structure with a tendency to develop large concentric cylindrical shells and aggregates. Later, wet chemical analysis established that cylindrite belongs to the sulfosalt mineral family and has a generic chemical formula Pb₃Sn₄FeSb₂S₁₄⁴¹. Artificial synthesis of the cylindrite structures with user-tailored composition has also been achieved^42,43. Recently, liquid-phase and mechanical exfoliation of cylindrite thin flakes down to a few nanometers thick have been reported, and it is shown that cylindrite flakes exhibit a narrow bandgap <0.85 eV and are stable under ambient conditions²⁴. Figure 1a is a simplified representation of the cylindrite crystal structure, indicating alternative three-atom-thick SnS₂-like layer (H-layer) and two-atom-thick PbS-like layer (Q-layer) bonded together via vdW forces. The optical image in Fig. 1b shows the characteristic appearance of a cylindrite mineral rock containing clusters of long cylinders with gray-black color. Figure1c is a scanning electron microscope (SEM) image of the cross-section of an individual cylindrite cylinder, which shows that the crystals are formed by concentric layers of materials held by vdW interactions. The zoomed-in view in Fig. 1d illustrating a region around the revolution axis of the cylinder justifies its cylindrical shape. The following notation is adopted to describe the crystal structure of cylindrite, as [100] (a-axis) is considered to be the out-of-plane direction which is the direction of layer stacking, while [010] (b-axis) and [001] (c-axis) are considered for the two in-plane directions. Also, α, β, and γ are the angles between the a, b, and c axis. The unit cell parameters of the H-layer are $a_{{{\mathrm{H}}}} = 11.709$ $\mathop {{{{\mathrm{A}}}}}\limits^ \circ$, $b_{{{\mathrm{H}}}} = 3.670$$\mathop {{{{\mathrm{A}}}}}\limits^ \circ$, $c_{{{\mathrm{H}}}} = 6.320$$\mathop {{{{\mathrm{A}}}}}\limits^ \circ$, $\alpha _{{{\mathrm{H}}}} = 90^\circ$, $\beta _{{{\mathrm{H}}}} = 92.85^\circ$, and $\gamma _{{{\mathrm{H}}}} = 90.58^\circ$, while the lattice parameters for the Q-layer are $a_{{{\mathrm{Q}}}} = 11.733$$\mathop {{{{\mathrm{A}}}}}\limits^ \circ$, $b_{{{\mathrm{Q}}}} = 5.790$$\mathop {{{{\mathrm{A}}}}}\limits^ \circ$, $c_{{{\mathrm{Q}}}} = 5.810$$\mathop {{{{\mathrm{A}}}}}\limits^ \circ$, $\alpha _{{{\mathrm{Q}}}} = 90^\circ$, $\beta _{{{\mathrm{Q}}}} = 92.38^\circ$, and γ_Q = 93.87°⁴¹. Due to the small mismatch of the stacking vectors ($a_{{{\mathrm{H}}}}$ and $a_{{{\mathrm{Q}}}}$) and their corresponding $\beta$ and $\gamma$ values, the coincidence cell or a supercell contains 12 layer pairs with a = 140 $\mathop {{{{\mathrm{A}}}}}\limits^ \circ$. Furthermore, the lattice mismatch in the b and c directions indicates that the supercell is made of 30:19 and 12:13 unit cells of H- and Q-type unit cells in these two directions, respectively. The inhomogeneous elastic deformation and local atomic alignment caused to force the commensuration between incommensurate H- and Q-layers will result in the interlayer moiŕe patterns in vdW heterostructures^24,25,27,28. As the crystal relaxes to minimize the total energy, strong in-plane strain is developed to form small out-of-plane rippling modulated perpendicular to the moiŕe pattern, resulting in a strong in-plane structural anisotropy in the crystal. Figure 1a schematically illustrates the rippling in the cylindrite crystal lattice. The TEM image of a mechanically exfoliated cylindrite thin flake transferred onto a TEM copper grid (details of the TEM sample preparation is included in the Methods section) in Fig. 1e clearly shows an in-plane striped pattern of periodic darker and lighter areas arising from the rippling on the surface of cylindrite flake along the c-direction. The high-resolution (HR) TEM image in Fig. 1f shows the atomic arrangements of H- and Q-layers in a cylindrite flake. The fringes occurring from the rippling along the c-direction are spaced with a period of 3.6 nm. Although the constituent H- and Q-lattices are comprised of isotropic vdW materials, the incommensurability between them introduces a strain-induced deformation in the crystals and hence induces the structural anisotropy in the heterostructure superlattice. Therefore, the interlayer moiŕe patterns observed in the TEM images confirm the anisotropic heterostructure of cylindrite with strain-induced deformation. The selected area diffraction pattern (SAED) from the crystal normal to the [100] crystal zone axis is shown in Fig. 1g. The SAED exhibits the reflections from the corresponding H- and Q-layers, which are marked in magenta and yellow colors, respectively. Two distinct sets of reflection spots in rectangular arrangements indicate different b and c unit cell parameters for the constituent H- and Q-layers, which manifests the heterostructure superlattice is composed of two alternatively stacked incommensurate 2D lattices with strain-induced deformation.

Fig. 1: Morphology and crystal structure of cylindrite.

a Schematic diagram of the crystal structure of cylindrite composed of alternating SnS₂-like pseudo-hexagonal layer (H-layer) and PbS-like pseudo-tetragonal layer (Q-layer). b Optical image of the cylindrite mineral rock. The scale bar is 10 mm. c SEM image of the cross-section of an individual cylindrite cylinder showing the crystals are formed by concentric layers of materials held by vdW interactions. The scale bar is 100 µm. d Zoomed-in SEM image of a region around the revolution axis showing the concentric layers. The scale bar is 5 µm. e TEM image of a mechanically exfoliated ultrathin cylindrite flake showing the characteristic striped patterns along the c-direction of the crystal, caused by the out-of-plane rippling due to the forced commensuration of the constituent H- and Q-layers. The scale bar is 10 nm. f Representative HRTEM image showing the atomic arrangements. The scale bar is 2 nm. g Corresponding SAED pattern consisting with diffraction spots related to H-layer (magenta) and Q-layer (yellow). The scale bar is 2 nm⁻¹.

EDXS analysis is performed to analyze the chemical composition of the exfoliated cylindrite flakes. The averaged EDXS spectrum in Fig. 2a acquired from a thin cylindrite flake shows the presence of main elements Pb, Sn, Sb, Fe, and S. Other elements such as C, O, and Si are also detected, which is expected for naturally occurring minerals where other elements can also coexist. However, the signal of Cu is present due to the copper TEM grid. Most of the signal of C, O, and Si may also come from the underlying carbon film support of the TEM grid, and the accumulation of carbon-based adsorbates during the transfer of the flakes on the grid (see Supplementary Figure 1). Figure 2b–f plot the TEM-EDXS elemental maps of the main elements Pb, Sn, Sb, Fe, and S from the area of a cylindrite crystal shown in the bright-field (BF)-TEM image in the inset of Fig. 2a. These maps show that all the elements are homogeneously distributed in the crystal when viewed from the top along [100] direction of the crystal. The quantification of each element in the crystal is listed in Table 1, from which the approximate chemical formula of the cylindrite flake under study is determined to be Pb_2.9Sn_4.0Fe_0.9Sb_1.9S₁₄. This approximate chemical formula is not completely charge-balanced which is mainly caused by the low sensitivity of EDXS, the overlap between S K_α and Pb M_α peaks, as well as the overlap between Sn L_β and Sb L_α peaks.

Fig. 2: Bulk and surface chemical composition of cylindrite.

a Averaged EDXS spectrum from a thin cylindrite flake. The peaks of the five main elements Pb, Sn, Sb, Fe, and S labeled in orange color are used for the quantification in Table 1. b–f TEM-EDXS maps of the main elements Pb, Sn, Sb, Fe, and S from the area of a cylindrite crystal shown in the BF-TEM image in the inset of (a). The scale bar is 50 nm. g Averaged XPS spectrum of cylindrite obtained from multiple measurements. The peaks labeled in orange color are used for the quantification in Table 2. h–l High-resolution spectra around the binding energy regions of Pb 4 f, Sn 3d, Sb 3d, Fe 2p, and S 2p, respectively.

Table 1 EDXS quantification of cylindrite.

Surface chemical composition of cylindrite

The surface chemical composition including the effect of oxidation on the surface of air-aged cylindrite is explored using XPS since the penetration depth of the XPS technique is limited to a few nanometers. The cylindrite crystal is kept under ambient conditions for three weeks after cleaving it out of the cylindrite rock. Figure 2g shows the XPS spectrum of air-aged cylindrite, which is used for the quantification of the main elements in Table 2. The approximate surface chemical formula of cylindrite obtained from the compositional stoichiometry analysis is ${{{\mathrm{Pb}}}}_{3.5}{{{\mathrm{Sn}}}}_{4.0}{{{\mathrm{Fe}}}}_{0.9}{{{\mathrm{Sb}}}}_{1.9}{{{\mathrm{S}}}}_{14}$. In comparison to EDXS, a slightly larger concentration of Pb is found in the derived stoichiometry from the XPS survey spectra. Although the obtained stoichiometry is very close to the ideal formula of cylindrite, there is the possibility of slight degradation or shear deformation of the surface while performing the exfoliation. Further, the oxidation states of individual elements are determined by the peak fitting in the high-resolution spectra at Pb 4 f, Sn 3d, Sb 3d, Fe 2p, and S 2p binding energy regions, as presented in Fig. 2h-l. As shown in Fig. 2h, lead is present in the material as Pb²⁺ (88%) and Pb⁴⁺ (12%) oxidation states. Tin (Fig. 2i) is present as Sn²⁺ (76%) and Sn⁴⁺ (24%), which justifies the presence of tin in both Q- and H-layers. Antimony (Fig. 2j) is present as Sb³⁺ (87%) and Sb⁵⁺ (13%), whereas iron (Fig. 2k) is present as Fe²⁺. Figure 2j also shows a very weak O 1 s peak which suggests there is only weak oxidation at the surface so that cylindrite is stable under ambient conditions. Sulfur is exclusively present in S²⁻ state as shown in Fig. 2l. The slight underrepresentation of Sb and Fe in the derived formula may also be attributed to the overlap between the Sb 3d peak with O 1 s peak, and Fe 2p peak with Sn 3p peak. The comparison between the XPS spectra of the pristine cylindrite crystal and the air-aged cylindrite crystal is shown in Supplementary Figure 2. It is observed that the XPS survey spectra of the pristine and air-aged crystals are just about identical and the surface stoichiometry of cylindrite remains almost intact. The amount of oxygen is further quantified by the peak fitting in the high-resolution spectra at the Sb 3d and O 1 s binding energy region. Over three weeks of air-aging, the atomic concentration ratio of Sb and O is only changed from 1:0.13 of the pristine crystal into 1:0.16. Therefore, it can be inferred that the effect of surface oxidation is very weak in cylindrite.

Table 2 XPS quantification of cylindrite surface.

Revealing structural anisotropy of cylindrite crystal

The cylindrite thin flakes are mechanically exfoliated using Nitto tape onto a quartz substrate (details in the Methods section). Figure 3a shows a reflection optical microscope image of an exfoliated thin cylindrite crystal. The topography of the crystal surface is analyzed using atomic force microscopy (AFM) in Fig. 3b. From the height profile of the AFM image (inset of Fig. 3b), the thickness is determined to be 10 nm approximately which indicates that the flake consists of only 9 layer pairs. First, the in-plane structural anisotropy of the crystal is studied by angle-resolved polarized Raman spectroscopy using a linearly polarized 632.8 nm He–Ne laser. Details of the experimental setup are described in the Methods section. Figure 3c plots the parallel polarization component of the Raman spectra as a function of the linear polarization angle of the excitation laser. The Raman spectrum shows six distinct bands in the 50–350 cm⁻¹ spectral region. The observed Raman peaks are located at 77, 90, 141, 182, 245, and 312 cm⁻¹. The Raman modes of cylindrite can be assigned according to the previously studied Raman spectra of the constituent components of PbS⁴⁴, SnS^45,46, SnS₂^47,48, FeS₂⁴⁹, and Sb₂S₃^50,51. The Raman peak at 77 cm⁻¹ is attributed to the A_g mode of Sb₂S₃ (73 cm⁻¹). The peak at 90 cm⁻¹ is assigned as the A_g mode of SnS (95 cm⁻¹). The 141 cm⁻¹ peak corresponds to a superposition of phonon modes of both SnS₂ lattices (the second-order effects at 140 cm⁻¹) and PbS lattices (transverse acoustic and optical modes at 154 cm⁻¹). The peak at 182 cm⁻¹ is most likely related to a combination of the A_g mode of SnS (192 cm⁻¹) and the A_g mode of Sb₂S₃ (191 cm⁻¹). The peak at 245 cm⁻¹ represents the vibrational mode of Sb₂S₃ (238 cm⁻¹). In addition, the peak at 312 cm⁻¹ is assigned as a combination of the A_g mode of Sb₂S₃ (312 cm⁻¹) and the A_1g mode of SnS₂ (315 cm⁻¹). It is noted that the Raman spectra of cylindrite are dominated by the contributions from Sb₂S₃, SnS₂, SnS, and PbS. The low energy and broad width of the Raman peaks are attributed to the superposition of the closely placed modes in the constituent lattices and the heavy damping present in the vdW heterostructure. Raman spectra of cylindrite flakes with other thicknesses are shown in Supplementary Figure 3. Almost no Raman shift is observed depending on the flake thickness due to the vdW stacking of the H- and Q-layers. From Fig. 3c, the intensity of each Raman mode in the parallel polarization configuration exhibits a strong dependence on the linear polarization angle of the incident beam. Such angular dependence of Raman intensity can be used to reveal the crystal structural anisotropy of cylindrite flakes and determine the crystal axes. Although the alternatively stacked SnS₂-like and PbS-like layers are isotropic, the incommensurability between these constituent 2D lattices introduces the strain-induced deformation in the crystals, resulting in the anisotropic heterostructure superlattice. The Raman analysis reveals the structural anisotropy of such heterostructure due to the strain-induced deformation. Figure 3d–g present the polar plots of the intensity variations of parallel polarization components for the Raman modes at 141, 182, 245, and 312 cm⁻¹ as a function of the incident linear polarization angle $\theta$. The triclinic crystal of cylindrite only supports A_g modes, where the Raman tensor of a given A_g mode is defined as $R = \left( {\begin{array}{*{20}{c}} u & w \\ w & v \end{array}} \right)$, where u, v, and w are the Raman tensor components. The parallel polarization component of the Raman A_g modes can be expressed as⁵²,

$$I_{{{{\mathrm{A}}}}_{{{\mathrm{g}}}}}^{||}\left( \theta \right) \propto \left( {u\cos ^2\theta + 2w\sin \theta \cos \theta + v\sin ^2\theta } \right)^2$$

(1)

Fig. 3: In-plane structural anisotropy measured through Raman spectroscopy.

a Reflection optical microscope image of an exfoliated cylindrite flake. x- and y-axes signify reference axes. b AFM image of the cylindrite flake with the inset showing that the flake thickness is around 10 nm. The scale bars are 2 µm. c Raman spectra acquired in the parallel polarization configuration from the same flake for different linear polarization angles between 0° and 180° with respect to the x-axis. d–g Polar plots of the intensity variations of parallel polarization components for Raman modes at 141, 182, 245, and 312 cm⁻¹. The data points are the experimentally measured values and the solid curves are the corresponding theoretical fits.

It is shown that the experimental measurements match well with the theoretical fits from Eq. (1). The anisotropic two-fold rotational symmetry of intensity is observed for each Raman mode with two unequal maxima as the incident linear polarization is along $0^\circ$ (x-direction) and $90^\circ$ (y-direction), where the x- and y-directions corresponding to the crystal orientations are shown in Fig. 2a. The higher maximum at $\theta = 0^\circ$ indicates that the x-direction corresponds to the rippling direction along the c-axis in the cylindrite crystal, whereas the y-direction is along the b-axis. Therefore, it can be inferred that the angle-resolved Raman mode analysis is a useful tool to show the crystal structural anisotropy and determine the crystal lattice orientations in cylindrite flakes.

Linear dichroism in cylindrite thin flake

To further explore the effect of the strain-induced lattice deformation on the anisotropic optical properties of the material, polarization-dependent DRS is performed on the cylindrite thin flakes, which is defined as,

$${{{\mathrm{DRS}}}} = \frac{{R - R_0}}{{R_0}}$$

(2)

where the change in the reflectance $R$ from the thin flake is measured with respect to the reflectance $R_0$ from the bare quartz substrate and normalized by $R_0$. By definition, the magnitude of DRS is proportional to the absorption of the sample, and it can be easily linearized in the case of ultrathin flakes with thickness d << λ, with $\lambda$ being the wavelength of incident light⁵³. For the normal incidence of light and transparent substrate such as quartz, DRS can be expressed as,

$${{{\mathrm{DRS}}}} = \frac{{8{\uppi}d}}{\lambda } \times \frac{{\epsilon_{{{{\mathrm{flake}}}}}^{\prime\prime} }}{{{\it{\epsilon }}_{{{{\mathrm{substrate}}}}}^\prime - 1}}$$

(3)

where ${\it{\epsilon }} = {\it{\epsilon ^{\prime}}} + {{{\mathrm{i}}}}{\it{\epsilon ^{\prime\prime}}} = (n + {{{\mathrm{i}}}}k)^2$ is the complex dielectric constant of the material with ${\it{\epsilon ^{\prime}}}$ and ${\it{\epsilon ^{\prime\prime}}}$ denoting the real and imaginary parts respectively, while $(n + {{{\mathrm{i}}}}k)$ represents the complex refractive index. Therefore, DRS is directly proportional to the imaginary part $\epsilon_{{{{\mathrm{flake}}}}}^{\prime\prime}$ of the cylindrite flake. Figure 4a plots the acquired DRS spectra from the same 10 nm thick cylindrite flake in Fig. 3a for different incident linear polarization angles with respect to the rippling direction (x-axis) of the crystal in the excitation wavelength range of 450–750 nm. Details of the experimental setup are described in the Methods section. The polar plots in Fig. 4b and c show the variations of the measured DRS and the extracted $\epsilon_{{{{\mathrm{flake}}}}}^{\prime\prime}$ as a function of the incident linear polarization angle at three different wavelengths of 520, 600, and 720 nm. A dipolar dependence is observed in all the cases. The DRS and $\epsilon_{{{{\mathrm{flake}}}}}^{\prime\prime}$ data are further fitted with a sinusoidal function of the form α_xcos² θ + α_ysin² θ, where $\alpha _x$ and $\alpha _y$ are constants proportional to the absorbance values as the incident white light is linearly polarized along x- and y-axes respectively. It is indicated that the absorbance reaches the maximum when the incident linear polarization is along the rippling direction of the crystal (along the c-axis), while it is the minimum when the incident polarization is perpendicular to the rippling direction (along the b-axis). Therefore, it is evident that strong in-plane structural anisotropy resulted from the lattice deformation in cylindrite crystal gives rise to linear dichroism, even though the individual H- and Q-layers are not expected to exhibit any in-plane anisotropy. In Fig. 4d, the linear dichroism is quantified by plotting the variation of the ratio $\epsilon^{\prime\prime}_{{{{\mathrm{flake,}}}}x} /\epsilon^{\prime\prime}_{{{{\mathrm{flake,}}}}y}$ between the values of $\epsilon_{{{{\mathrm{flake}}}}}^{\prime\prime}$ as the incident polarization is along the rippling direction of the crystal and its perpendicular direction. It is observed that the anisotropic linear optical response increases gradually with the excitation wavelength. Figure 4e plots the measured DRS spectra of four individual cylindrite flakes of different thicknesses 10, 45, 65, and 80 nm, when the incident white light is polarized along the rippling direction of each crystal. The crystal axes of each crystal are determined through angle-resolved polarized Raman spectroscopy measurements, which are shown in Supplementary Figure 4. A higher magnitude in DRS for a thicker crystal indicates that the absorption is more for a thicker sample. In addition, the redshift in the DRS spectrum as a function of thickness justifies a more yellowish glow in the reflection optical image of the thicker crystal compared to the greenish-yellow glow of a thinner one.

Fig. 4: Linear dichroism of cylindrite.

a DRS spectra from the 10 nm thick cylindrite flake as a function of the incident polarization angle between $0^\circ$ and $180^\circ$. b, c Evolution of the measured DRS and the extracted $\epsilon_{{{{\mathrm{flake}}}}}^{\prime\prime}$ as a function of the incident polarization angle for three different wavelengths of 520 nm (dark yellow), 600 nm (orange), and 700 nm (magenta). The measured data points are fitted with a sinusoidal function (solid curves). d $\epsilon^{\prime\prime}_{{{{\mathrm{flake,}}}}x} /\epsilon^{\prime\prime}_{{{{\mathrm{flake,}}}}y}$ response as a function of the excitation wavelength. e DRS spectra from cylindrite flakes of four different thicknesses 10, 45, 65, and 80 nm.

Anisotropic THG response from cylindrite flakes

The effect of in-plane structural anisotropy on the nonlinear optical properties of cylindrite crystal is investigated through character polarization-dependent anisotropic THG response. The details of the experimental setup are provided in the Methods section. When the c, b, and a axes of the triclinic cylindrite crystal is oriented along x, y, and z directions, the contracted form of its third-order nonlinear susceptibility tensor (χ⁽³⁾) can be expressed as⁵⁴,

$$\chi ^{\left( 3 \right)} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\chi _{11}} & {\chi _{12}} & {\chi _{13}} \\ {\chi _{21}} & {\chi _{22}} & {\chi _{23}} \\ {\chi _{31}} & {\chi _{32}} & {\chi _{33}} \end{array}} & {\begin{array}{*{20}{c}} {\chi _{14}} & {\chi _{15}} & {\chi _{16}} \\ {\chi _{24}} & {\chi _{25}} & {\chi _{26}} \\ {\chi _{34}} & {\chi _{35}} & {\chi _{36}} \end{array}} & {\begin{array}{*{20}{c}} {\chi _{17}} & {\chi _{18}} & {\chi _{19}} \\ {\chi _{27}} & {\chi _{28}} & {\chi _{29}} \\ {\chi _{37}} & {\chi _{38}} & {\chi _{39}} \end{array}} \end{array}} & {\begin{array}{*{20}{c}} {\chi _{10}} \\ {\chi _{20}} \\ {\chi _{30}} \end{array}} \end{array}} \right]$$

(4)

where the first subscript 1, 2, 3 represents x, y, z respectively, and the second subscript refers to the following combination of three components as,

$$\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {xxx} & {yyy} & {zzz} \\ 1 & 2 & 3 \end{array}} & {\begin{array}{*{20}{c}} {yzz} & {yyz} & {xzz} \\ 4 & 5 & 6 \end{array}} & {\begin{array}{*{20}{c}} {xxz} & {xyy} & {xxy} \\ 7 & 8 & 9 \end{array}} \end{array}} & {\begin{array}{*{20}{c}} {xyz} \\ 0 \end{array}} \end{array}$$

The electric field of the pump beam at frequency $\omega$ which is linearly polarized at an angle $\theta$ with respect to the x-axis can be written as ${{{\vec{\mathbf E}}}}^{\left( \omega \right)} = E_0\left( {\cos \theta \hat x + \sin \theta \hat y} \right)$, where $\hat x$ and $\hat y$ are the unit vectors along the x- and y-axis. Since the polarization of the excitation electric field always remains in the x-y plane in the current experimental configuration, no contribution from the components of χ⁽³⁾ containing z terms will be observed in the THG signal, so that there are only eight nonzero elements of χ⁽³⁾ as χ₁₁, χ₁₂, χ₁₈, χ₁₉, χ₂₁, χ₂₂, χ₂₈ and χ₂₉. Since the electric field of THG is proportional to χ⁽³⁾, the x- and y-polarized components of the THG output intensity can be expressed as,

$$\begin{array}{l}I_x^{3(\omega )} \propto \left( {\chi _{11}\,{{{\mathrm{cos}}}}^3\theta + \chi _{12}\,{{{\mathrm{sin}}}}^3\theta + 3\chi _{18}\cos \theta \,{{{\mathrm{sin}}}}^2\theta + 3\chi _{19}\sin \theta \,{{{\mathrm{cos}}}}^2\theta } \right)^2\\ I_y^{3(\omega )} \propto \left( {\chi _{21}\,{{{\mathrm{cos}}}}^3\theta + \chi _{22}\,{{{\mathrm{sin}}}}^3\theta + 3\chi _{28}\cos \theta \,{{{\mathrm{sin}}}}^2\theta + 3\chi _{29}\sin \theta \,{{{\mathrm{cos}}}}^2\theta } \right)^2\end{array}$$

(5)

The inset of Fig. 5a shows the transmission optical microscope image of the same 10 nm thick flake of Fig. 3a, where the green THG signal from the flaking area is observed as it is excited by a 1560 nm pump laser with the spot size of 1.5 µm. The THG spectrum in Fig. 5a shows the expected peak emission wavelength at 520 nm, which is exactly one-third of the fundamental wavelength. Figure 5b shows the log-scale plot of the THG power as a function of the pump power, where the cubic power-law dependence further confirms the THG process. For the average pump power of 4 mW corresponding to a peak pump irradiance of 27.2 GW cm⁻², the obtained THG conversion efficiency is around $2 \times 10^{ - 9}$ for the 10 nm thick flake. Next, the effect of the incident linear polarization angle on the THG emission power from the cylindrite crystal is characterized. The desired linear polarization of the pump beam is obtained by placing a linear polarizer oriented along the x-axis (parallel to the rippling direction of the flake) and a rotating half-wave plate. Figure 5c–f plot the angular dependence of the THG emission power as a function of incident linear polarization angle with respect to the x-axis for four cylindrite flakes with thicknesses of 10, 45, 80, and 150 nm. The black and magenta data points are the measured x- and y-polarized components of the THG power, while the dark cyan data points show the measured total THG power. The solid lines in Fig. 5c are theoretical fits of the experimental data using Eq. (5), showing a good agreement between the theoretical prediction and the experimental measurements. The optical images, the angle-resolved polarized Raman characterization for determining the crystal orientation, and the corresponding anisotropic THG characterization for the cylindrite flakes of different thicknesses are included in Supplementary Figure 4. It is evident that the polarization-dependent THG emission from the cylindrite crystal exhibits an anisotropic two-fold rotational symmetry, with the maximum THG emission occurring at the incident linear polarization along 0° (parallel to the rippling direction along the c-axis) and the second maximum THG emission along 90° (perpendicular to the rippling direction along the b-axis). The evolution of the anisotropy ratio of the THG signal $I_x^{(3\omega )}(\theta = 0)/I_y^{(3\omega )}(\theta = 90^ \circ )$ between the incident polarizations parallel and perpendicular to the rippling direction as a function of the flake thickness is plotted in Fig. 5g. It is shown that for all the cylindrite flakes, the anisotropy ratio almost keeps as a constant with an average value of 2.33 ± 0.19. The structural anisotropy in the cylindrite crystal is further reflected in the extracted relative magnitudes of the χ⁽³⁾ tensor. The relative magnitudes of the anisotropic χ⁽³⁾ tensor elements can be extracted from Eq. (5). It is found that the values of χ⁽³⁾tensor elements stay almost unchanged for cylindrite flakes of different thicknesses, with the average relative magnitudes as $\chi _{11}:\chi _{12}:\chi _{18}:\chi _{19}:\chi _{21}:\chi _{22}:\chi _{28}:\chi _{29} =$ 1 : 0.012 : 0.121 : 0.016 : 0.003 : 0.663 : 0.002 : 0.107, which further reveals the intrinsic nature of the anisotropic nonlinear optical response of the cylindrite crystal.

Fig. 5: Anisotropic THG response from cylindrite flakes.

a The measured spectrum of THG emission with a peak wavelength at 520 nm from the 10 nm thick cylindrite flake, which is exactly one-third of the fundamental wavelength of the pump beam at 1560 nm. Inset shows a transmission microscope image of THG emission from the flake. The scale bar is 2 µm. b The double log-scale plot of the measured THG emission power as a function of the pump power. c–f Angular dependence of the THG emission power as a function of incident linear polarization angle of the pump beam for four cylindrite flakes with the thicknesses of 10, 45, 80, and 150 nm. 0° corresponds to the x-axis along the rippling direction of the crystal. g Measured THG anisotropy ratio as a function of the flake thickness. h Evolution of the average THG emission power versus the flake thickness for incident linear polarization along the x-axis (black) and y-axis (magenta).

The air-stability of the 10 nm thick cylindrite flake is investigated by Raman spectroscopy and THG characterization after keeping the flake under ambient conditions for four weeks. The acquired Raman spectrum from the air-aged cylindrite flake is plotted in Supplementary Figure 5, which is almost exactly the same as that of the freshly exfoliated flake in Fig. 3c, indicating all the Raman modes are observed to be intact after air-aging. The comparison between the THG characterization of the freshly exfoliated flake and that of the air-aged one is presented in Supplementary Figure 6. It is observed that the THG emission power and THG conversion efficiency remain at the same level as the freshly exfoliated flake, even after four weeks of air-aging. In addition, the shape of the polarization-dependent anisotropic THG emission pattern remains intact after air-aging, with almost the same anisotropy ratio of the THG signal. Hence, it can be inferred that the exfoliated cylindrite flakes are stable in air.

Moreover, the dependence of the THG response as a function of the cylindrite flake thickness is investigated. It is possible to extract the magnitude of χ⁽³⁾ from the following equation for the average power of THG emission $P^{\left( {3\omega } \right)}$, which is derived by solving nonlinear Maxwell’s equations as³⁴,

$$P^{(3\omega )}\left( d \right) = \frac{{9\omega ^2d^2}}{{16\sqrt {n_3^2 + k_3^2} n_1^3{\it{\epsilon }}_0^2c^4}}\left| {\chi ^{(3)}} \right|^2\frac{{P^{(\omega )3}}}{{f_{{{{\mathrm{rep}}}}}^2W^4\tau ^2\left[ {\frac{{\uppi }}{{4{{{\mathrm{ln}}}}2}}} \right]^3}}\left( {\frac{{e^{ - \frac{{4{\uppi}k_3d}}{{\lambda _3}}} - 2\cos \left( {{{\Delta }}kd} \right)e^{ - \frac{{2{\uppi}k_3d}}{{\lambda _3}}} + 1}}{{d^2\left( {\frac{{4{\uppi}^2k_3^2}}{{\lambda _3^2}} + {{\Delta }}k^2} \right)}}} \right)e^{ - \frac{{4{\uppi}k_3d}}{{\lambda _3}}}$$

(6)

where $d$ is the thickness of the cylindrite flake, $n_1$ and $n_3$ are the real part of the refractive index at the fundamental frequency ($\omega$) and THG frequency ($3\omega$), $k_3$ is the imaginary part of the refractive index at $3\omega$, and ${{\Delta }}k = \frac{{6{\uppi}}}{{\lambda _1}}\left( {n_1 - n_3} \right)$ is the phase mismatch between the forward propagating fundamental beam and the third-harmonic emission beam in the transmission microscope arrangement. The parameters of the Gaussian fundamental pulsed laser include average pump power $P^{\left( \omega \right)}$, beam spot size $W$, pulse width $\tau$, and repetition rate $f_{{{{\mathrm{rep}}}}}$.

The black data points in Fig. 5h show the measured average THG power as a function of the cylindrite flake thickness at the pump power of 1.5 mW when the incident linear polarization is parallel to the rippling direction (x-axis). It is clearly shown that within the range of thickness less than 60 nm, the THG power increases quadratically up to around 8 pW as the flake gets thicker, giving the maximum conversion efficiency of 5.33 × 10⁻⁹. However, for the thickness greater than 60 nm, exponential decay of the THG signal is observed. This is attributed to the significant attenuation of the THG signal as it propagates through the flake due to the strong absorption of the THG signal in the cylindrite flake especially with large thickness. Therefore, by fitting the measured data with an exponentially decaying function (black curve) of the form $P^{\left( {3\omega } \right)}\left( d \right) = Cd^2\exp \left( { - \frac{{4{\uppi}k_3d}}{{\lambda _3}}} \right)$ with C as a constant, the imaginary part of refractive index at $\lambda _3$ = 520 nm for the incident linear polarization along the x-axis is extracted as $k_{3x} = 1.45$. The real part of the refractive index along the x-axis can then be estimated from the measured $\epsilon^{\prime\prime}_{{{{\mathrm{flake,}}}}x} = 2n_{3x}k_{3x}$ (Fig. 4c) as $n_{3x} = 3.12$. These extracted refractive index values for cylindrite are consistent with the measured ones for franckeite. Since cylindrite is having a similar chemical composition and layered structure as franckeite, it is expected to have a refractive index within the same range of franckeite. Next, the $\chi _{11}^{\left( 3 \right)}$ corresponding to the magnitude of χ⁽³⁾ for the x-polarized pump beam is extracted from the fitting in Eq. (6) by using the estimated values of $n_{3x}$ and $k_{3x}$. Taking into account the pump laser parameters of $P^{\left( \omega \right)}$ = 1.5 mW, $\tau$ = 90 fs, $f_{{{{\mathrm{rep}}}}}$ = 80 MHz, and W = 1.5 µm at $\lambda _1$ = 1560 nm with $n_1 = 4$, the magnitude of $\chi _{11}^{\left( 3 \right)}$ is estimated to be $3.06 \times 10^{ - 19}$ m² V⁻². In Fig. 5h, the magenta data points are the evolution of the measured THG power with flake thickness for the incident linear polarization perpendicular to the rippling direction (y-axis). Again, by fitting the THG power data with an exponentially decaying function and using the measured $\epsilon^{\prime\prime}_{{{{\mathrm{flake,}}}}y} = 2n_{3y}k_{3y}$, the refractive index values $k_{3y}$ and $n_{3y}$ along y-axis can be estimated as 1.38 and 2.79, respectively. The extracted magnitude of $\chi _{22}^{\left( 3 \right)}$ representing $\chi ^{(3)}$ for y-polarized pump beam is $2.33 \times 10^{ - 19}$ m² V⁻², which gives the $\chi _{11}^{\left( 3 \right)}/\chi _{22}^{\left( 3 \right)}$ ratio of 1.31 for the cylindrite crystal. Finally, to assess the THG response in cylindrite, the measured THG anisotropy ratio and χ⁽³⁾ value of the cylindrite crystal are compared with the already reported values of other existing anisotropic nonlinear 2D materials including BP, ReS₂, GeSe, GeAs and SiP, which are tabulated in Table 3. It shows that cylindrite has a strong nonlinear optical response with the third-order nonlinear susceptibility almost double that of BP.

Table 3 Comparison of THG response in cylindrite with other anisotropic nonlinear 2D materials.

Discussion

In conclusion, we have demonstrated in-plane structural anisotropy and highly anisotropic optical responses of the cylindrite crystal, where two individual SnS₂-like and PbS-like isotropic layers are stacked to create the natural vdW heterostructure. The vdW interaction, local atomic alignment, and lattice deformation while stacking of the constituent SnS₂-like and PbS-like layers create small out-of-plane rippling on the crystal surface, resulting in strong in-plane structural anisotropy in the crystal. It is shown that the direction of the ripples and hence the in-plane structural anisotropy in the crystal can be probed by angle-resolved polarized Raman spectroscopy. The effect of the strong structural anisotropy on the linear and nonlinear optical properties of the cylindrite crystal is further studied. Through polarization-dependent DRS measurements, it is revealed that cylindrite thin flakes exhibit strong linear dichroism, which is not only important in the context of understanding the origin of optical anisotropy in natural vdW heterostructures but also can be harnessed to design future polarization-sensitive photodetectors and transistors. Moreover, polarization-dependent anisotropic THG responses from cylindrite thin flakes is investigated. It is shown that the THG emission power as the pump beam is linearly polarized along the rippling direction of the crystal (c-axis) is more than two times compared to that along the perpendicular direction (b-axis). The relative magnitudes of the eight nonzero components of the χ⁽³⁾ tensor are extracted as an estimate of the nonlinear optical anisotropy in the crystal. Lastly, it is shown that the THG emission power is highly dependent on the flake thickness, where the THG response is maximum around the flake thickness of 60 nm. The estimated χ⁽³⁾ of cylindrite crystal is found to be $3.06 \times 10^{ - 19}$ m² V⁻², indicating that cylindrite has a strong nonlinear optical response as other anisotropic nonlinear 2D materials such as BP. Additionally, unlike BP, cylindrite is stable in ambient conditions, providing an added advantage for practical use in building future on-chip anisotropic nonlinear optical devices for optical information processing and communication.

Methods

Sample preparation

The quartz substrate is sonicated in deionized water, acetone, and isopropanol one after the other. The cylindrite thin flakes are mechanically exfoliated from the bulk natural cylindrite mineral (from Poopó town, Oruro, Bolivia) using Nitto tape (SPV 224). Then the flakes are transferred directly onto the quartz substrate by placing the tape containing the flakes in contact to the quartz surface and peeling it off fast. Now, flakes of different thicknesses on the quartz substrate can be readily identified by monitoring their optical transparency under a transmission optical microscope. Flakes with lower thickness exhibit higher optical transparency. Finally, the thickness of each identified flake is confirmed by AFM measurement.

For the TEM measurement, exfoliated cylindrite flakes are transferred from the quartz substrate to TEM copper grids following polymethyl methacrylate (PMMA)-assisted wet transfer method. First, the quartz substrate containing the thin flakes is spin-coated (2000 rpm, 60 s) with PMMA (950 kDa) and baked at 120 °C for 2 min to facilitate the adhesion between cylindrite and PMMA. Next, the sample is immersed in 3 M potassium chloride (KOH) solution and kept at 50 °C for 1.5–2 h. After quartz etching by KOH, the PMMA film containing the cylindrite crystals is washed with deionized (DI) water two times for 10 min at each step to remove the residual KOH. Finally, the PMMA film is fished out of the DI water on a TEM grid. The TEM grid covered with the PMMA film is left uncovered for drying out naturally. Finally, the PMMA is washed out by dissolving it in acetone while the crystals are transferred on the TEM grid.

Experimental setup

For analyzing the Raman spectrum, a 632.8 nm He–Ne laser beam is passed through a linear polarizer and a half-wave plate and then focused on the cylindrite flake using a 40× objective lens (NA = 0.65). The back‐reflected light is collected by the same objective lens and directed towards a spectrometer (Horiba, iHR 520) using a beam splitter. The excitation laser is blocked by placing a Rayleigh rejection filter (Semrock, LP02‐633RE‐25) in front of the spectrometer. The parallel polarization component of the Raman spectrum is analyzed using another linear polarizer before the spectrometer.

For the polarization‐resolved DRS measurement, light from a broadband white light source (Thorlabs, SLS201L, 360–2600 nm) is passed through a linear polarizer and focused on the cylindrite flake using a 100× objective lens (NA = 0.7). The reflection spectrum (R) is obtained by collecting the back‐reflected light from the sample using the same objective lens and directing it towards the spectrometer with a beam splitter. An iris is used to spatially filter out a small area of the flake. Then the flake is removed from the collection area and the reflection spectrum from the quartz substrate (R₀) is measured. Finally, DRS is calculated as ${{{\mathrm{DRS}}}} = \frac{{R - R_0}}{{R_0}}$.

For the THG measurements, the pump laser beam at the fundamental wavelength of 1560 nm (Calmar fiber laser, pulse width 90 fs, repetition rate 80 MHz) is transmitted through a linear polarizer and a rotating half-wave plate to set the incident linear polarization angle, and then focused onto the cylindrite flake using a 40× objective lens (NA = 0.65). The transmitted THG emission from the cylindrite flake is collected by a 100× objective lens (NA = 0.7), filtered spectrally with a short pass filter (Thorlabs, FESH 900) to remove the transmitted fundamental pump beam, and then focused onto a color charge-coupled device (CCD) camera. For the spectral characterization of the THG emission, a spectrometer (Horiba, iHR 520) is used instead of the CCD camera.