Second harmonic spectroscopy to optically detect valley polarization in 2D materials

TL;DR

This paper proposes a theoretical method using second harmonic spectroscopy to detect and quantify valley polarization in 2D materials, enabling contactless valley transport monitoring.

Contribution

It introduces a theoretical framework for using second harmonic generation to measure valley polarization in inversion-symmetric 2D materials, distinguishing intrinsic and valley contributions.

Findings

SH response depends on valley polarization, polarization, temperature, and electron density.

Proposed experimental setup allows direct quantitative measurement of valley polarization.

Method enables remote, contactless valley polarization mapping.

Abstract

Valley polarization (VP), an induced imbalance in the populations of a multi-valley electronic system, allows emission of second harmonic (SH) light even in centrosymmetric crystals such as graphene. Whereas in systems such as MoS$\mathrm{_2}$ or BN this adds to their intrinsic quadratic response, SH generation in a multi-valley inversion-symmetric crystal can provide a direct measure of valley polarization. By computing the nonlinear response and characterizing theoretically the respective SH as a function of polarization, temperature, electron density, and degree of VP, we demonstrate the possibility of disentangling and individually quantifying the intrinsic and valley contributions to the SH. A specific experimental setup is proposed to obtain direct quantitative information about the degree of VP and allow its remote mapping. This approach could prove useful for direct, contactless, real-space monitoring of valley injection and other applications of valley transport and valleytronics.

Equations

$$j_1^{(2)}({\omega }_1,{\omega }_2)$$

$$j_2^{(2)}({\omega }_1,{\omega }_2)$$

$$f_1$$

$$f_2$$

$$I/I_0$$

$$-8i({\bar{\sigma }}_1{\bar{\sigma }}_2^{\ast }-{\bar{\sigma }}_1^{\ast }{\bar{\sigma }}_2)\sin (2\phi ),$$

$$I_{\parallel }/I_0$$

$$\displaystyle+4|\bar{\sigma}_{2}|^{2}\big{[}\!\cos^{2}(3\zeta\!+\!2\varphi)\sin^{2}(2\varphi)\!+\!\sin^{2}(3\zeta\!+\!2\varphi)\big{]}$$

$$-2({\bar{\sigma }}_1{\bar{\sigma }}_2^{\ast }+{\bar{\sigma }}_1^{\ast }{\bar{\sigma }}_2)\sin (6\zeta +4\phi ){\cos }^2(2\phi )$$

$$-4i({\bar{\sigma }}_1{\bar{\sigma }}_2^{\ast }-{\bar{\sigma }}_1^{\ast }{\bar{\sigma }}_2)\sin (2\phi ),$$

Full text

I Introduction

Interactions of light with matter beyond linear response are a rich source of fundamentally interesting phenomena as well as of many-fold opportunities for applications Jackson (1999); Born and Wolf (1986); Shen (2002); Boyd (2008). In particular, the use of non-linear optical spectroscopy for characterizing the electronic properties of crystalline materials has emerged as a fruitful, simple and important technique because, among other advantages, it allows fast, non-invasive probing of electronic systems and is sensitive to intermediate coherent electronic transitions Shen (2002); Boyd (2008). Compared to linear optical absorption, for example, in non-linear optical spectroscopy there is a larger freedom in utilizing the expanded set of selection rules and conditions involving the polarization dependence or polarization state, in order to extract more microscopic details with the same type of measurement Ivchenko (2005); McIver et al. (2011). We can understand this in the simplest and most general way by recalling that the $n$–th order response is governed by a $(n+1)$–rank tensor and that, for a given crystalline symmetry, the number of independent optical constants increases with the order Haussühl (2007). Therefore, since a single frequency measurement at higher order of response can capture a larger number of independent quantities of the system, it more strongly constrains its microscopic details (e.g. the modeling of its bandstructure) while, at the same time, becomes a more versatile approach that is capable of probing a richer set of phenomena. This has a large potential for characterization and applications and, consequently, is of high interest.

