Raman Activity of Multilayer Phosphorene under Strain.

Using computational tools, we study the behavior of activities of lattice vibrational Raman modes in few-layered phosphorene of up to four layers subjected to a uniaxial strain of -2 to +6% applied in the armchair and zigzag directions. We study both high- and low-frequency modes and find very appreciable frequency shifts in response to the applied strain of up to ≈20 cm-1. The Raman activities are characterized by Ag 2/Ag 1 activity ratios, which provide very meaningful characteristics of functionalization via layer- and strain-engineering. The ratios exhibit a pronounced vibrational anisotropy, namely a linear increase with the applied armchair strain and a highly nonlinear behavior with a strong drop of the ratio with the strain applied along the zigzag direction. For the low-frequency modes, which are Raman active exclusively in few-layered systems, we find the breathing interlayer modes of primary importance due to their strong activities. For few-layered structures with a thickness ≥4, a splitting of the breathing modes into a pair of modes with complementary activities is found, with the lower frequency mode being strain activated. Our calculated database of results contains full angular information on activities of both low- and high-frequency Raman modes. These results, free of experimental complexities, such as dielectric embedding, defects, and size and orientation of the flakes, provide a convenient benchmark for experiments. Combined with high-spatial-resolution Raman scattering experiments, our calculated results will aid in the understanding of the complicated inhomogeneous strain distributions in few-layered phosphorene or the manufacture of materials with desired electronic properties via strain- or layer-engineering.

Introduction

Few-layered black phosphorus (BP) (Figure ) is a two-dimensional (2D) material that has recently received renewed attention.[1−9] This interest was fueled by its intrinsic strongly thickness-dependent direct bandgap ranging from 0.3 for bulk[10] to ≈2.0 eV for a monolayer,[1,6,8,9] thus bridging bandgaps between graphene and transition-metal dichalcogenides. The strong thickness dependence of electronic gaps was attributed to interlayer coupling, which is appreciably stronger in BP than in other 2D materials, such as graphene or transition-metal dichalcogenides.[11,12] The room-temperature free-carrier mobility approaches ≈1000 cm2 V–1 s–1,[3] which exceeds that of other typical 2D semiconductors, such as molybdenum sulfate (≈200 cm2 V–1 s–1 [1,13]). Combined with a high on–off ratio,[3] these properties hold prospects for electronic and optoelectronic applications.[1,2,14,15]

Figure 1

Top (a) and side (b) views of the atomic structure of few-layered phosphorene for the case of four-layered BP. (c,d) Examples of phonon polarization vectors illustrated for four-layered BP. (c) High-frequency Raman modes (Ag1 out-of-plane, B2g along ZZ, and Ag2 along AC). (d) Strongest-intensity low-frequency Raman breathing modes (B1(Ag) and B2(Ag)). The other LF modes have intensities that are several orders of magnitude lower.

In addition to layer-engineering, the unique electronic properties of few-layered phosphorene (FLP) can among others further be tuned by dielectric embedding and applied strain. Dielectric embedding refers to naturally entering a substrate for anchoring a 2D material and often, specifically for phosphorene, also as a capping material to avoid degradation due to oxidation. Dielectric embedding is expected to affect appreciably the dielectric function, exciton binding energies, quasiparticle gaps,[9] and vibrational properties.

The importance of strain-engineering, or straintronics,[16] in 2D materials is due to their ability to sustain much larger strains in excess of ≈10%, both compressive/tensile in-plane and normal compressive,[17] compared to just ≈2% by their three-dimensional counterparts.[16,18,19] For instance, single-layered MoS2 has a huge strain sensitivity of the bandgap from 1.8 to 0 eV in comparison with the poor window of tunability of only 0.25 eV achieved for strained silicon at 1.5% biaxial strain.[16] In addition to the in-plane strain, the effect of ripples and wrinkles with significant out-of-plane deformations has also been studied[20,21] and their strongly anisotropic effect on electronic properties and the chemical activity of phosphorene is predicted.[21] A wealth of properties under strain has been reported in FLP such as topological phase transition from normal to the quantum-spin Hall state,[22] onset of superconductivity,[23] a complex dependence of the bandgap ranging from metallic to semiconducting and featuring both direct and indirect bandgaps in single- and few-layered phosphorene,[24−26] a negative Poisson’s ratio,[27] a strain-engineered ultrahigh carrier mobility,[28] strain-modulated thermal conductivity,[29] or funneling of excitons in MoS2 toward the point of applied strain upon illumination.[30]

