Strain-tunable van der Waals interactions in few-layer black phosphorus

Interlayer interactions in 2D materials, also known as van der Waals (vdWs) interactions, play a critical role in the physical properties of layered materials. It is fascinating to manipulate the vdWs interaction, and hence to “redefine” the material properties. Here, we demonstrate that in-plane biaxial strain can effectively tune the vdWs interaction of few-layer black phosphorus with thickness of 2-10 layers, using infrared spectroscopy. Surprisingly, our results reveal that in-plane tensile strain efficiently weakens the interlayer coupling, even though the sample shrinks in the vertical direction due to the Poisson effect, in sharp contrast to one’s intuition. Moreover, density functional theory (DFT) calculations further confirm our observations and indicate a dominant role of the puckered lattice structure. Our study highlights the important role played by vdWs interactions in 2D materials during external physical perturbations.

R ecently, layered materials such as graphene, transition metal dichalcogenides (TMDCs), and black phosphorus (BP) have attracted increasing attention. The layers in these materials are held together by van der Waals (vdWs) interactions, which play a critical role in the material property, and hence device performance. One of the most fascinating aspects is that different two-dimensional (2D) layers can stick together through weak vdWs interactions without the constraint of lattice match, holding great promise for building hetero-structure-type materials with entirely new properties 1,2 . A variety of novel quantum phenomena have been reported in vdWs heterostructures, such as fractional quantum Hall effect 3,4 , gap opening in gapless graphene 5 , unconventional superconductivity 6 , tunable Mott insulator 7,8 , long life time Moire excitons 9 , and so on. Therefore, it provides us a new avenue to tailor the properties of layered materials by manipulating vdWs interactions. The strength of vdWs interaction is closely related to the interlayer distance, which can be shrunk or expanded by external perturbations, such as hydrostatic pressure 10,11 and ion intercalation 12,13 . Only recently, researches showed that in-plane strain can couple with interlayer interactions in bilayer MoS 2 14 and MoS 2 /WS 2 heterostructure 15 through Raman and photoluminescence spectroscopies.
BP is an emerging layered material with puckered lattice structure and thickness-dependent direct bandgap [16][17][18][19][20] . The vdWs interaction in few-layer BP splits a band into multiple subbands and causes a series of optical resonances in the infrared (IR) spectra 17,19 , which provides a unique opportunity to monitor the interlayer coupling strength. Experiments have shown that few-layer BP is very sensitive to in-plane strain, with a large strain tunability in the band gap and higher energy optical transitions 17,21,22 . However, there is no study on the strain effect on the interlayer coupling up to date. Typically, in-plane strain tunes the in-plane bond strength and how it affects the interlayer coupling is not clear. In this paper, we systematically investigate the evolution of band structures of few-layer BP under in-plane biaxial strain, using IR absorption spectroscopy. We separate the interlayer and intralayer contributions to the evolution of the bandstructure under in-plane strain, and gain quantitative insights into the van der Waals coupling effect in the strain engineering process. Surprisingly, we found that in-plane straining is a very efficient scheme to tune the vdWs interaction in BP, which opens a promising avenue for strain engineering of vdWs interactions.

