Soon after the advent of the laser, nonlinear optical effects were experimentally demonstrated, in which the second harmonic generation (SHG) was the first phenomenon to be observed (Franken et al., 1961). Since then, many applications have been accomplished in a diversity of research and technology areas. For instance, nonlinear optical phenomena are the major source of new wavelength range pulses obtained from changing the frequency of pulsed femtosecond laser by interactions with nonlinear materials (Manzoni and Cerullo, 2016). In the last two decades, the development of user-friendly tabletop femtosecond Ti:Sapphire laser systems has led to the implementation of nonlinear optical microscopy by point scanning the pulsed laser beam over a material area and measuring the intensity map of the generated nonlinear signal. Nonlinear microscopy by harmonic generation, especially the second and the third harmonic generation (SHG and THG), four-wave mixing (FWM) and fluorescence by two-photon excitation (2PEF) have allowed to optically obtain the structural characterization of 2D materials (Malard et al., 2013; Oliveira et al., 2015) and biological materials (Zipfel et al., 2003; Ouellette et al., 2021; Cunha et al., 2021; Gomes et al., 2023). The overwhelming amount of publications on nonlinear properties of 2D nanomaterials make it hard to provide a comprehensive review of the literature. Thus, in this mini-review, we present and discuss how different nonlinear imaging techniques was employed to investigate and characterize 2D materials. The discussion is separated into sections of second and third-order nonlinear optical imaging.

The continuous unveiling of novel properties in 2D materials with great appeal for future technologies is followed by a significant effort to improve growth techniques to provide large-area flakes in a scalable production (Zhang et al., 2019; Aras et al., 2022). Hence, methods to rapidly and easily characterize the crystal quality of these grown samples can highly contribute to developing the synthesis of 2D materials. In particular, synthesized large-area flakes ordinarily exhibit polycrystalline domains (Lin et al., 2016), highlighting the necessity of mapping the crystallographic orientation over the sample. For instance, transmission electron microscopy (TEM) measurements precisely determine the sample’s crystallographic orientation (Huang et al., 2011; Yu et al., 2011; Ly et al., 2014). However, the mandatory use of a TEM grid and the limited measured area make such an experiment time-consuming, hampering its broad use. The SHG, which is a second-order nonlinear optical effect in which two incident fields with the same frequency ω generate a third field with frequency 2 ω (see Figure 1A), is also sensitive to the crystalline properties of a material, enabling the ascertainment of its orientation (Malard et al., 2013; Li et al., 2013). Hence, as the SHG is a non-invasive technique that can be carried out for samples in standard substrates and can rapidly provide images over large areas, it has been extensively used for identifying crystallographic orientation. The first results were obtained for a few layers MoS 2 samples (Malard et al., 2013) where it has been shown that polarization-resolved SHG measurements provide crystallographic information of MoS 2 monolayer and trilayer exfoliated flakes, Figures 1B–D. Figure 1B shows a six-fold pattern of the SH intensity plotted as a function of the sample rotation angle, where the sketch indicates a top view showing the angle of the MoS 2 crystallographic orientation regarding the incident laser polarization direction ( e ̂ ω ) . Figure 1C shows SH images for the monolayer and trilayer MoS 2 samples together with the crystallographic orientation of the crystal lattice obtained from the polarization data. The AFM image of the sample indicating the number of layers and a SH intensity profile taken along the flakes are shown in Figure 1D. Note that due to the symmetry dependence, the SH signal is observed only for the odd layer positions in the sample (the χ ( 2 ) = 0 for the even layers). The difference in intensity observed for the monolayer and trilayer is due to resonance effects with the excitonic transitions (Malard et al., 2013). The crystallographic orientation is obtained by describing the electric field of the generated SH light (E (2 ω )) along a given direction ( e ̂ 2 ω ) in terms of the second-order susceptibility tensor ( χ ( 2 ) ) and input light polarization vector ( e ̂ ω ) as (Boyd, 2008; Shen, 2003): E 2 ω ⋅ e ̂ 2 ω = C e ̂ 2 ω ⋅ χ ( 2 ) : e ̂ ω e ̂ ω (1) where ω is the laser frequency, 2 ω is the SH frequency and C is a proportionality constant that contains local field factors determined by the local dielectric environment. Equation 1 describes the polarization dependence of the SHG and its specific form will depend on the χ ( 2 ) tensor. For the TMDs, such as MoS 2 , the odd-layered samples have D 3 h point group symmetry, resulting in a second-order susceptibility tensor with a single non-zero element (Boyd, 2008; Shen, 2003): χ M o S 2 ( 2 ) ≡ χ x x x ( 2 ) = - χ xyy ( 2 ) = - χ yyx ( 2 ) = - χ yxy ( 2 ) , where x corresponds to the armchair direction since it has a mirror plane symmetry (Boyd, 2008; Shen, 2003), and y is the zigzag direction, as indicated in Figure 1B. Thus, the generated SH electric field as a function of the sample angle for a pump laser polarization ( e ̂ ω ) parallel to the analyzer ( e ̂ 2 ω ) , is (Boyd, 2008; Shen, 2003): E 2 ω = C χ M o S 2 2 ⁡ cos 3 ϕ + ϕ 0 (2) where ϕ is the angle between the input laser polarization and the x direction, and ϕ 0 is the initial crystallographic orientation of MoS 2 sample. As the SH intensity is phase insensitive, there is an arbitrariness of π / 3 in the definition of the x -axis.

