Enhanced nonlinear optical response of graphene-based nanoflake van der Waals heterostructures.

The nonlinear optical properties of van der Waals bilayer heterostructures composed of graphene/h-BN and graphene/phosphorene nanoflakes are investigated using time-dependent density functional theory. Our calculated results show a significant enhancement of the first-hyperpolarizability value, β in heterostructures relative to the pristine nanoflakes at λ = 1064 nm. The calculated enhancement in optical nonlinearity mainly results from in-plane anisotropy induced by the interlayer electronic coupling between the adjacent nanoflake layers; a higher degree of anisotropy is induced by puckered phosphorene compared to atomically flat h-BN yielding χ (2) value corresponding to the second harmonic generation of ∼50 pm V-1 in the zigzag graphene/phosphorene bilayer heterostructure. The calculated results clearly show that graphene-based nanoflake heterostructures giving large NLO coefficients together with high electron mobility of these materials offer new opportunities as candidate materials of choice for next-generation photonics and integrated quantum technologies. This journal is © The Royal Society of Chemistry.

Introduction

In recent years, two-dimensional (2D) materials including graphene have attracted extensive research interest due to their novel properties and important applications in emerging electronics, photonics, and optoelectronics technologies. Specifically, 2D materials have been found to have large optical nonlinearities and ultra-fast optical responses and are, therefore, proposed as candidate materials for ultrafast laser systems, optical imaging, and nanophotonics.[1-5] Graphene, the most prominent member of the family, exhibits ultrafast NLO response resulting from a sub pico to femtosecond carrier dynamics.[6] The nonlinear effects observed in graphene now include saturable absorption,[7] third-order generation,[8] high-harmonic generation,[9] and nonlinear Kerr effect.[10]

In its purest, defect-free form, graphene lacks a second-order NLO response due to the inversion symmetry of its hexagonal lattice. However, a weak second harmonic generation (SHG) was observed from graphene supported on the oxidized silicon (001) due to substrate-induced symmetry breaking.[11] Furthermore, calculations have shown that the application of an in-plane electric field can induce tunable second-order nonlinear susceptibility in the AB-stacked bilayer graphene.[12] In trilayers, the stacking-induced SHG has been reported in ABA-stacked configuration, while no such response is observed in ABC-stacking, which retains inversion symmetry.[13] In-plane anisotropy can also be induced in a hexagonal lattice of graphene by forming bilayer heterostructures of graphene with the hexagonal (h)-BN monolayer. In such a heterostructure, the interlayer interaction leads to the inequivalent carbon sublattices that, in turn, results in a non-centrosymmetric configuration.[14-16]

The NLO properties of the h-BN monolayer, which is an isoelectronic system with graphene, but an insulator with a large bandgap of ≈6 eV in contrasts to the gapless graphene, has been the subject of several recent theoretical[17-19] and experimental[20-25] investigations. A pristine h-BN sheet exhibits NLO properties.[23] Theoretical calculations have emphasized the role of electron–hole-interaction,[17] as well as inversion symmetry, breaking in the lattice.[19] In addition to the symmetry breaking due to orientation-dependent stacking, experimental measurements have also revealed the effects of stacking order (odd–even number of layers),[21] presence of mechanical stress-induced wrinkles,[22] and microscopic defects,[25] on enhanced SHG response of h-BN nanosheets. In fact, as a natural hyperbolic metamaterial,[26] h-BN is considered an excellent material for a wide range of integrated photonic device applications.[27]

Among the elemental 2D materials, the recently exfoliated phosphorene, which is stable in a puckered configuration due to the preference of the hybridized sp3 orbitals over in-plane sp2 orbitals at the surface,[28,29] has been shown to possess a broadband linear[30] and NLO[31] response from the visible to the mid-infrared region of the spectrum. Recent experiments have also observed high NLO susceptibility corresponding to four-wave mixing[32] and saturation absorption[33] in phosphorene.

