Introduction

Ferroelectricity and bulk photovoltaic effect (BPVE)1,2 are two basic physical phenomena that emerge in condensed matter materials due to symmetry breaking. Generally speaking, BPVE, an intrinsic optical rectification phenomenon, appears in materials without inversion symmetry which can convert light to electricity as solar cells, hence it is naturally fulfilled the symmetry requirement in ferroelectric materials. As a result, the ferroelectric materials with switchable electronic polarization are one of the most studied BPVE materials1,2,3, and the corresponding energy conversion efficiency can exceed the Shockley-Queisser limit3,4,5.

It is known to us that the BPVE usually switches with the ferroelectric order. Under the inversion operation, both ferroelectric order and BPVE in all directions switch signs, i.e., BPVE switches synchronously with the ferroelectric order. Therefore, the ferroelectric polarization is often taken as a handle to manipulate the direction of BPVE photocurrent in, e.g., perovskite oxides BiFeO36,7,8, Pt/BiFeO3/SrRuO39, 2D (two-dimensional) CuInP2Se610, charge-transfer complex11, SnTe monolayer12, MX [M = (Ge, Sn) and X = (Se, S)]13 etc. However, from the aspect of symmetry, the ferroelectric polarization is a polar vector, while the BPVE coefficients form a rank-three tensor. They are hence subject to distinct symmetry transformation rules, and should not always synchronize. There should exist components of the BPVE tensor that do not switch with the ferroelectric order, which was visioned14, but has not been comprehensively studied before. For example, a mirror reflection perpendicular to the polarization will reverse the polarization just as the inversion operation does, while it cannot flip all the components of the BPVE. Thus, a precise understanding of the relationship between the ferroelectric order and the BPVE is of tremendous significance in both fundamental research and the design of photoelectric devices.

In this work, we think the exceptional case that non-synchronous BPVE with ferroelectric order can be found in 2D interlayer-sliding ferroelectric materials15. As demonstrated in recent experiments, the out-of-plane ferroelectricity can be achieved in bilayer WTe216,17,18,19, BN20,21, InSe22,23, MoS224,25, etc. by interlayer sliding, which effectively extends the family of 2D ferroelectric materials26,27,28,29. This ferroelectric flip mechanism may lead to the BPVE property beyond the switchable scenario.

Focusing on the four most common stacking cases of 2D bilayers with interlayer sliding, we use an abstract bilayer crystal model to locate the available configurations and symmetries with nontrivial ferroelectric order. We find two cases (Case 1b and Case 2a) can have the possibility to achieve the interlayer-sliding ferroelectricity, moreover the opposite ferroelectric states in the two cases are all related by the mirror operator. Coincidentally, the interlayer-sliding ferroelectric materials found in the experiments are attributed to these two cases. Using bilayer MoS2 and WTe2 as representative examples, we performed first-principles calculations and found that the in-plane BPVE does not change the sign with ferroelectric order, which is different from conventional ferroelectrics. In contrast to in-plane BPVE, out-of-plane BPVE switches with the reverse of the ferroelectric order. This non-synchronous character will lead to distinct photovoltaic phenomena in experiments.

Results

Polarization flipping via interlayer sliding