Due to the sensitivity of SHG to the material’s crystalline properties, one-dimensional defects such as edges, grain boundaries, and wrinkles are examples of common features presented by 2D materials that can be probed by SHG imaging (Oliveira et al., 2015; Yin et al., 2014; Karvonen et al., 2017; Carvalho et al., 2019). For instance, polarization-dependent SHG results have shown the formation of crystallographically-oriented origami-type wrinkles in annealed hexagonal boron nitride (h-BN) layers (Oliveira et al., 2015), Figures 1E, F. Additionally, Cunha et al. (2020) revealed the important role of defects in the increased efficiency of the SHG in h-BN flakes. SHG has also been used to reveal the crystalline details of grain boundaries. While the single crystalline grains exhibit uniform SHG intensities, edges and grain boundaries present a suppressed SH emission due to their translational symmetry breaking, which results in destructive interference of the SH fields and allows rapid visualization of these defective regions (Yin et al., 2014). Dark-field SHG imaging can also probe edges and grain boundaries in TMD monolayers (Carvalho et al., 2019). This method consists of blocking the central spot of the SH signal to collect only the emission at high angles, leading to an enhanced SH intensity at edges and grain boundaries, as displayed in Figures 1G–J. While there is destructive interference of the SH fields at small angles for these defective regions, large angles compensate for the fields’ phase difference and result in constructive interference of the SHG at grain boundaries and edges (Carvalho et al., 2019). These angle-dependent constructive and destructive interferences are confirmed by the spatial mode of the SH emission of a MoSe 2 monolayer (Carvalho et al., 2019), as shown in Figures 1G, H. Additionally, grain boundaries of TMD monolayers can also be imaged by polarized SHG experiments (Van Der Zande et al., 2013; Cheng et al., 2015; David et al., 2015; Karvonen et al., 2017; Rosa et al., 2022; Sousa et al., 2024b). As the polarized SH intensity of these materials depends on the crystallographic orientation (Malard et al., 2013), grains with distinct orientations might display different polarized SH intensities, highlighting the boundaries between them (Van Der Zande et al., 2013; Cheng et al., 2015; David et al., 2015; Karvonen et al., 2017; Psilodimitrakopoulos et al., 2018; Rosa et al., 2022; Sousa et al., 2024b). Nonetheless, it is worth underlining that while TMD monolayers present a three-fold rotational symmetry, their polarized SHG exhibits a six-fold pattern, resulting in similar SH emissions for grains with anti-parallel crystallographic orientation.

In addition to the rich information directly given by the SHG imaging, it is also possible to employ data processing of the measured images to map distinct properties of the sample. For instance, a map of the crystallographic orientations of a TMD monolayer can be generated from polarization-resolved SHG imaging by fitting the angle-dependent SH expression (Equation 2) to the SH data of each measured pixel (David et al., 2015; Psilodimitrakopoulos et al., 2018). Besides, it was also reported a strain mapping of a TMD monolayer from polarization-resolved SHG imaging (Mennel et al., 2018; Li et al., 2019). As strain breaks the crystal symmetry of the material, it induces asymmetric intensities in the polarized SH six-fold pattern due to modifications in the second-order susceptibility tensor (Liang et al., 2017; Mennel et al., 2018; Mennel et al., 2019); Li et al., 2019). Hence, a photoelastic tensor was introduced to account for the role of the strain tensor on the nonlinear susceptibility tensor, allowing the evaluation of the strain field on the sample from polarization-resolved SHG measurements (Mennel et al., 2018; Mennel et al., 2019). Therefore, strain fields over a TMD monolayer can be mapped using this approach for each measured pixel of the polarized SHG images (Mennel et al., 2018).