In view of their exotic NLO properties, graphene (semimetal), phosphorene (tunable semiconductor), and h-BN (insulator) are potential candidates for tunable, broadband photonic devices, especially for power limiting and higher harmonic generation applications. Therefore, it is of interest to explore the NLO properties of these materials are brought together, i.e., when the two of them are stacked to form a heterostructure. This is what we have addressed in this paper performing a time-dependent density functional theory (TD-DFT) study of graphene/h-BN and graphene/phosphorene nanoflakes with an aim to gain an understanding of their NLO responses at λ = 1064 nm.

Graphene nanoflakes are well characterized by both experimental and theoretical methods.[34-38] Depending upon the edge structures, they can be classified as zigzag or armchair nanoflakes. Note that graphene being zero-gap material exhibits zero optical emission, while the graphene nanoflakes are semiconducting with finite gaps. Several applications of graphene nanoflakes have been proposed, including its use in the spin-photovoltaic device,[39] and tumor-targeting theranostic drug-delivery vehicles.[40] Furthermore, the zigzag and armchair nanoflakes of h-BN and phosphorene were reported to be stable with finite band gaps.[41-47] It is worth mentioning that both h-BN monolayer and phosphorene form stable heterostructures with graphene.[48,49] Moreover, the use of h-BN instead of SiO2 as a substrate for the graphene-based devices was reported to provide better mobility and carrier inhomogeneity.[50]

Fig. 1 shows the hexagonal-shaped nanoflakes of graphene, and h-BN in both the zigzag and armchair edge structures considered for NLO calculations. For graphene, a zigzag C150H60 nanoflake has a width of ∼22 Å, and an armchair C180H72 flake has a width of ∼24 Å. These nanoflakes are significantly larger than the C54H18 nanoflake with a width of ∼12 Å, which was predicted to be stable by the semi-empirical PM3 molecular orbital calculations.[51] For the case of h-BN, the zigzag and armchair nanoflakes considered are B75N75H60 (∼22 Å) and B90N90H72 (∼24 Å), respectively. The zigzag and armchair phosphorene nanoflakes considered are of P96H24 (∼22 Å) and P84H24 (∼24 Å), respectively (Fig. S1, ESI†). Note that the edge atoms are passivated with H atoms, which has been suggested to minimize the reconstruction of its hexagonal-shaped configuration.[52]

Fig. 1

Top views of (a) zigzag graphene nanoflake (b) armchair graphene nanoflake, (c) zigzag BN nanoflake, and (d) armchair BN nanoflake. N and w, respectively, are the total number of atoms and the width of the nanoflake. Color code: C: green, B: purple, N: orange, H: blue.

The static and dynamic linear polarizability (α) and first- and second-hyperpolarizability, β, and γ, respectively, were calculated using TD-DFT theory implementation of the coupled-perturbed Hartree–Fock approach[53,54] in the Gaussian 09 quantum chemistry package. Table 1 lists α calculated as β calculated as corresponding to the static β (0; 0, 0) effect, electro-optic Pockels effect, β (−ω; ω, 0), and the second harmonic generation (SHG), β (−2ω; ω, w). The second hyperpolarizability is calculated for the static case, γ (0; 0, 0, 0), optical Kerr effect (OKE) γ (−ω; ω, 0, 0), and electric field induced second harmonic generation (EFISHG) γ (−2ω; ω, ω, 0). The values of the dynamic α, β, and γ values were reported for a fundamental wavelength (λ) of 1064 nm.

Zigzag graphene nanoflake with a width of ∼22 Å. The polarizabilities calculated using the PBE, PBE+D2, PBE0+GD3, B3LYP, wB97XD, and M062X functionals together with a 6-31G basis set; μ is dipole moment. α, β, and γ are the polarizability, the first-hyperpolarizability, and the second hyperpolarizability, respectively. The frequency-dependent (hyper) polarizabilities are calculated at λ = 1064 nm

