Effect of Interlayer Coupling and Symmetry on High-Order Harmonic Generation from Monolayer and Bilayer Hexagonal Boron Nitride

Abstract

High-order harmonic generation (HHG) is a fundamental process which can be simplified as the production of high energetic photons from a material subjected to a strong driving laser field. This highly nonlinear optical process contains rich information concerning the electron structure and dynamics of matter, for instance, gases, solids and liquids. Moreover, the HHG from solids has recently attracted the attention of both attosecond science and condensed matter physicists, since the HHG spectra can carry information of electron-hole dynamics in bands and inter- and intra-band current dynamics. In this paper, we study the effect of interlayer coupling and symmetry in two-dimensional (2D) material by analyzing high-order harmonic generation from monolayer and two differently stacked bilayer hexagonal boron nitrides (hBNs). These simulations reveal that high-order harmonic emission patterns strongly depend on crystal inversion symmetry (IS), rotation symmetry and interlayer coupling.

1. Introduction

Thanks to the significant advances in laser technology within the visible, infrared and terahertz regions, there has been a remarkable interest in understanding the nonlinear process in materials under the influence of a strong laser field. This is important not only for strong-field physics in condensed matter [1,2,3,4,5,6] but also for new technologies such as ultrafast optoelectronic and photonic applications.

One of the most representative nonlinear processes in solid state subjected to a strong laser field is high-order harmonic generation (HHG). HHG has been extensively studied and used in gases and molecules for table-top XUV sources as well as a tool to explore ultrafast electron dynamics [7]. Recently, research on HHG has been extended to solid-state physics, specifically on materials such as ZnO, MgO, GaSe, graphene, GaP meta surfaces, transition metal dichalcogenides (TMDCs), Cd${}_3$As${}_2$, hexagonal boron nitride (hBN) and topological materials [2,8,9,10,11,12,13,14,15,16,17,18]. HHG in solid-state materials may provide more compact and efficient XUV sources as well as information on electron structures and ultrafast electron-hole dynamics [19,20].

In recent years, two-dimensional (2D) materials have been studied for potential applications to next-generation nano-photonics and nano-electronic devices [21,22,23,24,25]. It is notable that 2D materials are unique in terms of their potential to combine various atomically thin layers to control material properties or even create new materials, such as so-called van der Waals heterostructures [26]. Since van der Waals interactions between layers is rather weak, 2D materials with different compositions can be readily assembled, providing additional knobs to manipulate material properties, such as conductivity or semi-metal characteristics [27,28,29].

Nowadays, the stacking of atomically thin layers with a high degree of control has become a reality [30,31,32,33]. Layer spacing can even be controlled by pressure [34] or doping [35]. It is therefore timely and demanding to understand how stacking scheme and interlayer coupling affect the optical properties of HHG in the context of strong-field physics. The identification of clear fingerprints related to a particular stacking or layer coupling could lead to HHG spectroscopy for monolayer, bilayer or few-layer systems [12].

Because of its high damage threshold, hBN is a promising candidate for light emission in the far UV region via HHG [36,37]. Theoretical investigations on the various characteristics of HHG from hBN, such as rotation anisotropy, time-frequency behavior and wavelength scaling of harmonic yields, have been pursued recently [16,38,39]. Improvements in growth technology have enabled hBN to be stacked in various ways [40,41,42]. Pioneering research has indicated that among five high symmetry stacking orders, the AB and AA${}^{\prime }$ stacking orders are the most stable configurations in either bilayer systems or in bulk [43,44,45].

In the present paper, using well-established semiconductor Bloch equations (SBEs) with a tight-binding approach to model the effect of interlayer coupling on the HHG process, we explore the HHG features of the monolayer and bilayer hBNs driven by linearly, circularly or elliptically polarized laser fields.

This paper is organized as follows: the methodology is discussed in Section 2. In Section 3, we present the calculated HHG spectra in the monolayer and AB and AA${}^{\prime }$ bilayers by varying the laser polarization angle, the ellipticity of driving laser field and interlayer coupling. The HHG spectra are analyzed in terms of symmetry and interlayer coupling. Finally, we summarize our findings and analysis in Section 4.

2. Methods

