Letter pubs.acs.org/NanoLett

Giant Optical Second Harmonic Generation in Two-Dimensional Multiferroics Hua Wang† and Xiaofeng Qian*,† †

Department of Materials Science and Engineering, College of Engineering and College of Science, Texas A&M University, College Station, Texas 77843, United States S Supporting Information *

ABSTRACT: Nonlinear optical properties of materials such as second and higher order harmonic generation and electro-optic effect play pivotal roles in lasers, frequency conversion, electro-optic modulators, switches, and so forth. The strength of nonlinear optical responses highly depends on intrinsic crystal symmetry, transition dipole moments, specific optical excitation, and local environment. Using first-principles electronic structure theory, here we predict giant second harmonic generation (SHG) in recently discovered two-dimensional (2D) ferroelectric−ferroelastic multiferroics−group IV monochalcogenides (i.e., GeSe, GeS, SnSe, and SnS). Remarkably, the strength of SHG susceptibility in GeSe and SnSe monolayers is more than 1 order of magnitude higher than that in monolayer MoS2, and 2 orders of magnitude higher than that in monolayer hexagonal BN. Their extraordinary SHG is dominated by the large residual of two opposite intraband contributions in the SHG susceptibility. More importantly, the SHG polarization anisotropy is strongly correlated with the intrinsic ferroelastic and ferroelectric orders in group IV monochalcogenide monolayers. Our present findings provide a microscopic understanding of the large SHG susceptibility in 2D group IV monochalcogenide multiferroics from first-principles theory and open up a variety of new avenues for 2D ferroelectrics, multiferroics, and nonlinear optoelectronics, for example, realizing active electrical/optical/mechanical switching of ferroic orders in 2D multiferroics and in situ ultrafast optical characterization of local atomistic and electronic structures using noncontact noninvasive optical SHG techniques. KEYWORDS: 2D materials, multiferroics, second harmonic generation, first-principles theory, group IV monochalcogenides, nonlinear optical properties

M

noncentrosymmetric point group D3h (6̅m2) with effectively only one independent SHG susceptibility tensor element. It has also been observed along the edge of MoS2.20 The SHG responses can be significantly enhanced by electric control21 and resonant excitonic excitation.22 Despite their monolayer or few layer nature, the effective nonlinear SHG susceptibility of hBN is close to the well-developed lithium niobate nonlinear crystals widely used in integrated optical waveguide and in MoS2 monolayer it is even one order magnitude higher than that of h-BN and lithium niobate.6 These discoveries have enkindled further investigations of 2D nonlinear optical materials that are free of phase-matching bottleneck and

aterials under intense optical excitation may generate large nonlinear optical responses,1−3 depending on the intrinsic crystalline symmetry, microscopic transition dipole matrix, and specific frequency and orientation of optical field applied. In centrosymmetric crystals, all even order electric susceptibility tensors vanish, as both electric field and electric polarization are polar vectors with odd parity. Noncentrosymmetric materials with large second-order electric susceptibilities such as optical second harmonic generation (SHG),4 linear electro-optic or Pockels effect, optical rectification, and sum frequency generation are highly valuable for a variety of applications such as lasers, frequency conversion, electro-optic modulators and switches, and so forth. Large SHG and third harmonic generation were recently discovered in a number of two-dimensional (2D) materials, including monolayer/multilayer MoS2,5−9 MoSe2,10−12 WS2,13−15 WSe2,13,16 hexagonal BN (h-BN),6,17 GaSe,18 and InSe.19 All of them belong to the © 2017 American Chemical Society

Received: May 30, 2017 Revised: June 30, 2017 Published: July 3, 2017 5027

DOI: 10.1021/acs.nanolett.7b02268 Nano Lett. 2017, 17, 5027−5034

Letter

Nano Letters

imation (GGA)38 and a Monkhorst−Pack k-point sampling for the Brillouin zone (BZ) integration. An energy cutoff of 400 eV for the plane-wave basis and a Monkhorst−Pack k-point sampling of 10 × 10 × 1 were applied. Dynamical stability of the crystal structures was confirmed by phonon dispersion from density functional perturbation theory calculations.39 A separate stand-alone package has been developed in the group to compute the second-order SHG susceptibility tensor40−42 interfaced with first-principles packages. A scissor operator was applied to account for the underestimated band gap in the DFT-GGA calculations. The Kramers−Kronig relation of the calculated SHG susceptibility tensors was verified, satisfying the causality condition. Benchmark calculations were performed to test the convergence with respect to number of bands, k-point sampling, and so forth. A dense k-point sampling of 72 × 72 × 1, 40 electronic bands, and a total of 1000 frequency grids between the energy range of [−6 eV, 6 eV] are enough for achieving converged SHG susceptibility tensor, except that 2000 frequency grids between [−12 eV, 12 eV] were used for h-BN due to the large optical gap. Only positive frequency range is presented because of the reality condition. The fundamental frequency ω in the denominator of susceptibility tensor carries a small imaginary smearing factor δ: ω → ω + iδ, and δ = 0.05 eV was used in this work. All symmetry-allowed and symmetry-disallowed SHG tensor elements have been extensively checked for noncentrosymmetric point groups C2v (mm2) and D3h (6̅m2). That is, the 23 SHG susceptibility elements for D3h and 20 elements for C2v are zero at all frequencies. Finally, as the thickness of 2D materials is not well-defined, sheet SHG susceptibility tensors (with the unit of pm2/V) are reported in this work. This is because every term of χSHG (−2ω; ω, ω) shown in the Supporting Information requires an integration in the first BZ. In practice, this is achieved by the discrete sum with the dense k-point sampling in the first BZ:

promising for applications in 2D nonlinear optics, for example, ultrathin nonlinear optical devices and spectroscopies.21 Very recently, group IV monochalcogenides (i.e., GeSe, GeS, SnSe, and SnS, denoted by MX with M = Ge, Sn and X = Se, S) were predicted to be 2D ferroelastic-ferroelectric multiferroics23−26 with giant spontaneous electric polarization23,24,27 and spontaneous lattice strain23−26 as well as sizable piezoelectricity.24,28 Ferroelectricity was experimentally observed by Chang et al.29 in an MX cousin, SnTe monolayer, with a similar type of spontaneous lattice distortion and ferroelectric polarization. Moreover, these 2D multiferroic materials possess coupled structure-optical properties owing to their correlated noncentrosymmetry and ferroelastic order.24,30−33 The presence of ferroelectricity necessitates the breaking of inversion symmetry, therefore in principle these 2D multiferroics can have finite second-order nonlinear electric susceptibilities, while specific nontrivial elements depend on the underlying point group symmetry. These MX monolayers belong to the noncentrosymmetric point group C2v (mm2), and have up to five independent SHG susceptibility tensor elements. Quite interestingly, their bulk holds a centrosymmetric point group D2h (mmm), hence has vanishing secondorder susceptibilities. In fact, MX with odd number of layers can have finite SHG response; conversely MX with an even number of layers have a vanishing SHG response like their “infinite” bulk. Nevertheless, neither theoretical nor experimental SHG has been reported for these 2D multiferroic materials. Here, using first-principles electronic structure theory we show that these group IV monochalcogenides possess giant SHG susceptibility, much higher than that in MoS2 and h-BN monolayers. The large SHG polarization anisotropy is pertinent to the ferroelastic order and noncentrosymmetric point group, and their extraordinary SHG responses are dominated by the large residual of two opposite intraband contributions in the SHG susceptibility governed by interband and intraband Berry connections. Remarkably, the SHG susceptibility in GeSe and SnSe monolayer is 1 order of magnitude higher than that in MoS2 monolayer, and 2 orders of magnitude higher than that in h-BN monolayer. In contrast, MoS2 and h-BN monolayers are largely determined by a single intraband term different from 2D MX. The unique ferroelasticity and ferroelectricity in these 2D multiferroics allow direct mechanical switching of lattice orientation and electrical switching of electrical polarization, accompanied by instantaneous direction and phase switching in the SHG susceptibility tensor, respectively. Our present findings not only predict colossal SHG nonlinear optical responses in group IV monochalcogenide monolayers and provide a microscopic understanding from first-principles theory but also open up a variety of new avenues for 2D ferroelectrics, multiferroics, and nonlinear optoelectronics, for example, achieving active control of their ferroic orders with electrical and/or mechanical stimuli and in situ ultrafast optical sensing of local electronic structures and domain evolution in these 2D multiferroics with noncontact noninvasive linear and nonlinear optical spectroscopy/imaging. Computational Methods. Ground-state crystal structures of group IV monochalcogenide monolayers were calculated by first-principles density-functional theory (DFT)34 implemented in the Vienna Ab initio Simulation Package35 with a plane-wave basis and the projector-augmented wave method.36 Here we used the Perdew−Burke−Ernzerhof (PBE) form37 of exchangecorrelation functional within the generalized-gradient approx-

3

∫ (2dπk)3 → ∑k Ω1 , where Ω is the volume of unit cell. For 2D

materials, Ω becomes ill-defined for arbitrary simulation cell length along plane norm (Lz assuming z is parallel to plane 1 normal). A more meaningful definition will be to use ∑k A , where A is the area of the 2D plane in the unit cell, and correspondingly the k-point sampling is on the 2D plane of the first BZ. It essentially leads to the following formula for sheet SHG susceptibility tensors we reported in this work, which is similar to the sheet conductance/resistance used for thin-film SHG materials and 2D materials: χSHG sheet (−2ω; ω, ω) = χbulk (−2ω; ω, ω) × Lz. Physically it means the ability of SHG response coming from a single 2D material. When comparing with other nonlinear optical bulk materials, one can define an effective Lz. This is often done by including the thickness of 2D material (d2D) and the van der Waals thickness on both sides of 2D material (approximated by ∼3.4 Å on each side), that is, Lz = 3.4 × 2 + d2D. However, it should be noted that this is certainly a rough estimation where the effect of symmetry and interlayer coupling has not taken into account. More detailed discussion can be found in Supporting Information. Additional SHG calculations were performed to check the impact of spin−orbit coupling which turns out to be small for the materials studied in this work, and the results are included in Supporting Information as well. Results and Discussion. The crystal structure of group IV monochalcogenide monolayers consists of two atomic layers

5028

DOI: 10.1021/acs.nanolett.7b02268 Nano Lett. 2017, 17, 5027−5034

Letter

Nano Letters