Two-dimensional crystals such as graphene, boron nitride, and transition-metal dichalcogenides (TMD) have been shown to have a particularly strong non-linear optical response, especially given their atomically thin character. In the latter, second harmonic generation (SHG) is particularly robust Ivchenko (2005); McIver et al. (2011); Wu et al. (2012) and is routinely used for simple characterization tasks such as identifying crystal orientation, uniformity, or layer number Shen (1989); Dean and van Driel (2009, 2010); McIver et al. (2012); Malard et al. (2013); Kumar et al. (2013); Li et al. (2013); Hsu et al. (2014); Zhou et al. (2015). With a setup that allows for translation of the beam along the sample, it becomes possible to spatially map the SHG by probing the sample in scanning mode with resolution limited by the spot size Malard et al. (2013); Kumar et al. (2013); Zhou et al. (2015). In graphene, on the other hand, SHG is forbidden in equilibrium by its $D_{6h}$ point group symmetry (PGS). As discussed below, the vanishing quadratic response in graphene arises at the microscopic level from the exact cancellation of finite contributions from the $\mathbf{K}$ and $\mathbf{K}^{\prime}$ valleys.

Previous theoretical calculations, indicate that disrupting the valley cancellation by population imbalance can lead to a finite SH with estimated magnitudes on par with conventional nonlinear crystals Golub and Tarasenko (2014); Wehling et al. (2015), which means that VP is expected to generate a very strong non-linear signal. Golub and Tarasenko Golub and Tarasenko (2014) calculated the frequency dependent SH susceptibility from an explicit integration of the time-dependent density matrix using an effective Dirac description of the electronic structure of graphene; trigonal warping is explicitly included both in the effective Hamiltonian and in the coupling to light. Wheling et al. Wehling et al. (2015) computed the SH optical conductivity from a diagrammatic expansion of the current response to second order, thus expressing all quantities in terms of Green’s functions. Their description is based on a tight-binding (TB) formulation of the electronic problem in graphene, including the TB derivation of the generalized coupling to light and the velocity operator in higher orders from a Peierls substitution in the hoppings. The different methods and approximations used in these references lead to not entirely compatible results. Moreover the scope of these calculations is limited by other strict approximations, e.g. zero temperature, very small VP and, above all, pertain to (gapless) graphene only.

Therefore, it becomes important to address the more general and rich problem of valley-induced SHG in 2D materials beyond graphene using an approach capable of addressing a more general set of conditions, such as finite temperature, variable carrier density, and polarization of the incoming and outgoing radiation. From the technical point of view, we address these characteristics (see §III) within the length gauge formalism Aversa and Sipe (1995); Hipolito et al. (2016). To be encompassing and allow a controlled breaking of inversion symmetry, we study the quadratic optical response to light of a generic two-band electronic system on a honeycomb lattice. This allows us to quantify the interplay between intrinsic contributions and those induced by a finite VP as a function of frequency and polarization. At the qualitative level this choice allows us to interpolate between the behavior of graphene (gapless) and that of semiconducting TMD (gapped).

Since much effort is currently invested to theoretically and experimentally develop methods and concepts to harness the valley degree of freedom in these and related systems for valleytronic applications Rycerz et al. (2007); Zeng et al. (2012); Xiao et al. (2012); Golub et al. (2011); Mak et al. (2012); Jiang et al. (2013); Mak et al. (2014); Yu et al. (2014); Jiang et al. (2014); Lee et al. (2016), it is important to establish practical, versatile and reliable probes able to quantify and track the degree of valley polarization (VP), just as, in spintronics, it is crucial to have probes capable of quantifying spin polarization, injection, relaxation, etc.

Light is a demonstrably effective means of inducing a VP in these materials Mak et al. (2012); Zeng et al. (2012); Cao et al. (2012) and, here, we discuss a specific proposal of its utilization as an effective qualitative and quantitative probe as well through SH spectroscopy. Even though SHG in otherwise SH-dark graphene provides direct access to the degree of an induced VP, our discussion extends from this case to systems with intrinsic SHG where the two effects can be present and contribute to independent second-order optical constants. To be encompassing and allow a controlled breaking of inversion symmetry, we do that by studying the quadratic optical response to light of a generic two-band electronic system on a honeycomb lattice. This allows us to quantify the interplay between intrinsic contributions and those induced by a finite VP as a function of frequency and polarization. As the key underlying physics is not dependent on specific microscopic details other than the crystal symmetry, we begin by discussing a specific experimental procedure that should allow one to use SHG as a useful probe in valleytronics and, subsequently, analyze the microscopic details of SHG in the framework discussed above.

