Polarization-resolved analysis of high-order harmonic generation in monolayer MoS2

We employ a theoretical model based on the density-matrix equation in the velocity gauge to calculate high harmonic generation from monolayer MoS2. This approach incorporates the tight-binding model, enabling the full consideration of both crystal symmetry and multiple band effects. In additional to the usual odd harmonics, even harmonics are also presented in the case of observing two different polarization components, which are parallel and perpendicular to the polarization of linearly driving pulses. We detailedly analyze the crystal orientation dependence for the parallel and perpendicular components of both odd- and even-harmonics. It is found that they exhibits different modulation behavior with rotating the crystal orientation. The simulation results capture all important orientation-dependent features observed in the recent experiment, thus demonstrating that the Berry curvature of MoS2 has been appropriately considered in our proposed model. In order to facilitate analysis of the underlying mechanism, we examine the channel current in terms of the contribution from different density-matrix elements, and identify their role in the orientation modulation of high harmonics. We further use simplified one-dimensional integral model to explain the appearance of perpendicular components of even harmonics. Our analysis shows that the multi-band coupling effect is the origin of the parallel even harmonics, while the broken inversion symmetry of phase difference of momentum matrix elements along two orthogonal directions determines the perpendicular even harmonic generation. Additionally, the relationship between the concept of Berry curvature and our theoretical framework is discussed. These demonstrations show that polarization-resolved high harmonics might provide an all-optical way for imaging material’s Berry curvature.