Figure 1. Atomic configuration and electronic structures of group IV monochalcogenide MX monolayers. (a) Top and side views of 2D MX. Dashed orange lines indicate the primitive cell with spontaneous polarization along the y-direction. (b) The 2D surface plot of electronic band structure in GeSe monolayer in the first Brillouin zone: CBM (upper) and VBM (lower). Λ denotes the location of energy gap in the reciprocal unit cell, Λ(kx, ky) = [0, ±0.39] in the fractional coordinate unit. (c) Quasiparticle band structure of GeSe monolayer with the quasiparticle band gap, the lowest exciton, and the exciton binding energy labeled in the plot.

Figure 2. Sheet SHG susceptibilities in GeSe, MoS2, and h-BN monolayers with their magnitude, imaginary, and real component. (a) GeSe monolayer which has seven nonzero susceptibility tensor elements with five independent ones due to its point group C2v. (b) MoS2 monolayer which has four SHG elements with only one independent element due to its point group D3h. (c) h-BN monolayer with only one independent element due to the same point group. Black dots indicate the experimental values.7

357, 260, and 181 pC/m for GeS, GeSe, SnS, and SnSe, respectively.24 In addition, the puckering also leads to highly anisotropic electronic band structure, as displayed in Figure 1b for the two bands near the Fermi level, namely conduction band minimum (CBM) and valence band maximum (VBM). The calculated DFT-GGA band gap is 1.13 eV at two k-points Λ(kx, ky) = [0, ±0.39] in the fractional coordinate unit (marked in Figure 1b,c). As GeSe monolayer is a 2D semiconductor, the reduced dimensionality enhances electron−hole Coulomb interaction and increases the density of states at the band edge, and consequently, large excitonic optical absorption is expected.24,30,43 Figure 1c displays the quasiparticle band structure of GeSe monolayer using many-body quasiparticle approach44 and first-principles quasiatomic orbital method45,46 with the exciton level marked by a red line. The computed exciton energy at the Λ point is 1.26 eV, which is only 0.13 eV higher than the DFT-GGA gap. The small difference is due to

that are puckered along either x- or y-direction, as their nondistorted parent structure has 4-fold rotation symmetry and four mirror planes. A representative 2D GeSe crystal structure is shown in Figure 1a. Here, the z-axis is oriented along the plane normal of 2D MX. Our previous calculation indicates large spontaneous lattice strain η in MX monolayers depicted by a transformation strain based on the Green-Lagrange strain tensor.24 Spontaneous tensile strain varies from 1% to 14% and the associated spontaneous compressive strain ranges from −1% to −7%. For the GeSe monolayer shown in Figure 1a, its distortion is along the y-direction, and consequently ηyy = 4.1% and ηxx = −2.7%. Such large spontaneous lattice strain has a nontrivial consequence in the electronic structure and crystal symmetry. Particularly, the puckering causes inversion symmetry breaking and results in the noncentrosymmetric polar point group C2v (mm2), consequently these MX monolayers exhibit large spontaneous polarization, i.e. 484, 5029

DOI: 10.1021/acs.nanolett.7b02268 Nano Lett. 2017, 17, 5027−5034

Letter

Nano Letters

Figure 3. Polarization anisotropy of sheet SHG susceptibilities in (a) GeSe, (b) MoS2, and (c) BN. Red (blue) line indicates the polarization component of the SHG response parallel (perpendicular) to the polarization E(ω) of the incident electric field. θ is the rotation angle between E(ω) and the crystal lattice (i.e., the x-axis here for all three 2D materials). Their maximum values are indicated by the green circles in the polar plot and listed at the bottom.

monolayer MoS2, and 6.38 × 104 pm2/V at 3.22 eV h-BN, demonstrating much higher SHG response in GeSe than that in MoS2 and h-BN. Moreover, GeSe has another independent (2) SHG tensor element χ(2) xyx = χxxy with a substantial magnitude of 6 2 more than 4 × 10 pm /V, and all the other three independent elements are much smaller. The computed frequency dependent SHG susceptibility satisfies the Kramers−Kronig relation as shown in Figure S3. It is worth noticing that although the energy gaps are corrected with respect to exciton energies, the computed SHG response is sensitive to the smearing factor, and the independent-particle approximation may also affect the accuracy. Consequently, the absolute magnitude may differ from experiment. Indeed, the calculated SHG response is higher than the measured SHG values MoS2 shown as dots in Figure 2b. Nonetheless, our results on MoS2 and h-BN monolayers are in good agreement with recent works within the independent-particle approximation.47−49 Meanwhile, the same smearing factor was used for all calculations of the six 2D materials studied to keep them on an equal footing. Very interestingly, the fundamental frequency for the first large peak of SHG susceptibility χ(2) yxx covers a wide range of optical spectrum, that is, 0.66, 1.18, 0.56, and 1.00 eV for GeSe, GeS, SnSe, and SnS, respectively, as shown in Figure S2. Nonlinear SHG responses can be characterized by shedding linearly polarized laser beam onto the materials and measuring different polarization component of the outgoing SHG response through an analyzer. By inspecting its angular dependence, for example, by rotating the samples, their crystallographic orientation and SHG polarization anisotropy and intensity can be determined.5−7 The SHG response in an arbitrarily oriented sample stems from all symmetry-allowed SHG nonlinear susceptibility tensor elements taking into account the relative orientation between the incident beam and the crystal lattice via proper Euler angles. Here, we