II Fingerprint of valley polarization in SHG

Threefold rotational symmetry severely restricts the in-plane components of the quadratic conductivity that obey $\sigma_{\alpha\beta\beta}{\,=\,}\sigma_{\beta\alpha\beta}{\,=\,}\sigma_{\beta\beta\alpha}{\,=\,}{\,-\,}\sigma_{\alpha\alpha\alpha}$, where $\alpha,\beta{\,\in\,}\{1,2\}$ Haussühl (2007). Alone, this symmetry reduces the number of independent components to just $\sigma_{111}$ and $\sigma_{222}$ which, for simplicity, we shall replace by the dimensionless counterparts $\bar{\sigma}_{i}{\,\equiv\,}\sigma_{iii}/\sigma_{0}$, where $\sigma_{0}{\,\equiv\,}e^{3}a/4\gamma_{0}\hbar$ sets the natural scale of the second order 2D conductivity (see below) Hipolito et al. (2016). Furthermore, in a honeycomb lattice in equilibrium whose mirror plane $\mathbf{e}_{m}$ is parallel to $\mathbf{e}_{2}$ [fig. 1(b)], only $\bar{\sigma}_{2}$ survives which defines the intrinsic quadratic response of the system. Since at frequencies much smaller that the bandwidth ($\omega{\,\ll\,}\gamma_{0}$) electronic processes are governed by states in the vicinity of the two inequivalent points $\mathbf{K}$ and $\mathbf{K}^{\prime}$ in the Brillouin zone (BZ), that intrinsic response is the combination of the contributions from each of these two valleys, which contribute independently (additively) in a translationally-invariant system. A crucial aspect, though, is that the PGS of $\mathbf{K}/\mathbf{K}^{\prime}$ is still $D_{3h}$ but with the mirror plane perpendicular to that of the real space lattice ($\mathbf{e}_{m}$). In other words, if taken independently, each valley contributes $\bar{\sigma}_{1}{\,\neq\,}0$, $\bar{\sigma}_{2}{\,=\,}0$ and it is their sum that yields an overall $\bar{\sigma}_{1}{\,=\,}0$, $\bar{\sigma}_{2}{\,\neq\,}0$ expected on symmetry grounds at equilibrium (in particular, if inversion symmetry is furthermore present as in graphene, the two valleys exactly cancel each other and $\bar{\sigma}_{1}{\,=\,}\bar{\sigma}_{2}{\,=\,}0$) Golub et al. (2011). However, if there exists an imbalance in the population of the valleys, there is no symmetry constraint to quench either $\bar{\sigma}_{1}$ or $\bar{\sigma}_{2}$; in particular, SHG can arise through $\bar{\sigma}_{1}{\,\neq\,}0$ in a lattice with inversion symmetry, which immediately suggests the detection of this valley-induced SHG response as a direct optical probe of VP. Generically, in a VP crystal without inversion both components will be present with a direct impact on the polarization dependence of the SHG that we explore here to disentangle and independently quantify the intrinsic ($\bar{\sigma}_{2}$) and valley-induced ($\bar{\sigma}_{1}$) SHG. We first derive the dependence of the SHG in the polarization state of the excitation field in general terms to establish the procedure for the individual component extraction, and afterwards analyze the frequency dependence of $\bar{\sigma}_{1}$ and $\bar{\sigma}_{2}$ in a microscopic model for graphene and for gapped graphene that applies qualitatively to the response in TMD.