Figure 1a shows a top view of the monolayer hBN, which has unit cell parameter a = b = 2.5 Å [46] in real space. The monolayer hBN belongs to the $D_{3h}$ point group. It has both 3-fold symmetry and broken inversion symmetry (IS). Figure 1b,c are the real space representations of the AB and AA${}^{\prime }$ stackings, respectively. Both stacking types share the same interlayer distance, 3.3 Å [46]. The corresponding Brillouin zone (BZ) is visualized in Figure 1d. We choose the $\Gamma -K$ direction as the x direction and the $\Gamma -M$ direction as the y direction. Figure 1e,f are the energy dispersion relations of the AB and AA${}^{\prime }$ stacking in comparison with the monolayer case. We choose the same stacking notation (AB or AA${}^{\prime }$) as done in the work by Ribeiro and Peres [43]. These two different stacking schemes are chosen to emphasize the symmetry change between the monolayer and bilayers—the symmetry of the AB bilayer is almost the same as the symmetry of the monolayer (belonging to point group $C_{3v}$), while the AA${}^{\prime }$ bilayer exhibits a huge difference (point group $D_{3d}$). In detail, the monolayer has the x-y plane as a mirror plane, but the AB bilayer does not have such a mirror plane, i.e., $H\left(z\right)\ne H(-z)$. This is the only difference between the monolayer and AB bilayer, and it is related to symmetry in z-direction. In this paper, we study 2D material; hence, ${\mathit{k}}_z$ is not a good quantum number. Moreover, we propagate a laser into the material in the surface normal direction, and do not consider applying an out-of-plane electric field. Hence, we can assume that the monolayer and AB bilayer have the same symmetry configuration in this paper.

Unlike the other two cases, the AA${}^{\prime }$ bilayer recovers IS, and this makes a large difference between the AA${}^{\prime }$ bilayer and the other two cases. It is also worth noting that the AA${}^{\prime }$ bilayer has degeneracy at K point, while the AB bilayer has band splitting, as shown in the inset of Figure 1e,f.

2.1. Hamiltonian of the hBN

We use a simple Hamiltonian that considers only one $p_z$ orbital (or hybridized orbital) from the B atom and N atom. Only considering nearest neighbor (NN) atoms, a monolayer Hamiltonian is written as

$$\begin{array}{c}\hfill H_0=\left[\begin{array}{cc}{M}_0& \gamma \\ {\gamma }^{*}& -M_0\end{array}\right].\end{array}$$

(1)

On-site energy difference is written as $M_0=({\epsilon }_B-{\epsilon }_N)/2$, where ${\epsilon }_B$ is the on-site energy in the B atom and ${\epsilon }_N$ is the on-site energy in the N atom. The interaction between B and N atoms is

$$\begin{array}{c}\hfill \gamma =t_1\sum _{i=1}^3{e}^{i\mathit{k}·{\mathit{a}}_i},\end{array}$$

(2)

where $t_1$ is the NN hopping parameter, and ${\mathit{a}}_i$ is the NN vector from the B atom to the N atom [47]. The first row or column corresponds to the $p_z$ orbital at the B atom ($p_z^B$), and the second refers to the $p_z$ orbital at the N atom ($p_z^N$). A bilayer Hamiltonian can be constructed by ‘stacking’ the monolayer Hamiltonian, $H_0$. Before stacking, it is mandatory to carefully consider the lattice orientation. For example, in the AA${}^{\prime }$ bilayer in Figure 1c, the top layer is rotated 180${}^{\circ }$ with respect to the bottom layer. This rotation makes ${\mathit{a}}_i\stackrel{}{\to }-{\mathit{a}}_i$, i.e., $H_0\stackrel{}{\to }{H}_0^{*}$ or $M_0\stackrel{}{\to }-M_0$ by exchanging the B with the N atom.

We can directly stack the monolayer Hamiltonian, $H_0$, to obtain the AB bilayer Hamiltonian because there is no rotation of the layers in AB stacking, as shown in Figure 1b. Then, the AB bilayer Hamiltonian basis can be referred to as $\{p_z^{B,\left(\mathrm{t}\right)},p_z^{N,\left(\mathrm{t}\right)},p_z^{B,\left(\mathrm{b}\right)},p_z^{N,\left(\mathrm{b}\right)}\}$, where (t) is a superscript for the top layer and (b) refers to the bottom layer. For interlayer direction, we only include interlayer NN interaction between two vertically overlapped atoms - N atoms of the top layer and B atoms of the bottom layer (Figure 1b). Then, the Hamiltonian of the AB bilayer is

