Natural layered mercury antimony sulfosalt livingstonite with anisotropic optical properties

Naturally occurring layered mineral livingstonite is identified as a new type of van der Waals (vdW) heterostructure based 2D material, consisting of two commensurately modulated alternating layers of HgSb2S4 and Sb2S4. The heterostructures of livingstonite crystal are prepared as thin flakes via mechanical exfoliation method. The prepared livingstonite crystals are further investigated in the context of vibrational, linear, and nonlinear optical properties, including anisotropic Raman scattering, wavelength-dependent linear dichroism (LD) transition effect, birefringence, and anisotropic third-harmonic generation (THG). Owing to the monoclinic crystal structure, livingstonite crystals exhibit strong anisotropic vibrational and optical responses. In contrast to conventional vdW heterostructures, the anomalous LD transition effect and the evolution of butterfly-shaped THG emission pattern in livingstonite crystals are demonstrated. Furthermore, the third-order nonlinear susceptibility is estimated for livingstonite crystal using the thickness-dependent THG emission response. Overall, the discussed outcomes establish livingstonite as a new type of naturally grown vdW heterostructure based 2D material and offer insights in tailoring linear and nonlinear light-matter interactions in such vdW heterostructures, which may find further relevance in polarized optical applications and on-chip integrated photonic circuits. © 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement


Introduction
Advances in nanoscale manipulation and deterministic placement of van der Waals (vdW) bonded materials have garnered significant attention in the context of tailoring physical properties, enabling tunable functionalities, and realizing applied device designs [1][2][3]. Since the successful mechanical isolation of graphene [4], extensive efforts on the exfoliation of other vdW materials, the manipulation and deterministic placement of vdW heterostructures, and the exploration of their novel physical properties have ushered many new avenues of fundamental scientific studies and technological modernization [1][2][3]. Among these, the manipulation of vdW heterostructures has recently gained substantial consideration in the context of reassembling individual atomic layer of different materials to design vdW heterostructures with distinct stacking symmetry and functionalities [1,2,5]. In contrast to conventional 2D materials, the properties of these engineered vdW heterostructures can be efficiently controlled via tailoring the interactions between different quantum states facilitated by individual atomic layers [1][2][3]5,6]. As a result, the physical properties such as charge carrier mobility, favorable band structure, phonon vibration, optical and magnetic responses can be effectively customized, which are very crucial to many applications in thermoelectrics, optoelectronics, nanophotonics, and optical communication [1,2,[5][6][7][8].
Based on these advances, it can be learned that the physical properties of these vdW heterostructures are customized by varying their constituent layer materials, layer stacking sequences, layer-to-layer interactions, component lattice modulations, and crystallographic alignment. In particular, franckeite (Pb 5 Sn 2 FeSb 2 S 14 ) [18], cylindrite (Pb 3 Sn 4 FeSb 2 S 14 ) [15] and lengenbachite (Pb 6 (Ag,Cu) 2 As 4 S 13 ) [16] exhibit incommensurate crystal lattice stacking with in-plane rippling patterns, whereas nagyágite [Pb 2 (Pb,Sb) 2 [20,21] have commensurate crystal lattice stacking. Therefore, comprehensive understanding of light-matter interactions and structural dependent responses in these vdW quantum materials become inevitable to enable their full-fledged engagement in future applications. Recently, several important findings are intended in this regard with incommensurate vdW heterostructures [14][15][16], however very few findings are directed towards commensurate lattice structures [17], which are in its infancy stage and several aspects remain elusive.
Motivated with this, herein we identify livingstonite as a new commensurate vdW heterostructure with anisotropic physical properties. Livingstonite is a mercury antimony sulfosalt with the idealized chemical formula of HgSb 4 S 8 , which was first described in 1874 for an occurrence in Huitzuco de los Figueroa, Guerrero, Mexico and named in honor of the Scottish explorer and missionary in Africa, David Livingstone [22]. It occurs in low-temperature hydrothermal veins usually associated with cinnabar, stibnite, sulfur, and gypsum. Livingstonite has the morphology in elongated needles, columnar to fibrous crystals, or globular masses of interlaced needles showing blackish-gray color and red streak. It is translucent in thin fragments with deep red internal reflections and has an elastic tenacity with a Mohs scale hardness of 2. Importantly, the recent demonstration of electronic properties and optical band gap study using computational approaches [23] has ushered an avenue of its engagement for integrated nanophotonic and optoelectronic devices. However, the vibrational, linear, and nonlinear optical properties of livingstonite crystal remain unexplored yet. Here, we systematically explore anisotropic vibrational, linear, and nonlinear optical responses of mechanically exfoliated commensurate vdW heterostructures of livingstonite crystal. First, we determine the crystal structure, chemical composition and empirical formulation using high-resolution transmission electron microscopy (HRTEM), and energy dispersive X-ray spectroscopy (EDXS) techniques. Afterwards, the evolution of anisotropic vibrational modes and their thickness dependence in livingstonite crystals are investigated. Subsequently, wavelength-dependent optical polarity switching along with the birefringent nature of livingstonite crystal are explored using polarization-resolved optical absorption spectroscopy and cross-polarized optical microscopy, respectively. Finally, the anisotropic nonlinear optical response with respect to the incident linear polarization angle is probed and the third-order nonlinear susceptibility of livingstonite crystal is estimated. We envisage that these outcomes and related discussions will be relevant in advancing future miniaturized optoelectronic and nanophotonic devices such as optical sensors, nonlinear signal processors, optical switches, and integrated photonic circuits.

Characterization of crystal structure and chemical composition of livingstonite
Livingstonite belongs to the monoclinic crystal with space group A2/a and the unit cell parameters of a = 30.154 Å, b = 3.995 Å, c = 21.426 Å, α = γ = 90°and β = 104.265° [24][25][26]. As illustrated in Fig. 1(a), the crystal structure of livingstonite is built up of two alternating layers of HgSb 2 S 4 and Sb 2 S 4 which run parallel to the c-axis. The top picture in Fig. 1 [29,30] which belong to the pavonite homologous series. The structure of grumiplucite consists of an AB stacked layers built up by Bi 2 S 4 chains, which are connected through Hg atoms as the similar Sb 2 S 4 chains of livingstonite. It has been found that isolated bilayers of HgSb 4 S 8 are weakly bonded by vdW interactions with an interlayer cohesive energy of 32 meV/Å 2 , which is much less than the approximate maximum interlayer cohesive energy of potentially exfoliable layered materials in 150 meV/Å 2 [23] indicating the feasibility of mechanical exfoliation of the bulk livingstonite mineral. Here the bulk livingstonite mineral (from Khaydarkan, Fergana Valley, Alai Range, Osh Oblast, Kyrgyzstan) as shown in Fig. 1(b) and 1(c) is mechanically exfoliated with Nitto tape (SPV 224) to obtain thin flakes with different thicknesses on glass substrate ranging from 26 nm to 192 nm. Figure 1(d) shows the AFM image of the 26 nm-thick livingstonite crystal, corresponding to the thickness of 24 layer pairs. In addition, the optical reflection and transmission microscope images with the corresponding AFM images of other exfoliated livingstonite crystals are also shown in Fig. 1(e)-(p), which emphasize the obtained flake size, flake thickness, and surface smoothness of the prepared samples. Previously, the lattice orientation of graphene hexagonal grain has been identified from its layer edges [31]. Here, it is observed that the exfoliated livingstonite flakes have sharp edges, which follow the orientations of the crystal's a-axis or b-axis.
Next, the livingstonite crystals are characterized using HRTEM and EDXS techniques to understand the atomic arrangement and the constituent element composition. Figure 2(a) shows the captured HRTEM image of a thin flake, displaying the Sb 2 S 4 double chain structures joined together by Hg atoms running along the b-axis, where the individual atoms are clearly visible. The brighter spots correspond to the atomic positions of Hg which has a higher atomic number than Sb. It also shows the measured lattice spacing of ∼ 7.31 Å and 2.00 Å along the a-axis and b-axis, respectively. The intersectional angle between these two sets of planes is found ∼ 90.2°. These lattice parameters are found consistent with the [400] and [020] sets of planes for the monoclinic livingstonite crystal. Figure 2(b) illustrates the indexed fast Fourier transform (FFT) of the HRTEM image shown in Fig. 2(a). The obtained reflections of the crystal lattice are systematically marked, which emphasize the single crystalline nature of the probed sample. To quantify the chemical composition and element distribution in the crystal, we further preform the EDXS measurement. Figure 2(c) shows the average EDXS spectrum collected from the livingstonite crystal. The recorded spectrum confirms the presence of prime constituent elements of mercury (Hg), stibnite (Sb), and sulfur (S) and provides a quantification of chemical stoichiometry of the individual heterolayers. The ratio between Hg, Sb and S is estimated as 0.94:4.91:8.00.  Table 1 summarizes the compositional stoichiometry analysis, resulting in an empirical formula of Hg 0.94 Sb 4.91 S 8.00 . It can be observed that the obtained formula is not the excellent match with the generic livingstonite formula HgSb 4 S 8 . This deviation between these two formulas can be attributed to the geological origins and the presence of small impurities of carbon (C), oxygen (O) and silicon (Si) in these naturally grown minerals. The observed copper (Cu) peaks in the EDXS spectrum are accounted from the TEM grid. The inset in Fig. 2(c) shows the dark-field TEM image of the crystal with the region of interest marked by the orange square box to probe the distribution uniformity and homogeneity of prime constituent elements in the livingstonite crystal. Figure 2(d)-(g) shows the magnified view of the dark-field TEM image in the orange square box as well as the corresponding TEM-EDXS elemental maps, emphasizing that mercury (Hg), antimony (Sb) and sulfur (S) are the prime components which are homogenously distributed in livingstonite crystal.

Anisotropic Raman scattering and crystal axis determination of livingstonite flakes
Further, the exfoliated thin livingstonite crystals are investigated to understand the effect of prevalent structural anisotropy on the evolution of vibrational modes using polarization-resolved Raman scattering in two different polarization configurations. Figure 3(a) shows the collected Raman spectrum of the 127 nm-thick flake using a 632.8 nm He-Ne laser. The excitation beam is focused onto the flake using a 40× objective lens (NA = 0.65) and the back reflected signal is collected via the same objective lens. The signal is filtered out by an edge filter (Semrock, LP02-633RE-25) and coupled into a spectrometer (Horiba, iHR 550) for recording the Raman spectrum. The Raman spectrum of livingstonite crystal displays multiple distinct Raman peaks in the 65-500 cm −1 range, which are located at 78, 102, 120, 136, 183, 239, 288, 314, 337, 432, and 495 cm −1 . The Raman modes of livingstonite can be correlated to the Raman spectra of its component binary sulfides including cinnabar HgS [32][33][34] and stibnite Sb 2 S 3 [35][36][37].
The dominant feature in the Raman spectra of livingstonite is a series of spectral bands that corresponds to the stretching and bending vibrations of Sb-S bonds and Hg-S bonds. The 78 cm −1 peak is assigned as the A g mode of Sb 2 S 3 (73 cm −1 ). The peak at 102 cm −1 is attributed to a combination of the E(TO) mode of HgS (103 cm −1 ) and the A g mode of Sb 2 S 3 (101 cm −1 ). The peaks at 120, 136 and 183 cm −1 are related to the S-Sb-S bending vibrations (140 and 184 cm −1 ). The 239 cm −1 peak represents the B 1g /B 3g mode of Sb 2 S 3 (238 cm −1 ). The 288 cm −1 peak belongs to the E(TO) mode of HgS (283 cm −1 ) and the A g mode of Sb 2 S 3 with the Sb-S stretching vibrations (283 cm −1 ). The peak of 314 cm −1 is assigned as the A g mode of Sb 2 S 3 with the Sb-S stretching vibrations (312 cm −1 ). The 337 cm −1 peak is attributed to the E(TO) mode of HgS corresponding to the Hg-S stretching vibrations (344 cm −1 ). The peaks of 432 and 495 cm −1 can be assigned to the Sb-S stretching vibrations (around 443 cm −1 ). In addition, all the identified Raman vibrational modes are summarized in Table 2. Subsequently, we quantity the Raman intensity variations of these vibrational modes by resolving in both parallel and perpendicular polarization configurations, where a linear polarization analyzer in the collection path is set to be in the parallel or perpendicular direction with respect to the input linear polarization. Figure 3(b) and (c) show the color maps of angle-resolved Raman spectra for the 127 nm-thick livingstonite crystal as a function of the incident linear polarization angle in the parallel and perpendicular polarization configurations, which confirm the presence of several A g modes at 78, 102, 120, 239, 288, 314, 432, 495 cm −1 and one B g mode at 183 cm −1 . It can be observed that the Raman modes are highly anisotropic in nature and exhibit periodic variations in the measured Raman intensity. Theoretically, the intensities of observed Raman modes are understood in the context of the incident (e i ) and scattered (e s ) beam polarization, and the corresponding Raman tensor (R) of each mode for the crystal system. Hence, the mathematical formulation of Raman intensity can be expressed as I ∝ |e i .R.e s | 2 . Under parallel polarization configuration, the associated unit polarization vectors would be same and can be expressed as e i = e s = (cosθ, sinθ, 0), whereas the scattered unit vector in perpendicular polarization configuration can be represented as e s = (−sinθ, cosθ, 0). Livingstonite crystal belongs to the monoclinic crystal family with space group A2/a, therefore the corresponding Raman tensor for A g and B g modes can be expressed as [38] R where a, b, c, d, e, and f are the Raman tensor elements. Considering all, the resultant expression for A g and B g Raman modes for parallel and perpendicular polarization configurations can be written as where // and ⊥ denote the parallel and perpendicular polarization configurations and ∅ ba = ∅ b − ∅ a denote the associated phase difference between a and b.  6) are represented with red solid lines. It can be observed that Raman A g modes show four-lobe patterns with two unequal maxima perpendicular to each other. The principal maxima for 78, 102, 120 cm −1 are aligned along x-axis (θ = 0°and 180°), while the principal maxima for 239, 288, 314, 432, and 495 cm −1 are along y-axis (θ = 90°a nd 270°). Based on the previous understanding of such crystal system, it can be inferred that θ = 0°and 90°should be either of the crystal axes of the investigated flake, which can be further utilized in probing the linear and nonlinear optical responses of livingstonite crystal. On the other hand, the Raman B g mode at 183 cm −1 (Fig. 3(g)) exhibits the characteristic isotropic four lobes pattern with the maxima at 45°, 135°, 225°and 315°. Furthermore, both A g and B g modes under perpendicular polarization configuration show the typical anisotropic four-lobe patterns with a period of 90° (Fig. 3(m)-(u)). The Raman intensity maxima of A g modes are observed at 45°, 135°, 225°and 315°, which match with the intensity minima of the B g mode at 183 cm −1 , emphasizing that the Raman intensity maxima of A g and B g modes cannot emerge at the same angle in two different polarization configurations.
Also, we investigate the thickness dependent effect on the observed Raman shifts and the anisotropy ratios of Raman A g modes. Figure 4(a) shows the recorded Raman spectra between 65-500 cm −1 for livingstonite crystals with five different thicknesses of 26, 38, 62, 82 and 127 nm. It shows there is no significant variations in the observed Raman shifts as a function of crystal thickness, which further validates the presence of vdW heterostructures in livingstonite crystal with constituent alternating layers of HgSb 2 S 4 and Sb 2 S 4 . Afterwards, we investigate the effect of flake thickness on the anisotropy ratios of Raman Ag modes, as shown in Fig. 4(b). It is worth noting that the observed anisotropy nature is consistent in all the flakes of different thicknesses, however the value of anisotropy ratio significantly depends on the flake thickness. The anisotropy ratio systematically increases up to 82 nm flake thickness and subsequently decreases, underlining that the anisotropy in phonon vibration in livingstonite crystal is sensitive to the number of heterostructure layers which influences the interactions between the lattice phonon vibrations and the flake interfaces.

Polarization-resolved linear optical properties and linear dichroism transition effect
Next, the effect of strong in-plane structural anisotropy on the appearance of intrinsic LD properties of livingstonite crystal is understood using polarization-resolved optical absorption microscopy in the visible range from 400 nm to 800 nm. The investigated crystal is illuminated with a broadband white light beam (Thorlabs, SLS201L) passing through a linear polarizer and a half-wave plate to set the incident light polarization. The incident beam is focused on the sample using an 80× objective lens (NA = 0.5). To measure the reflectance (R) spectra, the back-reflected light is collected using the same objective lens and routed to the spectrometer using a beam splitter. The reflected signal from the sample is normalized with the light source spectrum reflected from a silver mirror at the sample stage. In transmission mode, the transmitted light is collected with another 100× objective lens (NA = 0.7) and routed to the spectrometer. To measure the transmittance (T) spectra, the transmitted signal through the sample is normalized with the light source spectrum by removing the sample from the collection path. The absorbance (A) is then obtained by using A = 1 -R -T. Similarly, the reflectance spectra also show the initial declining trend with a dip at 427 ± 7 nm and upsurge afterwards. Nevertheless, in contrast to the transmittance spectra, the upsurge trend in the reflectance spectra makes a polarity switching after certain wavelength around 547 nm ± 3 nm. This observation strongly influences the resultant absorbance spectra of the crystal. Consequently, the absorbance spectra initially increase up to 445 nm ± 4 nm and afterwards show declining trend until 800 nm. The absorbance spectra sharply decline up to 582 nm ± 5 nm followed by gradual decreases with the reversed polarity in absorbance. To gain further insight about the anisotropic features and the wavelength-dependent polarity reversibility in the probed crystal, the angular plots are shown for transmittance, reflectance, and absorbance in Fig. 5(d)-(f), 5(g)-(i), and 5(j)-(l) respectively, at three different wavelengths, which are chosen in three different regions (before, at, and after the polarity transition). Furthermore, these experimentally observed values are fitted with the sinusoidal function of the form α (θ) = α x cos 2 (θ) + α y sin 2 (θ), where α x and α y are the magnitudes along x-axis and y-axis, θ is the polarization angle between the incident linear polarization and the crystal axis. The experimental data are shown with black squares, while the theoretical fits are denoted with red solid line. The anisotropy ratios of transmittance, reflectance, and absorbance at two different wavelengths of 520 nm (750 nm) are obtained as T x /T y = 1.31 (1.12), R x /R y = 1.04 (0.88), and A x /A y = 0.80 (1.5), respectively. Noticeably, the angular plots for transmittance at 520, 582, and 750 nm in Fig. 5(d)-(f) show the maxima along the direction of the crystal's a-axis (x-axis), whereas the angular plots for reflectance and absorbance exhibit systematic polarity switching as a function of the wavelength. The angular plot for reflectance at 520 nm in Fig. 5(g) exhibits the maxima along the crystal's a-axis (x-axis), but at 750 nm in Fig. 5(i) the orientation is rotated by 90°along the crystal's b-axis (y-axis). Similarly, the angular plots for absorbance in Fig. 5(j) and (l) show the polarity switching with the 90°rotation in two different regimes. The angular plots in Fig. 5(h) and (k) demonstrate the isotropic nature, which validates our perception of the transition point. We can infer that livingstonite crystals anisotropically absorb incident photon energy along the directions of a-axis (x-axis) and b-axis (y-axis). Furthermore, this intriguing observation of LD conversion effect facilitates additional insight related to such commensurately stacked vdW heterostructures, where the absorbance properties significantly depend on not only the incident linear polarization angles, but also the interacting beam wavelength. Hence, the LD reverses its polarity based on the excitation energy due to the presence of singularities in joint density of states along both a-axis and b-axis directions of the probed livingstonite crystal. Importantly, such wavelength-dependent LD conversion effect is reported for the first time in vdW heterostructures, however it has been recently observed in the context of complex vdW materials and explained in terms of the forceful localization rules at different parallel energy bands presented in the materials [38][39][40]. Owing to the insight, this process is a result of the strong localization of absorption peaks along two principal crystal axis directions in different parallel energy bands at different wavelengths in the livingstonite crystal. Nevertheless, this process is still elusive and further investigation is vital for the insightful comprehension about commensurate and incommensurate vdW heterostructures.

Anisotropic optical refraction of livingstonite crystal
The presence of strong linear dichroism further endorses the occurrence of birefringent phenomena in the probed livingstonite crystal. We investigate the birefringence effect in the 127 nm-thick livingstonite flake using cross-polarized optical microscopy. The flake is illuminated with a broadband white light source (Thorlabs, SLS201L) using a 20× objective lens (NA = 0.42) after passing through a linear polarization polarizer and a half-wave plate to control the incident linear polarization. The transmitted light through the flake is collected using another 50× objective lens (NA = 0.42). The signal is then passed through a linear polarization analyzer and routed towards an imaging camera to capture the optical transmission images from the flake. The polarization of the incident white light is rotated using the half-wave plate while maintaining the cross-polarization configuration by rotating the analyzer upholding the 90°difference. Other parameters such as exposure time and incident beam intensity are kept same while recording the images. The linear polarization angle of incident beam is varied from 0°to 180°with respect to the crystal's a-axis (x-axis). Figure 6(a) shows the captured optical transmission images from the flake at the incident polarization angles of 45°and 90°, showing that the higher brightness through the flake is obtained at 45°compared to that at 90°. We further observe that the brightness of the crystal exhibits systematic variation as a function of incident polarization angle and goes to the maximum when the incident polarization is aligned at 45°with respect to the crystal's a-axis and b-axis. On the contrary, the flake brightness shows the minimum when the incident polarization is aligned along any of the crystal axis. To further visualize the anisotropic optical refraction in the crystal, the recorded images are processed to filter out the background contribution by normalizing the background intensity. The relative brightness contrast between the probed livingstonite flake and the glass substrate is extracted at different incident polarizations. Finally, the normalized intensity (I − I min )/(I max − I min ) is plotted as a function of the incident linear polarization angle in Fig. 6(b) under the cross-polarization illumination configuration. The experimentally extracted values are indicated with black squares. The observed trend can be well fitted with the equation T = t sin 2 (2θ), where t denotes the transmittance along the crystal's axis. The fitted trend is shown with red solid line, which emphasizes that the measured flake brightness exhibits a four-lobe patterns with a period of 90°, reaching the maxima at 45°and 135°. This periodic variation in flake brightness for livingstonite crystal is ascribed as birefringence effect, which can be explained in the context of the outcoupled polarization state of transmitted light [41,42]. When the incident light polarization is aligned along any crystal axis (a-axis or b-axis), the outcoupled polarization state of transmitted light remains unchanged, where the polarization of transmitted light is same as the input beam polarization. Therefore, the transmitted light cannot pass through the linear polarization analyzer with cross-polarization configuration engaged in the collection path. As a result, the transmitted optical image appears dark. In contrast, when the incident beam polarization is away from the crystal axes, the transmitted light encounters the phase difference between x-and y-components. Consequently, the outcoupled polarization state of transmitted light is transformed to elliptical polarization so that the brightness contrast is visible in the transmitted image of the flake. The brightness contrast shows the maxima only along 45°and 135°, as demonstrated in our measurements.

Anisotropic nonlinear optical properties and determination of third-order nonlinear susceptibility of livingstonite crystal
The prevalent in-plane structural anisotropy and strong LD effect also endorse the presence of highly anisotropic nonlinear optical response. The anisotropic THG emission in livingstonite crystals are probed with a pulsed laser source at the wavelength of 1560 nm with a pulse width of 90 fs and a repetition rate of 80 MHz. The beam spot size is 1.5 µm and the average pump power is 1.1 mW with 8.65 GW/cm 2 peak irradiance. The laser beam is focused on the flake using a 40× objective lens (NA = 0.65). The transmitted THG signal from the flake is collected using another 100× objective lens (NA = 0.7) and routed to the spectrometer (Horiba, iHR 550) for recording the THG spectrum. A shortpass filter is engaged in the collection path to block the transmitted laser light. Figure 7(a) shows the collected THG emission spectrum from the 127 nm-thick livingstonite flake. The peak emission is recorded at 520 nm, which resides exactly at one third of the excitation beam wavelength of 1560 nm. The inset shows the captured optical image in transmission mode with the green THG emission. The crystal axis is indicated with black solid arrows. To further ascertain the third-harmonic generation process, the THG emission power is plotted as a function of the input pump power in logarithmic scale and fitted with the cubic power law. Figure 7(b) shows the log-scale plot and the retrieved slope endorses the THG process. Next, the in-plane anisotropic THG emission is probed through recording the THG emission intensity by varying the input pump polarization. The incident beam polarization is controlled via engaging a combination of a linear polarizer and a rotating half-wave plate in the illumination path, while the collected THG emission intensity is resolved along parallel (θ = 0°; I (3ω) x ) and perpendicular (θ = 90°; I (3ω) y ) to the crystal's a-axis (x-axis) by introducing a linear analyzer in the collection channel. First the crystal axes are identified using angle-resolved polarized Raman measurements for all the crystals, where the Raman A g mode at 78 cm −1 provides the information regarding the crystal axes. Hence, the direction of linear analyzer is aligned accordingly with θ = 0°for the crystal's a-axis in performing the anisotropic THG emission measurements. Figure 7(c) shows the measured angular dependence of THG emission power with respect to the input pump beam polarization for the 127 nm-thick livingstonite crystal. Noticeably, the THG emission is highly anisotropic in nature, and the resultant THG emission has a butterfly-shaped pattern in contrast to the conventional four-lobe pattern. Furthermore, the maximum THG power is not along the crystal axes. To further ascertain the consistency in measurements, several other livingstonite flakes with different thicknesses of 26, 38, 62, and 82 nm are characterized in Fig. 7(d)-(g). The experimental data is shown with red squares (x-component), blue dots (y-component), and black upper triangles (total power). It can be inferred that the resultant butterfly-shaped pattern of THG emission is consistent in all thickness range.
Next, to gain further insight regarding such unique butterfly-shaped THG emission pattern and other salient features, the experimental observations are we corroborated with a theoretical model of the third-order nonlinear susceptibility χ (3) . Since the livingstonite crystal belongs to the monoclinic crystal family, the contracted form of χ (3) can be written as follows [43,44] where the first term in subscript 1, 2 and 3 denotes x, y, and z respectively, and the second term in subscript refers to the combination of three components as xxx yyy zzz Considering the excitation electric field at frequency ω with linear polarization angle θ relative to x-axis of the crystal, the general form of electric field can be represented as ⃗ E (ω) = x(|E 0 |cosθ) +ŷ(|E 0 |sinθ). Since the excitation polarization is restricted in x-y plane, tensor elements containing z-axis contribution can be safely neglected, which results as only four non-zero elements χ 11 , χ 18 , χ 22 , and χ 29 contributing to the in-plane THG emission process. The resultant THG intensity of x-and y-components can be expressed as x ∝ ( χ 11 cos 3 θ + 3 χ 18 cosθ sin 2 θ) 2 (8) These equations are further used to fit the experimental data in the angular dependence of THG emission. The theoretically fitted curves are shown as the solid curves with the corresponding colors in Fig. 7(c)-(g), underlining a good agreement with the measured data. Also, it is evident that the y-component of THG emission consists of four lobes and the resultant maxima of THG emission is not aligned with the crystal axis. It can be understood in terms of the relative contributions of χ (3) tensor elements [45]. The shape of THG emission pattern substantially transforms as per the change in the values of χ (3) tensor elements. Based on the extracted values from the x-component of THG emission, χ 11 is significantly higher than χ 18 . However, according to the y-component of THG emission, χ 22 is comparable with χ 29 . Hence, the substantial contribution of χ 29 in the THG y-component leads to the four-fold polarization dependent patterns in all the flakes, which results in the butterfly-shaped THG emission patterns. Figure 7(h) shows the plot of the extracted χ (3) tensor elements as a function of the flake thickness, underlining the consistency of our comprehension for the crystals. It can be observed that the χ (3) tensor elements are almost unchanged with different flake thicknesses, and the average relative magnitudes of χ (3) tensor elements are χ 11 : χ 18 : χ 22 : χ 29 = 1 : 0.45 : 0.49 : 0.38, which also endorses the intrinsic nonlinear optical anisotropy in livingstonite crystals. Moreover, the | χ 11 | 2 /| χ 22 | 2 ratio denotes the THG anisotropy ratio along two crystal axes as I (3ω) In last, the third-order nonlinear susceptibility of livingstonite crystals is estimated by investigating the dependence of THG emission on the flake thickness [46]. Figure 7(i) shows the recorded THG emission power as a function of the flake thickness ranging from 26 to 192 nm. The THG emission power is recorded at 1.1 mW with 8.65 GW/cm 2 peak irradiance. It can be observed that the recorded THG emission power systematically increases up to 0.61 pW for the 82 nm-thick flake and afterwards exponentially decreases as the flake thickness further increases. Such thickness-dependent variation in THG emission power can be explained in the context of two ongoing competitive processes of optical gain and loss. At lower flake thickness, the THG emission power is directly proportional to the square of the flake thickness and the contribution of optical absorption remains negligible. Hence, the THG emission exhibits the increasing trend up to 82 nm. Nevertheless, as the flake thickness increases, the strong optical absorption plays very crucial role in attenuating the THG emission power propagation through the flake. This trend is observed in the exponential decay of the THG signal for the flake thickness larger than 82 nm. Further, such exponential decay trend in THG emission power is fitted with the equation of P (3ω) (d) = Ad 2 exp (︂ − 4πk 3 d λ 3 )︂ , where A is a constant and d is the flake thickness, to extract the imaginary part of the refractive index k 3 at λ 3 = 520 nm for livingstonite crystal. The extracted value is found as k 3 = 0.99. In addition, χ (3) can be estimated using the following relation expressed as [46] | χ ( where n 1 and n 3 is the real part of the refractive index of livingstonite crystal at λ 1 = 1560 nm and λ 3 = 520 nm, and the other terms are associated with the excitation laser system such as repetition rate f rep = 80 MHz, beam spot size W= 1.5 µm, pulse width τ = 90 fs. Based on the reported literature, we consider the refractive index values of livingstonite crystal are around n 1 = n 3 = 3 [47]. Taking into consideration of all these parameters, the estimated value of χ (3) of livingstonite crystal is 5.47×10 −20 m 2 /V 2 , which has the similar order of magnitude as those of recently explored vdW materials and heterostructures.

Conclusion
We have prepared the livingstonite thin flakes through mechanical exfoliation method. The prepared samples are extensively characterized in the context of flake geometry, atomic arrangement, crystal structural anisotropy, and chemical composition. The estimated empirical formula of livingstonite crystal is Hg 0.94 Sb 4.91 S 8.00 . We have further studied the vibrational, linear and nonlinear optical properties of the prepared livingstonite crystals, including anisotropic Raman scattering, linear dichroism transition, birefringence effect, and anisotropic third-harmonic generation. We have shown that the Raman vibrational modes are anisotropic in nature, attributed to the reduced in-plane symmetry in livingstonite crystal. These observations are further validated with the presence of strong linear dichroism and birefringence effect. Interestingly, in contrast to conventional vdW heterostructures, livingstonite crystals show unique wavelength-dependent polarity switching termed as linear dichroism transition effect. Moreover, we have demonstrated highly anisotropic THG emission with butterfly-shaped patterns in these crystals. We understood the evolution of such THG emission patterns in the context of the reduced in-plane crystal symmetry and the unconventional contribution from the off-axis nonlinear susceptibility tensor element, which can further pave a way to precisely design the resultant THG emission patterns in vdW materials. With this hindsight, we anticipate that livingstonite crystal will be a valuable addition to the vdW heterostructure fraternity. Furthermore, the interesting characteristics demonstrated by this commensurate vdW heterostructure will not only be relevant in understanding other naturally and artificially grown heterostructures, but also facilitate a rich testbed for novel polarized optoelectronic devices, frequency converters, wavelength-controllable photodetectors, multi-spectral imaging technology, wavelength -division multiplexing, encrypted optical signal processing, and optical communication.