We begin with the generic parameterization of an incoming monochromatic field $E_{0}$ normal to the sample, as illustrated in fig. 1(a), that is initially p-polarized before transmitting through a $\lambda/4$ plate with fast axis at an angle $\varphi$ with the plane of propagation. This permits the selection of any incoming polarization state, including linear polarization. For a general orientation of the propagation plane ($\zeta$) the (complex) amplitude of the electric field reaching the sample, $\mathbf{E}_{\omega}$, reads \mathbf{E}_{\omega}{\,=\,}E_{0}\,\big{(}a\sin\zeta+b\cos\zeta,\,-a\cos\zeta+b\sin\zeta,\,0\big{)} where $a{\,\equiv\,}i\sin(2\varphi)/\sqrt{2}$, $b{\,\equiv\,}[1{\,-\,}i\cos(2\varphi)]/\sqrt{2}$ and $\mathbf{e}_{1}$ is aligned with the lattice zigzag direction. The second order two-dimensional current density, $j_{i}^{(2)}(\omega_{1},\omega_{2})=\sum_{jk}\sigma_{ijk}^{(2)}E_{\omega_{1}}^{j}E_{\omega_{2}}^{k}$, can hence be written as

[TABLE]

where the auxiliary functions $f_{1}$ and $f_{2}$ read

[TABLE]

Even though we will be focusing on SHG arising from a single monochromatic source ($\omega_{1}{\,=\,}\omega_{2}{\,=\,}\omega$), we explicitly distinguish $\omega_{1}$ and $\omega_{2}$ to underline that our analysis applies to any second-order process. The sheet current eq. 1a radiates, in turn, an electromagnetic field with a flux density $I=\mu_{0}c\,|j^{(2)}(\omega_{1},\omega_{2})|^{2}/8$, or

[TABLE]

where $I_{0}{\,=\,}\mu_{0}c\,\sigma_{0}^{2}\;|E_{\omega_{1}}|^{2}|E_{\omega_{2}}|^{2}/32{\,=\,}(\mu_{0}c)^{3}\sigma_{0}^{2}I_{\omega_{1}}I_{\omega_{2}}/8$ $[\mathrm{W/m^{2}}]$. Whereas this shows that the total SHG intensity cannot discriminate the relative magnitudes of $\bar{\sigma}_{1}$ and $\bar{\sigma}_{2}$, that can be achieved by filtering the SH field with a linear polarizer parallel to the sample and rotated by $\xi$ with respect to $\mathbf{e}_{1}$ so that the electric field at the detector reads $\mathbf{E}_{\xi}{\,=\,}\mathbf{E}_{2\omega}{\,\cdot\,}(\cos\xi,\sin\xi,0)$. If the incoming light is linearly polarized parallel to the analizer ($\xi{\,=\,}\zeta$) the SH intensity at the detector reads

[TABLE]

while at cross orientation ($\xi{\,=\,}\zeta{\,+\,}90^{\circ}$) it is given by eq. II with the replacement $3\zeta{\,\to\,}3\zeta+90^{\circ}$.

Eq. II or any of its variants can thus be used to directly obtain $\sigma_{1}$ and $\sigma_{2}$ as well as the orientation of the lattice by fitting experimental SH intensities as a function of polarization. This is a concept similar to the usage of SHG as a remote, non-invasive probe of lattice orientation, layer number and other properties in recent applications of two-dimensional crystals having intrinsic SHG Shen (1989); Dean and van Driel (2009, 2010); McIver et al. (2012); Malard et al. (2013); Kumar et al. (2013); Li et al. (2013); Kumar et al. (2013); Hsu et al. (2014); Zhou et al. (2015). What we now propose and explicitly show is that, in addition, it follows from the general features of eq. II that the same concept can be applied to monitor and quantify the presence of VP, which is of clear interest for applications envisaged in the realm of valleytronics, somewhat similarly to the uses of the Kerr rotation to monitor and map spin accumulation in spintronics Kato (2004).