Now we introduce the mechanism of the polarization flipping via interlayer sliding. Figure 1a shows the out-of-plane polarization Pz in bilayer vdW (van der Waals) material without inversion nor horizontal mirror symmetry. Because the in-plane polarization is usually annihilated by in-plane rotational symmetry in real materials, the +Pz and –Pz polarization are denoted as +P and –P for short. The out-of-plane polarization can be flipped in two different ways: (i) intralayer movement, i.e., atomic relative displacements in each monolayer [Fig. 1b], and (ii) interlayer sliding15 where the crystal structure of each monolayer is invariant [Fig. 1c]. Although polarization flipping can be achieved in both situations, the underlying symmetries of the switch of ferroelectric states are different. For situation (i), the +P and –P ferroelectric states are correlated by the inversion operator, as indicated in unit cells in Fig. 1a, b. While for situation (ii), the +P and –P ferroelectric states are correlated by the mirror symmetry \(\hat M_{xy}\), as seen in unit cells in Fig. 1a, c (detailed analyses are presented in the Symmetry analysis section for real materials). The interlayer sliding mechanism is feasible experimentally and has been observed in bilayer vdW materials16,17,18,19,20,21,22,23,24,25 as stated above, which further extends the family of 2D ferroelectric materials26,27,28,29 effectively. Figure 1d–f shows the basis vector transformation under the inversion symmetry and mirror symmetry. The BPVEs in these interlayer-sliding ferroelectric materials do not always change synchronously with the ferroelectric order.

Fig. 1: Polarization flipping in a 2D vdW bilayer material via intralayer movement or interlayer sliding, and related symmetries.
figure 1

a A 2D system with positive out-of-plane polarization +P. b The polarization is reversed to –P by atom displacements in each monolayer, and the black arrow means the movement of the cations relative to the anions. c The polarization is reversed by interlayer sliding and the black arrows mean the interlayer sliding direction. d Three basis vectors and their transformations under e inversion symmetry f horizontal mirror symmetry. The black rectangles denote the bilayer unit cells in (a)-(c), and the gray rectangles denote the monolayer unit cells.

Full size image

Symmetry analysis

Next, we analyze how to design the bilayer vdW materials with the out-of-plane polarization via interlayer sliding in Fig. 1a, and analyze how positive and negative polarizations are related by symmetry operator. Most of the discovered monolayer materials are non-polar and non-ferroelectric (such as BN, transition-metal dichalcogenides, metal trihalide, post-transition-metal chalcogenides), and their monolayer crystal structures belong to the following two different cases, as shown in Supplementary Figs. 1, 2 in the Supplementary Information30,31:

  1. 1.

    Monolayer has the inversion symmetry \(\widehat I\), but has no horizontal mirror symmetry \(\widehat M_{xy}\), such as monolayer 1T-PtS2, 1T’-WTe2, BiI3, etc.;

  2. 2.

    Monolayer has no inversion symmetry \(\widehat I\) but has horizontal mirror symmetry \(\widehat M_{xy}\), such as monolayer BN, 1H-MoS2, 1H-WSe2, GaSe, InSe, etc.

Besides, the bilayer can be further divided into two subcases according to the stacking manner of two layers:

  1. (a)

    A/A stacking, i.e., two adjacent layers are identical where the A represents the crystal structure of the monolayer.

  2. (b)

    A/B stacking, i.e., two adjacent layers are relatively rotated by 180° with \({{{\mathrm{B}}}} = \hat C_{2z}{{{\mathrm{A}}}}\) (A and B mean the crystal types).

The above considerations can form four cases for 2D bilayer materials with interlayer sliding, and we can analyze their symmetries with an abstract bilayer crystal model which only considers the above crystal symmetries and ignores the specific crystal structures. The complete symmetry analysis for each case is summarized in Supplementary Note 1, and the corresponding results are shown in Table 1. As can be seen, Case 1b and Case 2a host the out-of-plane polarization, and the two opposite interlayer-sliding states are correlated by \(\widehat M_{xy}\) operator, instead of inversion \(\widehat I\) operator. In contrast, the crystal structures in Case 1a and Case 2b possess the inversion symmetry whatever the interlayer-sliding vector, which forbids the occurrence of ferroelectricity and BPVE.

Table 1 Symmetries of interlayer-sliding bilayer materials.
Full size table