$$\begin{array}{c}\hfill H_{\mathrm{AB}}=\left[\begin{array}{cccc}{M}_0& \gamma & 0& 0\\ {\gamma }^{*}& -M_0& t_{11}& 0\\ 0& t_{11}& M_0& \gamma \\ 0& 0& {\gamma }^{*}& -M_0\end{array}\right],\end{array}$$

(3)

where $t_{11}$ is the interlayer hopping parameter. It is worth noting that the interlayer NN vector only has a z-component; hence, there is no phase factor in interlayer terms.

The same procedure can be applied to the AA${}^{\prime }$ bilayer. For the case of the AA${}^{\prime }$ bilayer, the top layer is inverted with respect to the bottom layer. We choose $\{p_z^{B,\left(\mathrm{b}\right)},p_z^{N,\left(\mathrm{b}\right)},p_z^{N,\left(\mathrm{t}\right)},p_z^{B,\left(\mathrm{t}\right)}\}$ as the set of atomic orbitals corresponding to the AA${}^{\prime }$ bilayer Hamiltonian basis. The interlayer NN interaction is in-between different atomic species that are directly stacked on top of each other (see Figure 1c). The Hamiltonian of the AA${}^{\prime }$ bilayer is written as

$$\begin{array}{c}\hfill H_{{\mathrm{AA}}^{\prime }}=\left[\begin{array}{cccc}{M}_0& \gamma & t_{11}& 0\\ {\gamma }^{*}& -M_0& 0& t_{11}\\ t_{11}& 0& -M_0& \gamma \\ 0& t_{11}& {\gamma }^{*}& M_0\end{array}\right].\end{array}$$

(4)

The tight binding parameters from the study of Ribeiro and Peres [43] are given in

Table 1. Tight-binding parameters for the monolayer, AB bilayer and AA${}^{\prime }$ bilayer hBN from [43]. Parameters are fitted to the valence band calculated by DFT.

Stacking Type	${\mathit{M}}_0$ (eV)	${\mathit{t}}_1$ (eV)	${\mathit{t}}_{11}$ (eV)
Monolayer	1.96	2.33	-
AB	2.095	2.37	0.60
AA${}^{\prime }$	2.04	2.36	0.32

. The $M_0$ and $t_1$ values change slightly depending on the stacking type, even if the monolayer structure remains the same. This happens because the tight-binding Hamiltonian is constructed by projecting the result of the DFT calculation. Some effects, such as interlayer next-nearest-neighbor (NNN) coupling or interlayer third-nearest-neighbor (NNNN) coupling, are neglected in the model of this paper. These neglected effects are projected into the monolayer Hamiltonian and interlayer NN hopping parameter. This indicates that the interlayer coupling partially contributes to the monolayer Hamiltonian.

2.2. Semiconductor Bloch Equations on a Tight-Binding Basis

From the excellent work of Virk and Sipe [48], SBEs in length gauge (LG) are frequently used as a simulation framework to study the HHG process in solids. However, there is a problem related to the phase of Hamiltonian eigenstate [48,49]. An eigenstate of Hamiltonian has the freedom to choose its phase, $|{\psi }_{n,\mathit{k}}\rangle \to e^{i\varphi }|{\psi }_{n,\mathit{k}}\rangle$, without changing the property of the material. Even if this gauge transformation does not change the material characteristics, it affects the simulation result as a numerical error. Constructing a continuous, or good enough, wavefunction gauge to minimize numerical noise has been another big issue. To overcome this problem, we solve the SBEs on a tight-binding basis [50,51] instead of eigenstate basis as in other papers [47,48,49,52,53]. To distinguish the two different bases, we use superscript (W) to indicate the tight-binding basis gauge, and superscript (H) to indicate the eigenstate basis gauge as in [54].