Zigzag graphene nanoflake (C150H60)		Exchange and correlational functional
		PBE	PBE+GD2	B3LYP	wB97XD	M062X	PBE0+GD3
μ (D)		0	0	0	0	0	0
〈α〉 (×10−22 esu)	(0)	4.5	4.5	4.3	3.7	4.0	4.0
	(−ω)	7.3	7.3	5.5	4.1	4.8	4.8
β (×10−33 esu)	(0)	1.0	3.6	0.8	0	0	0
	(−ω)	27.2	148.4	2.7	0	0	0
	(−2ω)	295.8	1184.8	3.0	0	0	0
γ xxxx /γyyyy (×10−36 esu)	(0)	2831.7	2822.9	2401.0	1285.0	1630.0	2134.9
	(−ω)	68 419.5	98 045.0	13 520.0	2341.7	3438.5	7120.1
	(−2ω)	−208 343.3	−221 298.3	4050.0	35 856.6	−114 384	13 966.5
γ yyxx (×10−36 esu)	(0)	948.0	940.4	840.0	428.3	545.3	715.0
	(−ω)	127 715.3	97 624.8	9859.0	1161.7	1814.1	4741.8
	(−2ω)	−266 351.7	−255 805.0	2432.0	21 836.7	−78 158	10 910.3

Computational method

In DFT calculations, the range-separated hybrid wB97XD exchange and correlation functional[55] that incorporates a version of Grimme's dispersion term[56] and long-range correction terms[57] was employed. In general, the addition of the D2 dispersion term is necessary for an accurate description of the interlayer interactions in the graphene-based heterostructures.[15] The 6-31G basis set as given in the Gaussian 09 program package[58] was used for the constituent atoms of the nanoflake heterostructures. The convergence criteria for the RMS density matrix and total energy were set at 10−8 and 10−6 eV, respectively. For the geometry optimizations, the maximum force convergence criterion was set at 10−4 eV Å−1.

To establish the suitability of the basis set for NLO calculations, we calculated the hyperpolarizability values of a graphene nanoflake using the 6-31G and 6-31G+(d,p) basis sets within the wB97XD implantation of DFT. Note that the 6-31G basis set is the spilt-valence basis set, whereas the 6-31G+(d,p) basis set is supplemented by a diffuse function together with d and p-type polarization functions on H and other heavy atoms, respectively. As listed in Table S1 (ESI†), the calculated β values for zigzag graphene nanoflake are nearly zero due to the inversion symmetry of its hexagonal lattice. The negligibly small values of β may be a result of the edge passivation of nanoflakes. The calculated γ values the correct order, γ (0) < γ (−ω) < γ (−2ω), and only show a relatively small variation in going from 6-31G to 6-31+G(d,p) basis set. Thus considering the higher computational cost involved with the latter, the split-valence 6-31G basis set appears to be adequate to provide reliable values for NLO coefficients of graphene-based nanoflakes considered in this study. It is also worthwhile to note that the split-valence basis set has been shown to provide reasonable β values for heterocyclic organic molecules including pyrrole.[59]

Next, we investigated the appropriateness of the exchange and correlational functionals for NLO calculations by considering additional functional forms, namely the generalized gradient approximation (PBE),[60] hybrid B3LYP,[61] M062X,[62] and PBE0.[63] In general, the PBE functional form has been shown to provide reliable values of NLO properties of small molecules[64] and crystalline solids.[65] However, we find that this is not the case for the zigzag graphene nanoflake (Table 1). Both the PBE and PBE+GD2 functionals yield a large β for the pristine graphene nanoflake, though it should be nearly zero by symmetry. On the other hand, the hybrid B3LYP wB97XD, PBE0+GD3, and M062X functionals yield a nearly zero value of β.