With the abstract bilayer crystal model, the interlayer-sliding ferroelectric materials are further classified into two kinds (Case 1b and Cases 2a). Such classification outlines a potential roadmap to search for more 2D ferroelectric materials, which could greatly facilitate and advance the research and application of BPVE in 2D ferroelectric materials. Table 2 summarizes several bilayer candidates for the above four cases that can be experimentally realized by either simply mechanical exfoliation or tear-and-stack methods. In particular, the out-of-plane ferroelectricity of the bilayer WTe216,17,18,19 with Case 1b and BN20,21, InSe22,23, bilayer-TMD24,25 with Case 2a have been discovered in experiments recently. The abstract bilayer model is convenient and efficient in symmetry analysis for vdW materials with interlayer sliding because it does not rely on the specific crystal structure. The symmetry of trilayer vdW materials with interlayer sliding can also be obtained using a similar method, and the results are shown in Supplementary Note 7, where the transformations of ferroelectricity and BPVE are also discussed there.

Table 2 Candidate bilayer materials and corresponding symmetries of bulk parents.
Full size table

Now we analyze the characters of in-plane and out-of-plane BPVE of the interlayer-sliding ferroelectric materials in Case 1b and Case 2a under the normal incidence of light. The in-plane BPVE current density is

$$J_a = \mathop {\sum}\limits_{bc} {\sigma _{bc}^aE_b\left( \omega \right)E_c\left( { - \omega } \right)} ,$$
(1)

where \(a,b,c \in \left\{ {x,y} \right\}\), Eb and Ec are electric fields of light along b and c direction, ω is the frequency of light, and \(\sigma _{bc}^a\) is the BPVE coefficient. In contrast to x and y directions, the photo-excited shift of the wave package along the out-of-plane direction induces a static electric polarization instead of an electric current along the z-direction because the z-direction of 2D materials is discontinuous. The out-of-plane BPVE polarization is32

$$p_z = \mathop {\sum}\limits_{bc} {\sigma _{bc}^{\prime z}E_b\left( \omega \right)E_c\left( { - \omega } \right)} ,$$
(2)

where \(p_z\) is z-component polarizability.

Similar to other nonlinear optics tensors, the BPVE tensor \(\left[ {\sigma _{ab}^c} \right]_{3 \times 3 \times 3}\) in Eq. (1) [or Eq. (2)] obeys

$$\sigma _{jk}^i = \mathop {\sum}\limits_{abc} {R_{ia}R_{jb}R_{kc}\sigma _{bc}^a} ,$$
(3)

where \([R_{ia}]\) is the 3 × 3 matrix of the symmetry operator \(\widehat R\). As shown in Fig. 1f, \(\widehat M_{xy}\) will result in the reverse of the out-of-plane vector with invariant in-plane vectors, and the corresponding operator matrix is

$$\left[ {M_{xy}} \right] = \left[ {\begin{array}{*{20}{c}} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & { - 1} \end{array}} \right].$$
(4)

The opposite states are connected by the mirror operator [\(\widehat M_{xy}( + P) = - P\)] in interlayer-sliding ferroelectric materials. Therefore, the in-plane BPVE coefficient is invariant with the change of ferroelectric orders, i.e., \(\sigma _{ab}^c( + P) = \sigma _{ab}^c( - P),a,b,c \in \left\{ {x,y} \right\}\), according to Eq. (3). This is strikingly different from the out-of-plane case where the BPVE coefficient is reversed with the switch of ferroelectric orders, as indicated by \(\sigma _{bc}^{\prime z}\left( { + P} \right) = - \sigma _{bc}^{\prime z}\left( { - P} \right)\) (\({\kern 1pt} b,c \in \left\{ {x,y} \right\}\)). The calculated BPVE characters in Fig. 1a, c by a 1D effective model are also consistent with our symmetry analysis (see Supplementary Note 2).

In contrast, for the ferroelectric states correlated by inversion symmetry \(\hat I( + P) = - P\), the inversion operator \(\hat I\) will lead to the reverse of all the vectors (Fig. 1d), i.e., \(\sigma _{bc}^a( + P) = - \sigma _{bc}^a( - P)\) (\(a,b,c \in \left\{ {x,y,z} \right\}\)) according to Eq. (3). It is the reason why the BPVE coefficients reverse with ferroelectric orders in the conventional ferroelectric materials.