Exploring the 2D materials field, we find that enormous efforts are also being placed in producing and investigating 2D heterostructures (Geim and Grigorieva, 2013; Novoselov et al., 2016; Castellanos-Gomez et al., 2022). For vertically stacked 2D materials, for example, novel and promising physical phenomena were revealed to be highly dependent on the twist angle between the forming layers (Castellanos-Gomez et al., 2022). Hence, SHG imaging emerges as a powerful non-invasive technique to determine the relative orientation between layers in 2D TMD heterostructures as well. This crystallographic orientation monitoring is simpler when the fabricated vertical heterostructure presents individual monolayers with non-overlapping areas (Hsu et al., 2014; Kim et al., 2021; Yuan et al., 2023). In such cases, it is possible to determine the crystallographic orientation of each monolayer from their isolated regions. However, as mentioned before, samples with θ or π / 3 − θ crystallographic orientations might exhibit the same polarized SH response, thus the SH signal from the heterostructure should also be analyzed to ascertain the monolayers’ relative orientation (Hsu et al., 2014; Psilodimitrakopoulos et al., 2020; Kim et al., 2021; Yuan et al., 2023; Palekar et al., 2024). The SH emission of a twisted TMD bilayer, for example, is the coherent superposition of the SH emission from each monolayer, in which the bilayer SH intensity ( I T ) is given by I T = I 1 + I 2 + 2 I 1 I 2 cos ( 3 ( θ 1 − θ 2 ) ) , with I 1 and I 2 the SH intensities from the individual monolayers and θ 1 and θ 2 their orientations (Hsu et al., 2014). Therefore, Δ θ and π / 3 − Δ θ relative orientations can be differentiated from the type of SH interference presented at the twisted bilayer region (Hsu et al., 2014), as displayed in Figure 1K.

Twist angle analysis through SHG imaging experiments is not straightforward when the heterostructure contains monolayers without isolated regions. For a twisted TMD homobilayer in which one monolayer presents an isolated region and the other is entirely overlapped, it is elementary to obtain I 1 , θ 1 , and I T . Since it is a homobilayer, it is expected that I 1 = I 2 ; thus θ 2 and the twist angle can be determined (Psilodimitrakopoulos et al., 2019; Xu et al., 2024). Note, however, that this approach cannot be applied to heterobilayers. In turn, the crystallographic orientations of the forming monolayers of an entirely overlapped TMD heterobilayer can be ascertained using polarization-resolved resonant SH imaging (Paradisanos et al., 2022). When one of the virtual states of the SHG process is in resonance with a real electronic or excitonic state, the SH intensity is highly increased (Malard et al., 2013; Wang et al., 2015b; Seyler et al., 2015; Wang et al., 2015a; Zhao et al., 2016; Lafeta et al., 2021; Shree et al., 2021; Sousa et al., 2024c). Thereby, if the excitation energy (or twice the excitation energy) is in resonance with an excitonic state of one of the TMD monolayers of the twisted bilayer, the SH emission of the heterostructure will be dominated by the signal from this specific monolayer and its crystallographic orientation can be accessed (Paradisanos et al., 2022). The twist angle of a TMD heterobilayer can be determined thus by performing polarization-resolved SHG imaging in resonance with each monolayer (Paradisanos et al., 2022). A real-time measurement of the stacking angle in a TMD heterobilayer was also reported by imaging two orthogonal components of the SH signal directly on the monolayers’ overlapped region (Psilodimitrakopoulos et al., 2020), as Δ θ can be written as a function of these perpendicular quantities from the interference SHG equation. Beyond the vertical heterostructures, SHG imaging has also been employed to probe coherent superposition effects in TMD lateral heterostructures (Sousa et al., 2021). The reported enhanced SH emission at the sharp TMD lateral heterojunctions due to a constructive interference revealed the capability of SHG imaging to probe these sharp interfaces (Sousa et al., 2021), as shown in Figure 1L. Moreover, such an experiment also allows the exploration of phase differences arising from the distinct second-order susceptibilities between the TMD monolayers (Kim et al., 2020; Sousa et al., 2021).

