Infrared fingerprints of few-layer black phosphorus

Black phosphorus is an infrared layered material. Its bandgap complements other widely studied two-dimensional materials: zero-gap graphene and visible/near-infrared gap transition metal dichalcogenides. Although highly desirable, a comprehensive infrared characterization is still lacking. Here we report a systematic infrared study of mechanically exfoliated few-layer black phosphorus, with thickness ranging from 2 to 15 layers and photon energy spanning from 0.25 to 1.36 eV. Each few-layer black phosphorus exhibits a thickness-dependent unique infrared spectrum with a series of absorption resonances, which reveals the underlying electronic structure evolution and serves as its infrared fingerprints. Surprisingly, unexpected absorption features, which are associated with the forbidden optical transitions, have been observed. Furthermore, we unambiguously demonstrate that controllable uniaxial strain can be used as a convenient and effective approach to tune the electronic structure of few-layer black phosphorus. Our study paves the way for black phosphorus applications in infrared photonics and optoelectronics.

S ince the isolation of graphene in 2004 (ref. 1), tremendous attention has been paid to the family of two-dimensional (2D) materials. Recently, black phosphorus (BP) was reintroduced as a new 2D material [2][3][4][5] , exhibiting many intriguing properties such as highly tunable bandgap 6,7 , anisotropy 5 and relatively high carrier mobility 2 . It has been predicted that the bandgap of BP is always direct regardless of layer (L) number and ranges from 0.3 to 2 eV 6,8 , bridging the gap between zero-gap graphene and large-gap transition metal dichalcogenides 9 . Moreover, BP has a puckered hexagonal structure with two non-equivalent directions in the layer plane: armchair (AC) and zigzag (ZZ) (Fig. 1a). Arising from the structural anisotropy, strongly anisotropic mechanical 10 , thermal 11 , electrical 3,5 and optical 5,6,8,[12][13][14][15] properties have been highlighted in recent studies, opening up possibilities for conceptually new devices.
Compared with the bulk counterpart, one of the most intriguing distinctions for single or few-layer 2D materials is the highly tunable physical properties, through various techniques. This tunability is typically associated with the modification of the electronic band structure. A non-destructive and accurate characterization technique to monitor few-layer BP band structures is highly desirable, given its predicted strong dependence on layer thickness 6 , stacking order 16 , strain 7 and doping [17][18][19] . Previous photoluminescence (PL) 14 and differential reflectance 20 studies of BP are limited to the visible and near-infrared range, and are available only for thin BP layers with layer number o5. With the majority of the optical transitions expected in the mid-to near-infrared frequency range for few-layer BP, Fourier transform infrared spectrometer (FTIR)-based infrared spectroscopy is believed to be the superior characterization tool. However, up to date, such infrared study for mechanically exfoliated few-layer BP (o15L), with frequency ranging from the mid-to the near-infrared, is still lacking.
Here we systematically investigate the evolution of electronic structures in few-layer BP with layer number ranging from 2 up to 15, and report the experimental demonstration of highly tunable electronic structures in few-layer BP via controllable uniaxial strain [21][22][23][24][25] , using polarized infrared spectroscopy. For each few-layer BP, the infrared spectrum typically exhibits layer-dependent multiple optical resonances and can be readily served as its fingerprints. The infrared absorption shows strong polarization dependence, with strong optical resonances showing up in the AC direction. This dependence provides us a reliable way to determine the crystallographic orientation, which complements polarized Raman spectroscopy. For the latter, however, excitation wavelength and BP thickness complicate the polarization behaviour 15,26 . A simple tight-binding model, with only two fitting parameters, can well describe the major optical transitions for all of the measured BP layers. In addition to the main transitions, we also observed unexpected weak absorption features right in the middle of the adjacent main peaks, which we attribute to the forbidden optical transitions. These transitions, expected to be inactive in symmetric BP quantum wells (QWs), are made possible by unintentional doping. The physical origin of these weak features has been further confirmed by controlled chemical doping. By collecting the infrared spectra for different thickness BP layers, we provide a spectra database, which can be used to determine the layer thickness and identify unusual stacking order. More importantly, we also show that few-layer BP is highly sensitive to strain. The bandgap can be modulated by 420% with 1% uniaxial strain for a 6L BP. Surprisingly, this effect shows little dependence on the strain direction. Our study demonstrates that infrared spectroscopy is an ideal scheme for investigations of band structure engineering through mechanical strain, hydrostatic pressure 27 , electrical field 18 , magnetic field 28   doping. The rich band structures of few-layer BP and their potential large tunability promise a wide range of applications in infrared photonics and optoelectronics 29 such as polarizationsensitive photodetectors 12,30 , modulators, strain sensors and infrared lasers.