Further, as listed in Table 1, the results obtained with B3LYP do not follow the expected correct order of γ (0) < γ (−ω) < γ (−2ω). Also, the in-plane principal components, γ (−2ω) and γ (−2ω) predicted to be negative for PBE, PBE+GD2, M062X functionals for zigzag graphene nanoflake and the PBE0+GD3 functional form yields a negative γ (−2ω) value for armchair graphene nanoflake (Table S2 and Fig. S2, ESI†). Such discrepancies associated with the choice of the exchange and correlation functional form have been observed in previous studies[66] as well. Note that the negative γ value is associated with the self-defocusing effect of incident light.[67] Overall, the wB97XD functional form appears to be the appropriate choice to calculate the hyperpolarizability values of graphene-based nanoflakes.

para-Nitroaniline (p-NA) has been used in the literature[68-71] to benchmark the accuracy and appropriateness of basis sets and different methods. Note that p-nitroaniline consists of a benzene ring in which a nitro group is at the para position to an amino group [Fig. S3, ESI†]. Our calculated static and frequency-dependent β values for p-NA are slightly larger than those obtained using the TDHF method[68] but smaller than the corresponding values obtained using either the DFT (LDA)[69] or the DFT (B3LYP)[70] method. Also, it has been shown that the wB97XD functional form with the 6-311+G(d,p) basis set predicts the value of the β‖ (−2w; w, w) (i.e. component parallel to plane) to be 1032 a.u., which is in excellent agreement with the corresponding experimental value of 1072 a.u.[81] Moreover, we also performed additional calculations using the M062X, and PBE0+GD3 functional forms. The results are listed in Table S3 (ESI†) showing that the β(M062X), β(wB97XD), and β(PBE0+GD3) values are nearly similar, with β(PBE0+GD3) being slightly closer to the experimental value.

Overall, the DFT values for β (λ = 1064 nm) are predicted to be in the range of 11–17 (×10−30 esu), showing good agreement with the corresponding experimental value of 15.4 (×10−30 esu) measured in the gas phase (Table S3, ESI†).[71]

Results and discussion

Pristine nanoflakes

We begin with an investigation of the stability of the pristine nanoflakes which is defined in terms of the binding energy (EB) as where Eflake is the total energy of the flake, n is the number of i atoms in the flake, and E is the energy of the isolated i atom (e.g. C, H, B, N, or P). For the (∼22 Å) zigzag and (∼24 Å) armchair graphene nanoflakes, EB values are 8.03 and 8.05 eV per atom, respectively. For both the edge structures, the C–C and C–H bond lengths are 1.43 and 1.02 Å, respectively. A finite gap between the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) is found to be sensitive to the edges of the nanoflakes;[72] a higher gap for the armchair nanoflake (∼4.4 eV) is predicted as compared to that for the zigzag nanoflake (∼3.5 eV) which displays the localized edge states [Fig. S2, ESI†].

For h-BN nanoflakes, both the (∼22 Å) zigzag and (∼24 Å), armchair configurations have nearly the same B–N bond length of 1.45 Å, and the binding energy of 7.05 eV. The HOMO–LUMO gap of about ∼10.3 eV is found to be insensitive to the edge structure, which is in agreement with, previously reported results.[73] This is also the case with (∼24 Å) zigzag and (∼22 Å) armchair phosphorene nanoflakes who do not show edge-dependent properties; the calculated P–P bond length and the binding energy values, and HOMO–LUMO gaps are 2.37 Å, ∼3.70 eV, ∼6.6 eV, respectively. Overall, the calculated values of the structural and electronic properties of graphene, BN, and phosphorene nanoflakes are in excellent agreement with previously reported values.