Fig. 2 illustrates the two simplest and extreme cases of entirely intrinsic (red) and entirely VP-induced SHG (black). For definiteness, we consider $\varphi{\,=\,}0$ so that the incoming light reaches the sample linearly polarized. As exactly anticipated from the earlier discussion on the differences between the PGS of the crystal and that of each valley, the angular pattern of $I_{\parallel}$ is rotated by 30*∘* between the two cases: VP in an inversion symmetric crystal ($\bar{\sigma}_{2}{\,=\,}0$) leads to SHG whose intensity is a direct measure of the degree of polarization $\mu_{\mathbf{K}}-\mu_{\mathbf{K}^{\prime}}$ (more details below), and its flower-shaped pattern directly reveals the PGS of the $\mathbf{K}$ points. The more general case of a system already having an intrinsic SHG ($\bar{\sigma}_{2}{\,\neq\,}0$) is illustrated by the blue and gray curves. They reveal that an emerging VP is signaled by three distinct features: (i) the progressive rotation of the flower pattern away from the principal directions set by the lattice orientation; (ii) the increase in intensity as the contributions arising from $\bar{\sigma}_{1}$ add to the intrinsic SH intensity, as per eq. II; (iii) the minimal intensity is no longer zero. Since the zero of the intrinsic SH response is usually well resolved experimentally Hsu et al. (2014), any of these effects can be used for qualitative monitoring of the degree of VP in the system, a fast alternative to fitting the angular patterns to eq. II when the actual magnitudes of $\bar{\sigma}_{1,2}$ are not required. Note that $\bar{\sigma}_{1,2}$ are complex quantities and, hence, the orientation and shape of the pattern is determined not just by their relative magnitudes but also the relative phase, as easily seen in the case $\bar{\sigma}_{1}{\,=\,}\bar{\sigma}_{2}e^{i\delta}$: $I_{\parallel}/I_{0}{\,=\,}2|\bar{\sigma}_{1}|^{2}[1+\cos\delta\,\sin(6\zeta)]$ (notably, when $\delta=\pm 90^{\circ}$ the six-fold pattern vanishes and becomes isotropic).

An obvious but important implication of these modifications is that the orientation of the SH intensity pattern does not correlate directly anymore with the lattice orientation in the presence of both intrinsic and VP-induced SHG. But this same fact can be utilized to detect and quantify both intrinsic and valley-induced conductivities. As $\bar{\sigma}_{1}$ is to leading order linear in $\delta\mu$, a reversal of the VP ($\delta\mu{\,\to\,}{\,-\,}\delta\mu$) changes its sign. On account of the cross-term in eq. II, this translates into a rotation of the pattern by 30*∘, which is equivalent to a reflection about the principal directions set by the lattice orientation, as shown in fig. 2. Consequently, the intersection of two patterns associated with opposite VP defines the orientation of the lattice modulo 30∘*. While this still doesn’t uniquely distinguish ZZ and AC directions, we note that a unique identification is possible whenever $|\bar{\sigma}_{1}|\neq|\bar{\sigma}_{2}|$ (which comprises essentially all cases) because the two non-equivalent intersections will then occur at different SH intensities, and it follows from eq. II that the intersection at lower (higher) intensity, highlighted with black (red) markers in the plot, occurs along the direction $\mathbf{e}_{1}{\,\Leftrightarrow\,}$ZZ ($\mathbf{e}_{2}{\,\Leftrightarrow\,}$AC) when $|\sigma_{2}|{\,>\,}|\sigma_{1}|$. Conversely, when $|\sigma_{2}|{\,<\,}|\sigma_{1}|$ the lower (higher) intersection occurs along $\mathbf{e}_{2}{\,\Leftrightarrow\,}$AC ($\mathbf{e}_{1}{\,\Leftrightarrow\,}$ZZ). This is clearly seen in fig. 2 where the reversal of $\delta\mu$ allows the immediate conclusion that the direction $\mathbf{e}_{1}$ corresponds to $\zeta=0$ because the two curves intersect there with the lowest intensity. These considerations are relevant not just because they illustrate how to use all the available information for a facile and expedite characterization of the nonlinear optical constants, but also because the success of a full nonlinear fit of an experimental trace of $I_{\parallel}$ vs $\zeta$ to eq. II can depend strongly on the assumed alignment of the lattice. Finally, it is clear from eq. II that, if the lattice orientation is known, measuring $I_{\parallel}/I_{0}$ at three non-equivalent orientations such as $\zeta{\,=\,}0^{\circ},15^{\circ},30^{\circ}$ suffices to uniquely determine the magnitudes of $\sigma_{1}$ and $\sigma_{2}$, as well as their relative phase.

