Naturally occurring van der Waals heterostructure lengenbachite with strong in-plane structural and optical anisotropy

Lengenbachite is a naturally occurring layered mineral formed with alternating stacks of two constituent PbS-like and M2S3-like two-dimensional (2D) material layers due to the phase segregation process during the formation. Here, we demonstrate to achieve van der Waals (vdW) heterostructures of lengenbachite down to a few layer-pair thickness by mechanical exfoliation of bulk lengenbachite mineral. The incommensurability between the constituent isotropic 2D material layers makes the formed vdW heterostructure exhibit strong in-plane structural anisotropy, which leads to highly anisotropic optical responses in lengenbachite thin flakes, including anisotropic Raman scattering, linear dichroism, and anisotropic third-harmonic generation. Moreover, we exploit the nonlinear optical anisotropy for polarization-dependent intensity modulation of the converted third-harmonic optical vortices. Our study establishes lengenbachite as a new natural vdW heterostructure-based 2D material with unique optical properties for realizing anisotropic optical devices for photonic integrated circuits and optical information processing.


INTRODUCTION
Structural anisotropy in van der Waals (vdW) materials enables the variations in their physical properties when probed along different spatial directions. Therefore, in-plane anisotropy provides another degree of freedom for tailoring the vibrational, optical and electrical responses of vdW materials when excited with an external perturbation. Black phosphorus (BP), ReS 2 , ReSe 2 , and group-IV monochalcogenides such as GeSe are a few examples of such anisotropic vdW materials, which have been demonstrated for achieving new functionalities that are not possible with isotropic vdW materials, such as polarization-sensitive photodetectors 1-4 , modulators 5 , synaptic devices 6 , digital inverters 7 , and anisotropic resistors 8 for photonic integrated circuits. In addition to linear optical properties, anisotropy in these vdW materials also gives rise to exciting nonlinear optical effects 9-13 which can be exploited for developing high-performance nanoscale optical devices for all-optical signal processing in telecommunications. Previously, anisotropic third-harmonic generation (THG) has also been utilized for rapid optical probing and characterization of the crystal symmetries and layer thickness in BP 12 , ReS 2 9 , and GeSe 10 . On the other hand, vdW heterostructures created by combining different 2D materials have been widely studied to exhibit tailored electronic band structures, carrier mobility, optical and magnetic properties [14][15][16][17] . Apart from the nature of the vdW interaction, intercalation of ions or molecules, strain engineering, and introduction of a twist angle between stacked layers have also been explored as possible ways to modulate the electrical and optical properties of vdW heterostructures 18,19 . Recently, a capillary force-driven roll-up method for the realization of highorder vdW superlattices has been demonstrated, providing another way to control the nature of quantum confinement, dimensionality, and electronic band structures of these artificial materials 20 . These fabricated vdW heterostructures facilitate diverse functionalities beyond the reach of existing materials to realize nanoscale photonic and optoelectronic devices for optical communication and sensing 21 , transistors 22 , photodetectors 23,24 , and ultrafast lasers 25 . Until now, layer-by-layer mechanical restacking and sequential synthesis of individual layers of different vdW materials have been the most common practice for designing vdW heterostructures [25][26][27][28][29][30] . However, these procedures are cumbersome and susceptible to improper lattice alignment and the presence of undesirable interlayer adsorbates. Recent demonstration of the mechanical and liquid-phase exfoliation of natural vdW superlattices 31 including franckeite [32][33][34][35][36] and cylindrite 37,38 provides a new way for fabricating ultrathin air-stable and adsorbate-free vdW heterostructures with proper lattice orientations. Therefore, exploring the exfoliation of other natural vdW superlattices is important in the context of producing new types of vdW heterostructures with intriguing functionalities.
With that hindsight, here we study the in-plane structural and optical anisotropy of mechanically exfoliated ultrathin flakes of lengenbachite down to a few layer-pair thickness, where the superlattice is composed of alternating stacks of pseudotetragonal PbS-like layer and pseudo-hexagonal M 2 S 3 -like layer. Although the constituent 2D lattices are isotropic in nature, the incommensurability between them makes the superlattice structurally anisotropic. First, the anisotropic crystal structure and the chemical composition of lengenbachite are determined by transmission electron microscopy (TEM), energy-dispersive X-ray spectroscopy (EDXS), and X-ray photoelectron spectroscopy (XPS) analysis. Furthermore, the in-plane structural anisotropy of lengenbachite flakes is probed by angle-resolved polarized Raman spectroscopy. By performing polarization-dependent absorption measurements, we observe the optical anisotropy induced linear dichroism and a direct band gap of 2 eV in lengenbachite flakes. Moreover, we investigate the effect of crystal structural anisotropy on the optical nonlinearity by measuring polarization-dependent THG in lengenbachite flakes and estimating the third-order 1 nonlinear susceptibility. We also show that the anisotropic THG response can be further utilized for prompt determination of the lattice orientation, by characterizing the intensity profile of THG emission from the flake upon the excitation from a radially or azimuthally polarized vector beam. Finally, as a nonlinear optical device application, nonlinear optical vortex conversion from the flake is demonstrated with the polarization-dependent intensity modulation of TH optical vortices. Our results not only provide a deeper insight of the origination of structural anisotropy in natural vdW superlattice lengenbachite, but also offer a better understanding of the role of structural anisotropy in linear and nonlinear light-matter interactions in vdW materials. Furthermore, the demonstrated results are promising for the expansion of the current existing 2D material library through identification of new types of natural vdW heterostructures and can also be harnessed to realize polarization-sensitive anisotropic photonic and optoelectronic devices for future photonic integrated circuits and optical information processing.