the cancellation of quasiparticle energy correction to the DFTGGA gap (∼0.5 eV) and the large exciton binding energy (∼0.4 eV). To correct the underestimated band gap obtained from DFT-GGA calculations, we have used the quasiparticle GW gap and exciton binding energy to correct the optical gap of monolayer MX, MoS2, and BN. This is a reasonable remedy as the GW band structures are similar to DFT band structure for the 2D materials investigated in this work. With the major difference taken into account by the scissor shift in optical band gap, the dipole matrix elements after the correction should be reasonable. Future efforts shall be made to include quasiparticle GW electronic structure and exciton effect for monolayer MX in a computation and memory efficient way. The unique electronic structures and multiferroic nature of monolayer MX, particularly the absence of inversion symmetry, suggests a possibility of large nonlinear optical responses in group IV monochalcogenide monolayers. Under incident electric field E(ω) with fundamental frequency ω, the second-order nonlinear polarization P(2ω) at frequency 2ω is determined by a third-rank electric susceptibility tensor χ: Pi(2ω) = χijk(−2ω; ω, ω)Ej(ω)Ek(ω) which is further correlated with the emitting SHG field. Figure 2a−c shows the calculated sheet SHG susceptibility tensors in GeSe, MoS2, and h-BN monolayers, and the calculated SHG data for GeS, SnSe, and SnS are presented in Figures S1 and S2. Both MoS2 and BN belong to the D3h (6̅m2) point group with effectively only one independent nontrivial SHG susceptibility (2) (2) (2) tensor element that satisfies χ(2) yxx = χxyx = χxxy = −χyyy . Group IV monochalcogenides belong to the point group C2v (mm2) as mentioned earlier, thus they have five independent SHG (2) (2) (2) (2) (2) susceptibility tensor elements: χ(2) yxx, χyyy , χyzz , χxyx = χxxy, and χzzy (2) = χzyz . The magnitude of the calculated sheet SHG 6 2 susceptibility tensor element χ(2) yxx is 5.16 × 10 pm /V at 0.66 5 2 eV for monolayer GeSe, 3.02 × 10 pm /V at 1.67 eV for 5030

DOI: 10.1021/acs.nanolett.7b02268 Nano Lett. 2017, 17, 5027−5034

Letter

Nano Letters

χabc mod. These three contributions can be further classified into six terms, according to the location of the poles (either at ω or at v r abc 2ω): χabc intra ≡ χi(ω) + χi (2ω) + χi (2ω), χinter ≡ χe(ω) + χe(2ω), abc and χmod ≡ χm(ω). Detailed expressions of each term are given in the Supporting Information. Figure 4a,b presents the

consider a normal incidence geometry for the six 2D materials. As a result, the angular-dependent SHG susceptibilities for the point group C2v (mm2) and D3h (6̅m2) are given by (2) (2) (2) χ (2) (θ ; C2v) = (χxyx + χyxx )sin θ cos 2 θ + χyyy sin 3 θ (2) (2) (2) χ⊥(2) (θ ; C2v) = χyxx cos3 θ + (χyyy − χxyx )cos θ sin 2 θ (2) (2) χ (2) (θ ; D3h) = −χyxx sin 3θ , χ⊥(2) (θ ; D3h) = χyxx cos 3θ

where “∥” (“⊥”) stands for the polarization components of the SHG response parallel (perpendicular) to the polarization of incident electric field E(ω) at fundamental frequency ω. “θ” is the rotation angle between E(ω) and the crystal lattice (e.g., the x-axis here). For GeSe monolayer, (2) χyxx = (18.2 − 515.6·i) × 104

pm 2 V

(2) (2) χxyx = χxxy = ( −232.9 − 233.5·i) × 104

(2) χyyy = (9.0 − 11.2·i) × 104

(2) χyzz = ( −3.6 + 5.3·i) × 104

pm 2 V

pm 2 V pm 2 V

(2) (2) χzzy = χzyz = ( −1.1 + 0.1·i) × 104

pm 2 V

Figure 4. Frequency-dependent (a) imaginary and (b) real part of the interband (χe(ω), χe(2ω)), intraband (χi(ω), χi(2ω)), and modulation (χm(ω)) contributions to the sheet SHG susceptibility χ(2) yxx. ω and 2ω denote the poles of the specific terms. The blue and red dots indicate two major 2ω terms from the intraband contribution, χvi (2ω) and χri (2ω). Total imaginary and real part of the SHG susceptibility χ(2) yxx are marked by the black line.

and for ω = 0.66 eV. For MoS2 monolayer, (2) (2) (2) (2) χyxx = χxyx = χxxy = −χyyy = ( −1.7 − 30.2·i)

× 104

pm 2 V

for ω = 1.67 eV. For h-BN monolayer, (2) (2) (2) (2) χyxx = χxyx = χxxy = −χyyy = ( −2.8 − 5.7·i) × 104

pm V

(2) of GeSe frequency dependent SHG susceptibility χyxx monolayer and the contributions from the six terms with imaginary and real part, respectively. Among all the six terms, χvi (2ω) and χri (2ω) are two leading ones responsible for the giant SHG susceptibility in the GeSe monolayer governed by interband and intraband Berry connections. The explicit expressions are given below, which were originally derived by Sipe et al. using the generalized position operator for periodic crystals50,51

2