Table 2 lists the NLO polarizabilities of graphene, BN, and phosphorene nanoflakes. Due to the absence of inversion symmetry in the armchair graphene flake, we find that β (armchair) > β (zigzag) as expected. We also note that γ = γ due to the symmetry of its structure. Furthermore, γ values calculated for the zigzag and armchair graphene nanoflakes are consistent with previously reported values.[74,75] For example, the hybrid DFT calculations reported γ to be 210 (×10−36 esu) for a (C84H24) armchair nanoflake with a width of ∼16 Å (ref. 74) compared to our wB97XD value of γ of 1415 (×10−36 esu) for the (C180H72) armchair flake with a width of ∼24 Å (Table 2). For phosphorene, calculations predict relatively small values of γ, which may be associated with its puckered ground state configuration.[45] Note that the delocalized π-electron cloud in graphene yields significantly higher values of γ in contrast to those in the case of h-BN and phosphorene nanoflakes.

Dipole moment, HOMO–LUMO gap, and (hyper) polarizabilities of zigzag and armchair nanoflakes of graphene, BN monolayer, and phosphorene calculated using the wB97XD functional form and 6-31G basis set. μ is a dipole moment. α is the linear polarizability, β is the first-hyperpolarizability and γ is the second-hyperpolarizability. The frequency-dependent (hyper) polarizabilities are calculated at λ = 1064 nm

		Zigzag			Armchair
		Graphene (C150H60)	BN (B75N75H60)	Phosphorene (P96H24)	Graphene (C180H72)	BN (B90N90H72)	Phosphorene (P84H24)
Symmetry		C 2H	C S	C 1	C 1	C S	C 1
μ (D)		0	0.1	0.3	0	0.1	0
E (HOMO–LUMO gap) (eV)		3.9	10.3	6.6	4.2	10.3	6.7
Polarizabilities
〈α〉 (×10−22 esu)	(0)	3.8	1.9	5.1	4.2	2.3	4.4
	(−ω)	4.2	1.9	5.3	4.6	2.3	4.5
β (×10−33 esu)	(0)	∼0	39.6	9.6	6.4	35	11.1
	(−ω)	∼0	40.7	11.4	7.6	213.8	14.8
	(−2ω)	∼0	43.1	14.5	7.5	224.3	16.4
γ xxxx /γyyyy (×10−36 esu)	(0)	1285.0	45.0	95.0	1415.0	56.7	73.3
	(−ω)	2341.7	48.3	128.3	2301.7	60.0	98.3
	(−2ω)	3585.0	55.0	175.0	11 665.0	68.3	135.0
γ yyxx (×10−36 esu)	(0)	428.3	14.9	31.7	473.3	18.3	25.0
	(−ω)	1161.7	16.1	43.3	1013.3	20.0	33.3
	(−2ω)	21 836.7	18.5	61.7	5490.0	23.3	48.3

Nanoflake heterostructures

Graphene/BN nanoflake heterostructures

Fig. 2 and 3 display the AA- and AB-stacked graphene/BN heterostructures considered in the zigzag and armchair configurations, respectively. In the AA stacking, the atoms in the upper layer sit directly atop the atoms in the bottom layer. The AB stacking is formed by shifting one layer of the AA stacked configuration by half of the cell along with one of the in-plane directions.

Fig. 2

Top views of (a) AA- and (b) AB-stacked zigzag graphene/BN nanoflake heterostructures. Color code: C: green, B: purple, N: orange, H: blue.

Fig. 3

Top views of (a) AA- and (b) AB-stacked armchair graphene/BN nanoflake heterostructures. Color code: C: green, B: purple, N: orange, H: blue.

For the zigzag heterostructure, the calculated interlayer binding energy of the AA-stacked and AB-stacked configurations with respect to the constituent monolayers is 0.23 eV per atom and 0.19 eV per atom, respectively. The results suggest that both stacking configurations are nearly degenerate. This is also the case with the armchair nanoflake where the energy difference between AA- and AB-stacking configurations being 6 meV. In both, the AA- and AB-stacked heterostructures, as noted earlier, the intra-layer atomic bond-distance is 1.43 Å for graphene and is 1.45 Å for the h-BN monolayer. The calculated interlayer distance for the graphene/BN heterostructure is 3.4 Å and 3.2 Å for the AA- and AB-stacked configurations, respectively. The bandgap is calculated to be 3.44 and 3.86 eV for the AA- and AB-stacked zigzag heterostructures, respectively. These values are in good agreement with the previously reported results for graphene/BN bilayer heterostructure.[14]