Let us set localized Wannier functions $w_n(\mathit{r}-\mathit{R})=\langle \mathit{r}|\mathit{R}n\rangle$, where $\mathit{R}$ is the Bravais lattice vector. We assume that Wannier functions are orthogonal, $\langle {\mathit{R}}_m|{\mathit{R}}_n^{\prime }\rangle ={\delta }_{m,n}{\delta }_{\mathit{R},{\mathit{R}}^{\prime }}$. Tight-binding basis is represented from the Wannier functions as the Bloch-like basis

$$\begin{array}{c}\hfill |{\psi }^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\rangle =\frac{1}{\sqrt{N}}\sum _{\mathit{R}}{e}^{i\mathit{k}·(\mathit{R}+{\mathbf{\Delta }}_n)}|\mathit{R}n\rangle ,\end{array}$$

(5)

where N is the number of unit cells included as a normalization factor, and ${\mathbf{\Delta }}_n$ is the Wannier center (${\mathbf{\Delta }}_n=\langle \mathbf{0}n|\widehat{\mathit{r}}|\mathbf{0}n\rangle$). For the model used in this paper, Wannier centers indicate the locations of the atoms in an unit cell, as they are localized around an atom. The components of the Hamiltonian matrix on a tight-binding basis are

$$\begin{array}{cc}\hfill H_{mn}^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)& =\langle {\psi }_m^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\left|\widehat{H}\right|{\psi }_n^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\rangle \hfill \\ \hfill & =\sum _{\mathit{R}}{e}^{i\mathit{k}·(\mathit{R}+{\mathbf{\Delta }}_m-{\mathbf{\Delta }}_n)}\langle {\mathbf{0}}_m\left|{\widehat{H}}_0\right|{\mathit{R}}_n\rangle .\hfill \end{array}$$

(6)

By comparing Equations (2) and (6), it is clear that the Hamiltonians used in this paper are constructed from the basis defined in Equation (5). The NN vector ${\mathit{a}}_i$ in Equation (2) is the NN displacement between two atoms, which can be written as $\mathit{R}+{\mathbf{\Delta }}_N-{\mathbf{\Delta }}_B$.

It is worth noting that our basis choice is a bit different from previous research ($|{\psi }^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\rangle ={\sum }_{\mathit{R}}{e}^{i\mathit{k}·\mathit{R}}|\mathit{R}n\rangle$) [43,50,54]. We choose Equation (5) as a basis to handle dipole transition matrix elements (DTMEs) on a tight-binding basis. Well-localized Wannier functions can be approximated as $\langle \mathit{R}m|\widehat{\mathit{r}}|\mathbf{0}n\rangle ={\mathbf{\Delta }}_n{\delta }_{mn}$. With this so-called ‘diagonal’ approximation and our basis choice in Equation (5), we can eliminate DTMEs and the Berry connection on a tight-binding basis, $\langle u_m^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\left|{\nabla }_{\mathit{k}}\right|u_n^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\rangle =0$, where $|u_m^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\rangle =e^{-i\mathit{k}·\mathit{r}}|{\psi }_m^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\rangle$. This basis transformation only affects DTMEs or phase factor $e^{i\mathit{k}·(\mathit{R}+{\mathbf{\Delta }}_m-{\mathbf{\Delta }}_n)}$ in Equation (6), and do not change the terms like hopping parameters or on-site energies, $\langle {\mathbf{0}}_m\left|{\widehat{H}}_0\right|{\mathit{R}}_n\rangle$.

We calculate laser-matter interaction through the SBEs in LG with a tight-binding basis [50,51],

$$\begin{array}{c}\hfill i\frac{\partial }{\partial t}{\rho }^{\left(\mathrm{W}\right)}(\mathit{K},t)=[H_0^{\left(\mathrm{W}\right)}(\mathit{K}+\mathit{A}\left(t\right)),{\rho }^{\left(\mathrm{W}\right)}(\mathit{K},t)],\end{array}$$

(7)

where $\mathit{A}\left(t\right)$ is the vector potential of the laser, and ${\rho }^{\left(\mathrm{W}\right)}$ is the density matrix on a tight-binding basis whose components are ${\rho }_{mn}^{\left(\mathrm{W}\right)}=\langle {\psi }_m^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\left|\widehat{\rho }\right|{\psi }_n^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)\rangle$. $H_0^{\left(\mathrm{W}\right)}$ is a tight-binding Hamiltonian as we discussed in Equation (6). As in previous papers [47,53,55], we represent the dissipation term through a phenomenological constant, dephasing time $T_2$. Since we cannot apply the well-known decoherence term ${\Gamma }_{mn}=(1-{\delta }_{mn})/T_2$ directly to the tight-binding basis, we handle the dephasing term in the eigenstate basis gauge. In the eigenstate basis, SBEs, including the dissipation term, are written as