2 for ω = 3.22 eV. The SHG intensity is proportional to |χ(2) ∥ (θ)| 2 (2) and |χ⊥ (θ)| , hence their corresponding polar plots are presented in Figure 3a−c together with the corresponding frequencies and maximum values. The red and blue lines are 2 2 (2) |χ(2) ∥ (θ)| and |χ⊥ (θ)| , respectively. The results clearly show that the GeSe monolayer exhibits highly polarized colossal SHG response with its maximum value located on the perpendicular component at θ = 0 (corresponding to χ(2) yxx), which is also evident in the formula above. The polar plots of MoS2 monolayer in Figure 3b precisely reflect their D3h symmetry, which is in excellent agreement with several recent experiments.5−7 Because SHG susceptibility varies with fundamental frequency, the SHG polarization anisotropy will change accordingly, as illustrated in Figure S4. The giant SHG susceptibility and polarization anisotropy in MX monolayers can be observed by noncontact noninvasive SHG spectroscopy with lower pump beam intensity, thereby achieving accurate determination of their 2D crystal orientation. In general, the total SHG susceptibility χabc total (−2ω; ω, ω) contains an interband contribution (χabc ), a modification by inter intraband motion (χabc intra), and a modulation by interband 40,41 abc abc motion (χabc that is, χabc mod), total (−2ω; ω, ω) = χinter + χintra +

χi v (2ω) =

e3 ℏ2

⎛

∫ 4dπk3 ⎜⎜−8i ∑ ⎝

nm

a b c b c rnm ({vmm rmn} − {rmn vnn}) 2 ωmn

⎡ ⎤⎞ fnm ⎢ ⎥⎟⎟ ⎣ ωmn − 2ω ⎦⎠ χir (2ω) =

e3 ℏ2

⎛

∫ 4dπk3 ⎜⎜−2i ∑ ⎝

nml

a b c b c rnm ({vml rln} − {rml v ln}) 2 ωmn

⎡ ⎤⎞ fnm ⎢ ⎥⎟⎟ ⎣ ωmn − 2ω ⎦⎠ i(2π )3

* (r )∇k umk(r ) is the matrix element of Here, rnm = Ω ∫ dr 3unk position operator between the cell periodic part of eigenstates n 5031

DOI: 10.1021/acs.nanolett.7b02268 Nano Lett. 2017, 17, 5027−5034

Letter

Nano Letters

v r Figure 5. k-point dependent distribution of sheet SHG susceptibility χ(2) yxx and 2ω intraband terms χi (2ω) and χi (2ω) in the first BZ of GeSe monolayer. (a,b) Imaginary and real part of χ(2) yxx dominated by the contributions from two equivalent valleys at the two bandgap locations (marked by Λ). (c−f) Imaginary and real part of 2ω intraband terms χvi (2ω) and χri (2ω). (g) The distribution of velocity difference Δmn ≡ vmm − vnn between CBM and VBM in the first BZ. (h,i) Imaginary and real part of non velocity terms in χvi (2ω). (j,k) Imaginary and real part of χri (2ω) from three dominating sets of bands.

and m: umk (r) = e−ikrψmk(r). The superscripts “a, b, c” indicate the Cartesian directions. rnm is equivalent to interband Berry connection when n ≠ m, which is directly related to the interband velocity matrix element vml via vml = iωmlrml. ∂ω a = ∂km is the group Furthermore, ωmn = ωm − ωn, and vmm

Δmn ≡ vmm − vnn is the velocity difference between states m and n. Figure 5 panel g and panels h,i presents the vector field ⎡ −ir y r x ⎤ of Δmn and the imaginary and real part of ⎣⎢ ω 2 (ω nm −mn2ω) ⎦⎥, mn mn respectively. The velocity difference plot shows that close to the Λ points, Δxmn (δkx, Λy) = −Δxmn (−δkx, Λy), and the similar behavior can be seen for the other term owing to the mirror symmetry with a (001) mirror plane. Therefore, the product of the two vector fields for specific yxx component leads to the result shown in Figure 5c,d. Even though the numerator rynm rxmn is only constrained by the symmetry and strength of interband Berry connection, the denominator ω2mn (ωmn − 2ω) becomes relatively much larger away from the bandgap at Λ. Furthermore, Δmn vanishes at two Λ points as the velocity of both VBM and CBM vanishes. All these factors combined give rise to four spots concentrated around two valleys in the first BZ. For χri (2ω), a detailed analysis indicates that it primarily arises from n = VBM, m = CBM, and l = VBM − 2, CBM + 1, and CBM + 2, and the results are shown in Figure 5j,k accordingly. Similarly, the denominator ω2mn (ωmn − 2ω) largely dominate the response and concentrates it around the band gap at the two Λ points. In summary, our results identify that the giant SHG response in GeSe monolayer is controlled by two intraband contributions χvi (2ω) and χri (2ω) located near the two valleys. Despite the opposite contributions, the residual SHG susceptibility tensor elements χ(2) yxx remains extraordinary due to the much larger χri (2ω) term. It is worth to note that bulk MX belong to centrosymmetric point group mmm, hence their nature forbids SHG response. The difference between monolayer and bulk is therefore solely due to the stacking induced emerging inversion symmetry of individual noncentrosymmetric layers. As a result, multilayer MX materials shall exhibit alternating even−odd behavior. We have carried out first-principles calculations for multilayer GeSe

a