Results
Sample preparation and polarized infrared spectroscopy. Few-layer BP flakes were prepared by mechanical exfoliation of bulk crystals (HQ Graphene, Inc.), with areas varying from several 100 to 10,000 mm 2 . To minimize the influence of Fabry-Perot interference of the substrate, thick quartz substrates (B0.3 mm in thickness) were used to support BP flakes. Sample thickness was first determined from the optical contrast 31 under a microscope and then was further verified by infrared absorption spectra, as we will discuss later. According to our optical contrast analysis of many BP flakes (more than 100 few-layer flakes) with different thickness, we were able to quite accurately determine the thickness of BP flakes, especially for thinner ones (see Supplementary Fig. 1). To minimize the effect of sample degradation in air, infrared measurements were typically finished within 1 h after sample exfoliation. The infrared transmission (extinction) spectra with photon energy from 0.25 to 1.36 eV were obtained using a FTIR in conjunction with an infrared microscope at room temperature (see Methods). The lower bound cutoff photon energy is restricted by the quartz substrate. Figure 1b is an optical image of a representative 6L sample, with areas of B5,000 mm 2 , large enough for us to obtain accurate infrared extinction spectrum. For an atomically thin layer material sitting on a thick transparent substrate, the extinction (1 À T/T 0 ) is directly proportional to the real part of the optical conductivity s(:o) 32,33 , where T and T 0 denote the transmittance of samples on substrate and bare substrates, respectively. We measured polarization-resolved infrared absorption spectra of the 6L sample, with normal light incidence and polarization angles ranging from 0°to 360°in steps of 15°( Supplementary  Fig. 2). For clarity, only six spectra are shown in Fig. 1d. For such large area samples, the measured extinction can be directly translated into the real part of the optical conductivity according to the formula 32,33 Res(:o) ¼ (1 À T/T 0 ) Á (n s þ 1) Á c/8p, where n s is the refractive index of the quartz substrate and c is the speed of light. Two prominent peaks (labelled as E 11 and E 22 ) can be identified, revealing rich features of the underlying electronic structures. Owing to the confinement along the z direction, conduction bands and valence bands split into multiple subbands with well-defined quantum numbers n and optical transitions are allowed only from valence subbands to conduction subbands with the same index n for normal light incidence 17,18 . In low-dimensional materials, it is well known that the optical resonance peaks arise from excitons, due to significantly reduced screening of Coulomb interactions [34][35][36] . However, for clarity, we adopt the band to band transition picture, rather than the excitonic picture when we refer to the measured infrared resonance peaks. We assign the two peaks E 11 and E 22 to optical transitions v 1 -c 1 and v 2 -c 2 , respectively, illustrated in the left panel of Fig. 1c. Higher-order subband transitions are beyond the limit of our measurement range. Previous theory 6,8,17 and experiments 5,12 have demonstrated strong anisotropy in infrared absorption for thick BP layers: BP absorbs light polarized along the AC direction, with much less absorption for light with polarization along the ZZ direction, because mirror symmetry in the x-z plane forbids optical transitions for the ZZ polarization. Therefore, the optical conductivity of BP maximizes in the AC direction, as seen in Fig. 1e. By rotating the linear polarizer, the AC direction is determined to be B2°relative to the pre-selected x direction. This unique anisotropy makes infrared spectroscopy a very convenient and unambiguous method to determine the crystallographic orientation. For the peak E 22 , Re(s) reaches a maximum of 2.6s 0 (s 0 ¼ pe 2 /2h is the renowned universal optical conductivity of graphene 33,37 ). This conductivity gives an absorption of B6% at E 22 for a suspended 6L BP, indicating a 1% absorption per BP layer. The most interesting feature is the weak peak between E 11 and E 22 , labelled with an asterisk (*). We assign it to the transitions between v 1 and c 2 , and v 2 and c 1 , with transition energies denoted as E 12 and E 21 , respectively, as shown in the right panel of Fig. 1c. As these transitions involve bands with different quantum numbers 18,38 , we refer them as hybrid transitions. More detailed discussions will be given later.
The evolution of the electronic structures with layer number. We performed polarization-resolved infrared absorption measurements on few-layer BP from 2L to 15L, as shown in Fig. 2a-j. For comparison, we also measured the absorption spectrum for a bulk BP (Fig. 2k, thickness 4100 nm). All the experiments were performed at room temperature under ambient conditions. The optical features of monolayer BP lie in the visible frequency range, beyond the infrared measurement limit 3,14 . As seen in Supplementary Fig. 3, the bandgap PL emission peak is around 1.67 eV. Except for the 6L and a couple of other samples, the sample size may be smaller than the infrared beam size. As a result, the light extinction values are usually smaller than the true value and comparison of the absolute extinction values for different layers in Fig. 2 is not meaningful. Nevertheless, the peak positions are very informative. For the 2L sample, a salient peak around 1.14 eV is observed with incident light polarized along the AC direction (black line in Fig. 2a). We assign this peak to the lowest energy transition v 1 -c 1 , labelled as E 11 . As predicted by theory 17 , the spectrum in the ZZ direction is featureless (red line in Fig. 2a). With the increase of layer number, the peak E 11 undergoes a monotonous red shift from 1.14 eV (2L) to 0.38 eV (15L) and eventually reaches 0.34 eV in the bulk limit. The observed behaviour qualitatively agrees with previous theoretical predictions, in which strong interlayer interaction is responsible for the bandgap decrease 6,8,39 .
Interestingly, a new series of peaks appear in the absorption spectra for few-layer BP, in addition to the lowest energy peaks. All of the higher energy peaks exhibit similar polarization dependence to that of E 11 transitions. For the 5L sample, the transition E 22 (v 2 -c 2 ) peak is at 1.30 eV. Up to 15L, the peak E 22 red shifts to 0.48 eV, showing a similar layer dependence to that of the peak E 11 . For the 9L sample, the third peak (labelled as E 33 ) starts to emerge at around 1.22 eV, assigned to the third transition v 3 -c 3 . Again, E 33 shifts towards lower energies as the layer number increases. For the 13L and 15L samples, at least four peaks can be identified unambiguously. Apparently, each few-layer BP has a unique infrared spectrum and is readily distinguishable from each other, even for relatively thick samples such as 15L ones. For thinner BP flakes (typical layer number o5), sample thickness can be accurately determined by the E 11 peak position. Although the E 11 peak shows very little difference for thicker samples, the higher energy transition peaks (E 22 , E 33 and E 44 ) show more discernable differences. Therefore, it's quite accurate to determine the BP thickness by infrared spectroscopy with atomic-level precision, given the clear sequence of subband transitions. Supplementary Table 1 summarizes the peak energies in Fig. 2 and previous theoretical values for the quasiparticle bandgap are also shown for comparison.
For N-layer BP, the conduction and valence bands both split into N 2D subbands due to the layer-layer interactions 40 . These additional peaks originate from optical transitions between higher index subbands, revealing a QW-like band structure 17 . Following a tight binding model to take into account the nearestneighbour layer-layer interactions (see Supplementary Note 1) 40 , at G point of the 2D Brillioun zone, the transition energy between valence subband and conduction subband with the same index n (n ¼ 1,2,3,...,N) is given by where E g0 is the bandgap of monolayer BP, and g c and g v characterize the nearest interlayer couplings for the conduction and valence bands, respectively. Equation (1) describes the series of resonance peaks in infrared spectra. As seen in Fig. 3a, all of the four branches (E 11 , E 22 , E 33 and E 44 ) of optical transitions can be well fitted simultaneously, with fitting parameters E g0 ¼ 2.12 eV, g c À g v ¼ 0.88 eV. For completeness, we also put a PL data point (hollow square) for a monolayer BP in the figure. With increasing layer number N, the energy spacing between two adjacent subbands within the conduction or valence bands decreases monotonically. In the bulk limit, these subbands are so closely spaced that they evolve into quasi-continuous bands at last. Therefore, the spectrum for bulk BP in Fig. 2k has no resolvable peaks. Our results reveal the crossover of band structure from 2D to three dimensions. For the 13L and 15L BP, the relationship between transition energy E nn and subband index n (or quantum number) is well described by the QW based formula 17 where a and b are fitting parameters, as shown in Fig. 3b. This is a direct consequence of the QW confinement and can also be derived from equation (1), as shown in the Supplementary Note 1. The spectra set in Fig. 2 is a valuable fingerprint database for few-layer BP. Spectra with dramatic deviation from those typically signal large structural change. For instance, as predicated by theory, the stacking order can affect the band structure of BP layers 16 . Indeed, we observed a few anomalous spectra and attributed them to BP layers with different stacking order (see Supplementary Fig. 4).
In single and few-layer BP, it has been shown theoretically that the optical transitions may arise from excitons 6,41 . In the meantime, a theoretical calculation without taking into account excitonic effect also gives insightful results 17 . In our description of the optical resonance energy, we do not take into account electron-electron and electron-hole (excitonic effect) interactions, and adopt a simple tight binding model. Interestingly, the model is very consistent with our results. However, this   does not preclude the possible many-body effects in few-layer BP. Previously, for other low-dimensional materials, such as one-dimensional carbon nanotubes 34 , 2D transition metal dichalcogenides 36,42 and graphene 43 , it was reported that single particle band structure renormalization due to electronelectron interactions, which typically increases the optical transition energies, are largely cancelled out by the exciton binding energy. As a result, the experimentally observed optical transition energy can often be accounted for within the single particle description. Very likely, this is also the case for few-layer BP.
Strain engineering of the electronic structures. The electronic structure of mono-and few-layer BP has been predicted to be very sensitive to strain 7,[21][22][23][24]44 . Recently, optical spectroscopy in the visible frequency range showed strong spectral variation for wrinkled multilayer BP, indicating an inhomogeneous strain field 45 . Raman spectroscopy has been used to study the lattice vibrations of few-layer BP under uniaxial strain and strong anisotropy has been demonstrated 46,47 . In our experiment, we transferred few-layer BP on flexible polyethylene terephthalate substrate and employed a two-point bending apparatus to induce controllable uniaxial strain (e), as illustrated in Fig. 4a and Supplementary Fig. 5. We monitored the electronic structure change of the strained few-layer BP by infrared spectroscopy. The strain setup is detailed in Supplementary Note 2. Before applying strain, the crystallographic orientation of the BP flake was identified by polarized infrared spectroscopy, so that the uniaxial strain can be applied to the sample in the desired directions (AC and ZZ directions). To avoid the complication from sample slippage, the applied strain was typically kept below 1% during the entire strain process. Within such moderate strain range, the stretch process is reversible and repeatable. Figure 4c shows the infrared extinction spectra of a representative 6L BP sample under varying tensile strains up to e ¼ 0.92%, with strain applied along the AC and ZZ directions respectively. For clarity, the spectra are vertically shifted. As the tensile strain increases gradually, the two characteristic peaks of E 11 and E 22 both blueshift monotonically. As we know, BP is a highly anisotropic material, especially its mechanical properties exhibit strong crystallographic orientation dependence. For example, monolayer BP has been predicted to contract in the out of plane direction under tensile strain along the AC direction (positive Poisson's ratio), whereas it expands under ZZ tensile strain (negative Poisson's ratio) 48 . In addition, the Young's modulus of monolayer BP in the ZZ direction is predicted to be 3.8 times larger than that in the AC direction 10 . With these in mind, let us now examine the direction dependence of the uniaxial strain effect on the electronic structures. The peak positions in Fig. 4c are summarized in Fig. 4d as a function of tensile strain in both AC and ZZ directions. Surprisingly, little difference is observed for the two distinct strain directions. For each strain, the peak positions of the sample stretched in both directions are almost always the same. From the linearly fitted lines, we can extract that the blueshift rate for E 11 is 117 meV/% (AC direction) and 124 meV/% (ZZ direction), respectively, and 99 meV/% for E 22 in both directions. It is not a coincidence for this specific sample. We have performed similar measurements on multiple samples, no strain direction dependence beyond experimental uncertainty has been observed, in sharp contrast to the giant anisotropic Raman response to uniaxial strain 46,47 .
Within the strain range from e ¼ 0-0.92%, the optical bandgap (E 11 ) can be continuously tuned from 0.54 to 0.65 eV for this 6L sample. Such large tunability is highly desirable for the application in high-sensitivity strain sensors. From the extracted peak shifting rates, we conclude that 1% tensile strain leads to a 23% increase of the bandgap for a 6L BP. The fractional change of the bandgap is even greater for thicker BP samples, given that the bandgap is smaller and the shift rate is similar for all layer thickness, as discussed later on. These observations are consistent with previous first-principle calculations for monolayer and bulk BP [21][22][23][24] .
The observed electronic structure evolution can be understood within the tight-binding framework 25,39,44,49 . As strain is applied in the basal plane, it mostly affects the in-plane bonding. As a consequence, the strain effect (shift rates of the optical transitions) will be similar for mono-and few-layer BP, with little dependence on the layer thickness. Indeed, we observed a similar shift rate as that of the 6L BP sample for the E 11 transition of a 3L BP sample, as shown in Supplementary Fig. 6. In fact, such phenomenon can be inferred from equation (1). In-plane strain has strong effect on E g0 and smaller effect on g c À g v , which renders all optical transitions E N nn similar strain dependence. Therefore, we only need to consider the strain effect on monolayer BP, whose band structure around the band-edge can be quite satisfactorily captured by two hopping parameters t  in E g0 . The increase of the bandgap under uniaxial tensile strain for BP layers is in sharp contrast to mono-and few-layer MoS 2 , whose bandgap shrinks 51 .
More interestingly, the strain along the ZZ and AC directions have quantitatively the same effect on the optical transitions. This is very counter-intuitive, given the fact that BP is so anisotropic in terms of almost every property, in particular for the mechanical properties. To better understand this behaviour, we have to examine the strain effect more quantitatively within the tight-binding framework. With the hopping parameters obtained from DFT calculations 39 (at zero strain, t == 1 ¼ À 1:22 eV, t == 2 ¼ 3:665 eV) and, as a common practice, assume that they are inversely proportional to r 2 (r is the bond length), one can directly derive the bandgap dependence on the strain for monolayer BP 44,49 : with the unit of eV. In the equation, e x , e y and e z is the strain along x, y and z direction, respectively. It should be noted that for monolayer BP, e z reflects the change of the thickness of the single puckered sheet, that is, the vertical distance between atoms on the top and atoms at the bottom of the puckered sheet. In fact, in another 2D material-monolayer graphene, e z has no meaning since all atoms are on the same plane. According to equation (2), it is obvious that the pure ZZ strain (e y ) is more effective to tune the bandgap than the pure AC strain (e x ), which is consistent with our expectation. However, experimentally, due to the Poisson effect of the polyethylene terephthalate substrate, the strain is not truly uniaxial and deformations in other directions exist as well. BP layers stick to the underlying substrate and they deform in the same way as the substrate in the x-y plane. More specifically, the BP layer stretched in ZZ direction will slightly shrink in AC direction and vice versa. However, as the Poisson effect of the substrate has no orientation dependence, this effect cannot smear out the strain direction dependence shown in equation (2) and stretch along ZZ direction is still more efficient to change the electronic structure. On the other hand, in z direction, BP layers are free to shrink or expand and e z is non-zero. We attribute the observed lacking of orientation dependence for the strain effect to the change of the monolayer thickness under in-plane strain. More specifically, for a monolayer BP, the tensile strain along the AC direction leads to compression in z direction, which gives rise to additional bandgap enhancement (see equation (2), noticing that the coefficient of e z is negative). On the contrary, tensile strain along the ZZ direction results in expansion in the z direction due to a predicted negative Poisson's ratio (n zy o0) 48 , which partially cancels out the ZZ strain effect. More quantitatively, through equation (2), we find out that a Poisson's ratio difference (n zx À n zy ) of B0.12 in z direction will smear out the difference of the effect for strain along the ZZ and AC directions. This is consistent with DFT calculations, which have obtained n zx and n zy , with their difference ranging from 0.07 to 0. 3 (refs 24,48). It should be noted that the different Poisson's ratios in z direction for strain along the ZZ and AC directions are direct consequences of the hinge-like structure of a BP layer 48 .
Now we see that equation (2) semi-quantitatively describes the strain dependence of the electronic transitions in mono-and few-layer BP. However, the derived shift rate is smaller than what we observed experimentally. This discrepancy requires more effort in exactly determining the hopping parameters t == 1 and t == 2 , and reassessment of the assumption that tp1/r 2 . If we keep t == 1 and t == 2 values provided by Rudenko et al. 39 and assume that tp1/r b , we find out that the index bB6 can quantitatively account for the measured shift rates 25 .