The in-plane BPVE coefficients can be calculated by the standard second-order Kubo formalism theory33,34, and the out-of-plane BPVE coefficients can be obtained by a modified second-order Kubo formalism theory32 that was proposed recently. BPVE coefficients are calculated by the Wannier function (see Method for details)35,36,37. Recently, theory works indicate that38,39 BPVE coefficients are highly related to the geometry theory of electron band. In the following part, we choose a representative material for each case: bilayer MoS224,25 (Case 2a) and bilayer WTe216,17,18,19 (Case 1b) to perform the numerical first-principles calculations (see Methods section for calculation details).

First-principles calculation results

Bilayer MoS2

Bilayer MoS2 (Case 2a) that is consisted of two identical non-ferroelectric monolayers shows out-of-plane ferroelectricity15, which has been confirmed in experiments recently24,25. Noting that conventional bilayer MoS2 is stacked in A/B form (Case 2b) with inversion symmetry, which is different from the crystal structure discussed here. As shown in Fig. 2a, the out-of-plane ferroelectricity is reversed by the in-plane interlayer motion along the armchair direction. The monolayer possesses the \(D_{3h}\) symmetry. However, the interlayer-sliding bilayer breaks the \(\widehat M_{xy}\) and in-plane \(\widehat C_2\) rotation symmetries, which results in the bilayer MoS2 only having \(C_{3v}\) symmetry (\(D_{3h} = C_{3v} + \widehat M_{xy}C_{3v}\)). Therefore, +P and –P states only have out-of-plane polarization. Moreover, the two ferroelectric states are connected by the \(\widehat M_{xy}\) symmetry24,25, which is consistent with the above symmetry analysis. The calculated band structures of the +P and –P states are shown in Fig. 2b, which are invariant with the ferroelectric orders.

Fig. 2: Calculated results of ferroelectric bilayer MoS2.
figure 2

a Crystal structures of ferroelectric states (upper panel: side view, lower panel: top view). The arrows denote the transformation of the +P and –P states by interlayer \(\widehat M_{xy}\) symmetry. The definitions of +P and –P states are adopted from ref. 24. b Band structure of +P and –P states. c In-plane and e out-of-plane BPVE coefficients. d In-plane BPVE photocurrent and f out-of-plane BPVE polarization (in arbitrary unit) with the direction of normal-incidence polarization of light in the +P states. θ in d and f is the azimuthal angle of the polarization of light relative to the x axis, and the dash lines mean the vertical mirror planes.

Full size image

As shown in Fig. 2c, the calculated in-plane BPVE coefficients with +P and –P states are equal. The in-plane BPVE is invariant with the change of the ferroelectric order as predicted. In addition, \(\sigma _{xy}^x = \sigma _{xx}^y = - \sigma _{yy}^y\), and \(\sigma _{xx}^x = \sigma _{yy}^x = \sigma _{xy}^y = 0\) for each ferroelectric state due to the C3v symmetry [see Supplementary Note 3]. It is noted that this unswtichable in-plane BPVE character of interlayer-sliding ferroelectric materials is different from those of conventional ferroelectric materials6,7,8,9,10,11,12,40 and antiferromagnets41,42 where the BPVE is switched with the ferroelectric/ferromagnetic order. The in-plane BPVE photocurrent as a function of the polarization direction of light (see Supplementary Note 4 for details) is shown in Fig. 2d. The in-plane BPVE photocurrent shows the C3v symmetry, and its magnitude is independent of the direction of light.