velocity of eigenstate m along Cartesian direction a. The curly brackets indicate the intrinsic permutation symmetry χabc (−2ω; ω, ω) = χacb (−2ω; ω, ω) with respect to the Cartesian 1 a b a b b a directions, that is, {vml rln} = 2 (vml rln + rml vln). Finally, f mn ≡ f n − f m, where f n is occupation number of state n according to the Fermi−Dirac distribution. The above two equations clearly show that χvi (2ω) represents the susceptibility modified by intraband motion (evident by vbmm and vcnn), where χri (2ω) contains the susceptibility modified by interband motion (evident by vbml and vcln). As shown in Figure 4a,b, these two terms have large but opposite contribution. However, their sum χi (2ω) ≡ χvi (2ω) + χri (2ω) remains very large, which accounts for the major part of total SHG susceptibility χ(2) yxx compared to other terms including interband χe(2ω), interband χe(ω), intraband χi(ω), and so forth. To further understand the origin of the extraordinary SHG susceptibility in MX monolayers, we examine the detailed distribution of the total susceptibility χ(2) yxx (Figure 5a,b) and the two dominant SHG susceptibility terms, χvi (2ω) (Figure 5c,d) and χri (2ω) (Figure 5e,f) in the first BZ. The total susceptibility χ(2) yxx is unambiguously concentrated around two valleys where the bandgap is located. One of the two major terms, χvi (2ω), is distributed at four spots in the first BZ close to but a bit deviated from the two valleys (Figure 5c,d), and they are largely governed by VBM and CBM near the Fermi level. The yxx component of χ iv (2ω) can be rearranged as follows ⎡ −ir y r x ⎤ x e3 dk χi,vyxx (2ω) = 2 ∫ 3 ∑nm 8fnm ⎢⎣ 2 nm mn ⎥⎦Δmn , where ℏ ωmn(ωmn − 2ω) 4π

{

}

5032

DOI: 10.1021/acs.nanolett.7b02268 Nano Lett. 2017, 17, 5027−5034

Letter

Nano Letters

Figure 6. Strong coupling among ferroelectric (P), ferroelastic (η) order parameters, and sheet SHG susceptibilities χ(2). (a) Four different multiferroic states in the GeSe monolayer: (−Px,ηxx, χyxx = 0, − χxyy ≠ 0), (Py,ηyy, χyxx ≠ 0, χxyy = 0), (−Py,ηyy, χyxx = 0, − χxyy ≠ 0), and (Px,ηxx, χyxx = 0, χxyy ≠ 0) marked by cyan, red, magenta, and green boxes, respectively. (b−e) Frequency-dependent sheet SHG susceptibilities χyxx (solid line) and χxyy (dashed line) of the above four different multiferroic states in the unit of 104 pm2/V.

described hereby. Moreover, as the ultimate nanometer thickness of 2D materials is much smaller than their coherent wavelength, these 2D nonlinear optical materials will not suffer from the crucial phase-matching problem generally encountered in 3D nonlinear optical crystals. The extraordinary nonlinear responses together with excellent phase-matching make 2D multiferroics particularly promising for nonlinear optoelectronic devices.53 Very excitingly, the unique ferroelasticity and ferroelectricity in these 2D multiferroics, as illustrated in Figure 6a, allow direct mechanical switching of lattice orientation and electrical switching of electrical polarization, accompanied by instantaneous direction and sign switching of the giant SHG susceptibilities. One can envisage, for example, to achieve active control of the ferroic orders in these 2D multiferroics with electrical, optical, and/or mechanical stimuli while in situ characterizing their local electronic structures and domain evolution using noncontact noninvasive ultrafast optical SHG techniques.