Phonons, in addition to electrons, are the most important collective excitations as their interactions with electrons and photons determine material’s properties such as the carrier mobility, light-emitting efficiency, etc. In FLP, phonons have been found to be exceptionally anharmonic and anisotropic[31,32] as evidenced by angle-resolved absorption and Raman spectroscopies, which also contribute to electron–photon and electron–phonon anisotropies.[33] Furthermore, these anisotropies show a nontrivial dependence on the sample thickness and phonon and photon energies.[33] Vibrational properties of 2D materials can be used to determine properties of materials prepared by layer-engineering, dielectric embedding, and straintronics. For instance, the relative layer thickness of FLP, in addition to optical contrast measurement, can be determined from vibrational analyses[19] based on Raman active (R) low-frequency (LF) breathing modes active only in FLP, which are absent from single-layered and bulk BP.

Strain in real materials is unusually complicated and in a 2D structure usually induces an inhomogeneous distribution. Hence, measuring and understanding strain and its spatial distribution is a primary prerequisite for strain-engineering. Phonons and Raman active modes, due to their high spatial resolution, play an important role also in extracting strain information. Indeed, substantial frequency shifts of Raman modes under various strain conditions have been reported in FLP experimentally[12,32,34] and in density functional theory (DFT) modeling,[26,32,34,35] including strain effects on the low- and ultralow-frequency modes,[12] paving the way to the measurement of strain and its inhomogeneous distributions.

Furthermore, due to the strong anisotropy of BP, the susceptibility of the vibrational modes is also very selective to the strain applied in the armchair (AC) and zigzag (ZZ) directions.[26,32,34,35]

It is expected that the vibrational properties will also be affected by the dielectric embedding via modification of the dielectric function. However, to focus on a single phenomenon, here, we study freestanding FLP, which should be kept in mind when comparing the simulations with the experimental results.