The discussion so far was done for a linearly polarized excitation field ($\varphi=0$). An alternative consists in analyzing the SH signal as a function of the polarization state of the excitation field determined by $\varphi$, and which can be tuned continuously with the rotation of a $\lambda/4$ plate Ganichev and Prettl (2003); Ivchenko (2005); McIver et al. (2011). Since the roles of $\varphi$ and $\zeta$ are very much equivalent in eq. II, an analysis analogous to the one above can be straightforwardly done in this case. For example, with a fixed analyzer at $\xi{\,=\,}\xi_{zz}{\,\equiv\,}0$ ($\xi_{ac}{\,\equiv\,}30^{\circ}$), the intensity can still be read from eq. II with the replacement $3\zeta{\,\to\,}2\zeta$ ($3\zeta{\,\to\,}2\zeta+90^{\circ}$). Since the description of the $\varphi$-dependence is similar to the one above, we omit it for brevity.

III Frequency dependence of $\bar{\sigma}_{1}$ and

$\bar{\sigma}_{2}$

In order to determine the typical dependence of both $\bar{\sigma}_{1}$ and $\bar{\sigma}_{2}$ on excitation frequency and chemical potential for representative cases, we focus the analysis now on graphene-based systems, where recent reports have demonstrated the possibility of generating valley-polarized currents with high valley relaxation lengths, both in mono and bilayers Gorbachev et al. (2014); Shimazaki et al. (2015). The electronic degrees of freedom of a graphene monolayer are extremely well described by a single-orbital tight-binding (TB) model for electrons in the honeycomb lattice of fig. 1(b). This is a single (hopping) parameter model for graphene, which can only have $\bar{\sigma}_{1}{\,\neq\,}0$. In addition, in order to study the characteristics and relative magnitudes of $\bar{\sigma}_{1}$ and $\bar{\sigma}_{2}$ in a non inversion-symmetric system, it is desirable to have a model where inversion symmetry can be broken in a controlled way. That is easily incorporated in the single-orbital tight-binding via a sublattice potential $\pm\Delta/2$ which explicitly breaks the sublattice symmetry, and allows one to study the effects of VP in a more general “gapped graphene” setting. Whereas the case of graphene should be captured with good accuracy within this framework, the case of “gapped graphene” is expected to convey the main qualitative features expected in gapped systems such as in doped MoS2 (doping suppresses excitonic effects, and renders a single-particle description of the optical response appropriate). The second order conductivity tensor is computed perturbatively for a translationally-invariant system treating the interaction with light via the direct coupling, $\mathbf{r\cdot E}$, in the dipole approximation as described in references Aversa and Sipe (1995); Hipolito et al. (2016). We consider only the clean limit, but account phenomenologically for disorder broadening of the conductivity. Each component $\bar{\sigma}^{(2)}_{\lambda\alpha\beta}(\omega_{1},\omega_{2})$ is obtained from the formal result (25) of reference Hipolito et al. (2016).

Our results show explicitly that, as expected, the $\mathbf{k}$-space integration for small photon energies $\hbar\omega\leq\gamma_{0}$, is dominated by the vicinity of the K points. This allows us to use the equilibrium results to compute the contribution of each separate valley at different chemical potential by restricting the momentum integration to either of the shaded regions in fig. 1(c), while still working with the full tight-binding. Being able to keep the full tight-binding bandstructure rather than a Dirac-type approximation is important because, on the one hand, this allows us to immediately accommodate any refinement of the bandstructure model or straightforwardly extend the analysis to a different material. On the other hand, SHG in the clean limit arises only when the trigonal warping of the bands is explicitly considered Golub and Tarasenko (2014), which is guaranteed in the TB scheme.