Morphology and chemical composition of lengenbachite
Lengenbachite is a sulfosalt mineral with the approximate chemical formula of Pb 6 (Ag,Cu) 2 As 4 S 13 , which was discovered in 1904 by R. H. Solly and named after its type locality of the Lengenbach quarry in Valais, Switzerland [39][40][41] . It belongs to the subclass of minerals called misfit compounds that are composed of stacks of alternating semi-commensurate layers of a four-atomthick PbS-like pseudo-tetragonal (Q) layer and a five-atom-thick M 2 S 3 -like pseudo-hexagonal (H) layer. Previous studies on the mineral indicate that both of these layers contain Pb and As atoms in cation sites, with stoichiometries for Q-layer being (Pb,As) 4 S 4 and H-layer being (As,Pb) 2 S 3 . Figure 1a represents the simplified crystal structure of lengenbachite, indicating alternative Q-layers and H-layers bonded together via vdW forces, where the a-axis [100] is considered to be the out-of-plane stacking direction, while b-axis [010] and c-axis [001] are considered to be the two in-plane directions. Also, α, β, and γ are the angles between a-axis, b-axis and c-axis. Lengenbachite has a triclinic crystal structure belonging to the P1 space group. The geometric crystallographic data suggests that the unit cell parameters for the Q-layers and H-layers are a Q ¼ a H ¼ 18: 33 . These data indicate that the unit cell parameter of the coincidence cell or a super cell consists of two-layer-pair repeat along the a-direction with α = 18.45 Å. In b-direction and c-direction, lattice parameters are semi- . Therefore, the supercell is made of 2:3 and 12:11 Q-type and H-type unit cells in b-direction and c-direction. Figure 1b shows the characteristic appearance of lengenbachite clusters formed on a big piece of crystalline dolomite rock which is acquired from Lengenbach quarry, Binn valley, Switzerland. Figure  1c is a zoomed-in view of one such cluster, where several tabular and blade-like plates (about 3-5 mm long) of steel-gray colored lengenbachite crystals with perfect [100] cleavage are visible to be contained in a vug in white dolomite. The lengenbachite thin flakes are mechanically exfoliated using Nitto tape onto a quartz substrate. Details of the mechanical exfoliation are included in the Methods section. Figure 1d shows a reflection optical microscope image of a few exfoliated ultrathin lengenbachite crystals. The topography of the crystal surface is analyzed using atomic force microscopy (AFM) in Fig. 1e. From the height profiles of the flakes in the AFM image (inset of Fig. 1e), the thicknesses of two flakes are determined to be 28 and 36 nm approximately. The 28 nm thickness indicates that the flake consists of only 15-layer pairs. Figure 1f shows the reflection microscope image of an isolated lengenbachite flake on quartz substrate, which is used for all of our experiments. The associated AFM image in Fig. 1g confirms that the flake is 44 nm thick.
First, we analyze the crystal structure of mechanically exfoliated lengenbachite flakes by TEM studies. The exfoliated ultrathin lengenbachite flakes on quartz substrate are transferred onto a TEM nickel grid (details of the TEM sample preparation is included in the Methods section). Figure 1h is the TEM image of a typical lengenbachite flake where an in-plane striped pattern of periodic darker and lighter areas is clearly visible, indicating an out-of-plane rippling modulated to its perpendicular direction on the surface. Such rippled structures known as interlayer moiré patterns are formed in vdW heterostructures due to the elastic deformation caused by the presence of periodic strain in the crystal lattice for forcing local atomic alignment between the constituent layers. Here, the pattern is arising to force the commensuration between the two incommensurate Q-type and H-type lattices. Previously, similar interlayer moiré patterns have also been observed in other natural vdW heterostructures of franckeite [32][33][34] and cylindrite 37 . The high-resolution (HR) TEM image in Fig. 1i is a zoomed-in view on the crystal surface showing the atomic arrangements of the component layers. The fringes occurring from the ripples are spaced with a period of 3.55 nm which is approximately half of the lattice constant c = 7.016 nm of the supercell. The selected area diffraction pattern (SAED) from the crystal, normal to the [100] crystal zone axis is shown in Fig. 1j. The SAED exhibits the reflections from the corresponding Q-layers and H-layers, which are marked in yellow color and red color, respectively. Many weak superlattice diffraction spots are also visible. The SAED further indicates that the rippling on the surface is modulated along the c-direction ([001] zone axis) of the crystal. Therefore, c-direction is assigned as the rippling direction.
The chemical composition of the exfoliated lengenbachite flakes is then analyzed by performing EDXS. The averaged EDXS spectrum in Fig. 2a acquired from a thin flake confirms the presence of main elements Pb, As, Ag, Cu, and S. Peaks from other elements such as C, O, and Si are present since the sample is prepared by exfoliating naturally occurring minerals. However, the signal of Ni is present due to the nickel TEM grid. The presence of the underlying carbon film support and the accumulation of carbon-based adsorbates during the flake transfer on the TEM grid may also contribute majorly to the detected signal of C, O, and Si. The compositional stoichiometry of the lengenbachite flake is deduced by the quantification of each element in the crystal, which is provided in Table 1. The approximate chemical formula is determined to be Pb 5.9 Cu 0.73 Ag 1.44 As 4.2 S 13 . Although this approximate chemical formula is close to that of the previously determined chemical formula of lengenbachite, it is not completely charge-balanced which is mainly caused by the low sensitivity of EDXS, the heavy overlap between S K α and Pb M α peaks, as well as the overlap between Pb L α and As K α peaks. TEM-EDXS elemental maps in Fig. 2b-f of the main elements Pb, Cu, Ag, As, and S are acquired from the area of a lengenbachite crystal shown in the dark-field (DF)-TEM image in the inset of Fig. 2a. These maps further indicate a homogeneous distribution of all the elements in the crystal when viewed from the top, along a-direction of the crystal.
Next, we perform XPS on air-aged lengenbachite flakes to determine the surface chemical composition and the effect of oxidation on the surface, due to the limited penetration depth of a few nanometers in XPS. Figure 2g is the averaged XPS spectrum of air-aged lengenbachite. The quantification of the main elements is listed in Table 2. The approximate surface chemical formula of lengenbachite deduced from the compositional stoichiometry analysis is Pb 6 Cu 0:8 Ag 1:3 As 4:3 S 13 . The observed C 1s peak in the XPS spectrum suggests the presence of carbon adsorbates on the crystal surface. However, it may also be occurring due to the underlying carbon tape. The presence of a faint O 1s peak may be A. Dasgupta et al. attributed to very weak oxidation at the surface of the flake but oxygen may also be present inside the carbon adsorbates. In addition, the oxidation states of individual elements are determined by the peak fitting in the high-resolution spectra at Pb 4f, Cu 2p, Ag 3d, As 3d and S 2p binding energy regions, as presented in Fig. 2h-l. The spectra in Fig. 2h indicate that lead is present in the material in Pb 2+ (91%) and Pb 4+ (9%) oxidation states, while the spectrum in Fig. 2k shows that arsenic is present as As 3+ (61%) and As 5+ (39%). Copper (Fig. 2i), silver (Fig. 2j), and sulfur ( Fig. 2l) are exclusively present as Cu + , Ag + , and S 2− , respectively. The slight over representation of As in the derived formula may be attributed to the small deformation of the crystal surface during the sheer-assisted mechanical exfoliation process.
Anisotropic Raman scattering from lengenbachite flakes Figure 3a shows the Raman spectrum collected from the 44 nmthick lengenbachite flake in Fig. 1f The oxidation states of the individual elements in lengenbachite determined from the XPS measurements suggest that the molecular bonds between the cationic and anionic elements in lengenbachite should be analogous to those present in the crystals of its constituent components such as PbS, As 2 S 3 , As 4 S 4 , Ag 2 S, and CuS. Therefore, the Raman modes of lengenbachite can be assigned according to the Raman spectra of the constituent components, including orpiment As 2 S 3 42-44 , realgar As 4 S 4 45,46 , galena PbS 47 , acanthite Ag 2 S 48-50 , and covellite CuS 51,52 . The 77 cm −1 peak is assigned as the Raman modes of As 2 S 3 (69 cm −1 ) and CuS (65 cm −1 ). The peak at 101 cm −1 is attributed to the silver lattice vibrational mode of Ag 2 S (93 cm −1 ). The 134 cm −1 peak corresponds to both the A g mode of As 2 S 3 (136 cm −1 ) and the mode of CuS (142 cm −1 ). The 158 cm −1 peak represents a combination of the A g mode of As 2 S 3 (154 cm −1 ), the transverse acoustic and optical phonon modes of PbS lattices (154 cm −1 ), and the silver lattice vibrational mode of Ag 2 S (147 cm −1 ). The peak at 182 cm −1 corresponds to the combination of the Raman modes of As 4 S 4 (183 cm −1 ) and As 2 S 3 (179 cm −1 ) as well as the Ag-S stretching mode of Ag 2 S (188 cm −1 ). The 208 cm −1 peak is attributed to a combination of the A g mode of As 2 S 3 (203 cm −1 ) h TEM image of a lengenbachite flake displaying the characteristic striped patterns along the c-axis, caused by the out-of-plane rippling due to the forced commensuration of Q-layers and H-layers. Scale bar is 40 nm. i HRTEM image illustrating the atomic arrangements. Scale bar is 5 nm. j Corresponding SAED pattern consisting with diffraction spots associated with Q-layer (yellow) and H-layer (red). Scale bar is 1 nm −1 .