Discussion
In addition to the main peaks in Fig. 2, there are weak absorption peaks between E 11 and E 22 transitions, indicated by asterisks (*). For 4,5,6-layer BP, the weak feature shows one bump, whereas for 7,8,9-layer case, it splits into two peaks, as indicated by double asterisks. We assign these weak peaks to the hybrid transitions from the valence subband to the conduction subband with different quantum numbers: v 1 to c 2 (E 12 transition) and v 2 to c 1 (E 21 transition), as illustrated in the right panel of Fig. 1c. From the energy diagram in Fig. 1c, it is straightforward to have (E 11 þ E 22 )/2 ¼ (E 12 þ E 21 )/2, which means that the average peak position of the E 12 and E 21 transitions is right in the middle of the E 11 and E 22 transitions. Indeed, this is exactly what we observed. Figure 5a plots the average peak position of E 12 and E 21 as a function of (E 11 þ E 22 )/2. For the not well-split peaks below 6L, (E 12 þ E 21 )/2 is just the peak position of the single bump. In addition to the data points extracted from the spectra in Fig. 2, other data points for additional 7L or 8L samples are also shown. The solid line is a straight line with (E 11 þ E 22 )/2 ¼ (E 12 þ E 21 )/2. We see that data points all follow the line, which is consistent with hybrid transitions.
The origin of these hybrid transitions deserve some discussion. Ideally, these hybrid transitions are forbidden and the optical oscillator strength should be zero for a symmetric QW with a normal-incident infrared beam 38 . We believe that the unintentional doping from the substrate and/or air can break the symmetry of the BP QWs and relax the selection rules 38,52 . This is consistent with a theoretical study, which shows doping indeed can activate such transitions 18 . To verify this conjecture, we intentionally doped a 9L sample with nitric acid vapour. Such doping scheme has been proved very efficient to increase hole carrier density in graphene 53 . Figure 5b shows the extinction spectra before and after doping. Clearly, the hybrid transition intensity becomes almost two times as strong after treatment, as shown in Fig. 5c. Meanwhile, E 11 slightly redshifts and its oscillator strength decreases dramatically, which can be understood as due to reduced wavefunction overlap for valence and conduction states 18 . These behaviours are all consistent with doping effect 18 . This controlled doping experiment suggests that doping can relax the optical transition selection rules and make forbidden transitions possible. It is natural to ask whether sample degradation makes hybrid transitions observable, as BP flakes are sensitive to oxygen and water in air and degrade at ambient condition 54 . To examine the degradation effect, we monitored the spectra for samples as a function of time. We find that the degradation is a slow process and freshly cleaved few-layer samples (L43) can keep the same spectra for several hours. However, thinner samples (Lo3) degrade much faster and spectra can change in 2 h after cleavage. The spectra shown in Fig. 2 were typically measured within 1 h after cleavage and we believe the degradation effect on the spectra is minimal. Supplementary Fig. 7 presents two spectra for the 4L sample right after cleavage and after 1 day. A blueshift of the main peak E 11 can be observed, presumably due to defects introduced by the degradation 55 . In addition, the hybrid peak (marked with *) disappears, suggesting that degradation is not the origin of the observed hybrid transitions. It should be noted that for some cases, even freshly cleaved BP samples do not show the hybrid peaks, presumably due to a low doping concentration.
Although with small oscillator strength, the hybrid transitions are very informative. With both the main peaks and hybrid peaks showing up, the energy levels of v 1 , v 2 , c 1 , c 2 (labelled as E v1 , E v2 , E c1 , E c2 ) with respect to a common reference energy can in principle all be determined. For 4,5,6-L BPs, the splitting between E 12 and E 21 are tiny and they are nearly degenerate, which indicates that the energy spacing between v 1 and v 2 (E v1 À E v2 ), c 1 and c 2 (E c2 À E c1 ) are almost identical, and electron and hole are nearly symmetric in z direction 17 . For those BP layers with split hybrid transitions (L46), the energy spacings E v1 À E v2 and E c2 À E c1 are different, which indicates an electron-hole asymmetry in the z direction. Previous cyclotron resonance experiments for bulk BP showed that the effective mass in z direction for electrons and holes are very different, with hole mass almost two times of electron mass 28 . This is consistent with our observations that the BP layers with larger thickness have split hybrid peaks and shows more pronounced electron-hole asymmetry. The reason why thin BP layers (Lo6) have almost degenerate E 12 and E 21 transitions merits further theoretical and experimental investigations. Meanwhile, the linewidths and relative intensities of E 12 and E 21 transitions, as well as their doping dependence deserve careful examinations as well. These studies will give valuable information concerning the band structure, especially electron-hole asymmetry in few-layer BP samples.