Results
Sample characterization and the scheme for biaxial strain. We studied exfoliated BP samples on polypropylene (PP) substrate (see the section "Methods" for sample preparation). The IR extinction (1 − T/T 0 ) spectrum of few-layer BP was obtained using Fourier transform IR (FTIR) spectrometer, where T/T 0 is the ratio of light transmission with/without BP flakes (see details in the section "Methods") 17 . From the polarized IR spectrum, we can directly identify the layer number and crystal orientation 17 . Figure 1b is a typical IR extinction spectrum of a 6L BP. Two   Fig. 1 Band structure engineering by biaxial strain. a Schematic illustration of the experiment setup used for applying in-plane biaxial strain by heating or cooling the PP substrate. b A typical IR extinction spectrum (1 − T/T 0 ) for a 6L BP on PP substrate under zero strain, with the incident light polarized along the AC direction. Top inset is a schematic illustration for optical transitions between different subbands. Bottom inset is an optical image of this 6L BP sample. Scale bar: 20 μm. c IR extinction spectra (1 − T/T 0 ) for the 6L BP under varying biaxial tensile (>0) and compressive (<0) strains. For clarity, the spectra are vertically offset. Dashed lines are guides to the eye. d The E 11 and E 22 peak energies as a function of biaxial strain. The solid lines are linear fits to the data, which give shift rates for E 11  Transition index-dependent strain effect. Figure 1c shows the typical IR extinction spectra of a 6L BP sample under various biaxial tensile and compressive strains. Both of the characteristic peaks (E 11 and E 22 ) exhibit blueshift with tensile strain and redshift with compressive strain. Interestingly, the E 11 peak shifts faster than the E 22 peak, with shift rates of 222 and 167 meV/%, respectively, as extracted from linear fittings in Fig. 1d. The inset of Fig. 1d shows that the energy difference between E 11 and E 22 decreases as tensile strain is applied. Although optical resonances are dominated by excitons in low-dimensional semiconductors 19,20 , the shift is mainly attributed to the change of the electronic band structure and the modification of exciton binding energy due to such small strain is negelected 23 . In addition, the observed shift of resonances might be a mixed result of strain and pure temperature effect. To account for this possibility, controlled experiments were carried out. We transferred few-layer BP samples to quartz substrates, which have a much smaller thermal expansion coefficient than PP. Therefore, the strain effect of quartz can be ignored. As shown in Supplementary  Fig. 3, the pure temperature effect in the range of 25-75°C is almost negligible compared to the strain effect. Thus, in our analysis, we simply neglect the pure temperature-induced shift. The fact that E 11 shifts faster than E 22 is not limited to biaxially strained 6L and can be observed in all strained few-layer BP samples, even uniaxially strained samples (see more data in Supplementary Figs. 4 and 5). Generally, we find that optical resonances originated from higher-index subbands exhibit smaller shift rate under strain. This transition index-dependent strain effect indicates that interlayer vdWs interaction changes with in-plane strain, since the splitting of subbands directly reflects the strength of interlayer coupling. Quantitative description will be provided in the section "Discussion".
Layer-dependent strain effect. For E 11 transitions of samples with different thickness, the shift rates are different as well. We can obtain BP samples with different thickness in one big flake, which are presumably under the same strain condition all the time. An example is shown in Fig. 2a, where 3L and 4L BP are adjacent to each other. To confirm this, Raman measurements were performed. According to previous studies 24 , it is reasonable to assume that the strain effect on some of the Raman peaks has little sample thickness dependence. Therefore, if the shift rates of Raman peaks under strain are the same, we can guarantee the same strain magnitude for different BP layers all the time. Given that there is Davydov splitting of A 2 g mode 24,25 , and the intensity of B 2g mode is relatively weak under excitation with light polarized along armchair direction, A 1 g peak shift serves as an indicator for the strain magnitude.  Indeed, the shift rates of the A 1 g mode are almost the same, with values of −5.17 and −4.92 cm −1 /% for 3L and 4L, respectively (see Supplementary Fig. 8), which assures the same strain condition. Figure 2b shows the IR extinction spectra of this sample under different strains, with peak positions as a function of biaxial strain summarized in Fig. 2c. E 11 of 3L shifts more slowly than that of 4L, with shift rates of 158 and 185 meV/%, respectively. This result shows that the shift rate of E 11 transition also depends on the sample thickness. Indeed, this dependence is systematic, as evidenced by a large number of measured samples with thickness ranging from 2L to 10L. The averaged shift rates for each layer number are plotted in Fig. 3a, with the error bar indicating the shift rate range of multiple measured samples (for each layer number, at least three samples were measured to ensure the reproducibility of the data). To make sure the strain was accurately calibrated for different thickness samples, the Raman A 1 g peak shift induced by strain was served as a monitor. As we can see in Fig. 3a, the shift rates of E nn transitions are strongly subband index and layer number dependent, particularly for thin samples (layer number below 6L). Our measurements reveal that thinner samples and transitions associated with higher-index subbands are less sensitive to strain (smaller shift rates), which implies that inplane strain can tune the vdWs interaction, as detailed in the section "Discussion".
To further confirm this, we performed uniaxial strain experiments for BP with different layer number, with the strain direction along AC and ZZ, respectively. Results are presented in Supplementary Figs. 6 and 7. Notable layer number dependence of the strain effect is observed both for AC and ZZ strains, consistent with biaxial strain results. Note that the shift rates induced by uniaxial strain are almost half of those induced by biaxial strain, which is as expected, since the sample deforms in two directions under biaxial strain while only one direction under uniaxial strain. According to Fig. 3a, we can infer that the shift rate of bulk samples under biaxial strain is a little larger than 200 meV. This suggests a shift rate of~100 meV/% under uniaxial strain, which is in good agreement with previous results (99 meV/% for strain along AC and 109 meV/% for strain along ZZ in ref. 21 ).