$$\begin{array}{cc}\hfill i\frac{\partial }{\partial t}{\rho }_{mn}^{\left(\mathrm{H}\right)}(\mathit{K},t)& ={[H_0^{\left(\mathrm{H}\right)}(\mathit{K}+\mathit{A}\left(t\right)),{\rho }^{\left(\mathrm{H}\right)}(\mathit{K},t)]}_{mn}\hfill \\ \hfill & +\mathit{E}\left(t\right)·{[{\mathit{D}}^{\left(\mathrm{H}\right)}(\mathit{K}+\mathit{A}\left(t\right)),{\rho }^{\left(\mathrm{H}\right)}(\mathit{K},t)]}_{mn}-i\frac{1-{\delta }_{mn}}{T_2}{\rho }_{mn}^{\left(\mathrm{H}\right)}(\mathit{K},t).\hfill \end{array}$$

(8)

There are two additional terms in comparison with Equation (7). The first one is a term which contains laser electric field $\mathit{E}\left(t\right)$ and the dipole matrix in an eigenstate basis ${\mathit{D}}^{\left(\mathrm{H}\right)}$, whose diagonal components are the Berry connection and whose off-diagonal components are DTMEs. This term is omitted in Equation (7), since ${\mathit{D}}^{\left(\mathrm{W}\right)}=0$ in our approximation. The final term in Equation (2) is the dephasing term, which comes from electron-electron, electron-phonon scattering [53] and propagation effects [56,57]. Additionally, it is worth noting that $H_0^{\left(\mathrm{H}\right)}$ is the Hamiltonian matrix in the eigenstate basis, hence the diagonal matrix of band energy, ${\left(H_0^{\left(\mathrm{H}\right)}\right)}_{mn}={\epsilon }_m{\delta }_{m,n}$. We handle the dephasing time in the eigenstate basis and calculate other terms on a tight-binding basis. The gauge transformation is performed in every time step as in [50,51].

The microscopic electric current is defined with a density matrix that reads:

$$\begin{array}{cc}\hfill \mathit{j}(\mathit{K},t)& =-\mathrm{Tr}\left(\widehat{v}\widehat{\rho }\right)\hfill \\ \hfill & =-\sum _{m,n}{\rho }_{m,n}^{\left(\mathrm{W}\right)}(\mathit{K},t){\mathit{P}}_{nm}^{\left(\mathrm{W}\right)}(\mathit{K}+\mathit{A}\left(t\right)),\hfill \end{array}$$

(9)

where ${\mathit{P}}_{nm}^{\left(\mathrm{W}\right)}\left(\mathit{k}\right)=\langle {\psi }_n^{\left(\mathrm{W}\right)}\left|{\partial }_{\mathit{k}}{\widehat{H}}_0^{\left(\mathrm{W}\right)}\right|{\psi }_m^{\left(\mathrm{W}\right)}\rangle$ [49,50,58,59]. We integrate this over the first BZ, and calculate the macroscopic currents

$$\begin{array}{c}\hfill \mathit{J}\left(t\right)={\int }_{\mathrm{BZ}}\frac{d\mathit{k}}{{\left(2\pi \right)}^2}\mathit{j}(\mathit{K},t).\end{array}$$

(10)

Finally, the intensity of the HHG spectrum is calculated by Fourier transforming the current $\mathit{J}\left(t\right)$:

$$\begin{array}{c}\hfill I\left(w\right)={\left|\mathcal{F}\left[\frac{d}{dt}\mathit{J}\left(t\right)\right]\right|}^2.\end{array}$$

(11)

In our calculation, we use a laser with the vector potential

$$\begin{array}{c}\hfill \mathit{A}\left(t\right)=\frac{E_0}{w_0}f\left(t\right)\left(-\frac{1}{\sqrt{1+{ϵ}^2}}sin\left(w_0(t-t_0)-{\varphi }_0\right)\widehat{{\mathit{e}}_x}+\frac{ϵ}{\sqrt{1+{ϵ}^2}}cos\left(w_0(t-t_0)-{\varphi }_0\right)\widehat{{\mathit{e}}_y}\right),\end{array}$$