The calculated out-of-plane BPVE coefficients are shown in Fig. 2e. For each ferroelectric state, \(\sigma _{xx}^{\prime z} = \sigma _{yy}^{\prime z}\), and \(\sigma _{xy}^{\prime z} = \sigma _{yx}^{\prime z} = 0\). Moreover, the out-of-plane BPVE coefficients switch with the ferroelectric order following \(\sigma _{xx}^{\prime z}( + P) = \sigma _{yy}^{\prime z}( + P) = - \sigma _{xx}^{\prime z}( - P) = - \sigma _{yy}^{\prime z}( - P)\), which is consistent with our symmetry analysis. The out-of-plane BPVE is isotropic, and photo-induced polarization is independent of the polarization direction of light due to C3v symmetry, as shown in Fig. 2f.

Similarly, experimentally prepared bilayer BN20,21, InSe22,23, GaSe, etc., possess similar interlayer-sliding ferroelectricity (Case 2a) and the same symmetry. Therefore, similar in-plane and out-of-plane BPVE features are also expected to exhibit in these materials.

Bilayer WTe2

Bilayer WTe2 is composed of two 180°-rotated two monolayers (Case 1b), as shown in Fig. 3a. Even though each monolayer (with \(C_{2h} = \{ E,\widehat I,\widehat C_{2x},\widehat M_{yz}\}\) symmetry) has the inversion \(\widehat I\) symmetry, the bilayer breaks this symmetry. Besides, the interlayer-sliding vector along b axis further breaks \(\widehat M_{xy}\) and in-plane \(\widehat C_2\) symmetry. Thereby the bilayer has the \(C_{1v} = \{ E,\widehat M_{yz}\}\) symmetry (\(C_{2h} = C_{1v} + \widehat M_{xy}C_{1v}\)).

Fig. 3: Calculated results of ferroelectric bilayer WTe2.
figure 3

a Crystal structure. +P and –P states are adopted from ref. 43. Rectangle in a represents the unit cell, and the double-arrows line means the monolayer transformation under \(\widehat M_{xy}\) symmetry. The lower panel in a denotes the top view of bilayer WTe2 with +P state. b Band structures of +P and –P states. c In-plane and e out-of-plane BPVE coefficients. d In-plane BPVE photocurrent and f out-of-plane BPVE polarization of the +P states with the different directions of polarization of light. The length of arrows in d indicates the magnitude of in-plane BPVE response (in arbitrary unit), where \(\sigma _{xy}^x:\sigma _{xx}^y:\sigma _{yy}^y = 1:0.7: - 1.6\). The width and the color of the ribbon in f denote the magnitude and direction of out-of-plane BPVE response, and \(\sigma {\prime}_{xx}^z:\sigma {\prime}_{yy}^z = 1: - 0.6\).

Full size image

Interlayer sliding43,44 induced out-of-plane ferroelectricity in bilayer WTe2 has been verified experimentally16,17,18,19. The –P state is the mirror image of +P state, as shown in Fig. 3a. Figure 3b shows the band structures of +P and –P states, where the valence and conduction bands overlap in our calculation.

The calculated in-plane BPVE coefficients are shown in Fig. 3c. There are three independent tensor elements \(\sigma _{xy}^x\), \(\sigma _{xx}^y\), \(\sigma _{yy}^y\), and \(\sigma _{xx}^x = \sigma _{yy}^x = \sigma _{xy}^y = 0\) due to the C1v symmetry. The calculated in-plane BPVE coefficients are the same in two opposite ferroelectric orders: \(\sigma _{xy}^x( + P) = \sigma _{xy}^x( - P)\), \(\sigma _{xx}^y( + P) = \sigma _{xx}^y( - P)\), \(\sigma _{yy}^y( + P) = \sigma _{yy}^y( - P)\), which is consistent with the above symmetry analysis. The in-plane BPVE photocurrent with the direction of polarization of light is shown in Fig. 3d, which shows \(\widehat M_{yz}\) symmetry, and its magnitude is dependent on the light (see Supplementary Note 4 for details). Recently, it was shown that bilayer WTe2 exhibits unswitchable in-plane nonlinear anomalous Hall effect14,18.