Fig. 2
Chemical composition characterization of lengenbachite crystal. a Averaged EDXS spectrum from a lengenbachite thin flake. The peaks of main elements Pb, Cu, Ag, As, and S labeled in red color are used for the quantification in Table 1. b-f TEM-EDXS maps of main elements Pb, Cu, Ag, As, and S acquired from the crystal area marked in the DF-TEM image in the inset of a. Scale bar is 50 nm. g Averaged XPS spectrum of lengenbachite from multiple measurements. The peaks labeled in red color are used for the quantification in Table 2. h-l Highresolution XPS spectra in the binding energy regions of Pb 4f, Cu 2p, Ag 3d, As 3d, and S 2p, respectively.  and the longitudinal optical phonon modes of PbS lattices (204 cm −1 ). The 221 cm −1 peak is assigned as the Raman mode of As 4 S 4 (221 cm −1 ). The peak at 245 cm −1 belongs to the Ag-S stretching mode of Ag 2 S (243 cm −1 ). The 319 cm −1 peak corresponds to the A g mode of As 2 S 3 (311 cm −1 ). The 340 cm −1 peak is attributed to the mode of As 4 S 4 (342 cm −1 ). The peak at 355 cm −1 represents the combined Raman modes of As 2 S 3 (355, 359 cm −1 ) and As 4 S 4 (354 cm −1 ). The peak at 370 cm −1 is assigned as the Raman modes from As 4 S 4 (369, 375 cm −1 ). The 432 cm −1 peak belongs to the S-S stretching mode of Ag 2 S (430 cm −1 ). Finally, the 492 cm −1 peak is related to the S-S stretching mode of CuS (475 cm −1 ). The in-plane structural anisotropy of triclinic lengenbachite crystal can be simply revealed by angle-resolved polarized Raman spectroscopy. Previously the method has been used for determining the crystal axes of other naturally occurring vdW heterostructures such as franckeite and cylindrite. Therefore, we perform polarization-resolved Raman measurements for the parallel and perpendicular polarization components of the scattered light from the 44 nm-thick lengenbachite crystal. In Fig. 3b, c, we plot the evolution of parallel and perpendicular components of the Raman spectra with the linear polarization angle of incident light, where 0°indicates the linear polarization along x-axis as indicated in the inset of Fig. 1f. It is obvious that the intensities of the Raman modes vary periodically depending on the excitation polarization angle, and no identifiable shifts are observed. According to the classical Raman selection rule, Raman intensity can be expressed as I / e i Á R Á e s j j with e i and e s being the unit polarization vectors of the incident and scattered light, while R is the Raman tensor of each Raman mode. Under parallel configuration e i ¼ e s ¼ cos θ sin θ 0 ½ , and under perpendicular configuration e s ¼ Àsin θ cos θ 0 ½ . For triclinic lengenbachite crystal, the Raman tensor of the A g Raman modes can be expressed as, (1) where a, b, c, d, e, and f are the complex elements of the Raman tensor. Therefore, the Raman intensity of the A g modes under parallel and perpendicular configurations can be written as, and where ϕ ab , ϕ ad , ϕ bd represent the phase differences between a and b, a and d, b and d, respectively. The polar plots of parallel and perpendicular components of each Raman mode are plotted in Fig. 3d, e. The experimental data (black points) match well with the theoretical fits (red curves) with Eqs. (2) and (3). It is noted that the parallel components of the Raman peaks at 134, 245, 355, 432, and 490 cm −1 show an anisotropic two-fold symmetric pattern with a maximum at 0°and a secondary maximum at 90°. On the other hand, the Raman peaks at 158, 182, 208, and 319 cm −1 exhibit a similar anisotropic pattern but with a maximum at 90°a nd a secondary maximum at 0°. The polarization-resolved Raman scattering have been used to identify the crystal axes of anisotropic crystals, where the Raman intensity of parallel polarization component of the A g mode reaches a maximum or a secondary maximum value as the excitation linear polarization is oriented along any of the crystal axes. In order to validate these results and unambiguously determine the crystal axis, we also perform both polarization-resolved Raman scattering and TEM characterizations on one thin flake of lengenbachite transferred onto a nickel TEM grid. The experimental results are included in Supplementary Fig. 1 to show the one-to-one correspondence between the polar plots of anisotropic Raman modes and the rippling direction of the crystal. It is evident that the Raman modes at 245 and 355 cm −1 reach the maximum intensities along the rippling direction (c-axis) of the lengenbachite crystal, whereas the Raman modes at 208 and 319 cm −1 exhibit the secondary maximum intensities along the rippling direction. By comparing the anisotropic Raman modes of this flake with those of the 44 nm lengenbachite crystal shown in Fig. 3d, we infer that the rippling direction (c-axis) of the 44 nm crystal displayed in Fig. 1f is aligned along 0°(x-axis). The perpendicular components of most Raman peaks show an anisotropic four-lobe polar pattern, where the intensities reach a minimum value when the excitation linear polarization is along any of the crystal axes in 0°(x-axis) or 90°(y-axis).