The calculated values of the NLO polarizabilities for the graphene/BN nanoflake heterostructures are listed in Table 3. β values for the heterostructure, while significantly larger than those calculated for the pristine monolayers (Fig. 4), are higher in magnitude for the AB- over AA-stacked configurations. The calculated results lack a clear trend due to the flake edge difference (zigzag vs. armchair) in the heterostructures. For example, whereas β (0) and β (−ω) for armchair bilayer are considerably larger than their values in zigzag nanoflakes, the reverse is true for β (−2ω), which has much higher magnitudes in the latter arrangement than in the former. We also note that, while the expected order in the frequency-dependent value of β, i.e. β (−2ω) > β (−ω) > β (0), is generally maintained for zigzag (AA as well as AB stacking) bilayers, the armchair nanoflakes only show this for the AB stacking. The calculation fails to predict the correct order between β (−ω) and β (−2ω) for the armchair bilayer in the AA stacking. Judging from the correct qualitative and quantitative trends predicted by the current method for PNA molecule, as discussed in the previous section, it should be pointed out that the present calculations predict considerably high β (−ω) coefficients at ∼105 × 10−33 esu for the AB-stacked graphene/BN bilayer, regardless of the edge geometry. The experimentally measurable β, which is zero for the pristine nanoflakes, is calculated to be as high as ∼104 × 10−33 esu for the heterostructure.

Structure, interlayer distance, interlayer binding energy, HOMO–LUMO gap, dipole moment, and (hyper) polarizabilities of zigzag and armchair nanoflakes of graphene/BN bilayer calculated with wB97XD functional and 6-31G basis set. α, β, and γ have the same meaning as in Table 2. The frequency-dependent (hyper) polarizabilities are calculated at λ = 1064 nm

Graphene/BN nanoflake heterostructures
		Zigzag		Armchair
		AA-stacked	AB-stacked	AA-stacked	AB-stacked
Symmetry		C 3V	C 1	C 1	C 1
Interlayer distance, Å		3.4	3.2	3.4	3.2
Interlayer binding energy, eV		0.23	0.19	0.19	0.19
E (HOMO–LUMO gap), eV		3.4	3.9	4.2	3.9
Dipole moment, μ (D)		3.3	3.5	3.4	4.3
Polarizabilities
〈α〉 (×10−22 esu)	(0)	5.4	5.0	5.8	5.8
	(−ω)	6.1	5.4	6.1	6.2
β (×10−33 esu)	(0)	5620.8	23 521.8	6993.4	26 221.2
	(−ω)	7116.5	30 028.1	22 539.2	34 760.1
	(−2ω)	78 504.8	596 865.1	12 657.4	103 297.9
γ xxxx /γyyyy (×10−36 esu)	(0)	1630.0	1043.3	1108.3	1081.7
	(−ω)	4160.0	1988.3	1858.7	1796.7
	(−2ω)	−695.0	75 535.0	12 304.7	10 956.7
γ yyxx (×10−36 esu)	(0)	542.6	335.0	596.7	588.3
	(−ω)	2543.7	958.3	1083.3	1051.7
	(−2ω)	10 659.5	49 201.7	7276.7	4450.0

Fig. 4

The frequency-dependent β for the pristine and zigzag bilayers calculated at a fundamental wavelength, λ of 1064 nm. C is graphene, BN is boron nitride, and P is phosphorene.