On the contrary, the out-of-plane BPVE coefficients switch with the change of ferroelectric order, \(\sigma _{xx}^{\prime z}( + P) = - \sigma _{xx}^{\prime z}( - P)\), \(\sigma _{yy}^{\prime z}( + P) = - \sigma _{yy}^{\prime z}( - P)\) [Fig. 3e]. Unlike the bilayer MoS2, the out-of-plane BPVE coefficients show anisotropy \(\sigma _{xx}^{\prime z} \,\ne\, \sigma _{yy}^{\prime z}\) (\(\sigma _{xy}^{\prime z} = \sigma _{yx}^{\prime z} = 0\)) for each ferroelectric state due to the absence of in-plane C3 symmetry, leading to the out-of-plane BPVE polarization is dependent on the polarization of light as shown in Fig. 3f. Recently, ZrI2, a sister material of polar semimetals WTe2, is also demonstrated45,46,47 to have interlayer-sliding ferroelectricity. Similar BPVE behaviors with bilayer WTe2 are expected to be observed due to the same symmetry and interlayer-stacking way (Case 1b).

The calculated values of the in-plane BVPE coefficient of bilayer MoS2 and bilayer WTe2 are comparable to the reported monolayer MX [M = (Ge, Sn) and X = (Se, S)]13,48 and larger than 3D conventional ferroelectric materials (such as BiFeO349, BiTiO3, and PbTiO32). The out-of-plane BPVE coefficients are comparable to twisted bilayers graphene32, which induces voltage are detectable experimentally. We can distinguish the photo-induced polarization and intrinsic polarization of interlayer-sliding ferroelectric materials in two subsequent detections with/without light in experiments.

We can calculate the light-induced out-of-plane electric dipole moment and voltage caused by the out-of-plane BPVE, and estimate the possibility to flip the ferroelectric order. With light power of 2.65 × 1011 W/cm2 (equivalent electric field E=1 V/nm) and a typical value of \(\sigma _{xx}^{\prime z} = \sigma _{yy}^{\prime z}\)=1 e/V2 for the out-of-plane BPVE coefficients, the light-induced polarization pz = 1 e/nm2, meanwhile the light-induced voltage difference between the two layers (\(V = p_zd/\varepsilon _0\), d = 1 nm) is found to be 1.8 V. The photo-induced polarization is already bigger than intrinsic polarization of interlayer-sliding ferroelectric materials (about 0.01 e/nm2)15, and the photo-induced out-of-plane voltage is also bigger than the flipping voltage (about 0.1 V/nm) in experiments16,17,18,19,20,21,22,23,24,25. Therefore, we think the photo-induced polarization can flip the ferroelectric order when the laser intensity is larger enough. Recently, experiment50 shows that the stacking order in bulk WTe2 can be manipulated by ultrafast optical excitation, which indicates the possibility for the flip the ferroelectric order by optical irradiation. The kinetic process of ferroelectric phase transition under light is beyond the scope of current work, which needs to be further studied.

In our work, we only studied BPVE caused by the electronic shift current mechanism. While recently studies51,52,53 showed that excitons can enhance BPVE even above the band gap. The exciton effect on the BPVE in interlayer-sliding ferroelectric materials is deserved further study, however, the BPVE behavior with the ferroelectric order should not be changed because it is constrained by symmetry.

Discussion

According to the above two specific examples, the in-plane BPVE is invariant, while the out-of-plane BPVE reverses with the change of ferroelectric order in interlayer-sliding ferroelectric materials. The ferroelectric order and BPVE in bilayer interlayer sliding ferroelectric materials are subject to different symmetry transformation rules, and do not always synchronize. The study of non-synchronous BPVE with the ferroelectric order in mirror symmetry-related ferroelectric materials allows the comprehension of the relationship between them in another view.