Polarization-dependent absorption and linear dichroism
Next, we explore the evolution of the absorbance of lengenbachite crystals as a function of film thickness by monitoring the optical absorption spectra from thin flakes of three different thicknesses of 44, 100, and 300 nm within the wavelength (λ) range of 420-800 nm. Description of the experimental setup is included in the Methods section. Figure 4a, b plot the measured reflection (R) and transmission (T) spectra. The reflectance of the 44 nm flake remains flat within the wavelength range 420-650 nm and then decreases gradually, whereas the reflectance of the 100 nm flake gradually increases within 420-650 nm and remains flat afterwards. This justifies the yellowish and reddish-yellow glow of the 44 and 100 nm flakes in the reflection image. On the other hand, the 300 nm flake shows an interesting trend with a reflection dip around 660 nm, and therefore exhibits a greenish color under reflection microscope. The transmittance of the 44 and 100 nm flakes shows a similar trend with gradually increasing value within the wavelength range 420-800 nm range, while the 300 nm flake exhibits a very low transmission in 420-600 nm and a transmission peak around 680 nm. These trends justify the reddish-black color of the crystals in the transmission images. The absorption (A = 1 -R -T) spectra from the flakes of different thicknesses are plotted in Fig. 4c, showing a thicker crystal has a higher absorbance as expected. For the 44 and 100 nm crystals, an absorption peak is observed around 460 nm, whereas for the 300 nm crystal, this peak is shifted to 500 nm. The opposite oscillations in the reflection and transmission spectra around the wavelength of 660 nm is attributed to the thin film interference of light inside the 300 nm flake. However, these are almost canceled out in the absorption spectra.
The rippling caused from the periodic lattice strain can further result in anisotropic optical properties. Therefore, polarizationdependent absorption spectroscopy is performed to study the linear dichroism response of lengenbachite. Figure 4d shows a series of absorption spectra acquired from the 44 nm-thick flake as a function of the rotation angle of the linear polarizer. It is observed that the absorbance varies periodically with the linear polarization angle, exhibiting a maximum absorbance when the incident polarization is along 0°(x-axis) and a minimum along 90°( y-axis). Figure 4e and f are the polar plots of the evolution of absorbance as a function of the linear polarization angle at two different wavelengths of 460 and 520 nm for the 44 nm-thick flake. The experimental data is fitted with a periodic function of A θ ð Þ ¼ A x cos 2 θ þ A y sin 2 θ, where A x and A y are the absorbance values for incident linear polarization along x-axis and y-axis. It has been reported that for franckeite the highest absorbance is obtained as the incident polarization is along the rippling direction (c-axis) of the crystal 34 . Therefore, we identify x-axis as the c-direction and y-axis as the b-direction for the lengenbachite crystal. Another important inference drawn from Fig. 4e, f is that lengenbachite exhibits strong linear dichroism with the ratio A x /A y = 1.13 and 1.14 at λ = 460 nm and 520 nm, even though the constituent Q-layers and H-layers are isotropic. Linear dichroism is also observed in the 100 and 300 nm-thick flakes, where the obtained ratio at the wavelength of 460 nm is A x /A y = 1.20 and 1.18, respectively. The measured polarization-dependent absorbance spectra and the polar plots of the evolution of absorbance at the wavelength of 460 nm for the 100 and 300 nm flakes are included in Supplementary Fig. 2. Furthermore, the moiré patterns formed by the periodic in-plane strain in the vdW superlattice could also strongly modulate the electronic band structure and the band gap of the material. For the 44 nm-thick flake, we observe a 30 meV modulation in the direct allowed optical transition depending upon the linear polarization angle of the incident light (see Supplementary Fig. 3), showing the direct optical band gap varies from 1.96 to 1.99 eV as the linear polarization angle is rotated from the rippling direction (c-axis) of the flake to its perpendicular direction (b-axis). For lengenbachite crystal, the observed band structure anisotropy is not very strong because the strain-induced modulation of the electronic band structure happens along both the principal axes of the crystal. The optical band gap E g of the lengenbachite flakes are achieved by the Tauc plot based on the relation of αhν ð Þ n ¼ mðhν À E g Þ, where α is the absorption coefficient corresponding to the photon energy hν (h is the Plank's constant) with ν being the frequency of the incident photon and m is a constant. The value of n = 2 signifies an allowed direct transition, while n = 1/2 or n = 2/3 represents an allowed indirect or forbidden transition. The value of α can be extracted from the measured R and T plots from Fig. 4a, b with the flake thickness d as 53 , Figure 4g plots the (ahv) 2 as a function of the photon energy (hv) for the three flakes. A good linear fit in the Tauc plots is obtained for all the cases, which signifies an allowed direct transition and hence the presence of direct band gap. The extracted direct optical band gap from the Tauc plots for 44, 100, and 300 nm-thick lengenbachite flakes are estimated to be 1.96, 2.03, and 1.98 eV respectively. It is noteworthy that the band gap 2 eV of lengenbachite tends to lie between the band gaps of its component binary sulfides of As 2 S 3 at 2.37 eV and PbS at 0.37 eV, which is consistent with the previous studies in the optical band gaps for ternary and quaternary sulfide minerals 48 . In Supplementary Fig. 4, we further confirm the presence of the direct optical band gap in the lengenbachite crystals by measuring the photoluminescence spectra from all the three flakes with a 405 nm excitation laser. We observe a photoluminescence peak around the wavelength of 616 nm for each flake, which is consistent with the energy of the measured direct band gap of 2 eV estimated from the Tauc plot.

Anisotropic THG in lengenbachite crystals
To investigate the effect of structural anisotropy on anisotropic nonlinear optical properties, we explore the polarizationdependent THG from lengenbachite thin flakes. Details of the experimental setup are provided in the Methods section. Figure 5a shows the transmission optical microscope image of THG emission from the area of the 44 nm-thick flake of Fig. 1f, when it is illuminated with a 1560 nm wavelength pump laser of spot size 1.5 µm. The recorded THG spectrum in Fig. 5b shows that the peak at 520 nm is exactly one-third of the excitation wavelength. The cubic power-law dependence in the log-scale plot of the average THG power as a function of the average pump power in Fig. 5c further confirms the THG process. For the average pump power of 2 mW corresponding to a peak pump irradiance of 13.6 GW cm −2 , the achieved THG conversion efficiency is around 1.1 × 10 −8 for the 44 nm-thick flake. When the crystal axes c, b, and a are oriented along x, y, and z-directions, the contracted form of the third-order nonlinear susceptibility tensor (χ (3) ) for the triclinic crystal system of lengenbachite is expressed as 54 , where the first subscript 1, 2, 3 represents x, y, z respectively, and the second subscript refers the following combination of three components as, We consider the fundamental beam has the linearly polarized electric field along an angle θ to the x-axis at frequency ω as E ω ð Þ ¼ E 0 cos θx þ sin θŷ ð Þ , wherex andŷ are the unit vectors along x-axis and y-axis. Note that the polarization of the excitation electric field always remains in the x-y plane, there is no contribution from the χ (3) components containing z terms in the THG emission. Hence, only eight non-zero elements of χ (3) as χ 11 , χ 12 , χ 18 , χ 19 , χ 21 , χ 22 , χ 28 , and χ 29 will be contributing to the THG signal and the x-polarized and y-polarized components of the THG intensity can be expressed as, ð Þ x / χ 11 cos 3 θ þ χ 12 sin 3 θ þ 3χ 18 cos θsin 2 θ þ 3χ 19 sin θcos 2 θ À Á 2 I 3ω ð Þ y / χ 21 cos 3 θ þ χ 22 sin 3 θ þ 3χ 28 cos θsin 2 θ þ 3χ 29 sin θcos 2 θ À Á 2 (6) Figure 5d plots the dependence of the THG emission power from the 44 nm-thick lengenbachite crystal as function of incident linear polarization angle with respect to x-axis for fundamental beam power at 1.1 mW. The desired linear polarization of the pump beam is obtained by placing a linear polarizer oriented along x-axis (c-axis of the crystal) and a rotating half-wave plate. The black, red and blue data points represent the measured x-component, y-component and total THG power, while the solid curves are theoretical fits using Eq. (6), which agree well with the experimental data. It is clear that the polarization-dependent THG emission from the lengenbachite crystal displays an anisotropic two-fold symmetric pattern, with the maximum THG power at 0°a s the incident linear polarization is parallel to the rippling direction along the c-axis and the second maximum at 90°, which is perpendicular to the ripples along the b-axis. We further confirm this result by measuring the polarization-dependent THG emission from the thin flake on nickel TEM grid, which is characterized in Supplementary Fig. 1. According to the measured data in Supplementary Fig. 5, it is shown that the maximum THG power occurs as the incident polarization is orientated along the rippling direction (c-axis), which is determined directly by TEM imaging. Furthermore, the dependence of THG emission on the lengenbachite flake thickness is investigated. Figure 5h plots the evolution of the measured average THG power (P (3ω) ) as a function of the flake thickness (d) at the pump power of P (ω) = 1.1 mW. The black data points are the measured THG power from 12 different flakes of thicknesses ranging in 28-156 nm, when the incident linear polarization is parallel to the rippling direction (c-axis) of each flake. It is observed that the THG power increases quadratically up to around 4.63 pW as the flake thickness is increased to 85 nm, giving the maximum conversion efficiency of 4.21 × 10 −9 . However, the THG power exhibits an exponential decay for the thickness larger than 85 nm, where the significant attenuation of THG power is attributed to the strong absorption through thick lengenbachite flakes. Thereby, the measured thickness-dependent THG response can be fitted with an exponentially decaying function P 3ω ð Þ d ð Þ ¼ Cd 2 exp À 4πk3d λ3 with C as a constant and k 3 being the imaginary part of refractive index at λ 3 = 520 nm, where the fitting (red curve) shows k 3 = 1.02. Additionally, χ (3) value of lengenbachite crystal can be estimated by using the following relationship 11 , where n 1 and n 3 are the real part of the refractive index at fundamental frequency (ω) and THG frequency (3ω), and the other terms are the parameters of the Gaussian fundamental pulsed laser including beam spot size W = 1.5 µm, pulse width τ = 90 fs, and repetition rate f rep = 80 MHz. It is noted that the refractive index information of lengenbachite crystal is still unknown, but the refractive indices of its main constituent layers of PbS and As 2 S 3 are available, with n 1 = 4.24 and n 3 = 4.34 for PbS, n 1 = 2.44 and n 3 = 2.73 for As 2 S 3 55 . By assuming an averaged refractive index of n 1 = n 3 = 3.5 for lengenbachite crystal, the estimated χ (3) value is 2.18 × 10 −19 m 2 V −2 . Also as the refractive index varies from 2.5 to 4.5, the χ (3) value stays between 1.13 × 10 −19 and 3.56 × 10 −19 m 2 V −2 . Nonetheless, the χ (3) of lengenbachite crystal is of the same order as the widely studied anisotropic 2D material  Table 3 the measured THG anisotropy ratio and χ (3) value of lengenbachite crystal are compared with the reported values of other anisotropic nonlinear 2D materials such as franckeite, cylindrite, BP, ReS 2 , GeSe, GeAs, and SiP.

Prompt determination of crystal axes with vector beams
Radially or azimuthally polarized doughnut-shaped vector beams carry all the linear polarization components equally distributed in its spatial intensity profile. Since the THG response from lengenbachite flakes is directly related to the in-plane structural anisotropy of the crystal, one can directly probe the crystallographic information by imaging the THG intensity profile upon the excitation with a radially or azimuthally polarized vector beam. Figure 6a is the THG image from the 44 nm-thick lengenbachite flake in Fig. 1f when it is excited with a doughnut-shaped radially polarized vector beam. The rippling direction of the same flake has already determined to be along x-axis from the Raman and THG characterization. Figure 6b shows the same THG image when the background white light is turned off. It is evident that the THG emission has two unequal maxima along the x-direction and ydirection, indicating the anisotropic χ (3) tensor elements of the crystal. The more intense maxima along the x-direction occur due to the fact that the polarization of the pump vector beam is oriented along the x-axis corresponding to the crystal's rippling direction (c-axis) at these positions, where the largest χ 11 element in χ (3) tensor dominates. While the secondary maxima along the y-direction are observed as the pump vector beam at these locations is polarized perpendicular to the crystal's rippling direction (b-axis), where the χ 22 element in χ (3) tensor is present. Figure 6c plots the THG intensity profile as a function of the Table 3. Comparison of THG response in lengenbachite with other anisotropic nonlinear 2D materials.

Material
Crystal system Space group (bulk) Thickness (nm) Fundamental wavelength (nm) Ref.  azimuthal angle. The extracted intensity profile (black data points) is fitted (red curve) with Eq. (6). Figure 6d, e further show the THG images from the crystal when it is pumped with an azimuthally polarized vector beam. Opposite to the previous case, here the more intense maxima in the THG emission are observed along the y-direction while the secondary maxima are in the x-direction. The corresponding THG intensity profile depending on the azimuthal angle is plotted in Fig. 6f. Therefore, imaging the THG emission from the lengenbachite crystal upon the excitation with a doughnut-shaped radially or azimuthally polarized vector beam is useful for rapid characterization of the crystal axes and in-plane structural anisotropy, which can be further extended for other anisotropic vdW materials. In contrast to polarization-dependent Raman scattering, angle-resolved photocurrent measurement, and angle-resolved THG experiment, where the crystallographic characterization requires multiple measurements, the current technique directly provides the crystallographic information with only one single THG image collection.
Third-harmonic optical vortex conversion As a potential nanoscale device application, we demonstrate TH optical orbital angular momentum (OAM) conversion by exciting the lengenbachite flake with a fundament vortex beam. An optical vortex beam possesses a helical phase structure of exp(ilϕ) and thus carries the OAM of lℏ per photon and exhibit the doughnutshaped intensity profile, where l is the topological charge (TC) and ϕ is the azimuthal angle around the phase singularity. The unbounded dimension of TC of light has made it one of the most vital parameters to store, control and transport information in optical communication. Therefore, optical vortex beams have been harnessed for various applications such as large-scale optical data transmission 57 , high-dimensional quantum information processing 58 , quantum memory 59 , and nanoscale optical tweezers.
In the context of our work, the lengenbachite thin flakes provide a platform to realize nonlinear vortex conversion devices, where the TC tripling in THG process is demonstrated. The electric field profile of the fundamental vortex beam with TC = l 1 at frequency ω can be expressed as, , where r is the distance from the beam center, and w 1 is the radius for which the Gaussian term falls to 1/e of its on-axis value. Now that the electric field of the TH vortex beam at frequency 3ω is given , with a tripled TC l 3 = 3l 1 and w 3 reduced by a factor of ffiffi ffi 3 p from w 1 . Figure 7a is the recorded transmission image of the doughnut-shaped intensity pattern for the fundamental vortex beam with TC = 1. Details of the experimental setup for generation of the vortex beam are included in the Methods section. The associated TC of the fundamental vortex beam is confirmed by performing an astigmatic transformation of the recorded image using a cylindrical lens placed before the CCD camera where the high-TC optical vortex splits into its constituent elementary vortices with TC = 1 to form an tilted dark-striped pattern near the focal plane of the cylindrical lens. The number of dark stripes in the converted image indicates the TC of the vortex beam. Figure 7b confirms that the fundamental beam has TC = 1. Figure 7c is the recorded image of the converted TH vortex beam from the lengenbachite nanoflake. The corresponding cylindrical lens image in Fig. 7d is a confirmation that the TH vortex beam has TC = 3. The measured line profiles of the intensity of the fundamental (I fund ) and TH (I THG ) vortex beams are plotted in Fig.  7e, where the experimental data points are fitted by the calculated profiles of E 2 fund and E 2 THG (solid curves). It shows that the measured w3 w1 ratio is 0.58, which is close to the expected value of 1= ffiffi ffi 3 p . Further, we explore the effect of the in-plane structural anisotropy of the crystal on the intensity of the converted TH vortex beam by performing the polarization-resolved THG study upon excitation with the fundamental vortex beam. The intensity variation of the TH vortex beam as a function of the linear polarization angle of the fundamental vortex beam is plotted in Fig. 7f, which shows a similar polarization dependence to that of the anisotropic THG response in the case of fundamental Gaussian beam excitation in Fig. 5d. It proves that the nonlinear optical anisotropy of lengenbachite crystal remains intact under the excitation from a fundamental optical vortex beam. Figure 7g displays the THG images of the converted vortex beams for different fundamental beam polarization angles at 0°, 45°, 90°, 135°, and 180°, showing a clear intensity modulation depending on the incident polarization angle. Therefore, in-plane structural anisotropy in lengenbachite flakes can be exploited for polarization-multiplexing nonlinear vortex generation and OAM conversion in photonic integrated circuits and optical communication.

DISCUSSION
In summary, we have introduced lengenbachite as a new type of naturally occurring vdW heterostructure composed of alternating stacks of PbS-type and M 2 S 3 -type isotropic layers. It is demonstrated that lengenbachite can be mechanically exfoliated down to 15 layer-pair thickness. With further study and the use of other exfoliation methods such as liquid-phase exfoliation, it may be possible to exfoliate lengenbachite down to single unit cell thickness. From the TEM analysis of ultrathin lengenbachite flakes, we have observed the interlayer moiré patterns as the small outof-plane rippling in the crystal surface and thus a strong in-plane structural anisotropy of the crystal, which is caused by the vdW interaction and lattice deformation between the alternately stacked constituent layers. The chemical composition of bulk lengenbachite crystals and their surface chemical stoichiometry further prove the presence of the constituent PbS-like and M 2 S 3like layers containing the elements Pb, Cu, Ag, As, and S. Optical probing of the in-plane structural anisotropy and the rippling direction of ultrathin lengenbachite flakes is demonstrated through angle-resolved polarized Raman spectroscopy. These results provide further insight for understanding the origin of structural anisotropy in lengenbachite crystal and other naturally occurring vdW heterostructures, where the initial in-plane lattice symmetry of the individual constituent layer is broken. The effect of the structural anisotropy on the linear and nonlinear optical properties of the lengenbachite crystal is further demonstrated. The lengenbachite flakes are found to exhibit strong linear dichroism through polarization-dependent absorption spectroscopy, which may be harnessed for designing polarizationsensitive photodetectors and modulators. A direct optical band gap is also observed around 2 eV with the Tauc plot analysis. Furthermore, anisotropic THG response from lengenbachite thin flakes are investigated to show that the THG power as the linearly polarized pump beam is along the rippling direction of the crystal (c-axis) is around twice of that along the perpendicular direction (b-axis). The anisotropic χ (3) tensor elements are extracted and the χ (3) value of lengenbachite crystal is estimated to be 2.18 × 10 −19 m 2 V −2 , indicating that lengenbachite exhibits a strong nonlinear optical response as other anisotropic nonlinear 2D materials. In addition, it is demonstrated that the lattice orientation of lengenbachite flake can be promptly identified by characterizing the THG image from the flake upon the excitation with a radially or azimuthally polarized vector beam. The air-stable nature and strong anisotropic χ (3) tensor of ultrathin lengenbachite flakes make them to be a promising platform for realizing future on-chip anisotropic nonlinear optical devices for optical information processing. As one potential application, here we have demonstrated the nonlinear OAM conversion of optical vortices from lengenbachite thin flakes, where the TC of the converted TH vortex beam gets tripled from that of the fundamental vortex beam. The nonlinear optical anisotropy in lengenbachite allows the intensity modulation of the converted TH vortex by changing the polarization of the fundamental vortex beam. These results may be useful for harnessing natural vdW heterostructure-based nonlinear optical devices for frequency, polarization, and OAM multiplexing, data decoding and encoding in optical information processing and optical communication. Further study of anisotropic electrical and magnetic properties of lengenbachite thin flakes will also be useful for establishing this new type of vdW heterostructure for exploring other on-chip device applications.

Sample preparation
The quartz substrate is sonicated in deionized water, acetone, and isopropanol one after the other. The lengenbachite thin flakes are mechanically exfoliated from the bulk natural lengenbachite mineral acquired from Lengenbach quarry, Binn valley, Valais, Switzerland 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. For preparing the TEM sample, exfoliated lengenbachite flakes are transferred from the quartz substrate to TEM nickel 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 lengenbachite 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 lengenbachite 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
The Raman spectrum is acquired in a reflection microscope setup, where a 632.8 nm He-Ne laser beam is passed through a linear polarizer and a halfwave plate (HWP) and then focused on the lengenbachite flake using a 40× objective lens (NA = 0.65). The back-reflected light is collected by the same objective lens and directed to an optical spectrometer (Horiba, iHR 520) using a beam splitter. The Rayleigh scattered light of the laser is then filtered out using a longpass filter (Semrock, LP02-633RE-25) in front of the spectrometer. The parallel and perpendicular polarization components of the Raman spectrum are analyzed using another linear polarizer before the spectrometer.
For the polarization-dependent absorption measurement, light from a broadband white light source (Thorlabs, SLS201L, 360-2600 nm) is passed through a linear polarizer and focused on the lengenbachite flake using a 80× objective lens (NA = 0.5). The reflection spectrum is obtained by collecting the back-reflected light from the flake using the same objective lens and directing it to the spectrometer with a beam splitter, while the transmission spectrum is collected using a 100× objective lens (NA = 0.7). An iris is used to spatially filter out a small area of the flake in both the cases. After normalizing the reflection and transmission spectra with the source spectra, reflectance (R) and transmittance (T) are measured. Finally, the absorbance (A) spectrum is obtained by using the relationship A = 1 -R -T.
For the THG measurements, a femtosecond pulsed 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 HWP to set the incident linear polarization angle, and then focused onto the lengenbachite flake using a 40× objective lens (NA = 0.65). For the spectral characterization, the transmitted THG emission from the flake is collected by another 100× objective lens (NA = 0.7) and directed towards a spectrometer (Horiba iHR 520). The contribution of the fundamental beam is spectrally filtered out using a shortpass filter (Thorlabs, FESH 900). For the imaging collection, the THG emission is focused onto a color charge-coupled device (CCD) camera instead of the spectrometer. For the generation of radially and azimuthally polarized vector beams, the linearly polarized fundamental laser beam is passed through a zero-order vortex HWP (Thorlabs, WPV10L-1550) and then focused onto the lengenbachite flake. For the generation of the vortex beam with TC = 1, the linear polarization of the incident laser beam is converted into circular polarization by passing through a quarter-wave plate and a zero-order vortex HWP.