(12)

where $E_0$ is the electric field peak strength, $w_0$ is the laser angular frequency, ${\varphi }_0$ is the carrier-envelope phase (CEP), $ϵ$ is the ellipticity of the laser, and $f\left(t\right)$ is the laser pulse envelope. We use Gaussian envelope $f\left(t\right)=exp[-4log2\frac{{(t-t_0)}^2}{t_{\sigma }^2}]$ where $t_0$ is the center of the pulse and $t_{\sigma }$ is the full width at half maximum (FWHM) of the laser. The electric field is given by $\mathit{E}\left(t\right)=-{\partial }_t\mathit{A}\left(t\right)$. Equation (12) is constructed in a way that the electric field has its maximum when $t=t_0$ and ${\varphi }_0=0$ at least in a slow varying envelope. We use laser parameters of $E_0$ = 0.3 V/Å, $w_0$ = 0.02398 a.u. (corresponding to wavelength $\lambda$ = 1.9 $\mathsf{\mu }$m), and FWHM of 10 cycles with a Gaussian envelope.

3. Results and Discussion

We start by analyzing the effect of layer stacking and symmetry change on the HHG spectra in monolayer and bilayer hBNs. For a 10-cycle driving field of 1.2 $\times {10}^{12}$ W/${\mathrm{cm}}^2$ intensity at 1.9 $\mathsf{\mu }$m wavelength, the calculated HHG spectra are displayed for linearly polarized and circularly polarized driving laser fields in Figure 2. The intensities of HHG spectra are normalized to the number of layers in the system. In Figure 2a,c,e, the HHG spectra show symmetry effects as a function of layer stacking, i.e., monolayer, AB and AA${}^{\prime }$ bilayer. We observe remarkable inversion symmetry effects on the selection rules for all optical harmonic orders (HOs) between Figure 2a,c,e. The inversion symmetry of the AA${}^{\prime }$ bilayer allows us to observe only odd HOs independent of the driving laser field direction (Figure 2e). For the monolayer and AB bilayer, the $\Gamma -M$ direction is a mirror plane, and thus, even HOs vanish. When the polarization of the laser is along the $\Gamma -K$ direction, even HOs appear, and their directions are perpendicular to the driving laser field (see Figure 2a. In Figure 2b,d,f, we show the HHG spectra from a right circularly polarized (RCP) driving field. Here, the harmonic emissions are analyzed as a sum of the left circular polarization (LCP) and RCP components. Depending on the harmonic order, one of the RCP or LCP components dominates the harmonic emissions. The HHG signal even disappears in a certain HO. The monolayer and AB bilayer hBNs follow the 3n ± 1 selection rule. When HO is 3n−1, HHG have opposite circular polarization to the driving laser field (counter-rotating order), and when HO is 3n + 1, the HHG have the same circular polarization to the input beam (co-rotating order). There is no signal with the HO 3n. For the case of the AA${}^{\prime }$ bilayer hBN, even order harmonics are canceled due to the effect of non-vanishing inversion symmetry. As a result, it follows the 6n ± 1 selection rule, which is similar to the graphene case [60].

The basic selection rules come from the symmetry of the system, but there are some minor details that cannot be explained by the symmetry. If the interlayer coupling is not considered, the HHG signal from the bilayer system can be simply calculated as a sum of two HHG signals from each band. Then, the symmetry can lead constructive interference in the HHG signal of the AB bilayer as shown in Figure 2c,d. For the case of the AA${}^{\prime }$ bilayer, the odd order harmonics show constructive interference, while the even order harmonics are canceled by destructive interference (see Figure 2e,f). However, this simple constructive (or destructive) interference cannot perfectly explain the HHG behavior. There are slight differences between the HHG spectra from the monolayer, AB bilayer and AA’ bilayer. The interlayer coupling needs to be considered to explain these differences.

However, Figure 2 only shows the HHG signal induced by linearly (along the x-axis in Figure 1d) and circularly polarized driving fields. To have a better understanding, we analyze the HHG spectra with respect to the polarization angle, ellipticity of the driving laser, and interlayer coupling parameter in the following subsections.

3.1. Polarization Angle Analysis in HHG Spectrum

Figure 3 shows the dependence of the HHG spectrum on the laser polarization angle, which illustrates: (1) the effect of the symmetry and (2) the effect of the interlayer coupling by varying laser polarization angles. We analyze the parallel and perpendicular components of HHG with regards to the laser polarization angles of 0${}^{\circ }$ to 120${}^{\circ }$. Here, we set the laser polarization angle to be 0${}^{\circ }$ when it is along the $\Gamma -K$ direction in BZ (see Figure 1d). The results from three geometrically different types of hBN share many features. Every harmonic order has a periodicity of 60${}^{\circ }$, which comes from the 6-fold symmetry of the band structure, even though the crystal symmetry is 3-fold. The perpendicular components of the odd HOs spectra are much smaller than those of other HOs. The perpendicular components of the even HOs have maxima in the $\Gamma -K$ direction, while the parallel components have maxima in the $\Gamma -M$ direction for both even and odd HOs. Additionally, some parallel components of odd harmonics span over the whole angle, and there is no clear minima. All of these observations agree with previous works, such as the monolayer and AA${}^{\prime }$ bilayer hBN calculations [61] or monolayer MoS${}_2$ (the same symmetry as monolayer hBN) experiment [11].

Although HHG spectra from geometrically different hBNs have some similar tendencies, there are differences due to the change of symmetry and interlayer coupling. In the perpendicular components, there is no noticeable difference between the monolayer, AB bilayer and AA${}^{\prime }$ bilayer hBNs. The only difference is that the even HOs in the AA${}^{\prime }$ bilayer are missing due to IS. The parallel components show similar structures, but they are not strictly similar, as in the perpendicular case. For example, the parallel component of HO9 in the monolayer has its maximum intensity at 30${}^{\circ }$, but it has the minimum intensity in the AB and AA${}^{\prime }$ bilayers. It is interesting that the AB bilayer is closer to the AA${}^{\prime }$ bilayer rather than the monolayer in terms of patterns in parallel components of the HHG, even though the symmetries of the AB and AA${}^{\prime }$ bilayers are different. This result suggests that the perpendicular components heavily depend on crystal symmetry, while the parallel components are more affected by energy dispersion and other terms such as interlayer coupling.

3.2. Ellipticity Dependence of HHG Spectrum

The ellipticity dependence of HHG emission is shown in Figure 4 for monolayer and bilayer hBNs. We use the laser defined by Equation (12) and change the ellipticity of the driving laser from −1 (left circular polarization) to 0 (linear polarization) to 1 (right circular polarization). The major axis of the elliptical polarization or linear polarization direction is set to be in the $\Gamma -K$(x) direction. We analyze the intensity of the HHG spectrum in HO5 and HO7, and also show the ellipticity of the output HHG signal as a function of the ellipticity of the driving laser field. Since HHG emission from the AA${}^{\prime }$ bilayer hBN is limited by the 6n ± 1 selection rule for a circularly polarized driving field, we will focus on HO5 and HO7 for our study. We already observed that HO5 and HO7 have a constructive interference for both bilayers in Figure 2b,d,f. Figure 4 shows that they also have a constructive interference in HO5 and HO7 with the elliptically polarized laser. For the case of HO5, the ellipticity dependencies of the HHG spectra in monolayer and bilayer cases are almost similar both in intensity (see Figure 4a,c,e) and ellipticity analysis (see Figure 4b,d,f). They have a similar and typical pattern—maximized at the linear polarization and minimized at the circular polarization. It is impressive that even the tendencies of ellipticity also match each other very well (see Figure 4b,d,f). For the case of HO7, the total intensity is similar (see Figure 4g,i,k), but there are some differences in the ellipticity of the output beam (see Figure 4h,j,l). There is a local maximum of output ellipticity in the monolayer case. The AB bilayer shows the enhancement of the local peak, while the AA${}^{\prime }$ bilayer cancels the peak and flattens the response. It is interesting that although the strength of the local minimum is different, the location of the peak is the same in each case. HO11 and HO13 also behave similarly. It is clear that the output beam ellipticity response is very sensitive to changes in material structure. However, there are some preserved properties that might come from the lattice structure.

3.3. Effect of Interlayer Coupling on HHG Spectrum

Figure 5 presents the effect of interlayer coupling ($t_{11}$) on the HHG spectrum. We plot the HHG spectra as a function of three $t_{11}$ values ($t_{11}$ = 0 eV, 0.54 eV, 1.08 eV) for the AB (Figure 5a) and AA${}^{\prime }$(Figure 5b) bilayers. We observe that high HOs (HOs after the first plateau) get stronger as $t_{11}$ increases. Although there was an observation of amplification of HHG intensities in previous experiments for graphene/hBN hetero structure [34], their explanation cannot be applied here. In Ref. [34], the enhancement is only on the perpendicular components of the HHG signal in all HOs. In this study, only the cutoff region shows clear enhancement, while the plateau region remains complex. We divide the HHG spectrum into low-order and high-order regions. When the strong interlayer coupling is applied, the additional transition path become stronger and the band splitting becomes larger. Additionally, the energy gaps are also larger as $t_{11}$ increases. A larger maximum band gap and multiple paths lead to the enhancement of HHG signal in the high HOs.

To see the change in low HOs in detail, we scan through $t_{11}$ values in Figure 5c,d. The effect of $t_{11}$ is straight-forward in the AB bilayer. Although the variation is small, the HHG intensities have a monotonically decreasing tendency as the $t_{11}$ value increases. In the case of the AA${}^{\prime }$ bilayer, the change is more complex. When the $t_{11}$ value increases, the HHG intensities decrease in the beginning, but at some point, they start to increase. In the end, a stronger HHG spectrum is produced in the AA${}^{\prime }$ bilayer if $t_{11}$ is sufficiently increased. This can be understood through the analysis of dipole transition matrix elements (DTMEs). The absolute values of the DTMEs between the highest valence band (VB2) and two conduction bands (CBs) are plotted for the AB (Figure 6) and AA${}^{\prime }$ (Figure 7) bilayers.

In the AB bilayer, we observe that the DTMEs between CB2 and VB2 (${\mathit{d}}_{v_2,c_2}$) are nearly zero without interlayer coupling ($t_{11}\approx$ 0), as shown in Figure 6d). The bilayer can be treated as two independent layers without interlayer coupling. Therefore, the existing DTMEs in Figure 6a indicate that the VB2 band and CB1 belong to one layer, while VB1 and CB2 belong to the other layer. In addition, the DTMEs between VB2 and CB1 seem to be transferred to the DTMEs between VB2 and CB2 as the interlayer coupling gets stronger. This explains the monotonic change shown in Figure 5.