Also regarding 2D heterostructures, it has been recently demonstrated that the stacking of TMD monolayers can be used to boost the nonlinear optical gain and reach strong SHG signals (Trovatello et al., 2021), achieving optical parametric amplification in 2D TMDs. Furthermore, high nonlinear conversion efficiencies were obtained by controlling the phase-matching in multilayer 3R- MoS 2 with large thickness (Xu et al., 2022). Finally, SHG imaging was also employed to directly probe an incommensurate to commensurate phase transition in graphene/h-BN vertical heterostructures (Stepanov et al., 2020). While the centrosymmetric character of graphene is not affected by the h-BN bottom layer for an incommensurate alignment, the transition to a commensurate phase results in the breaking of graphene inversion symmetry and thus a finite SH emission (Stepanov et al., 2020).

TMD monolayers also exhibit singular valley phenomena (Xiao et al., 2012; Mak et al., 2012; Cao et al., 2012; Sousa et al., 2024a) due to their strong spin-orbit coupling and broken inversion symmetry, leading to a spin-valley locking effect at the inequivalent K and K′ valleys. This coupling between spin and valley degrees of freedom, which can be selectively accessed by circularly polarized light (Mak et al., 2012; Cao et al., 2012), sheds light on a novel field with promising potential for data storage, manipulation, and readout named valleytronics (Schaibley et al., 2016). Hence, control over valley polarization is essential in view of technological applications, resulting in a demand for methods to probe valley asymmetries in 2D materials. In this sense, a SHG imaging technique was recently used to map the valley imbalance in TMD monolayers (Mouchliadis et al., 2021). A valley polarization out of equilibrium impacts the crystal symmetry, thus inducing new terms in the material’s second-order susceptibility tensor and modifying the SH emission (Hipolito and Pereira, 2017). Particularly, the intrinsic ( χ int ( 2 ) ) and valley-polarization-induced ( χ vp ( 2 ) ) second-order susceptibility terms for TMD monolayers lead to orthogonal SHG polarizations, resulting in a rotation of the polarized SH six-fold pattern dependent on the valley imbalance (Ho et al., 2020; Mouchliadis et al., 2021; Herrmann et al., 2023). Therefore, by using an elliptical polarized excitation and detecting the polarization-resolved SH emission, it is possible to create a valley imbalance in TMD monolayers and indirectly track it through the rotation of the polarized SH pattern (Ho et al., 2020; Mouchliadis et al., 2021; Herrmann et al., 2023). As valley imbalance is mainly generated by a valley-polarized exciton population in TMD monolayers, resonant excitation also plays a major role in this phenomenon (Ho et al., 2020). However, it is worth stressing that further effects such as valley-exclusive optical Stark or Bloch-Siegert shifts can also induce valley polarization by breaking the time-reversal symmetry in TMD monolayers (Sie et al., 2017), which was demonstrated through SHG as well (Herrmann et al., 2023). Notably, these valley asymmetries in TMD monolayers are widely investigated by circularly polarized photoluminescence and absorption experiments (Mak et al., 2012; Cao et al., 2012). Nonetheless, these measurements cannot be used for gapless materials such as graphene, which also present valley phenomena (Rycerz et al., 2007; Xiao et al., 2007; Yao et al., 2008). In contrast, it is possible to measure valley population imbalances even in centrosymmetric crystals as graphene (Golub and Tarasenko, 2014), since this valley asymmetry lowers the crystal symmetry leading to non-zero terms in its second-order susceptibility.

Beyond the largely investigated 2D materials like graphene, h-BN, and TMDs, SHG has also been employed to probe magnetic symmetries in emergent layered magnets such as CrI 3 , CrSBr, and MnPS 3 (Wu et al., 2024). Although these materials are centrosymmetric in the paramagnetic phase, their antiferromagnetic ordering leads to an inversion symmetry breaking. Hence, an enhancement of the SH emission is observed below the Néel temperature, providing an all-optical and rapid method to probe magnetic phase transitions in these materials (Sun et al., 2019; Chu et al., 2020; Lee et al., 2021). In particular, SHG imaging was used to map this magnetic phase transition on a CrI 3 bilayer (Sun et al., 2019).