Introduction
The successful observation of solid high harmonic generation (HHG) beyond band gap energy opens up a new research frontier, which aims at investigating the highly nonlinear response of periodic condensed matters to the ultrashort laser field [1]. A wide range of driving wavelength from the near infrared [2] and midinfrared [3] to the terahertz [4] region has been used to excite HHG in solids under the condition that the laser field amplitude should be controlled to avoid material damage. This new coherent radiation can not only provide a potential pathway to a compact extreme ultraviolet light sources and highly efficient terahertz frequency synthesis [5], but also encode the detailed information on the electronic structure and the ultrafast carrier dynamics in target materials [6].
The microscopic mechanism of solid-state HHG is always interpreted in the reciprocal space [7][8][9] using the concept of the intraband current, which describes the electron Bloch motion within a single band, and the interband current, which considers the transition coupling among different energy bands. It has been demonstrated that the two mechanisms lead to the different dependence of HHG on driving pulse wavelength [10] and ellipticity [11,12]. The simple single-particle picture constitutes the basis of high-harmonic spectroscopy towards the all-optical reconstruction of band gap in ZnO [13] and the measurement of Berry curvature [14].
Besides the specific electronic energy band, the periodic lattice structure has the spatial symmetry of translation, rotation, inversion and reflection, which strongly affect HHG in solids. For the molecular system, it has been reported that the harmonic polarization is sensitive to the molecular geometry symmetry and nuclear dynamics [15]. Similarly, the crystal symmetry also allows to generate harmonics from solids with controlled polarization state [16], as well as the temporal and spectral properties [17]. Due to the lattice structure, the HHG also exhibits the specific dependence on the circular and elliptical polarization properties of the driving light [18][19][20][21][22]. By combining the time properties of driving fields, spatiotemporal symmetries can lead to various selection rules imposed to HHG in solids [23,24]. Even if a linearly polarized driving pulse is used, the complex spatiotemporal waveform of harmonic radiation can be also generated. Thus, the study of the polarization components of high harmonics along two orthogonal directions, parallel and perpendicular to driving field, would be helpful to get deep insight into the microscopic physical process. This kind of polarization-resolved measurement upon varying crystal orientation has been implemented experimentally in bulk semiconductor GaSe [25], and SiO 2 with different levels of crystallinity [26].
Recently, there has been a growing interest to extend HHG in two-dimensional (2D) materials [5,19,[27][28][29] due to their distinctive electronic and optical properties. As a prototypical material, the single-layer MoS 2 possesses the broken spatial inversion symmetry and the nonzero direct band gap. The HHG from MoS 2 exhibits an important characteristics that the even harmonics are predominately polarized perpendicular to the driving pulse polarization [30], which is attributed to the Berry curvature inducing anomalous current [31]. Although the measurement can be qualitatively explained with the single-band model, the quantitative calculation of HHG should consider the coupling transition among different bands. Thus, addressing the issue becomes one of the important goal of our work. To date, various numerical methods have been proposed to calculate HHG in solids [32]. The length-gauge semiconductor Bloch equation (SBE) using the field-free Bloch state [33][34][35][36] or the Houston state [10,37] as basis is one of the relatively simple method, which can easily incorporate real material parameters and consider relaxation processes phenomenologically. The physical quantity in the SBE describing the interband transition is the dipole matrix, whose diagonal element associated with the Berry connection [31] is always ignored in the calculation. Thus, the previous model is not applicable to the HHG process where the Berry curvature plays an important role. In addition, the previous model usually requires the continuous phase imposed to dipole element, which is always a difficult or time-consuming task. If more bands are included or there are two bands that are nearly degenerate, the previous model suffers from the numerical difficulty, mainly due to the fact that the dipole element is inversely proportional to band gap and has the near discontinuity at band edge.
In this work, for appropriately considering the Berry curvature effect on HHG from single-layer MoS 2 , we start from the equation of motion for the density matrix in the velocity gauge, which holds the legitimacy of Bloch's theorem. The model is applicable to multiple band structures without any numerical difficulty. The full Hamiltonian of MoS 2 is constructed from the tight-binding (TB) approximation involving up to the third-nearest-neighbor hoppings. This treatment leads to three-band energy structure that can cover the path interference effect of electron excitation beyond the two-band model. In this framework, the excitation coupling is described by the momentum matrix element, whose value can be determined by several TB parameters. The phase of momentum matrix can be chosen in an arbitrary way. It is found that the proposed model can reproduce almost all important crystal orientation-dependent features reported in the recent experiment [30]. We focus on the interesting physics underlying the polarization-resolved measurement of orientation-dependent HHG in single-layer MoS 2 . The origin of the observed modulation is analyzed in terms of the band structure and the relevant density-matrix element. Moreover, a simplified one-dimensional (1D) integral expression of current is obtained under the two-band model. According to the formula, we show the generation mechanism of the perpendicular even harmonics from a different perspective. We also discuss the relationship of this interpretation with the Berry curvature.

Theoretical method
The basic unit block of MoS 2 consists a trigonal prismatic structure where Mo atoms are in the middle of two-layer of S atoms [38]. The top view of its atomic structure is regarded as the analog of graphene, as shown in figure 1(a), where black-S atoms and blue-Mo atoms are alternatingly arranged in a 2D hexagonal honeycomb lattice structure. The two in-plane primitive lattice vectors are denoted by a 1 = a 2 (−1, where a is the nearest-neighbor in-plane S-S or Mo-Mo distances. Figure 1(b) shows the corresponding in-plane hexagonal Brillouin zone (BZ), together with high-symmetry points Γ = (0, 0), . To obtain the electronic band structure of MoS 2 , we consider a general TB eigenstate where k is the crystal momentum in BZ, and n is the band index. The basis states |B α (k)> are expressed by the linear combinations of atomic orbitals as where there are N lattice sites in the crystal, |φ α (R) > is an atomic orbit centered at real-space position R, and the label α represents the combined index of both atom species and its orbital symmetry.
The interaction of MoS 2 with an ultrashort pulse is described by density-matrix equations within the independent-particle approximation. This approach can take into account multiple bands and Pauli blocking of interband transitions [40]. The relaxation processes due to scattering effect beyond the mean-field approximation can also be included via phenomenologically introducing the dephasing term.
The temporal evolution of density operator ρ is described by Liouville-von Neumann equation The total time-dependent Hamiltonian in the dipole approximation and in the velocity gauge is given by where p is the canonical momentum operator, U(r) is the periodic potential, and A(t) is the time-dependent laser vector potential, related to the electric field F(t) = −dA(t)/dt. The Hamiltonian equation (14) still maintains its spatial periodicity in the presence of the high-intensity electric field, which ensures the validity of Bloch's theorem [41]. We can remove the spatially independent term A 2 (t)/2 in equation (14), since it only contributes to a time-dependent phase and does not affect the expectation value of the current and polarization induced by the laser field. Thus, the Hamiltonian given by equation (14) can be rewritten in the Coulomb gauge as where the field-free full Hamiltonian H 0 is constructed with TB eigenstates as Here we have considered such fact that there is no coupling between different k values in the velocity gauge. By substituting equations (15) and (16) into (13), and including the relaxation terms phenomenologically, the dynamic equation for the temporal evolution of density-matrix elements reads where theρ is the density matrix with its element defined as ρ mn (k, t) = m, k|ρ|n, k . The last two terms describe the transverse and longitudinal relaxation process with a constant dephasing time T 2 and T 1 , respectively. The f 0 m (k) is the initial distribution of the m-band electrons in the BZ. The momentum matrix p(k) consists of the element p mn (k) = m, k|p|n, k , which can be expressed as [42] The latter term in equation (18) represents the intra-atomic contribution, which cannot be obtained with the TB parameters due to the unknown φ α (0)|r|φ β (0) . It has been demonstrated that the intra-atomic term leads to an approximate k-independent shift of momentum matrix [42]. Consequently, p mn (k) can be modified as where we introduce a factor γ to rescale momentum matrix. In our simulation, we choose γ = 0.1 to make the amplitude of p mn (k) at K point be comparable to ab initio result, and therefore p mn (k) can be completely determined by the TB model.
After propagating density-matrix elements in equation (17), the time-dependent current density at a particular k is obtained by where Tr denotes trace, and N vb is the number of valence bands. The total current density is then calculated by the integral of j k (t) over the first BZ The two components of J(t), parallel and perpendicular to the driving laser polarization, are expressed as J // (t) = J(t) · e // (ξ) and J ⊥ (t) = J(t) · e ⊥ (ξ), respectively. Here the two unit vectors are defined as e // (ξ) = (cos ξ, sin ξ) and e ⊥ (ξ) = (−sin ξ, cos ξ), where ξ denotes the crystal orientation described by the polar angle of driving polarization, as shown in figure 2(b). For the polarization components of high harmonic radiation, parallel and perpendicular to the linearly polarized fundamental field F(t), the resulting harmonic spectrum is calculated by where T FF is the fast Fourier transform algorithm. For clarity, in the following text we call I // (ω) and I ⊥ (ω) as 'parallel harmonics' and 'perpendicular harmonics', respectively. Equation (17) is numerically solved for each independent k by the classical 4th order Runge-Kutta method. We sample the Brillouin zone with uniform grid along two non-orthogonal directions of in-plane reciprocal lattice vectors. The equal k-grid spacing is employed for the two directions. The grid interval of crystal momentum and evolution time is fixed to δk = 0.004a · u. and δt = 0.1a · u. We have examined and confirmed the simulation convergence by the comparison of results obtained with decreasing the δk and δt value.

Calculation of crystal-orientation dependence
We have firstly calculated the polarization-resolved harmonic spectrum upon varying the MoS 2 crystal orientation. The driving mid-infrared laser field is modelled by a Gaussian temporal envelope with 160 fs duration and 1.0 TW cm −2 peak intensity at a central photon energy of 0.3 eV. These laser parameters are the same as the ones used in the experiment [30], which is convenient for examining our theoretical model. The initial condition for solving equation (17) is that the valence band are fully populated across the entire BZ, leaving all conduction bands empty. The dephasing times are set as T 2 = 20 fs and T 1 = 500 fs. We have checked that the faster dephasing time T 2 only induces the decreased high harmonic intensity, and does not affect the conclusion. Figure 2(a) shows the total harmonic spectrum including the parallel and perpendicular component at a fixed crystal orientation angle ξ = 0. Both the clear odd and even harmonics are present above the direct bandgap of ∼1.8 eV, with the cut-off order extending to 19, reflecting highly nonlinear response process in HHG. Note that harmonic intensity is plotted with linear scale throughout our work. The previous study [41,43] shows that at high field strengths the cut-off energy is no longer limited by the band-gap energy. Instead, the momentum operator matrix element p cv between valence and conduction bands becomes an important factor that contributes to the cut-off position of harmonic spectrum. Our model gives the cut-off order larger than the experimental result [30]. This is probably due to the fact that the magnitude of k-dependent p cv is not exactly obtained by the simplified TB model. However, the symmetry of p cv in the k-space is correctly considered with this model, so that this will not affect the study of modulation behavior of harmonic spectrum.
We point out that the appearance of even harmonics with the intensity comparable to the odd harmonics are important features in the HHG of single-layer MoS 2 . The more information can be obtained by further analyzing the different polarization components of high harmonics, i.e., the parallel harmonic spectrum I // (ω) and the perpendicular harmonic spectrum I ⊥ (ω). To facilitate the comparison with experimental results [30], we choose from the 6th to 13th harmonic order (HH6-HH13) as an observable. Figures 2(b) and (c) show the calculated parallel and perpendicular high harmonic spectrum as a function of the crystal orientation angle, respectively.
We begin to describe main features in the parallel harmonics shown by figure 2(b), where the harmonic signal is dominated by the odd orders. Although the odd harmonics can be observed for all crystal orientations, their intensity exhibit 60 • modulation with changing the crystal orientation. In addition to the odd harmonics, figure 2(b) also shows the weak even harmonic signals (10th and 12th) in the linear false color scale, which are modulated with 60 • periodicity, and the modulation is in phase. However, different characteristics are presented for the perpendicular harmonics in figure 2(c), where the even harmonics (6th, 8th, 10th, 12th) dominate and exhibit in-phase intensity modulation with a period of 60 • . Although the parallel and perpendicular components of even harmonics have the same 60 • modulation, there is 30 • phase shift between them. Moreover, the modulation contrast of even harmonics in both components is close to unity, which represents the even harmonic signals cannot be produced in any crystal orientation. That is to say, at some particular crystal orientation, the parallel or perpendicular components of even harmonics are prohibited.
All these results about the orientation modulation phenomenon reported in figure 2 are in excellent agreement with the recent experimental observation [30], thus validating our theoretical model. Note that the definition of orientation angle is respect to the Γ-K direction, as shown in figure 1(b), and there is a 30 • shift in comparison with the definition in reference [30]. Another distinct feature in figure 2(c) is the appearance of rather weak odd harmonics (13th), whose intensity is modulated with a period of 60 • but the single peak is split to two sub-peaks. The split position where the emission is suppressed occurs at the orientation angle ξ = 60 • × n. This phenomena is not observed in the experiment [30], probably due to the low photon flux. It is worth pointing that this kind of modulation behavior for the perpendicular odd harmonics is observed experimentally in another crystal, GaSe [25], which has the same crystallographic point group as MoS 2 .
The modulation depth and period can be more clearly observed by calculating the normalized yield of each harmonics as a function of the crystal orientation. We focus on the 12th and 13th harmonics, since their parallel (//) and perpendicular (⊥) components are both visible in figures 2(b) and (c). Figure 3 shows the normalized yield versus the crystal orientation, calculated for the 12th parallel (H12 // , black-solid), 12th perpendicular (H12 ⊥ , blue-solid), 13th parallel (H13 // , red-solid), and 13th perpendicular (H13 ⊥ , green-solid) components. The yield is defined by integrating each harmonic peak and normalized to the corresponding maximum. We can see from figure 3 that only perpendicular odd harmonics (H13 ⊥ ) represents the different modulation behavior compared with other components. Except the parallel odd harmonics (H13 // ), the modulation depth of other harmonic components is nearly unity.

Interpretation of orientation dependence in terms of channel currents
To track back the physical origin on the different crystal-orientation dependence, we analyze the channel current associated with the contribution from different density-matrix elements. Instead of using the full three-band model, we remove the second conduction band (C 2 ) and only consider two-band model (V and C 1 ). This treatment can highlight the essence and help to reveal the multi-band effect if comparing with three-band results.
According to the general equation (17), the two-band density-matrix equation can be explicitly expressed as and where c denotes the first conduction band (C 1 ), v denotes the highest valence band (V), P(k, t) = ρ cv (k, t) represents the off-diagonal element, and f n (k, t) = ρ nn (k, t) represents the diagonal element, with n = c, v. The time-dependent band energyε captures the electron Bloch motion in k space, which can be seen by the relationshipε n (k, t) = E n [k + , t), and the Rabi frequency is Ω(k, t) = A(t) · p cv (k). The s = +1 is for valence band (v) and s = −1 is for conduction band (c). As shown by equations (20) and (21), the total current density can be separated into two channel parts (25) represents the polarization current associated with the off-diagonal element of density matrix, and represents the transport current associated with the diagonal element of density matrix. We have neglected the N vb A(t) term due to without affecting the final harmonic spectrum. This kind of decomposition for the current is different from the so-called interband and intraband current, which is always discussed in the length gauge and associated with the position operation matrix element [44]. The purpose of our decomposition method is to study how different parts of the density matrix element determine the harmonic spectrum. Similar to three-band case, the parallel and perpendicular components of two currents can be obtained by J p,// (t) = J p (t) · e // (ξ), J p,⊥ (t) = J p (t) · e ⊥ (ξ), J t,// (t) = J t (t) · e // (ξ), and J t,⊥ (t) = J t (t) · e ⊥ (ξ). From these currents, we can calculate the corresponding harmonic spectrum. Figure 4 shows the harmonic spectrum as a function of the crystal orientation, obtained from different current components: (a) J t,// (t), (b) J t,⊥ (t), (c) J p,// (t), and (d) J p,⊥ (t). Several important conclusions can be drawn: (i) the parallel components of harmonics, calculated from channel currents J t,// (t) and J p,// (t), have no even orders for all orientation angles, as shown in figures 4(a) and (c). This result is different from the one obtained by three-band model, see figure 2(b). Thus, we identifies the importance of indirect excitation pathways from multi-band coupling in generating the parallel components of even harmonics. This can be further interpreted by observing figure 1(c), which shows avoided cross points between C 1 and C 2 band at Γ-M direction. At these avoided crossing points, the transition from C 1 to C 2 has the large probability, so that the indirect excitation is enhanced. Consequently, the even harmonics would become evident. This interpretation is consistent with figure 2(b). The strongest parallel even harmonics occur at ξ = 30 • + 60 • × n, which is exactly the Γ-M direction; (ii) the comparison of figures 4(a) and (c) shows the parallel components of odd harmonics originated from J t,// (t) and J p,// (t) have the same modulation behavior, except that the contribution from J p,// (t) is dominated; (iii) for the perpendicular components shown in figures 4(b) and (d), odd and even harmonics are derived from different channel currents, exhibiting different modulation behavior. Additionally, the intensity of odd harmonics is nearly two orders of magnitude lower than the intensity of even harmonics. This poses challenge for experimental observation of perpendicular odd harmonics from MoS 2 .
We can see from figures 2 and 4 that the harmonic intensity is mainly modulated with two different characteristics. The 60 • modulation can be attributed to the combination of crystal and laser symmetry together [23]. However, the 60 • modulation including two sub-peaks from the perpendicular odd harmonics cannot be intuitively understood by rotation symmetry. We have demonstrated that the perpendicular transport current J t,⊥ (t) contributes to this kind of modulation. Consequently, the generating mechanism can be related with the band curvature, which is similar to results reported in GaSe [25] due to the same crystallographic point group.

Analytical investigation of perpendicular even harmonic generation
As shown in figure 4, polarization currents J p,// (t) and J p,⊥ (t) dominate the harmonic emission, and the J p,⊥ (t) is responsible for generating the perpendicular even harmonics, which is regarded as the important characteristics of HHG in MoS 2 . Therefore, we only focus on the channel current J p (t) and analyze the physical mechanism of generating even harmonics.
We first rewrite equation (22) in an integral form where D(T 2 , t − t ) = exp[−(t − t )/T 2 ] describes the decoherence relaxation process, and ε g (k, τ ) =ε c (k, τ ) −ε v (k, τ ) is the energy difference between C 1 and V band. Note that equation (27) is completely equivalent to equation (22) without making any approximation. We decompose momentum matrix element into two orthogonal components as By combining equations (25), (27), and (28), we can obtain two orthogonal components of channel current J p (t), parallel and perpendicular to fundamental pulse polarization, in the form of and where δφ(k ξ , ξ) = φ // (k ξ , ξ) − φ ⊥ (k ξ , ξ) represents the phase difference of momentum matrix element along two directions e // (ξ) and e ⊥ (ξ). Here we approximately consider the integration in one-dimension (1D) direction, which corresponds to the polarization of fundamental field. The 1D integral path C in equations (29) and (30) passes through Γ point (k ξ = 0) and is symmetry with respect to Γ point. The simplified 1D integration method has been used in the previous work [29,33]. We calculate the high harmonic spectrum as a function of crystal orientation, using the integral current expression equations (29) and (30), where w(k ξ , t ) is obtained in advance via numerically solving equations (22) and (23). The resulting spectrum from J p,// (t) and J p,⊥ (t) are reported in figures 5(a) and (b), respectively. In comparison with figures 4(c) and (d), one can see that the integral 1D formula can reproduce not only the correct origin of odd [from J p,// (t)] and even [from J p,⊥ (t)] harmonics, but also the modulation with a period of 60 • . In particularly, figure 5(b) shows that the modulation phase of even harmonics predicted by the 1D approximation agrees with the exact simulation. Thus, it is reasonable to use equation (30) for revealing the generation mechanism of perpendicular even harmonics.
The phase term in equation (29) and (30) is key to determine the harmonic emission. Compared with J p,// (t), J p,⊥ (t) includes an additional phase δφ(k ξ , ξ), depending on the crystal orientation. It is therefore expected that the symmetry of δφ(k ξ , ξ) plays an important role in generating even harmonics. The impact of δφ(k ξ , ξ) can be represented by its sinusoidal and cosine value, i.e., sin δφ and cos δφ. Figures 5(d) and (e) respectively shows the calculated sin δφ and cos δφ as a function of crystal orientation ξ and 1D crystal momentum k ξ . We can see that for all orientation angels cos δφ always maintains the inversion symmetry with respect to Γ point (k ξ = 0). However, sin δφ only has inversion symmetry with respect to Γ point at some specific crystal orientation denoted by the vertical black-dot line in figure 5. These orientation angles, given by ξ = 30 • + 60 • × n with n being an integer, corresponds to fundamental polarization parallel to the MoS 2 real-space mirror plane. The orientation-dependent yields for even harmonics shown by figure 5(b) is plotted in figure 5(c). One can easily see that even harmonics are strongly suppressed at these angles where both sin δφ and cos δφ have the inversion symmetry. For ξ = 60 • × n denoted by the vertical red-dashed line in figure 5, sin δφ represents the maximal asymmetry with respect to Γ point. In this case, even harmonics are the most efficiently emitted. The polarization direction of the perpendicular even harmonics corresponds to the mirror reflection symmetry axis of the real-space MoS 2 lattice. We emphasize that the similar phenomenon is also present for the H + 3 /H 2+ 3 molecular system [15] with the same three-fold symmetry as MoS 2 . An exact simulation of this molecular model shows that even or odd, or both harmonics appear for different laser orientation. In particular, only perpendicular even harmonics can occur due to the reflection symmetry of the molecule with respect to the perpendicular emission direction. Thus, this simple and important symmetry holds for molecules and solids.
It is worth pointing out that although equation (30) indicates that the δφ(k ξ , ξ) is responsible for even harmonics, this does not mean J p,⊥ (t) cannot induce odd harmonics. In order to confirm this point, we directly calculate high harmonic spectrum from J p (t) at a non-high symmetry crystal orientation ξ = 10 • , using the exact equation (25) instead of the 1D approximate model. Figure 6 shows the J p (t) induced parallel (blue-solid) and perpendicular (red-dashed) components of high harmonics, on the logarithmic scale. One can conclude that J p,// (t) only leads to odd harmonics, while J p,⊥ (t) can induce very weak odd harmonics with the intensity lower nearly two orders of magnitude than even harmonics. This explains why only perpendicular even harmonics are observed in the experiment.
In reference [30], Liu et al use the concept of Berry curvature to qualitatively understand the perpendicular even harmonics. We turn to simply discuss the relationship between our theoretical model and Berry curvature. With a Kubo-like formula [45], the Berry curvature of conduction band can be written as With a single-band model, the anomalous current induced by Berry curvature is expressed as [31] J where ρ(k, t) denotes the k-space electron distribution after excitation. By substituting equations (28) and (31) into (32), the anomalous current can be written as We find from equation (33) that some important factors, such as p // (k, ξ), p ⊥ (k, ξ), and δφ(k, ξ), are also presented in our theoretical formula equation (30). More importantly, equation (33) shows that the δφ(k, ξ) through its sinusoidal value sin δφ determines the J c (t), which is consistent with our above analysis, see figure 5(d). These agreements support the fact that our proposed theoretical frame can capture the Berry curvature and consequently reproduce the experimental result very well.

Conclusions
In summary, for the first time we have theoretically investigated a polarization property of high harmonics from monolayer MoS 2 . An alternative theoretical framework based on the combination of density-matrix equation in the velocity gauge and tight-binding electronic structure is proposed to calculate high harmonic spectrum as a function of the crystal orientation. The theoretical model can not only take into account the crystal symmetry and the multiple band coupling, but also appropriately include the material's Berry curvature. Thus, the excellent agreement between the simulation and experiment is obtained. We provide the detailed analysis of the polarization-resolved high harmonics with scanning the crystal orientation. The origin of parallel and perpendicular components of high harmonics is interpreted in terms of channel current associated with the different density-matrix element. In addition to reproducing the key experimental result, we clearly identify that the parallel components of even harmonics are generated as a consequence of multi-band effects, involving the interference of direct and indirect excitation pathways. We further emphasize the mechanism of generating perpendicular even harmonics by a 1D integral model. It is demonstrated that the broken inversion symmetry related with the phase difference of two orthogonal components of momentum matrix elements is the key to induce the perpendicular even harmonics. We also reveal that the perpendicular components of odd harmonics can in principle exist. However, the spectral intensity is estimated to be two orders of magnitude lower than the accompanying even harmonics, leading to the difficulty in experimental observation. The domination of even-order harmonics in the perpendicular component can be attributed to the sufficiently large berry curvature of specific materials. Finally, we qualitatively discuss why our model can capture the Berry curvature effect. Thus, our proposed model is able to quantitatively calculate high harmonic generation in a regime where the Berry curvature cannot be ignored. It might become a useful tool for the investigation of the material's Berry curvature through high harmonic spectroscopy.