In the case of the AA${}^{\prime }$ bilayer, different behavior is observed. As shown in Figure 7a,d corresponding to an extremely small $t_{11}$ value ($t_{11}$≈ 3.7 × 10${}^{-9}$ eV), the shape and strength of the DTMEs are already similar to those at equilibrium ($t_{11}$ = 0.32 eV), as seen in Figure 7b,e. This implies that bands from different layers are already mixed even if the $t_{11}$ is small. If the bands from two layers are already mixed, the increase in $t_{11}$ does not simply mean the increase in interaction strength between the bands from different layers, but the strength and shape of the interaction has a complex modification. This can also be predicted by the energy dispersion degeneracy at $\mathit{K}$ point (Figure 1f). This degeneracy cannot be avoided because it comes from the structure of the Hamiltonian, and directly indicates the mixing of the bands from two layers.

4. Conclusions

In summary, we study the effect of interlayer coupling and crystal symmetry with respect to laser polarization angle, driving laser ellipticity, and interlayer coupling $t_{11}$.

• The HHG spectrum from the AB or AA${}^{\prime }$ bilayer can be approximately understood as the constructive/destructive interference of two independent monolayer HHG spectra;
• The effect of interlayer coupling is well manifested in parallel odd components;
• Which aspects are affected by the interlayer coupling are specified;
• More importantly, the effect of interlayer coupling on HHG intensity is analyzed and explained via dipole analysis.

We believe this work will help in the understanding of the interlayer mechanics of bilayer materials, and offer a guide to possible experiments. The role of the interlayer coupling and changes in symmetry could be extended to investigate or design materials from a two-dimensional monolayer to a multiple layer system such as a van der Waals heterostructure.