Starting with the case of graphene, it is instructive to consider first the individual contribution of each valley. Since the point group symmetry of K has no inversion, each valley contributes a finite SHG (through $\bar{\sigma}_{1}{\,\neq\,}0$) but, in equilibrium ($\mu_{\mathbf{K}}=\mu_{\mathbf{K}^{\prime}}$), time-reversal symmetry forces an exact cancellation of each valley’s contribution and a system such as graphene has no intrinsic SHG (cf. dashed curves in the top panel of fig. 3). When a VP is induced as in fig. 3, there is no cancellation anymore and the overall effect at finite frequency is the appearance of two features at $\omega{\,=\,}\mu$ and $2\mu$ and a strong enhancement when $\omega\to 0$, consistent with two previous reports based on a related calculation in the Dirac approximation Golub and Tarasenko (2014); Wehling et al. (2015). The bottom panel shows the decomposition of $\bar{\sigma}_{1}$ in terms of the inter and intraband contributions defined in references Aversa and Sipe (1995); Hipolito et al. (2016) 111We follow the notation of Aversa and Sipe (1995); Hipolito et al. (2016), where $ee$ represents purely interband, $ie$ and $ei$ mixed inter-intraband processes, and $ii$ purely intraband processes. These results are in qualitative agreement in with Golub and Tarasenko (2014), but exhibit magnitudes ${\,\sim\,}500$-fold larger, as the intensity of the resonant features increases significantly at low temperature.. Whereas the behavior at $\hbar\omega{\,\sim\,}\mu$ and $2\mu$ is due to the resonant denominators coming from interband processes, the signal is much amplified towards the DC limit because the lack of inversion within each valley implies that when $\omega{\,\to\,}0$ the response is dominated by purely intraband transitions, since the conductivity terms $\sigma_{1}^{\mathrm{(ii)}}$ are now strictly finite (they cancel by symmetry in equilibrium Hipolito et al. (2016)). We note that, despite being at the same level of single-particle approximation, our results in fig. 3 disagree with the previous calculations of the SHG in valley-polarized graphene Golub and Tarasenko (2014); Wehling et al. (2015) (which, in turn, themselves disagree with each other 222With respect to Ref. Wehling et al. (2015), despite capturing the divergence ${\propto\,}1/\omega^{2}$ in the DC limit, our results are qualitatively different in the range $\mu\lesssim\hbar\omega\lesssim 2\mu$ in what regards the shape of the resonances at $\mu$ and $2\mu$, as well as the relative sign of $\bar{\sigma}_{1}$ in the two limits $\omega\to\,0,\,\infty$. Regarding the results of Ref. Golub and Tarasenko (2014), there is some qualitative agreement, but a detailed analysis shows that the results do not match quantitatively in the amplitude and shape of the resonances.). We attribute these differences to the long-standing problem of taking proper account of the intra-band contributions in the calculation of nonlinear response functions. Since a VP leads to explicitly finite intra-band terms even in the presence of a band-gap (see below), such contributions must be handled with care, which is addressed here in the framework of Aversa and Sipe that has been proven reliable in the DC limit Blount (1962); Sipe and Ghahramani (1993); Aversa and Sipe (1995).

The dependence of this valley-induced SHG on $T$ and $\mu$ is addressed in figure fig. 4. The features at $\omega{\,=\,}\mu$ and $2\mu$ are strongly temperature-dependent and disappear as soon as $k_{B}T{\,\gtrsim\,}\delta\mu$ (top panel) because, at this point, the temperature broadening of the Fermi-Dirac function whittles down the effective valley polarization. At fixed $T$, the response grows linearly with $\delta\mu$ Golub and Tarasenko (2014) when $\omega{\,\lesssim\,}\mu$, except near the resonances at $\omega{\,\sim\,}\mu$ and $2\mu$.

Setting $\Delta{\,\neq\,}0$ explicitly breaks the inversion symmetry of graphene and an intrinsic SHG ($\bar{\sigma}_{2}{\,\neq\,}0$) obtains. For definiteness, consider the case when $\Delta/2$ and $\mu$ are comparable, which we illustrate in fig. 5 for $\delta\mu{\,=\,}10$ meV at $T{\,=\,}50$ K. We see that lower frequencies $\omega\leq\mu{\,\sim\,}\Delta/2$ are dominated by the VP mechanism ($|\bar{\sigma}_{1}|{\,\gg\,}|\bar{\sigma}_{2}|$), while the intrinsic response dominates for $\omega{\,\gtrsim\,\mu}$. This happens because $\bar{\sigma}_{2}$ is Pauli-blocked at $\omega{\,\lesssim\,}\mu$ but, as seen above, $\bar{\sigma}_{1}$ is enhanced at lower frequencies and varies weakly with $\mu$ (except if $|\mu\pm\delta\mu|{\,\approx\,}\Delta/2$, cf. black curve, since then VP is strongly affected by small changes in $\mu$). As a result, even in a situation where $\bar{\sigma}_{1}$ and $\bar{\sigma}_{2}$ might have comparable maximum magnitudes, it is possible to separate the VP-dominated and intrinsic-dominated regimes by tuning the relative position of the excitation frequency and $\mu$, since the latter can be used to push up the spectral region for which $\mathrm{Re}\,\sigma_{2}\neq 0$. Furthermore, the rapid whittling of $\bar{\sigma}_{1}$ when $k_{B}T{\,\gtrsim\,}\delta\mu$, in contrast with the robustness of $\bar{\sigma}_{2}$ up to temperatures, significantly above room temperature Hipolito et al. (2016) results in a strong temperature dependence of the ratio $|\bar{\sigma}_{1}|/|\bar{\sigma}_{2}|$.