Despite the standard polarized SHG techniques presented above as consolidated and powerful tools to determine the orientation and symmetry properties of 2D materials, other innovative approaches using SHG improved and extended this second-order nonlinear technique to determine other properties and explore different materials. An example of this recent development is the Fourier space SH imaging, which provides a straightforward method to determine the crystallographic orientation and symmetry based on the images generated by the SH signal in Fourier space using an azimuthal laser mode to excite the sample (Lafeta et al., 2025). In addition, there is also a great interest in developing the potential of SHG techniques to improve their spatial resolution. Although considerable efforts are employed to achieve increased spatial resolutions for SHG imaging (Psilodimitrakopoulos et al., 2018), there is the intrinsic confocal optics resolution limitation of hundreds of nm due to the diffraction limit of light. Therefore, nonlinear near-field techniques emerge as a powerful alternative for imaging the nanoscale SH responses (Yao et al., 2022; Luo et al., 2023). For example, localized variations of the stacking order in a TMD homobilayer as well as excitonic resonances were recently probed by nano-SHG imaging with a spatial resolution down to 20 nm (Yao et al., 2022). Moreover, nano-SHG imaging could also detect localized variations in the symmetry of a TMD monolayer (Luo et al., 2023).

Similar to second-order optical effects, third-order nonlinear processes play an important role and significantly contribute to the study of 2D materials. The general third-order nonlinear optical phenomenon is known as four-wave mixing (FWM), which depends on the third-order electrical susceptibility ( χ ( 3 ) ) and can be generated by any combination of three electric fields of frequencies ω 1 , ω 2 , and ω 3 generating a fourth one of frequency ω FWM = | ± ω 1 ± ω 2 ± ω 3 | in a nonlinear medium. This phenomenon has exciting applications in 2D materials, mainly because it could be applied to centrosymmetric materials, contrasting with second-order nonlinear processes limited to noncentrosymmetric materials (Boyd, 2008). The FWM presents some important particular cases, such as the degenerated four-wave mixing (DFWM), in which two of the incident fields are frequency degenerate, for example, ω 1 = ω 3 , resulting in ω DFWM = 2 ω 1 − ω 2 (see Figure 2A) (Lafeta et al., 2021). Another special case commonly exploited in the 2D materials investigation is the third harmonic generation (THG), which occurs when the three incident fields have the same frequency ( ω 1 = ω 2 = ω 3 ) ; thus ω THG = 3 ω 1 (Wen et al., 2019).

Previous studies have explored the dependence of polarization in exceptional cases of FWM. In particular, Wang et al. (2014); Woodward et al. (2016) investigated the polarization dependence of THG, showing that third-harmonic signals from D 3 h symmetry materials align with the polarization of the three incident fields when they share the same linear polarization. This result was reproduced by Rosa et al. (2018), as illustrated in Figure 2B. The polarization dependence of DFWM, which exhibits additional degrees of freedom compared to THG, has also been reported for D 3 h symmetry materials. Li et al. (2016) examined DFWM signals for three polarization configurations: rotating only ω 1 incident field, rotating only ω 2 field, and rotating both ω 1 and ω 2 fields. However, they noted that when the experiments rotated the ω 2 field, the results were inconclusive due to the broadband nature of the ω 2 beam they used. Later, Balla et al. (2018) showed that when the polarizations of both incident beams are parallel and rotated with a fixed analyzer, the DFWM signal behaves like THG, aligning with the polarization of the incident fields. Dai et al. (2021) corroborated this, using lasers with parallel and fixed polarizations and rotating a polarizer in front of the detector. Furthermore, Dai et al. (2021) reproduced the first condition presented by Li et al. (2016) rotating ω 1 with ω 2 fixed, indicating that the rotation of the broadband ω 2 beams probably caused earlier inconsistencies. Dai et al. (2021) further demonstrated enhanced DFWM signals using hybrid MoS 2 -plasmonic structures of varying sizes, as shown in Figure 2F. They also examined DFWM polarization dependence by rotating both incident beams in pristine and hybrid MoS 2 -plasmonic structures. Maximum signals were observed when incident fields were aligned parallel to the longitudinal direction of the cavity.