Moreover, these BPVE characters in 2D interlayer-siding ferroelectric materials will lead to distinctive photovoltaic phenomena and applications in photoelectric devices. The in-plane BPVE photocurrents in a single crystal of interlayer-siding ferroelectric materials are invariant with two opposite ferroelectric domains, thus the photovoltaic currents between different domains superimpose rather than cancel, as shown in Fig. 4. Since ferroelectric domains are often unavoidable in real materials, especially for the interlayer-sliding ferroelectric materials fabricated via the tear-and-stack method20,21,25,54. Therefore, this unswitchable BPVE feature boosts the development of large-scale photoelectric devices that are immune to the polarization of ferroelectric domains. Nevertheless, the out-of-plane BPVE inverses between the opposite ferroelectric orders (Fig. 4), which can be utilized to construct switchable photoelectric devices.

Fig. 4: In-plane and out-of-plane BPVE response with opposite ferroelectric domains for interlayer-sliding ferroelectric materials.
figure 4

The horizontal and vertical arrows denote the in-plane BPVE current and out-of-plane BPVE polarization, respectively.

Full size image

In conclusion, we use the abstract bilayer crystal model to analyze the symmetries of the four most common cases of bilayer vdW materials with interlayer sliding and found two cases (Case 1b and Case 2a) can have the possibility to achieve the interlayer-sliding ferroelectricity, which the interlayer-sliding ferroelectric materials found in experiments fall into. We revealed that the opposite ferroelectric states in the interlayer-sliding ferroelectric materials are linked by the horizontal mirror symmetry, leading to the in-plane BPVE being invariant while the out-of-plane BPVE switches with opposite ferroelectric states. BPVE in all directions of interlayer-sliding ferroelectric materials are not flipped synchronously with the ferroelectric orders, which is different from the scenario occurs in traditional ferroelectric materials. Our theoretical study provides a strategy for a comprehensive understanding of the relationship between the ferroelectric order and BPVE in 2D interlayer-sliding ferroelectric materials. Moreover, the switchable out-of-plane BPVE can be used to design switchable photoelectric devices. On the other hand, the in-plane BPVE is robust against the ferroelectric order, and the unswitchable character is beneficial to construct larger-scale photoelectric devices.

Methods

In-plane and out-of-plane BPVE theory

According to the second-order Kubo formalism theory33,34, \(\sigma _{bc}^a\) in the Eq. (1) of the main text can be expressed as

$$\sigma _{bc}^a \equiv \frac{1}{2}{\rm{Re}} \left\{ {\chi _{bc}^a + \chi _{cb}^a} \right\},$$
(5)

where

$$\chi _{bc}^a(\omega ) = - \frac{{e^3}}{{(\hbar\omega)^2}}{\int} {\frac{{d{{{\bf{k}}}}}}{{(2\pi )^3}}} \mathop {\sum}\limits_{mnl} {\frac{{f_{lm}v_{lm}^b}}{{E_{ml} - \hbar \omega + {{{\mathrm{i}}}}\delta }}} \left( {\frac{{v_{mn}^av_{nl}^c}}{{E_{mn} + {{{\mathrm{i}}}}\delta }} - \frac{{v_{mn}^cv_{nl}^a}}{{E_{nl} + {{{\mathrm{i}}}}\delta }}} \right),$$
(6)

\(\delta = \hbar /\tau\), and τ is the lifetime, ω is the frequency of the light, \(f_{lm}\) and \(E_{ml}\) are the difference of occupation number and band energy between bands l and m.

Within a modified second-order Kubo formalism theory32, the out-of-plane BPVE coefficient in Eq. (2) of the main text can be expressed as

$$\sigma _{bc}^{\prime z} \equiv \frac{1}{2}{\rm{Re}} \left\{ {\chi _{bc}^{\prime z} + \chi _{cb}^{\prime z}} \right\},$$
(7)