IV Concluding remarks

We studied the generation and polarization dependence of SH in threefold symmetric 2D materials with a finite VP, and performed specific microscopic calculations of the SH conductivity for a model that applies (accurately) to graphene and (qualitatively) to semiconducting TMD such as MoS2. Our results show that VP and intrinsic (when present) quadratic response generate distinct contributions with contrasting symmetry properties, which can be disentangled by analyzing the dependence of SHG on the orientation of polarization plane [cf. eq. II] or on the state of polarization $\varphi$. To achieve this, the SHG signal $I_{\parallel}$ can be used to determine the orientation of the lattice by either reversing the sign of the VP, or by probing the dependence on $\zeta$ at photon energies above (below) $\mu$, in the regime dominated by intrinsic (VP) where the maxima indicate the armchair (zigzag) directions of the lattice. Knowledge of the lattice orientation (thus obtained or otherwise), permits a direct application of eq. II to extract the two independent nonlinear optical constants $\bar{\sigma}_{1}$ and $\bar{\sigma}_{2}$. Since $\bar{\sigma}_{1}$ is proportional to $\delta\mu$, the SH fingerprint of these systems can be used to directly identify and quantify an underlying imbalance between the populations in the valleys $\mathbf{K}$ and $\mathbf{K}^{\prime}$.

If performed with a small spot size in a scanning mode, such measurements provide a means to directly map valley polarization throughout a system, measure the spatial decay of valley currents, and investigate the possibility or efficacy of their injection across heterostructure junctions and interfaces.

Our data for $\sigma_{iii}$ is presented in units of $\sigma_{0}{\,=\,}2.88\times 10^{-15}$ S m/V. If converted to 3D quadratic susceptibilities using an effective graphene thickness of $d{\,\equiv\,}3.4$ Å, this corresponds to $\chi^{(2)}{\,\equiv\,}\sigma_{0}/(\omega\varepsilon_{0}d){\,\approx\,}6.3\,\mathrm{nm/V}$ at $\omega{\,=\,}0.1$ eV. As a reference, $\chi^{(2)}$ in a good non-linear bulk crystal has typical values of $0.01\,\mathrm{nm/V}$ (ZnO) Chan et al. (2006), $0.5\,\mathrm{nm/V}$ (GaAs, MoS2) Bergfeld and Daum (2003); Nastos and Sipe (2006); Wagoner et al. (1998), $2\,\mathrm{nm/V}$ (monolayer GaSe) Zhou et al. (2015). Hence, SHG due to valley polarization can exceed largely the typical non-linear response of bulk materials such as GaSe. The ability to vary the reference Fermi level in most atomically thin crystals through gating Novoselov et al. (2005), when combined with frequency dependent measurements, further expands the versatility of SH spectroscopy to assess valley-dependent properties, rendering it a valuable characterization tool in the nascent field of valleytronics.

Acknowledgements.

FH thanks A. H. Castro Neto and M. Milletarì for their support and discussions throughout this project. He was supported by the National Research Foundation Singapore under the CRP award NRF-CRP6-2010-05 and by the QUSCOPE center sponsored by the Villum Foundation. VMP was supported by the Singapore Ministry of Education through grant MOE2015-T2-2-059. Numerical computations were carried out at the HPC facilities of the NUS Centre for Advanced 2D Materials.