Unlike β, the calculated γ coefficients exhibit a rather higher degree of consistency in magnitude with respect to the stacking as well as edge configurations, except γ (−2ω) for AA-stacked zigzag bilayer, which is calculated to have a negative value that may be associated with the self-defocusing effect of incident light.[67] In general, the in-plane principal components (γ/γ) for the AA-stacked bilayers have larger values than their AB-stacked counterparts do, with the zigzag configuration giving higher values than the armchair. For the γ component, no clear pattern emerges. This may be due to the sensitivity of third harmonics to the details of the edge-dependent gap. On the other hand, β values are inversely proportional to the energy gap given by the HOMO–LUMO gap in the present case.

Graphene/phosphorene nanoflake heterostructures

Since the lattice constants of graphene and phosphorene are significantly different, we employ the “commensurate” supercell approach for calculations that minimizes the lattice mismatch in the hetero-bilayers. Under this constraint, the width of a phosphorene nanoflake consisting of the P6 rings was taken to be the same as that of the zigzag or armchair graphene nanoflake in a heterostructure configuration [Fig. 5].

Fig. 5

Top views of graphene/phosphorene heterostructures: (a) armchair and (b) zigzag. Color code: C: green, P: dark blue, H: light blue.

Table 4 lists some of the structural and electronic properties of the graphene/phosphorene nanoflake hetero-bilayers. The zigzag and armchair configurations have nearly the same interlayer binding energy of 0.2 eV per atom, and an interlayer distance of about 3.5 Å. The calculated HOMO–LUMO gap of 3.9 and 4.2 eV for the zigzag and armchair heterostructures, respectively, exhibit a slight edge-dependence. It is also noteworthy that the calculated dipole moments are significantly larger than the dipole moments calculated for the pristine monolayer nanoflakes, suggesting an interlayer coupling-induced large asymmetric charge distribution in the heterostructures. This, in turn, results in significant enhancement of β values, which are significantly larger in magnitude compared to those predicted for graphene/BN heterostructures.

Structure, interlayer separation, binding energy, dipole moment, and (hyper) polarizabilities of zigzag and armchair nanoflakes of graphene/phosphorene hetero-bilayer calculated with wB97XD functional and 6-31G basis set. α, β, and γ have the same meaning as in Table 2. All dynamic (hyper) polarizabilities are calculated at a fundamental wavelength, λ = 1064 nm

Graphene/phosphorene nanoflake heterostructures
		Zigzag	Armchair
Symmetry		C 1	C 1
Interlayer distance, Å		3.5	3.6
Interlayer binding energy, eV		0.21	0.20
E (HOMO–LUMO gap) eV		3.9	4.2
Dipole moment, μ (D)		8.6	8.2
Polarizabilities
〈α〉 (×10−22 esu)	(0)	7.3	7.1
	(−ω)	7.7	7.5
β (×10−33 esu)	(0)	102 595.5	80 663.0
	(−ω)	130 010.6	94 616.0
	(−2ω)	3 624 891.0	132 134.0
γ xxxx /γyyyy (×10−36 esu)	(0)	576.7	615.0
	(−ω)	1078.3	1006.7
	(−2ω)	−190 798.3	6452.6
γ yyxx (×10−36 esu)	(0)	190.00	205.0
	(−ω)	525.7	436.7
	(−2ω)	118 491.2	3068.3

It is worth mentioning that the β (−2ω) value of 3 624 891 × 10−33 esu calculated for the zigzag graphene/phosphorene hetero-bilayer structure is comparable to the β (−2ω) value of the widely used NLO materials; e.g. β333 at 1064 nm for potassium titanyl phosphate (KTP) is measured[76] to be 31 600 × 10−33 esu. In the organic crystals, one of the largest β is reported to be 2 480 000 × 10−33 esu associated with thiophene compounds.[77,78]