Methods
Sample preparation. Few-layer BP samples were prepared by a polydimethylsiloxane assisted mechanical exfoliation method. In brief, BP flakes were first cleaved on a polydimethylsiloxane substrate, then transferred onto a quartz substrate once thin flakes were identified under microscope. After the transfer, sample thickness was determined using a Nikon inverted microscope (Eclipse Ti-U), in combination with infrared spectroscopy, as mentioned in the main text.
Polarized infrared spectroscopy. Polarization-resolved infrared spectroscopy was performed using a Bruker FTIR spectrometer (Vertex 70v) integrated with a Hyperion 2000 microscope. A combination of tungsten halogen lamp and globar was used as light sources to cover the wide energy range from mid-to near-infrared (0.1-1.36 eV). The incident light was focused on BP samples with a Â 15 infrared objective, the polarization was controlled by a broadband ZnSe grid polarizer. Infrared radiation was collected by a liquid nitrogen cooled MCT detector. In general, aperture size was set to be larger than the sample size, to improve signal to noise. All the measurements were conducted at room temperature in ambient conditions.
Chemical doping of few-layer BP samples. The intentional doping for the 9L sample shown in Fig. 4 was done through HNO 3 vapour exposure. The sample was exposed for 5 s. Typically, longer time of exposure will gradually damage the sample.
Data availability. The data that support the findings of this study are available from the corresponding author upon request.