and

$$\chi _{bc}^{\prime z} = \frac{{e^2}}{{(\hbar\omega)^2}}{\int} {\frac{{d{{{\mathbf{k}}}}}}{{(2\pi )^3}}} \mathop {\sum}\limits_{mnl} {\frac{{f_{lm}v_{lm}^b}}{{E_{ml} - \hbar \omega + {{{\mathrm{i}}}}\delta }}} \left( {\frac{{(p_z)_{nm}v_{nl}^c}}{{E_{mn} + {{{\mathrm{i}}}}\delta }} - \frac{{v_{mn}^c(p_z)_{nl}}}{{E_{nl} + {{{\mathrm{i}}}}\delta }}} \right),$$
(8)

where \((p_z)_{nm} = \left\langle n \right|\hat p_z\left| m \right\rangle\) is the out-of-plane dipole matrix element, and \(\hat p_z\) is the out-of-plane dipole operator. The interatomic electric dipoles due to the overlapping of wave functions should be very small, which was ignored in our calculation. \(\hat p_z\) contributed by interatomic electric dipole only has nonzero diagonal components which can be expressed as

$$\hat p_z = - e\left[ {\begin{array}{*{20}{c}} {r_{1,z}} & 0 & \cdots \\ 0 & {r_{2,z}} & \cdots \\ \vdots & \vdots & \ddots \end{array}} \right]$$
(9)

where \(r_{i,z}\) is the z-component position of the i-th atom. We define the middle of the 2D materials as the origin: \(\mathop {\sum}\limits_i {r_{i,z}} = 0\). Equation (6) is very similar to Eq. (8), except for the position operator \(\hat p_{i,z}\) instead of the velocity operator.

First-principles calculations

The first-principles calculations based on density functional theory (DFT) are performed by using the VASP package. General gradient approximation (GGA) according to the Perdew–Burke–Ernzerhof (PBE) functional is used. The energy cutoff of the plane wave basis is set to 400 eV. The Brillouin zone is sampled with a 12 × 12 × 1 (12 × 8 × 1) mesh of k-points for MoS2 (WTe2). To simulate the monolayers, vacuum layers (~15 Å) are introduced. The vdW force with DFT-D2 correction is considered in the main text. The calculation results and corresponding discussions using vdW force with DFT-D3 correction are shown in Supplementary Note 6. Spin-orbital coupling effects are considered for bilayer MoS2 and WTe2 in the band structure and BPVE calculations.

The DFT Bloch wave functions are iteratively transformed into maximally localized Wannier functions by the Wannier90 code55,56. Mo-d and S-p (W-d and Te-p) orbitals are used to construct the Wannier functions for MoS2 (WTe2). The in-plane and out-of-plane BPVE coefficients are calculated by our own program WNLOP (Wannier Nonlinear Optics Package) based on effective tight-binding (TB) Hamiltonian. A convergence test of k-mesh is performed, and 500 × 500 × 1 (600 × 300 × 1) k-mesh is sufficient in BPVE calculations of MoS2 (WTe2). We adopt δ = 0.02 eV for Eq. (6) and Eq. (8) in our calculation to take into account various relaxation processes.

The 3D-like BPVE coefficients are obtained assuming an active single-layer with a thickness of \(L_{{{{\mathrm{active}}}}}\)13

$$\sigma _{3D} = \frac{{L_{{{{\mathrm{slab}}}}}}}{{L_{{{{\mathrm{active}}}}}}}\sigma _{{{{\mathrm{slab}}}}},$$
(10)

where \(\sigma _{{{{\mathrm{slab}}}}}\) is the calculated BPVE coefficient, and \(L_{{{{\mathrm{slab}}}}}\) (\(L_{{{{\mathrm{active}}}}} \,<\, L_{{{{\mathrm{slab}}}}}\)) is the effective thickness slab.