The calculated results clearly show that the NLO coefficients for the graphene/phosphorene heterostructures have significantly higher values than those for graphene/BN. It can be associated with the anisotropic nature of phosphorene leading to distinctly different structural and electronic properties compared to those exhibited by the atomically flat BN nanoflakes. In phosphorene, for example, HOMO is associated with the localized (lone pair) p orbitals, and LUMO is associated with delocalized sp3 hybridized orbitals. This is not the case with BN nanoflakes, where the frontier molecular orbitals (HOMO, LUMO) describe the delocalized multi-center bonding feature. Furthermore, we show a simplified model of the chemical bonding in the pristine graphene, h-BN monolayer, and phosphorene in Fig. S4 (ESI†). In the case of phosphorene, the lone pair electrons lead to a highly polarizable delocalized electron charge density. This is reflected in the charge density plots [Fig. S5 (ESI†)], clearly showing the delocalized nature of electron density in graphene/phosphorene heterostructure.

In general, an occurrence of the light-induced second harmonic depends on the dipole-allowed transition between frontier orbitals, facilitated by a non-centrosymmetric structure. Its magnitude depends on the transition's oscillator strength, which is determined by the extent of overlap of the electronic wave functions associated with valence and conduction bands of the material and the electronic bandgap. The degree of overlap, in a bilayer configuration, is associated with the strength of the interlayer electronic coupling between the constituent monolayers.[79] This can be also verified by performing calculations of the second harmonic coefficient by varying the interlayer distance. A decrease in the interlayer distance corresponds to an increased overlap between the wave functions of the individual monolayers (i.e. strength of the interlayer electronic coupling) as reflected in the increased value of the calculated dipole moment of graphene/BN (Table S4, ESI†). The calculated results for the zigzag AA-stacked graphene/BN show that a mere 10% decrease in the interlayer distance leads to about 50% enhancement in the second harmonic coefficient (Fig. S6, ESI†).

To facilitate a comparison of the calculated β values with experiments, we now estimate the SHG intensity, defined in terms of the macroscopic second-order susceptibility, χ(2), using the following expression:[80]where V is the volume of a heterostructure configuration which can be approximated as (tlayer × Rinterlayer). Here t is layer thickness and Rinterlayer is the interlayer separation. The parameters employed in calculations are listed in Table S5 (ESI†).

Fig. 6 displays the SHG intensity that follows the trend exhibited by β values in showing significant enhancement in heterostructures relative to pristine nanoflakes. For the zigzag case, the estimated SHG values are 2.2, 17.1, and 50 pm V−1 for AA-stacked graphene/BN, AB-stacked graphene/BN, and graphene/phosphorene heterostructures. On the other hand, the estimated SHG values are 0.3, 2.4, and 1.5 pm V−1 for AA-stacked graphene/BN, AB-stacked graphene/BN, and graphene/phosphorene armchair heterostructures. Therefore, it is worth stressing that the zigzag nanoflake heterostructures exhibit a strong SHG response, which is comparable to that exhibited by the widely used inorganic crystals, e.g. KTP.[74]

Fig. 6

The calculated SHG values in terms of χ(2) for the graphene-based heterostructures at λ = 1064 nm.

Summary

The NLO response of the graphene-based nanoflake heterostructures was investigated using time-dependent density functional theory. Since a graphene bilayer configuration is centrosymmetric, in-plane asymmetry is induced by forming the bilayer heterostructures with h-BN or phosphorene in either zigzag or armchair configurations. Test calculations on PNA and comparison with literature suggest that the wB97XD functional together with the 6-31G basis set are adequate to provide reliable values microscopic NLO coefficients.

The results show that the interlayer electronic coupling of graphene with h-BN or phosphorene significantly enhances the first-order NLO response in a given nanoflake hetero-bilayer structure. Moreover, the delocalized π-electron cloud in graphene leads to significantly higher values of the second-order NLO coefficients, γ, in these heterostructures. Our study clearly shows that graphene-based nanoflake heterostructures giving large NLO coefficients together with high electron mobility of these materials offer new opportunities as candidate materials of choice for next-generation photonics and integrated quantum technologies.

Conflicts of interest

There are no conflicts to declare.