and SnSe and the results are shown in Figure S6 of Supporting Information, which demonstrates that monolayer and trilayer GeSe (SnSe) have very similar SHG strength and shape at low energy range, and they start to deviate more at relatively higher frequency range with additional oscillations and peaks, largely due to the weak van der Waals stacking of multiple layers. Group IV monochalcogenide monolayers are 2D ferroelasticferroelectric multiferroics with a few intrinsic order parameters, include the ferroelectric spontaneous polarization vector (P), the second-rank ferroelastic spontaneous strain tensor (η), and the SHG second-order nonlinear susceptibility (χ(2) abc (−2ω; ω, ω), a third-rank tensor). These order parameters are intrinsically coupled with each other. As illustrated in Figure 6(a), there are four possible configurations upon ferroelastic and/or ferroelectric transition. Each configuration has a peculiar set of order parameters that can be switched via ferroelastic/ ferroelectric phase transition. For example, the structure at the top of Figure 6a may experience ferroelastic transition to the one on the right with order parameters evolving from (Py, ηyy, χyxx ≠ 0, χxyy = 0) to (Px, ηxx, χyxx = 0, χxyy ≠ 0) where χxyy (right) = χyxx (top). Their SHG susceptibility tensor will change accordingly, displayed in Figure 6c,e. A subsequent ferroelectric transition from the right to the left will change it from (Px, ηxx, χyxx = 0, χxyy ≠ 0) to (−Px, ηxx, χyxx = 0, − χxyy ≠ 0), and its corresponding SHG susceptibility is shown in Figure 6b. The switching of all the above order parameters is simply due to its inherent symmetry including lattice, sites, and electric polarization. Very importantly, our results demonstrate that the SHG susceptibility in four ferroelectric−ferroelastic multiferroic states are distinctly different from each other, and hence the SHG intensity and phase information52 can be adopted as a unique noncontact measure for accessing multiferroic states. Although linear optical absorption may be used to determine the orientation of crystal lattice based on the anisotropic optical response,24 it cannot distinguish different ferroelectric polarization. The great sensitivity of SHG on the inversion symmetry breaking and the possibility of obtaining the SHG phase information are therefore highly advantageous. Conclusions. In summary, our present findings provide a theoretical prediction of extraordinary SHG nonlinear optical responses and a microscopic understanding in group IV monochalcogenide monolayers. More importantly, it will open up a variety of new opportunities for 2D ferroelectrics, multiferroics, and nonlinear optoelectronics. For example, it could be of fundamental interest to see the effect of photostriction31 on the first and higher-order optical responses

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.nanolett.7b02268. Calculation details, sheet SHG susceptibilities in all six monolayers, benchmark calculation of the Kramers− Kronig relation, SHG polarization anisotropy at different frequencies, the effect of spin−orbit coupling on SHG, and the results of monolayer and trilayer MX as well as detailed definition of sheet SHG susceptibility (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Xiaofeng Qian: 0000-0003-1627-288X Author Contributions

X.Q. conceived the idea and supervised the project. H.W. and X.Q. wrote the computational package for SHG nonlinear susceptibility tensors, performed the calculations, and analyzed the results. X.Q. wrote the initial draft of the manuscript with contributions from H.W. Notes

The authors declare no competing financial interest. 5033

DOI: 10.1021/acs.nanolett.7b02268 Nano Lett. 2017, 17, 5027−5034

Letter

Nano Letters

■

(28) Fei, R.; Li, W.; Li, J.; Yang, L. Appl. Phys. Lett. 2015, 107, 173104. (29) Chang, K.; Liu, J.; Lin, H.; Wang, N.; Zhao, K.; Zhang, A.; Jin, F.; Zhong, Y.; Hu, X.; Duan, W.; Zhang, Q.; Fu, L.; Xue, Q.-K.; Chen, X.; Ji, S.-H. Science 2016, 353, 274−278. (30) Shi, G.; Kioupakis, E. Nano Lett. 2015, 15, 6926−6931. (31) Haleoot, R.; Paillard, C.; Kaloni, T. P.; Mehboudi, M.; Xu, B.; Bellaiche, L.; Barraza-Lopez, S. Phys. Rev. Lett. 2017, 118, 227401. (32) Cook, A. M.; M. Fregoso, B.; de Juan, F.; Coh, S.; Moore, J. E. Nat. Commun. 2017, 8, 14176. (33) Rangel, T.; Fregoso, B. M.; Mendoza, B. S.; Morimoto, T.; Moore, J. E.; Neaton, J. B. 2016, arXiv preprint, arXiv:1610.06589 https://arxiv.org/abs/1610.06589, accessed October 20, 2016. (34) Kohn, W.; Sham, L. Phys. Rev. 1965, 140, A1133−A1138. (35) Kresse, G.; Furthmüller, J. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (36) Blöchl, P. E. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979. (37) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (38) Langreth, D. C.; Mehl, M. J. Phys. Rev. B: Condens. Matter Mater. Phys. 1983, 28, 1809−1834. (39) Baroni, S.; de Gironcoli, S.; Dal Corso, A.; Giannozzi, P. Rev. Mod. Phys. 2001, 73, 515−562. (40) Aversa, C.; Sipe, J. E. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 52, 14636−14645. (41) Hughes, J. L. P.; Sipe, J. E. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 53, 10751−10763. (42) Sharma, S.; Dewhurst, J. K.; Ambrosch-Draxl, C. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 67, 165332. (43) Gomes, L. C.; Trevisanutto, P. E.; Carvalho, A.; Rodin, A. S.; Castro Neto, A. H. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 94, 155428. (44) Hedin, L. Phys. Rev. 1965, 139, A796−A823. (45) Qian, X.; Li, J.; Qi, L.; Wang, C. Z.; Chan, T. L.; Yao, Y. X.; Ho, K. M.; Yip, S. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 245112. (46) Qian, X.; Li, J.; Yip, S. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 195442. (47) Guo, G. Y.; Lin, J. C. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 075416. (48) Grüning, M.; Attaccalite, C. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 081102. (49) Wang, C.-Y.; Guo, G.-Y. J. Phys. Chem. C 2015, 119, 13268− 13276. (50) Adams, E. N.; Blount, E. I. J. J. Phys. Chem. Solids 1959, 10, 286−303. (51) Sipe, J. E.; Ghahramani, E. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 48, 11705−11722. (52) Kemnitz, K.; Bhattacharyya, K.; Hicks, J. M.; Pinto, G. R.; Eisenthal, B.; Heinz, T. F. Chem. Phys. Lett. 1986, 131, 285−290. (53) Zhao, M.; Ye, Z.; Suzuki, R.; Ye, Y.; Zhu, H.; Xiao, J.; Wang, Y.; Iwasa, Y.; Zhang, X. Light: Sci. Appl. 2016, 5, e16131.

ACKNOWLEDGMENTS H.W. and X.Q. acknowledge the start-up funds from Texas A&M University. Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing.

■

REFERENCES

(1) Shen, Y.-R. The Principles of Nonlinear Optics; Wiley, 1984. (2) Butcher, P. N.; Cotter, D. The Elements of Nonlinear Optics; Cambridge University Press: 1991; Vol. 9. (3) Boyd, R. W. Nonlinear Optics, 3rd ed.; Elsevier Inc.: 2008. (4) Franken, P.; Hill, A. E.; Peters, C.; Weinreich, G. Phys. Rev. Lett. 1961, 7, 118. (5) Kumar, N.; Najmaei, S.; Cui, Q. N.; Ceballos, F.; Ajayan, P. M.; Lou, J.; Zhao, H. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 161403. (6) Li, Y. L.; Rao, Y.; Mak, K. F.; You, Y. M.; Wang, S. Y.; Dean, C. R.; Heinz, T. F. Nano Lett. 2013, 13, 3329−3333. (7) Malard, L. M.; Alencar, T. V.; Barboza, A. P. M.; Mak, K. F.; de Paula, A. M. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 201401. (8) Wang, R.; Chien, H. C.; Kumar, J.; Kumar, N.; Chiu, H. Y.; Zhao, H. ACS Appl. Mater. Interfaces 2014, 6, 314−318. (9) Jiang, T.; Liu, H.; Huang, D.; Zhang, S.; Li, Y.; Gong, X.; Shen, Y.-R.; Liu, W.-T.; Wu, S. Nat. Nanotechnol. 2014, 9, 825−829. (10) Gong, Z. R.; Liu, G. B.; Yu, H. Y.; Xiao, D.; Cui, X. D.; Xu, X. D.; Yao, W. Nat. Commun. 2013, 4, 2053. (11) Kim, D. H.; Lim, D. J. Korean Phys. Soc. 2015, 66, 816−820. (12) Wang, G.; Gerber, I. C.; Bouet, L.; Lagarde, D.; Balocchi, A.; Vidal, M.; Amand, T.; Marie, X.; Urbaszek, B. 2D Mater. 2015, 2, 045005. (13) Zeng, H. L.; Liu, G. B.; Dai, J. F.; Yan, Y. J.; Zhu, B. R.; He, R. C.; Xie, L.; Xu, S. J.; Chen, X. H.; Yao, W.; Cui, X. D. Sci. Rep. 2013, 3, 1608. (14) Janisch, C.; Wang, Y. X.; Ma, D.; Mehta, N.; Elias, A. L.; PereaLopez, N.; Terrones, M.; Crespi, V.; Liu, Z. W. Sci. Rep. 2015, 4, 5530. (15) Ye, Z. L.; Cao, T.; O’Brien, K.; Zhu, H. Y.; Yin, X. B.; Wang, Y.; Louie, S. G.; Zhang, X. Nature 2014, 513, 214−218. (16) Ribeiro-Soares, J.; Janisch, C.; Liu, Z.; Elias, A. L.; Dresselhaus, M. S.; Terrones, M.; Cancado, L. G.; Jorio, A. 2D Mater. 2015, 2, 045015. (17) Kim, C. J.; Brown, L.; Graham, M. W.; Hovden, R.; Havener, R. W.; McEuen, P. L.; Muller, D. A.; Park, J. Nano Lett. 2013, 13, 5660− 5665. (18) Karvonen, L.; Saynatjoki, A.; Mehravar, S.; Rodriguez, R. D.; Hartmann, S.; Zahn, D. R. T.; Honkanen, S.; Norwood, R. A.; Peyghambarian, N.; Kieu, K.; Lipsanen, H.; Riikonen, J. Sci. Rep. 2015, 5, 10334. (19) Deckoff-Jones, S.; Zhang, J.; Petoukhoff, C. E.; Man, M. K. L.; Lei, S.; Vajtai, R.; Ajayan, P. M.; Talbayev, D.; Madéo, J.; Dani, K. M. Sci. Rep. 2016, 6, 22620. (20) Yin, X. B.; Ye, Z. L.; Chenet, D. A.; Ye, Y.; O’Brien, K.; Hone, J. C.; Zhang, X. Science 2014, 344, 488−490. (21) Seyler, K. L.; Schaibley, J. R.; Gong, P.; Rivera, P.; Jones, A. M.; Wu, S. F.; Yan, J. Q.; Mandrus, D. G.; Yao, W.; Xu, X. D. Nat. Nanotechnol. 2015, 10, 407−411. (22) Wang, G.; Marie, X.; Gerber, I.; Amand, T.; Lagarde, D.; Bouet, L.; Vidal, M.; Balocchi, A.; Urbaszek, B. Phys. Rev. Lett. 2015, 114, 097403. (23) Wu, M.; Zeng, X. C. Nano Lett. 2016, 16, 3236−3241. (24) Wang, H.; Qian, X. 2D Mater. 2017, 4, 015042. (25) Mehboudi, M.; Dorio, A. M.; Zhu, W.; van der Zande, A.; Churchill, H. O.; Pacheco-Sanjuan, A. A.; Harriss, E. O.; Kumar, P.; Barraza-Lopez, S. Nano Lett. 2016, 16, 1704−1712. (26) Mehboudi, M.; Fregoso, B. M.; Yang, Y.; Zhu, W.; van der Zande, A.; Ferrer, J.; Bellaiche, L.; Kumar, P.; Barraza-Lopez, S. Phys. Rev. Lett. 2016, 117, 246802. (27) Fei, R.; Kang, W.; Yang, L. Phys. Rev. Lett. 2016, 117, 097601. 5034

DOI: 10.1021/acs.nanolett.7b02268 Nano Lett. 2017, 17, 5027−5034