Discussion
For unstrained few-layer BP, a 1D tight-binding model accounting for the interlayer coupling can well describe the optical transitions, with only two relevant parameters: the monolayer band gap E g and the interlayer coupling strength 17,19 . It is apparent that E g depends on the in-plane strain, since it is merely determined by the intralayer bond. The dependence of interlayer vdW interaction on the in-plane strain is not clear yet. An intuition may imply that the interlayer coupling would be enhanced with in-plane tensile strain, since the interlayer distance typically decreases due to the Poisson effect. However, our experiment gives opposite result, as detailed in the following analysis. From 1D tight-binding model (see details   Fig. 3 Layer and subband index dependence of the biaxial strain effect. a Averaged shift rates of E 11 , E 22 , and E 33 peaks as a function of layer number in 2-10L BP. The solid curves are fitted to the data using the tight-binding model shown in the text. The error bar is defined from the data spread of multiple samples. For each layer, at least three samples were measured. b DFT calculated shift rates for 1L, 2L, 3L, and 8L BP induced by in-plane biaxial strain. c Illustration of two in-plane hoping parameters (t 1 jj and t 2 jj ) and one out-of-plane hopping parameter (t ⊥ ) in a 2L BP. d is the height of an individual layer and D is the gap between two layers. When biaxial in-plane tensile stain is applied, the average distance between two layers (D + d) decreases due to Poisson effect, while the gap (D) increases accompanied by a stronger decrease of d (see Supplementary Fig. 9). d Schematic illustration of the band structure evolution of a bilayer BP under tensile and compressive strain. The orange dashed curves are the bands for a monolayer BP. The change of subband splittings causes the shift rate of E 22 smaller than that of E 11 . Layered materials governed by van der Waals (vdW) interactions offer opportunities for interlayer tuning of the materials' properties. Here, the authors demonstrate that in-plane tensile strain can effectively tune the vdW interactions of fewlayered black phosphorus and weaken its interlayer coupling even though the sample shrinks in the vertical direction ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10483-8 in Supplementary Information), the shift rate of transition energy due to strain is where ε is the applied biaxial strain, N is the layer number, n is the index of subband, h and k are fitting parameters describing the change rates of monolayer band-gap and vdWs interaction strength under strain, respectively. Based on this equation, we can conclude that the shift rate of transition energy is layer (N) dependent and transition index (n) dependent for k ≠ 0 . To extract the values of h and k, we used Eq. (1) to globally fit all of the experiment data, shown as the solid curves in Fig. 3a.
The fitting values of h and k are 66 and −86 meV/%, respectively. The overall trend of the fitting curves is in reasonable agreement with the experiment data. The deviation of data points from the fitting curves is mainly attributed to the experimental uncertainty in determining the strain value. Of course, Eq. (1) does not take into account excitonic effects and other manybody interactions, which may compromise its accuracy. Nevertheless, the basic behavior is well captured by the model. Since h and k are in the same order of magnitude, the modification of vdWs interaction under in-plane strain contributes as significantly as that of intralyer bonding does to the tuning of the band structure. Previously, we experimentally determined the transition energy E nn = E g0 − 2Δγ cos(nπ/(N + 1)), with Δγ = 880 meV, for unstrained few-layer BP 17 . We can see that 1% in-plane strain changes the interlayer coupling (Δγ) almost 10%, which is quite remarkable.
To further confirm our experimental results, we performed DFT calculations of the strain effect. The details of the calculation are presented in the section "Methods". Figure 3b shows the DFT calculated shift rates for 1L, 2L, 3L, and 8L BP, in good agreement with our experiments (see Fig. 3a). The infinite lattice periodicity in the ab plane (in-plane) gives rise to dispersive bands, which have the band gap at Γ point. For those few-layer BPs, the interlayer vdWs interaction along the c-axis (out-ofplane) lifts the energy degeneracy at Γ, thus forming different groups of subbands. For example, we plot in Fig. 3c the lattice structure of a 2L BP, and in Fig. 3d the illustrative band structure of 2L BP and its evolution against strain. When an in-plane tensile strain is applied (see the right part of Fig. 3d), the increasing intralayer atomic distance and thus the decreasing intralayer interaction reduces the band width and enlarges the band gap. Moreover, our calculation also indicates a decrease of the interlayer interaction, since the subband splitting gets smaller. As a result, for the lower-index subband transitions (E nn , n < N/2) in the finite N-layer BP, both the decreasing intralayer and interlayer interactions contribute positively to the increasing E nn , giving a positive energy shift rate against the strain. However, for the higher-index transitions (E nn , n > N/2), the decreasing intralayer (interlayer) interaction yields a positive (negative) contribution to the E nn , making its shift rate smaller than the lower-index E nn . All these mechanisms are well captured by the DFT calculations: for a given few-layer BP, the energy shift rate of E nn gets smaller with the increasing index, and E nn can even become negative, e.g., E 33 for 3L BP and E nn (n ≥ 5) for 8L BP. There is overall a quite good agreement between the DFT calculations and the experiments, as shown in Fig. 3a, b. It should be emphasized that the fitting parameter k in Fig. 3a is negative. A negative k means a weakened interlayer coupling under tensile strain. This is rather astonishing, given that in-plane tensile strain typically compresses the sample in the out-of-plane direction due to the Poisson effect. Now we see that a shorter distance between layers gives a weaker vdWs coupling between them. How could this happen? In order to make sense out of it, we have to revisit the puckered lattice structure of fewlayer BP.
According to the tight-binding model, the band structure in the vicinity of Γ point of the Brillouin zone in few-layer BP is mainly determined by three hopping parameters: two intralayer ones (t 1 jj and t 2 jj ) and one interlayer one (t ⊥ ) 26 , as illustrated in Fig. 3c. The hopping parameter t is generally related to the bond length r in the form of t / 1 r 2 27 . It should be noted that the vdWs interaction between layers is a complicated and sometimes vague concept. Here we just capture the major contributing factor t ⊥ for our discussion, and treat it on the same footing as inplane bonds. For monolayer BP, the bandgap can be expressed as E g = 4t 1 jj + 2t 2 jj . While for 2L (or few-layer) BP, interlayer interaction (t ⊥ ) also contributes to the band structure, and leads to the splitting of valence and conduction bands, resulting in two subband transitions: E 11 and E 22 . The energy difference between E 11 and E 22 (Δ = E 22 − E 11 ) is proportional to t ⊥ 26 .
Based on the above analysis, one can directly understand how the band structure of 2L BP evolves with the in-plane biaxial strain. Under biaxial tensile strain, the absolute values of t 1 jj and t 2 jj both decrease. But it is still not straightforward to predict how E g changes, since the two hopping parameters have different signs (t 1 jj < 0, t 2 jj > 0) 26 . However, the experiment results show that both E 11 and E 22 blueshift, indicating biaxial tensile strain induces an increase of E g . Moreover, when tensile strain is applied, atoms connected by t ⊥ are split apart laterally, since the two atoms are not in registry (one right above the other). More importantly, those two atoms are split apart vertically as well, even though the Poisson effect compresses the sample in the vertical direction, as shown by our DFT calculations (see Supplementary Fig. 9).
This unusual scenario originates from the puckered structure within each BP layer, which is constituted of two sub-layers. When biaxial in-plane tensile stain is applied, the average distance between two adjacent layers (D + d, shown in Fig. 3c) decreases due to the Poisson effect. According to the DFT calculations ( Supplementary Fig. 9), when 1% biaxial tensile strain is applied, the interlayer distance (D + d) between two phosphorene layers decreases by 0.031 Å, corresponding to an out-of-plane Poisson's ratio ν ¼ À dε z dε xy~0 .5. However, at the same time, the distance between two sub-layers (d) within each layer decreases even more thanks to the puckered structure (0.087 Å, with 1% biaxial tensile strain, see Supplementary Fig. 9), causing an increased gap (D) between two layers (0.056 Å, with 1% biaxial tensile strain). As a result, the distance between the two atoms responsible for t ⊥ increases, leading to the decrease of t ⊥ and hence the decrease of interlayer interactions. The puckered structure of BP efficiently facilitates this counterintuitive phenomenon. Therefore, when the interlayer interaction decreases, the energy difference between E 11 and E 22 also decreases. Meanwhile, tensile strain leads to the blueshift of both E 11 and E 22 , hence the shift rate of E 11 is larger than that of E 22 , as illustrated in Fig. 3d. With above argument, we conclude that the higher the transition index is, the smaller the shift rate is. As demonstrated by DFT calculation in Fig. 3b, the shift rate even changes to negative sign when the subband index n is >5 for the 8L BP. The split subbands of few-layer BP provide us a unique platform to monitor the vdWs interactions under strain or other physical perturbations.
In summary, we systematically investigate the evolution of band structures in few-layer BP under in-plane biaxial strain. Our results show that strain effect of few-layer BP exhibits layer number and subband index dependence. This dependence is closely related to the tuning of interlayer coupling. Surprisingly, in-plane biaxial tensile strain weakens such coupling, which is counterintuitive. Indeed, it is the puckered structure of BP which facilitates this phenomenon. Our study paves the way for strain engineering of vdWs interactions in few-layer BP and other 2D materials.

Methods
Sample preparation and the strain setup. Few-layer BP flakes were mechanically exfoliated from bulk crystals (HQ graphene Inc.) onto PDMS substrate, and then transferred to PP substrates by a dry-transfer method. To apply biaxial strain, the sample was loaded into a heating and cooling chamber (Linkam, FTIR600). To avoid sample degradation, all sample preparation procedures were conducted in a glove box filled with N 2 gas (O 2 and H 2 O < 1 ppm). For optical measurements, which were outside of the glove box, the chamber was flushed with N 2 gas to protect BP samples.
Optical measurements. The IR extinction (1 − T/T 0 ) spectrum of few-layer BP was recorded using a FTIR spectrometer (Bruker Vertex 70V) in conjunction with an IR microscope (Hyperion 2000). By using tungsten halogen lamp as IR source, quartz beam splitter and liquid nitrogen cooled MCT detector, wide range extinction spectrum from 0.25 to 1.36 eV could be obtained. The incident light was polarized along the armchair (AC) direction through a broadband ZnSe grid polarizer. The visible reflection (ΔR/R 0 ) measurements of few-layer BP were performed using an Andor grating spectrometer (SR550i) in conjunction with a Nikon inverted microscope (Elipse Ti-U). By using a broadband LED light source and a CCD detector, the visible reflection spectrum ranging from 1.55 to 3.14 eV could be covered. Raman measurements were conducted using a commercial Raman system (Horriba, Jobin Yvon HR-Evolution 2). 532 and 633 nm lasers were available for excitation.
DFT calculations. All the DFT calculations were performed using the Vienna ab initio simulation package (VASP) 28,29 . The generalized gradient approximation of Perdew-Burke-Ernzerhof (PBE) 30 type was employed for the exchange-correlation energy. The interlayer interactions, vdWs interactions, were included by using the optB88 vdW density functional 31,32 . To describe the band gaps more accurately, the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional 33 was also used. Wave functions were expanded in terms of plane waves, for which the kinetic energy cutoff was 500 eV, and the energy integration was performed on a 10 × 12 × 1 Monkhorst-Pack k-point grid. All the structures in our calculations were fully optimized with a force criteria of 0.01 eV/Å. For the calculations of all the few-layer BPs, slab models were employed, in which a 20 Å vacuum was added along the c direction.

Data availability
The data that support the findings of this study are available from the corresponding author upon request.