In FLP, two very different types of Raman modes are present: high-frequency (HF, >350 cm–1) intralayer and low-frequency (LF, <100 cm–1) interlayer modes.[19] While the HF modes involve vibrations from intralayer chemical bonds, the LF modes probe the relative motion of the layers and, hence, the predominantly van der Waals (vdW) interlayer interactions. In an unstrained FLP, Ag1, B2g, and Ag2 modes are normally experimentally observed at frequencies of ≈365, ≈440, and ≈470 cm–1,[19] respectively (see Figure . Following the vibration-mode symmetry and back-scattering setup, the intensity of the Ag modes changes with the variation period of 180°, but the maximum intensity occurs at different polarization angles, depending on the sample thickness.[19,32] The intensity of the Ag1 mode appears to be especially thickness-sensitive.[32] The intensity of the B2g mode has a variation period of 90° and the minimum intensity corresponding to the armchair and zigzag polarizations.[19,32] Because of the crystal symmetry and corresponding selection rules, the LF modes are not active in single-layered and bulk BP. They appear in FLP in the form of three types of LF R modes, two shear and one breathing (B1g, B3g, and Ag),[19,26] experimentally observed at ≈26.2, ≈75.6, and 85.6 cm–1.[19] While all three R LF modes have intensities significantly weaker than the R HF modes, the breathing Ag mode again stands out as having an intensity significantly higher than those of the other two shear modes. In particular, starting from four layers, the breathing Ag mode splits into two modes, B1 and B2, as shown in Figure .

While several experimental[12,32,34] and simulation studies[26,32,34,35] of strain effects on vibrational modes in FLP exist, to the best of our knowledge, no computational study of strain effects on intensities of Raman vibrational modes in FLP exists. Based on a computational density functional theory (DFT)[36] study, the aim of the present paper is to fill the gaps in understanding vibrations in FLP under strain, including the R intensities of both HF and LF modes. This is done by focusing on the strongly anisotropic effects of uniaxial AC and ZZ strains on vibrational properties of FLP of up to four layers, a thickness at which the vibrational properties are already reminiscent of those of the bulk.[37]

Results and Discussion

We now discuss the effect induced on the Raman modes by the different strains (AC and ZZ) applied to mono-, bi-, tri-, and four-layered BF.

Layer Dependence

To set the stage, we show in Figure the dependence of the activity of R HF modes on the number of layers L for unstrained FLP for L = 1–4. One can clearly see that the Ag1 and Ag2 modes have activities appreciably larger than that of the B2g mode, which has a negligible activity for L = 1, 2, and 3 and a fairly weak activity only for L = 4 of the B2g mode. The other trend clearly visible is the increase in the activity of the Ag1 mode relative to the Ag2 with the number of layers.

Figure 2

Activity of HF Raman modes in FLP with different numbers of layers.

Qualitatively, both trends are supported by experiments, which show that the intensity of the Ag1 mode is extremely sensitive to and increases with the sample thickness.[32]

The experimental intensity of the B2g mode also increases with the number of layers,[37] and hence this trend is visible in our four-layer results. The intensity of the B2g mode also shows an extreme angular dependence, and consequently, the experimental intensities vary considerably with the polarization angle (see Figure S4). The activity of the B2g mode is significantly smaller than that of the Ag modes, and consequently, the B2g mode is barely visible in our angularly averaged activities dominated by the A modes. It should be noted that the published ratios of intensities of the R HF A and B modes vary appreciably between different experiments.[5,34,37,38] We attribute this to the angular dependence of these modes (see Figure S4).

Effect of Strain on High-Frequency Raman Modes

The effect of strain on activities of HF R modes was studied for FLP in the range of L = 1–4. An example of results for a bilayer is shown in Figure ; for the rest of the results, see the Supporting Information (SI). In analogy with experiments,[32,34] we first analyze the strain-induced shifts of R peaks. For strains in the window of ⟨−2, +2%⟩, the AC and ZZ strains cause linear shifts of the peak positions. More importantly, the AC and ZZ strains have completely opposite effects on the highest activity of R HF modes. ZZ strain causes red shifts of Ag2 and B2g modes, while the Ag1 mode exhibits only a negligible peak shift. In contrast, AC strain causes blue shifts of Ag2 and B2g modes, while the Ag1 mode exhibits a red shift. The largest frequency shifts we observe in the studied range of strains ⟨−2, +4%⟩ are very significant (≈20 cm–1). The B2g mode for strain applied in the ZZ direction exhibits an even larger response, albeit at very low activity. In general, the shifts in the frequencies of the R HF modes are larger for strains applied in the ZZ direction than in the AC direction. Similar trends apply also to the other L (see the SI).

Figure 3

Intensity of the key R HF modes under strain in bilayer BPs. Top, AC strain; bottom, ZZ strain. Dashed black lines are to guide the eyes.

Our results on stress-induced frequency shifts can be compared with experimental results corresponding to uniaxial strain applied separately approximately in the AC and ZZ directions albeit on significantly thicker phosphorene samples (2–9 layers[34] and 5–30 layers)[32] and to DFT calculations[32] on single- and double-layered phosphorene in a set-up (PBE xc-functional[39] and optB86b-vdW[40] to model the vdW interlayer interactions) very similar to ours (PBE xc-functional and DFT-D2 model to describe the vdW interactions). Our calculated strain-induced frequency shifts (Figures and S4) can in a limited range of strains of approx. ±2% be approximated by straight lines. Experimentally, a linear behavior is observed for a much wider range of strains of up to ≈7% on much thicker samples.[32] The calculated rates of strain-induced frequency change per strain-percentage are compared with experimental[32,34] and alternative simulation results[32] in Table . The strain response is significantly stronger in the ZZ than in the AC direction. Qualitative agreement with experiments is very good. Similar to the alternative simulation results,[32] our calculated trends exhibit a stronger strain response of the B2g and Ag2 modes in the ZZ direction than in the experiments. The weaker experimental response may be due to the different experimental sample thicknesses and possible sample misalignment. Perhaps surprisingly, even the experimental results exhibit quantitative differences of ≈50% or occasional opposite trends, which we speculate may result from sample misorientation, different thicknesses, or dielectric embedding. We have tested the influence of the interlayer electronic coupling by comparing selected vibrational frequencies in the presence and absence of the D2 correction[41] to PBE. The results in the SI show that the inclusion of the D2 correction invariably increases all Raman vibrational frequencies, both high and low frequency. This also means that even in the absence of interlayer electronic coupling and/or vdW interactions, purely classical long-range dipole or higher multipole interactions can play a role in influencing the vibrational properties.

Table 1

Experimentally Determined[32,34] Rates of Strain-Induced Frequency Change per Strain-Percentage of the Principal High-Frequency Raman Ag1, B2g, and Ag2 Modes Compared to Previous[32] and Present Simulations

mode	Ag1				B2g				Ag2
strain	exp.[34]	exp.[32]	th.[32]	th.	exp.[34]	exp.[32]	th.[32]	th.	exp.[34]	exp.[32]	th.[32]	th.
AC	–1.3	+1.8	–2.9	–2.4	+0.1	+1.3	+0.8	–0.9	+0.3	+0.5	+0.8	+0.5
ZZ	–0.1	0	+2.1	–0.3	–2.4	–3.6	–8.2	–8.6	–1.1	–1.6	–4.4	–3.4

Compared to frequency shifts, the intensities are more difficult to quantify as both experimental spectra[42] and calculated spectra may use different normalizations. Therefore, the ratio of intensities of the main peaks, Ag2/Ag1, is often used, making the ratio less susceptible to normalization. In addition, the intensity ratio is also used to determine the number of layers[37] or oxidation level in photo-assisted oxidation.[5] The Ag2/Ag1 peak ratio characteristics are based on the sensitivity of the Ag1 peak to the layer thickness.[32] In contrast, for a given L, the Ag2/Ag1 peak ratio has only a negligible angular dependence (see also Figure S4), making results of different experiments easily comparable.

As an approximation to the experimental Raman Ag2/Ag1 intensity ratio, we calculate Raman Ag2/Ag1 activity ratios in equilibrium and also as a function of strain applied along the AC and ZZ directions. The results are shown in Figures and 5. The different types of strains have vastly different effects on the Ag2/Ag1 activity ratios. While the application of an AC strain leads essentially to a uniform linear increase of the ratio in the studied range of strains, application of the ZZ strain exhibits a highly nonuniform behavior with an essentially linear increase from tensile to moderate compressive strains of ≈2% and then a sharp drop for larger compressive strains. We note that this behavior is unlikely to be induced by a direct to indirect gap transition or by a semiconductor to metal transition in response to the applied strain as the gap response in our DFT model is broadly similar for the strain applied in the AC and ZZ directions (see the SI).

Figure 4

Intensity ratio of the key high-frequency modes as a function of the applied strain (top) and the number of layers (bottom) for strain applied in the AC direction.

Figure 5

Intensity ratio of the key high-frequency modes as a function of the applied strain (top) and the number of layers (bottom) for strain applied in the ZZ direction. Note the deviation of the results for L = 3 from general trends likely caused by numerical factors. These results are not included in the bottom panel.

Figures and 5 also show the Ag2/Ag1 activity ratio plotted against the number of layers, L, for strain applied in the AC and ZZ directions, respectively, which is the usual form in the unstrained FLP.[5,37] The plot for AC strain shows the typical saturation around L ∼ 3 also found experimentally[37] at zero strain, which corresponds to the intensity saturation of the Ag1 peak with L (see also Figure ). For the strain applied in the ZZ direction, the plot is initially qualitatively similar to the AC strains of up to ≈+2% followed by a steep drop for strains ≥2%. We note that measured Ag2/Ag1 ratios of Lu et al.[38] and Zhang et al.[37] for unstrained FLP differ by a factor of 5, with our results being close to those of Lu et al.[38] The main quantitative difference between the AC and ZZ strains is that the applied compressive strain tends to flatten the curve and that the effect is stronger for the ZZ strain than for the AC strain up to ≈2% at which strain a significant change in the Ag2/Ag1 activity ratio is predicted in the ZZ direction, reconfirming the strongly anisotropic behavior of the Raman vibrational modes.

Effect of Strain on Low-Frequency Raman Modes

Despite the, in general, much lower intensities of the R LF modes, they play a very important role in the FLP.[12,19] Their behavior under applied AC and ZZ strains is shown in Figure for a four-layered BP. Four-layered BP is special in that, here, as shown in Figure , two breathing modes are found, B1 and B2. For thinner FLP, the peaks shrink into a single feature, while for much thicker systems, intensities of the breathing modes are expected to vanish.[19,26] From Figure , three trends are evident: (1) the frequency of both breathing modes decreases with increased compressive strain applied both along the AC and ZZ directions, except for the largest strain of +6% applied in the AC direction; (2) the shifts of the Raman frequencies are again quite significant, up to ≈15 cm–1, especially for the AC strain; (3) the activities of the B1/B2 modes are complementary, that is, where the activity of one increases in response to applied strain, the other decreases and vice versa. For clarity, the last feature is depicted in Figure where the complementary intensities are clearly visible. Interestingly, unlike the B2 mode, which shows appreciable activity at any strain load, including equilibrium, the B1 mode appears to be strain activated, equally by both AC and ZZ strains, and its predicted activity in the unstrained state vanishes.

Figure 6

Intensity of the key R LF modes under strain in four-layered BP. Top, AC strain; bottom, ZZ strain.

Figure 7

Raman intensity (normalized; I1 + I2 = 1) of the key low-frequency modes as a function of strain for AC (top panel) and ZZ strains (bottom panel).

The low-frequency vibrational modes have been experimentally studied.[12,19] Ultralow- (≈15–40 cm–1) and low-frequency (≈70–90 cm–1) modes[12] with completely different properties were reported. While the ultralow-frequency mode requires the presence of a substrate, the low-frequency mode, experimentally located at ≈80 cm–1 in four layers, corresponds to our B2 mode (≈90 cm–1 without strain). Importantly, no B1 mode (≈40 cm–1 from our simulation) was detected in the experiment without applied strain, which is in agreement with our model.

Conclusions

In conclusion, using density functional theory modeling, we have presented the lattice vibrational study of Raman modes of strained few-layered phosphorene, up to a thickness of four layers, focusing on Raman activities, which are closely related to the experimentally measured intensities. We find that the Raman modes exhibit a qualitatively complex strain response, which is quantitatively strong. Furthermore, in addition to strongly anisotropic excitons,[6] electron mobility,[28] and thermal conductivity,[29] strong anisotropy is also found in vibrational properties, which yields qualitatively very different responses to strains applied in the AC and ZZ directions. Indeed, our simulations predict a significant modification of HF vibrational properties for compressive strains in excess of 2% applied in the ZZ direction.

Our study also pays attention to the low-frequency modes, which, up to date, have received much less attention[12,19] despite their indisputable importance. We observed two LF R interlayer breathing modes (B1, B2) with complementary intensities. The B1 mode appears to be strain activated. The B2 mode is in good agreement with the experimentally detected LF mode.[12] We note that the strong interlayer interaction in BP is not due to the electronic hybridization of the lone electron pairs as suggested previously,[11] but, as follows from accurate quantum Monte Carlo (QMC) calculations of Shuleburger et al.,[43] is rather purely vdW. The strong interaction can naturally be explained by the fairly large c6 constant of phosphorus, compared to those of carbon (graphene) or sulfur (molybdenum sulfide).[44]

A comparison of our calculated predictions with experiments shows generally excellent qualitative agreement with experiments, where available. At the quantitative level, we find differences not only between calculated and experimental results but also between different experimental sets.[5,37,38] A comparison of experimental results for 2D materials is complicated by differences in dielectric embedding, presence of defects, and sizes of the flakes, which cause scatter in the experimental data. In addition to vibrational properties discussed here, scatter between different experimental sets is also observed in electronic properties, which causes a scatter of ≈1 eV in the measured optical band gap of single-layer phosphorene.[9] Hence, benchmark results free of such experimental complexities, such as ours, are required.

Our calculated results, combined with experimental high-spatial-resolution Raman scattering experiments, aid in a qualitative understanding of the complicated inhomogeneous and anisotropic strain distributions in few-layered phosphorene. Since strain is predicted to be an efficient tuning tool for bandgap[2,16] and carrier mobility,[28] we expect our results to be of use also for the manufacture of materials with desired electronic properties via strain- or layer-engineering.

Computational Methods

The structure of FLP (single to four layers) was represented by a 5 × 5 supercell with a vacuum of 17 Å included in the z-axis to avoid spurious interactions between periodic images. All structures were fully optimized using DFT techniques[36] as implemented in the VASP code.[45,46] In FLP, the van der Waals interlayer interaction is critical for the description of the interlayer modes. A comparison with ultra-accurate quantum Monte Carlo (QMC) calculations[43] shows that the interlayer interactions both in bulk and in double-layered BP are very well-described with the PBE-D2[39,41] approach, which we have also used in the present study. It should be noted that none of the DFT models tested generates charge densities in complete agreement with QMC.[43] The important message to be noted with regard to the QMC-calculated charge densities is that there is absolutely no interlayer hybridization and, hence, no interlayer chemical bonding as concluded from the DFT description.[11] Instead of hybridization, the interlayer vdW attraction causes Pauli repulsion between the BP layers.[43] A projector augmented-wave pseudopotential[47] with a plane-wave cutoff of 500 eV and a Monkhosrt–Pack mesh[48] of 12 × 8 × 1 was used.

In the next step, the Raman activities were determined. To reduce the complex phonon anisotropic behavior,[19] unless stated otherwise, the Raman activities Am of the ωm mode were computed as spherical averages of the Raman tensor of that mode, αm, over incoming and scattered light-beam polarizations,with the Raman tensor, αm, being the derivative of the macroscopic dielectric tensor with respect to the set of atomic displacements that correspond to the phonon eigenvector.[49] The dielectric tensor components were computed in the linear response theory (DFPT) as implemented in the VASP code.[45,46] The effect of the averaging is that the A modes are significantly stronger weighted than the B mode. This is because the average is dominated by the much higher activity of the A modes (see Figure S4). The B mode has higher activities at angles around 45° where the A modes have low activities, making activities of all three modes comparable at these angles (see Figure S4). Only vibrational spectra at the Γ point of the Brillouin zone were considered. The calculation of Raman activities followed that of Porezag and Pederson[50] as implemented in ref (51). The related Raman intensities, which take into account also the experimental temperature, the exciting frequency of the laser light, etc., could be calculated from the activity data.[52] In either case, a line shape, such as Lorentzian or Gaussian, has to be imposed, which affects the peak heights.

A compressive/tensile uniaxial in-plane strain was applied both in the AC and ZZ directions. All strains were applied in the range of −2 to ≈+6%. These limits were chosen so as to limit the strain effects on the closure of the HOMO–LUMO gaps in our DFT description (see the SI and ref (26). In this limited strain range, Dirac minicones are formed at compressive strains (< −2%) and indirect gaps are formed in the tensile strain regime. More general technical details can be found in the SI.

We note that despite limiting our presentation essentially only to averaged activities, our database of results contains full angular information, and any angle-resolved Raman intensity could be reconstructed and compared with experimental results (for an example, see Figure S4).