In general, TMDs present a strong FWM signal that increases with the number of layers (Li et al., 2016; Säynätjoki et al., 2017; Rosa et al., 2018; Autere et al., 2018; Balla et al., 2018; Lange et al., 2024), as shown in Figures 2C–E, H–J. This FWM response complements the second-order techniques that cannot probe centrosymmetric materials such as even-layered TMDs due to their D 3 d symmetry (Malard et al., 2013; Li et al., 2013). Furthermore, resonances with excitons were also observed by FWM experiments when incident fields lead to a ω FWM that matches the excitons frequency in TMDs (Lafeta et al., 2021), demonstrating the potential of this technique to characterize, for example, the impact of doping on the excitonic responses of TMD materials (Sousa et al., 2024a; Sousa et al., 2024c). Bauer et al. (2022) also studied the effect of excitonic resonances on the dynamics of second- and third-order effects on 2D TMDs, revealing that temporal delays between incident fields can be required to maximize the resonant nonlinear optical responses of TMDs. Besides, Karvonen et al. (2017) used THG imaging to investigate grain boundaries in synthesized MoS 2 monolayers, showing an enhanced contrast of this one-dimensional defect compared with SHG imaging; thus enabling its rapid visualization.

Moreover, FWM techniques using broadband lasers were also used to investigate TMDs. This method employs broadband incident laser pulses to generate a combination of third-order nonlinear responses over a large range of frequencies, which enables the investigation of several effects such as resonances with excitons in one single measurement (Ko et al., 2019) and the increase in the nonlinear signal produced by the control of the phase of the incident fields (Lange et al., 2024).

FWM imaging has also been employed to study other 2D materials such as graphene and h-BN, with a special interest in monolayer graphene because of its inversion symmetry and consequent absent SHG signal (Shan et al., 2018). Such investigations in these materials exploit a particular case of DFWM called coherent anti-Stokes Raman scattering (CARS), which occurs when the energy difference of the incident fields ( ω 1 and ω 2 ) is equal to the energy of a vibrational mode ( ω Ω ) of the material, i.e., ω Ω = ω 1 − ω 2 . This condition is a special type of stimulated Raman process that amplifies the anti-Stokes signal (Cheng and Xie, 2004). The CARS technique provided remarkable results in the studies of biological materials and more recently in carbon nanotubes (Duncan et al., 1982; Cheng and Xie, 2004; Potma and Xie, 2004; Polli et al., 2018; Li et al., 2020; Cunha et al., 2021; Paddubskaya et al., 2020; Gordeev et al., 2023). Since CARS is a label-free technique that uses just the Raman fingerprints as a marker, several works explore such experiments in biological applications. Nevertheless, CARS is also used to study 2D materials (Malard et al., 2021). For instance, an anomalous nonlinear behavior was reported for graphene when the DFWM satisfies the resonant CARS condition. Lafeta et al. (2017) showed an expected strong enhancement in the h-BN DFWM signal at the CARS resonance, while a decreased DFWM intensity at the CARS condition was observed for graphene, which was explained with a Fano resonance model (Fano, 1961). A posterior study confirmed this anomalous behavior and exploited the impact of applying a temporal delay between ω 1 and ω 2 incident fields, unveiling the temporal dynamics between the electronic transition and the vibrational resonant lifetimes (Virga et al., 2019), as shown in Figure 2G. Another work involving FWM and stimulated Raman scattering (SRS) imaging demonstrated the capacity of the latter to remove the nonresonant background in contrast to CARS in h-BN samples (Ling et al., 2019).

In this mini-review, we explored the progress of nonlinear optical imaging techniques for the study of nonlinear responses in 2D materials, focusing on second-order, specifically SHG, and third-order, FWM and its special cases. In particular, we discussed the potential of SHG experiments to probe symmetry information of h-BN and odd-layered TMD samples. This capability of SHG measurements to image the crystallographic orientations of these materials enables, for example, the study of one-dimensional defects such as wrinkles, edges, and grain boundaries, the determination of twist angles in 2D heterostructures, and the mapping of strain fields over the samples. Additionally, we discussed the importance of FWM experiments to investigate nonlinearities of centrosymmetric materials as graphene and even-layered TMDs, which present a negligible second-order nonlinear emission. Furthermore, we also addressed the potential of resonant SHG and FWM experiments to probe excitonic effects in TMDs. In summary, we highlighted that these non-invasive nonlinear techniques are powerful tools for the rapid imaging of 2D material properties, paving the way for the development of sample fabrication and the research of fundamental aspects of nonlinear optics. Moreover, we shed light on different possibilities to increase the potential of nonlinear imaging. For instance, the recent use of near-field SHG and FWM experiments allows the imaging of 2D materials with nanometric resolution, while the constant developments in nonlinear optics theory and data analysis enable to unveil notable properties as valley population imbalance in these samples.