s41699-023-00384-2.pdf

ARTICLE OPEN Unconventional conductivity increase in multilayer black phosphorus Maria Kole śnik-Gray1, Laura Meingast 1, Martin Siebert 1, Tim Unbehaun 1, Tobias Huf 1, Günter Ellrott 1, Gonzalo Abellán 2,3, Stefan Wild 2, Vicent Lloret 2, Udo Mundloch 2, Julian Schwarz 4, Michael Niebauer 4, Maximilian Szabo 4, Mathias Rommel 5, Andreas Hutzler 6, Frank Hauke 2, Andreas Hirsch 2 and Vojislav Krsti ć 1 ✉ Multilayers of so-called 2D van der Waals materials have gained considerable attention as active components of next-generation electronic and optoelectronic technologies, with semiconducting black phosphorus (BP) regarded as one of the most promising systems. The applicability and performance limits of BP in both stand-alone and heterostructure-based multilayer devices are determined by individual ﬂake charge transport properties, which synergistically depend on the number of layers and the strength of interlayer coupling between those. In this work, we study the DC electrical transport characteristics of high-quality BP ﬁeld-effect devices within a wide range of ﬂake thicknesses at room temperature. The experimental data show a non-trivial increase in conductivity and hole density with a reduced number of layers while maintaining constant ﬁeld-effect mobility due to the prevalence of electron –phonon scattering. Based on the solution of the 1D Schrödinger –Poisson equation, we ﬁnd that the observed phenomena are a direct consequence of non-negligible interlayer coupling, which in turn causes a local redistribution of free charge carriers towards the central layers. Our data show that due to the electrostatic conditions at the ﬂake surfaces, a naturally protected 2D hole gas can be encapsulated in ﬂakes as high as 10 nm, which preserves the bulk-like bandgap and effective carrier masses due to the electrostatic environment. npj 2D Materials and Applications (2023) 7:21 ; https://doi.org/10.1038/s41699-023-00384-2 INTRODUCTION Black phosphorus (BP) is a synthetically derived layered semi- conductor material1 representing a class of monoatomic so-called two-dimensional (2D) van der Waals materials 2,3. The individual layers which form BP — phosphorene— have a puckered hexago- nal structure4 and are bound together by both van der Waals and ionic forces 5,6. As-grown (undoped) BP crystals are hole con- ductors7,8 with distinctive anisotropy of electrical transport properties along the principal axis 8. The p-type character of the majority of charge carriers has been attributed to the presence of electron trapping sites within the crystal lattice 9,10. Monolayer phosphorene, in particular, has been identi ﬁed as one of the most promising 2D electronic materials 11–13. Its semiconducting properties, including a bandgap comparable to silicon14 and theoretically predicted hole mobility in the 10 3 cm2/ Vs range 15, is expected to signi ﬁcantly boost the electrical conduction and switching properties of associated phosphorene ﬁeld-effect transistor devices regardless of their actual (effective) carrier concentration or their spatial distribution within the system. On the other hand, few-layer components, which constitute the vast majority of experimental data reported in the literature, present a higher degree of structural — and thus electronic — complexity. In particular, the overlap of the electronic wavefunc- tions between individual phosphorene layers 5,6 is non-negligible, therefore, black phosphorus is not a pure van der Waals system. This is also re ﬂected by the three to ten times larger layer-to-layer binding energy compared to graphene and MoS 2, respectively (Supporting Information). Consequently, such non-negligible interaction also facilitates some charge carrier transfer between the neighbouring 2D sheets while simultaneously preserving the preferential motion of carriers within individual layers under applied in-plane external electric ﬁelds. Therefore, the question of how the thickness of a multi-layered BP system affects the spatial distribution of free charge carriers and, consequently, its overall electrical conduction properties becomes highly relevant, not only from the fundamental but also application viewpoint. In particular, the layer number dependent charge transport ef ﬁcacy will strongly impact the performance of BP-based photodetec- tors 16,17, photocatalysts 18, energy conversion and storage devices19,20, biomedical devices and sensors 21. In the present work, we demonstrate through a comprehensive experimental study that the room temperature electrical con- ductivity of BP outside of the band-gap opening layer range reveals an unexpected regime of unconventional increase as the number of layers is reduced. Considering BP in this layer range in electrical terms to be a simple small-sized bulk semiconductor, such an increase is contrary to the expectation — a continuous reduction of size should lead to a decrease of conductivity upon height (layer number) reduction. The here observed conductivity augmentation is intimately connected to an increase in the effective carrier density in the ﬂakes while charge-carrier mobility remains constant. Using a simple 3D-semiconductor description, we ﬁnd that the experimentally observed effects can be explained by local charge redistribution towards the central 1Department of Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Staudtstr. 7, 91058 Erlangen, Germany. 2Department of Chemistry and Pharmacy and Joint Institute of Advanced Materials and Processes (ZMP), Friedrich-Alexander University of Erlangen-Nürnberg (FAU), Nikolaus-Fiebiger-Str. 10, 91 058 Erlangen, Germany. 3Instituto de Ciencia Molecular, Universidad de Valencia, Catedra ́tico Jose ́Beltrán 2, 46980 Paterna, Spain. 4Electron Devices (LEB), Department of Electrical Engineering, Friedrich- Alexander-Universität Erlangen-Nürnberg (FAU), Cauerstraße 6, 91058 Erlangen, Germany. 5Fraunhofer Institute for Integrated Systems and Device Technology IISB, Schottkystraße 10, 91058 Erlangen, Germany. 6Forschungszentrum Jülich GmbH, Helmholtz Institute Erlangen-Nürnberg for Renewable Energy (IEK-11), Cauerstraße 1, 91058 Erlangen, Germany. ✉email: [email protected] www.nature.com/npj2dmaterials Published in partnership with FCT NOVA with the support of E-MRS 1234567890():,; layers due to surface electrostatic conditions. Our results also provide insight into resolving the in-literature experimentally observed mobility limit of ~300 cm 2/Vs in few-layer BP devices at room temperature 22–25. RESULTS AND DISCUSSION Sample contacting, Raman and AFM characterisation BP ﬂakes were exfoliated from bulk crystals on the substrate and immediately capped with PMMA, both under an inert atmosphere (cf. “Methods” section). Optical images of PMMA capped exemplary ﬂakes are shown in Fig. 1a. After identi ﬁcation of ﬂakes, Ti/Au electrodes were electron-beam-lithographically de ﬁned on top (cf. “Methods” section), resulting in devices as shown in Fig. 1b. Immediately after electrode de ﬁnition the devices were again covered by a protective PMMA layer. For further characterisation, Raman spectroscopy was carried out (cf. “Methods” section) on fully processed devices, revealing a spectrum as shown in Fig. 1c. The spectrum shows only the well-known intrinsic ﬁrst ðA 1 g; B2g and A2 gÞ and second order (ca. 790 –920 nm−1 range) Raman signatures 26,27, while there is no trace of lines associated with oxidation 28,29 or any other uncommon lines discernible. Combined, the presence of the two intrinsic Raman responses and the absence of any other signs indicative of oxidation or degradation shows that the BP processed into devices is free of any unintentional in ﬂuences. This is further supported by AFM characterisation (cf. “Meth- ods”) of our processed devices from which we also extracted the height of the samples (Fig. 1d). Regarding the latter, the height range of the electrically characterised ﬂakes ranged from ca. 6 nm up to 160 nm. That is, the thinnest ﬂakes characterised were still beyond the bandgap broadening limit 14. Also, our AFM imaging directly after electrical measurements (after PMMA removal) did not show any discernible blistering associated with BP ﬂake degradation (Fig. 1d inset), in agreement with the Raman ﬁndings. Speciﬁcally, the data analysis revealed root mean square (RMS) values consistent with the absence of any relevant surface oxidation or degradation in general (Supporting Information) in agreement with previous AFM studies 30. Layer-dependent charge transport properties Electrical measurements were carried out in a four-point con ﬁg- uration at room temperature and under DC conditions. The highly doped silicon substrate (300 nm SiO2 gate dielectric on top) served as a global back-gate for transfer characteristic measurements. In- plane-orientation selection of the ﬂakes was discarded as experi- mental evidence showed that the conductivity and mobility of BP ﬂakes in armchairs and in zigzag direction differ at room temperature by less than a factor of 1.6 11. All devices showed linear two-point current–voltage (I–V) characteristics demonstrating low contact resistance and ohmicity of the electrode/BP interface (cf. Supporting Information). In all devices studied, unipolar p-type ﬁeld switching behaviour was observed, with weak on/off switching for gate bias exceeding ±60 V (Fig. 2a). The hole mobilities, µ, extracted from the linear part of the subthreshold slope were found to be largely independent of the ﬂake height (as corroborated by Hall measurements, see Supporting Information), and to lie around 200 cm 2/Vs (Fig. 2b), which is in agreement with the majority of earlier experimental studies 22–25. Fig. 1 Mechanical exfoliation and electrical contacting of multilayer BP. Exemplary optical micrographs of a BP ﬂakes transferred onto the SiO2-coated silicon substrate. Inset: single BP ﬂake (scale bar 10 μm). b A ﬂake after electrical contact with Ti/Au electrodes in a 4-point conﬁguration. Large electrodes are current-carrying leads, side electrodes are 4-point voltage probes (cf. “Methods”). c Raman spectrum of two ﬂakes of different thicknesses after device processing. Only the common BP signature lines Ag 1, B2g and Ag 2 and the intrinsic BP second- order Raman response (ca. 790 –920 nm−1 range) are observed (Si line from the substrate as expected at 520 cm −1). No trace of other lines associated with oxidation or degradation is present. The expected positions of theoretically predicted prominent oxidation lines are indicated by the arrows 28,29. d Flake height proﬁle extracted from an atomic force micrograph (inset) after electrical characterisation. The extracted ﬂake height is (15.1 ± 1.8) nm corresponding to (29 ± 3) phosphorene layers. M. Kole śnik-Gray et al. 2 npj 2D Materials and Applications (2023) 21 Published in partnership with FCT NOVA with the support of E-MRS 1234567890():,; Figure 2c shows the electrical conductivity, σ, as a function of ﬂake height h. In the case of thicker ﬂakes, the electrical conductivity is relatively constant and of the order of a few tens of S/cm. Notably, a marked increase in σ is observed for ﬂakes with thicknesses below ca. 20 nm, reaching values in the range of 10 2 S/cm. This result is non-trivial, as in the conventional nanoscaled, but still bulk, semiconductor case; one would expect a decrease in conductivity towards smaller layer numbers due to the mere reduction of the cross-sectional area. Taking the experimentally measured µ values into account, we estimated the free mobile charge carrier (hole) density nC (Fig. 2d). It needs to be emphasised here that the hole densities of about 1×1 0 18 cm−3 for our thicker ﬂakes are consistent with recent experimental work on the density of p-type vacancies (about 60–63 ppm) in high-quality black phosphorus crystals 10. For a BP monolayer with an atomic density of 2.53 × 10 15 cm−2, it can be therefore estimated that there is only one acceptor-inducing electron trap per several tens of thousands of P atoms in the crystal lattice. This demonstrates the high quality of the BP material used in our study. Remarkably, we found that n C— being the effective density of mobile charge carriers contributing to the electrical conduction – also changes with the number of layers (height). Speci ﬁcally, nC ﬁrst remains rather constant for thick(er) ﬂakes, while from about 20 nm height on, it displays a continuous increase, reaching values in the range from 5 to 8 × 10 18 cm−3 for heights in the 10 nm range. Intriguingly, such volume carrier densities in this height range are commonly observed. For instance, a recent publication by Yan et al. 28 reports on experimental sheet free carrier densities of 5 × 1012 cm−2 in a 10 nm high ﬂake; in a similar height range of 8–15 nm, Xia et al. 11 reports on an average sheet carrier density of about 6.7 × 10 12 cm−2. Strikingly, these values result in a volume- free carrier density of 5 × 10 18 cm−3 and 4.5 –8.4 × 1018 cm−3, respectively, which falls quantitatively into our experimental values in the same height range. While our experimentally determined hole densities are consistent with other experimental studies, the continuous increase with reduction of height (equivalently layer number) is non-trivial. The experimental data also reveal that σ is directly propor- tional to n C (inset Fig. 2d), which allows estimating the average charge carrier mobility of (215 ± 49) cm 2/Vs in all devices. This indicates that a dominant scattering mechanism limiting the motion of the charge carriers exists, which is independent of the layer number and charge carrier density. In fact, recent theoretical 31 and experimental work 32 has shown that the hole mobility in BP ﬂakes is limited exclusively by electron-phonon scattering at room temperature. Speci ﬁcally, the theoretically expected mobility in monolayer phosphorene at room tem- perature is about 250 cm 2/Vs regardless of the carrier density. It should be noted that any ionise d charged impurity scattering becomes only dominant at temperatures below 100 K and is at room temperature negligible 32. Comparison of this value with our experimental (average) mobil ity suggests, therefore, that indeed, the main scattering contribution of the charge carriers is in-plane electron-phonon scattering; however, also that addi- tional scattering channels exist. These are naturally occurring Fig. 2 Height-dependent electronic properties of BP ﬂakes at room temperature. a Representative four-point ﬁeld-effect characteristic of a BP ﬂake showing unipolar p-type behaviour and weak on/off switch. The hole mobility was extracted from the linear part of the subthreshold slope (denoted in the graph by a dashed grey line). b Charge carrier mobility values are estimated based on the ﬁeld-effect characterisation. The mobility values do not show a height-dependent trend and show a few variations around 200 cm 2/Vs, which fall within the range of the weakly developed in-plane anisotropy in charge transport at room temperature 11. c Electrical conductivity and d charge carrier density show a steep increase for ﬂakes with a height below ca. 20 nm (dashed lines as a guide to the eye), whereas, in thicker ﬂakes, the values remain approximately constant. Inset d Linear dependence between conductivity and carrier density implies a height-independent average charge carrier mobility of ca. 215 cm 2/Vs. Error bars result from data acquisition and error propagation in data evaluation. M. Kole śnik-Gray et al. 3 Published in partnership with FCT NOVA with the support of E-MRS npj 2D Materials and Applications (2023) 21 defects, electronic trapping sites as well as layer-to-layer hopping of free charge carriers. Spatial distribution of charge carriers To elucidate the origin of the observed conductivity and carrier density increase, we ﬁrst note that due to the presence of a ﬁnite charge density between the individual layers 5,6, the electrical current carrying charges can move both within a single layer and in between neighbouring ones. In a simpli ﬁed picture, such an electronic system can therefore be treated as a 3D semiconductor with some degree of in-plane preference. As schematically shown in Fig. 3a, the BP devices in our experiment are bound on each side by insulating SiO 2 and PMMA, which act both as (rigid) insulating barriers for charge transfer. In nanoscaled systems, such barriers are well known to induce valence and conduction band bending and produce sub-surface space-charge (carrier depletion and accumulation) regions in semiconductor materials 33. Consequently, in the case of BP multilayers, the few layers close to the surfaces become depleted of charge carriers and thus do not participate in the overall charge transport 34. Since typical space- charge regions extend into standard semiconductor bodies by up to a few nm, the impact of surface electrostatics is non-negligible in the case of multilayer BP and, therefore, must be taken into account. In view of the above complexity, we carried out the self- consistent evaluation of the one-dimensional Schrödinger-Poisson equation 35 to estimate the extent of band bending at the interfaces of BP with SiO 2 and PMMA perpendicular to the layers (Supplementary Information). Within our approach, we excluded the contributions from the contact interfaces and ﬂake edges. The one-dimensional descrip- tion used here was suitable as the lateral dimensions of our ﬂakes are by orders of magnitude larger ( μm) than our experimental height range (nm). The simulations were carried out for layer numbers ranging from relatively thick, bulk-like >100-layer stacks down to ﬂakes consisting of only 10 layers ( ca. 6 nm height — thinnest ﬂakes measured in our study). We also note that we are here outside of the band-gap opening regime for BP. For hole density, we took the experimentally observed average carrier density in thick ﬂakes of about 2 × 10 18 cm−3 in the constant conductivity regime (details on parameters for the model in the Supporting Information). To visualise the impact of space charge formation on the free carrier landscape, we also estimated the spatial distribution of holes across individual ﬂakes. Although the theoretical approach used here is simplistic, it provides not only a qualitative but also a good quantitative description of the multilayer BP system, as shown below. Figure 3b shows the energetic positions of the valence band (VB) with respect to the Fermi level ( EF = 0) and the extent of band bending at the SiO 2/BP and BP/PMMA interfaces for ﬂakes with varying layer numbers (for completeness, bulk-VB is also plotted as a limit of very thick ﬂakes as grey dashed line). The VB bends upwards at both interfaces, indicating the presence of charge accumulation close to the surfaces. As the ﬂake height is reduced, the VB bottom (corresponding to the central parts of the ﬂakes) is shifted towards the Fermi level, and noﬂat-band (potential) landscape is observed within the ﬂake centre for ﬂakes with ca. 35–45 layers, which falls close to our experimentally observed onset Fig. 3 Space charge region formation and local redistribution of majority charge carriers in multilayer BP close to the ﬂake surfaces. a Schematic diagram depicting an electrically contacted BP multilayer on a SiO 2 substrate capped with a PMMA protective layer. b Valence band energies within BP with layer number varied from 30 to 10. The dashed grey line denotes the VB position expected in bulk BP as a reference. For clarity, the ﬂake heights are normalised as indicated by the double arrow. c Layer number dependence of the VB bottom energies (orange) and associated with it the augmentation of local carrier density within the central layers of BP ﬂakes. d Examples of local distribution of carriers within BP ﬂakes of different layer numbers ( ﬂake heights normalised as in ( b) and e the corresponding colour maps depicting local carrier densities within BP ﬂakes as a function of height/layer number. As the layer number is reduced, the majority of charge carriers become con ﬁned to the innermost layers. Already for ﬂakes with heights of 10 nm, the free holes are electrically con ﬁned into the innermost 1–2 layers, with the carrier-depleted layers not contributing to the charge transport. In all calculations, the average carried density of 2 × 10 18 cm−3 was used. M. Kole śnik-Gray et al. 4 npj 2D Materials and Applications (2023) 21 Published in partnership with FCT NOVA with the support of E-MRS of conductivity increase. With further reduction of layer number (below ca. 25 layers), also the VB edges at both surfaces move closer to the Fermi level, and for ﬂakes with 10–11 layers (corresponding to the thinnest ﬂakes measured in our study), the VB appears to cross the EF close to the ﬂake surfaces. In Fig. 3c, the distance of the VB bottom (orange line) from EF is plotted as a function of layer number. The energetic separation of the VB from the Fermi level is well-known to be intimately connected with the charge carrier density 33. Based on that, we estimated the effective augmentation of the charge carrier density, ncentre/nbulk, in the central layers of the ﬂakes as compared to the corresponding bulk value. As shown in Fig. 3c, we would expect up to a threefold increase in the charge carrier density within the central layers for the thinnest ﬂakes. Remarkably, this is of the same order as observed in our experiments, which further supports the validity of our approach. For completeness, we had a closer look at the spatial distribution of charge carriers a mong the layers in individual ﬂakes. Figure 3d shows the (local) spatial distribution, n loc,o f free charge carriers (holes) in ﬂakes with varied layer numbers. In all cases, there is a distinct carrier depletion region beneath both surfaces of the ﬂakes. For thicker ﬂakes (>75 layers), two narrow charge carrier accumulation layers are formed below the depleted regions, as seen in the ‘rabbit ear ’curve shapes, while t h eh o l ed e n s i t yi nt h ec e n t r a ll a y e r sr e m a i n sc o n s t a n t .W i t h height reduction, however, the free charge carriers are squeezed towards the central layers due to the electrostatic potential at the surface, increasing the hole density in the centre. Eventually, for layer numbers below 18, a delta-like dopant pro ﬁle is observed (Fig. 3e). Since within the depletion region, the electrical current ﬂow is by large negligible 34, the ﬁndings imply that the charge carriers in the thin BP ﬂakes are spatially limited to the innermost layers only. This effect becomes more prominent for ﬂake heights smaller than about 20 nm, remarkably consistent with our experimental ﬁndings that σ and n suddenly increase (Fig. 2c, d, respectively). In fact, the calculations suggest that for ﬂake heights as large as 10 nm, already only one to two central phosphorene layers are able to carry the electrical current since all mobile charges are conﬁned to these. This ﬁnding is crucial, as it reveals that an electrostatically protected, quasi-2D hole gas con ﬁned to only one or two monolayers of phosphorene in the middle of a multilayer ﬂake can be achieved in ﬂakes with about 10 nm height. In particular, it is then this 2D hole gas which carries (predominantly) the electrical current in the entire system. It is worthwhile noting that the con ﬁnement of the mobile charge carriers to only one or two layers within the ﬂake does not necessarily imply that these layers experience a bandgap opening, such as isolated mono- layer phosphorene 14. The layers here are still embedded within a larger structure, that is, are subject to the inherent interlayer interactions present in BP. Therefore, the bandgap of these layers to which the 2D hole gas is con ﬁn e di se x p e c t e dt os t i l l have its bulk-value 36,37. Remarkably, this (rather simplistic) theoretical divagation reﬂects the experimentally observed phenomena. Indeed, the development of space-charge regions close to the ﬂake surface can only be possible due to the quasi-3D nature of the BP system. In other words, the combined experimental and model data unequivocally prove the non-pure van der Waals character of interlayer bonding in BP. Summarising, our ﬁndings demonstrate that due to non- negligible interlayer coupling in multilayer BP, an unprecedented increase in the electrical conductivity with decreasing number of layers is found, where the ﬂow of charge carriers is not entirely limited to the constituting monolayers. Consequently, the presence of an insulating environment induces the formation of charge depletion regions extending beneath the ﬂake surfaces. The observed experimentally drastic increase of electrical con- ductivity and effective charge carrier density with a reduction of layer number is a natural consequence of the electrostatic conﬁnement of free holes towards the central BP layers. Most notably, already in BP ﬂakes as thick as 10 nm, a naturally encapsulated quasi-2D hole gas can be expected, as the free charge carriers are limited to one or two central layers. This gas should be, by large, protected from external in ﬂuences and preserve the bulk electronic properties of BP, such as a narrow bandgap. This is further corroborated by the limited layer number and charge-carrier density independent hole mobility at room temperature. Its limitation in this temperature range is associated with electron-phonon scattering 28,29, and, in particular, its magnitude coincides well with, according to theoretical predic- tions29. Remarkably, we demonstrate that this phenomenon holds at room temperature, which opens the perspective to re-evaluate the use of multilayered BP in particular and comparable multi- layered semiconductor materials for on-chip integrated applica- tions in general. METHODS Sample preparation BP ﬂakes were mechanically cleaved from high-purity single crystals (2D Semiconductors), transferred onto SiO 2(300 nm)/Si(n++)s u b - strates and then capped with ca. 300 nm PMMA. This entire process was carried out under Ar atmosphere within a glove box. Individual ﬂakes covered under PMMA were identi ﬁed via optical means and electrically contacted by standard electron-beam lithography. As electrode material, electron-beam-evaporated Ti/Au (5 nm/50 nm thick) was used. Immediately after fabrication, the devices were capped again by a PMMA layer (ca. 300 nm) and baked out. Electrical measurements Electrical measurements were carried out with device-protecting PMMA layer-protected devices on a standard probe station. Driving and gate voltage was applied by K2450 source units. Voltage drops and currents were measured with the aid of K2000 digital multimeters and SR530 current –voltage ampli ﬁers. Both conductivity and mobility were measured by 4-point probing independently. Conductivities were measured with zero gate bias in an open-circuit con ﬁguration. All devices were electrically contacted in standard Hall-bar con ﬁguration, 2 oppositely oriented source/drain electrodes for the driving voltage and 4 voltage probes in-between (cf. Fig. 1b and inset Fig. 1d). The voltage probes were used to measure the 4-point voltage drop V 4p along the sample when current I (measured simultaneously) was ﬂowing through the device. The intrinsic conductivity then reads σ ¼ G4p L wh where G4p ¼ I V4p , L is the channel length (distance between probing electrodes), w the channel width (distance between probing electrodes) and h is the ﬂake height as measured by AFM. The intrinsic 4-point mobility was extracted from the linear (subthreshold) slope of the G4p vs. backgate voltage Vg curve using the expression μ ¼ ∂G4p ∂Vg /C12/C12 /C12 /C12/C12 /C12 L w tox εr ε0 where tox = 300 nm is the thickness of the SiO 2, εr = 3.9 the dielectric constant of SiO 2 and ε0 the vacuum permittivity. The carrier density n was obtained via σ = enμ where e is the electron charge. By intrinsic, we refer here to the electronic transport properties as inherent to the real material under study under the given experimental conditions. AFM measurements The height of each ﬂake was determined by a Park Systems NX10 AFM directly after the electrical measurements after stripping the protective PMMA layer. M. Kole śnik-Gray et al. 5 Published in partnership with FCT NOVA with the support of E-MRS npj 2D Materials and Applications (2023) 21 Raman microscopy Raman spectroscopic characterisation was carried out on a Horiba Jobin Yvon LabRAM Aramis confocal Raman microscope (excitation wavelength: 532 nm) with a laser spot size of <1 μma t room temperature. To not compromise the Raman measure- ments, the spectra were taken without the removal of the device- protecting PMMA layer. No relative position adjustment of the ﬂake and laser beam to account for the in-plane anisotropy of BP was undertaken. Schrödinger– Poisson calculations We used a FreeWare programme available online ( https:// www3.nd.edu/gsnider/) which solves the one-dimensional Poisson and Schrödinger equations self-consistently. We assumed the thickness of both SiO 2 and PMMA to be 300 nm and varied the thickness of the BP sandwiched between them incrementally layer-by-layer. The physical parameters of the system are listed in Table 1 in the Supplementary Information. DATA AVAILABILITY The data supporting the ﬁndings of this study are available within the paper and its Supplementary Information ﬁle. Data are also available from the corresponding author upon reasonable request. CODE AVAILABILITY The FreeWare programme we used is available fully online at https://www3.nd.edu/ gsnider/. Received: 19 July 2022; Accepted: 10 March 2023; REFERENCES

Xia, F., Wang, H., Hwang, J. C. M., Castro Neto, A. H. & Yang, L. Black phosphorus and its isoelectronic materials. Nat. Rev. Phys. 1, 306 –317 (2019).
McCreary, A., Kazakova, O., Jariwala, D. & Al Balushi, Z. Y. An outlook into the ﬂat land of 2D materials beyond graphene: synthesis, properties and device appli- cations. 2D Mater. 8, 013001 (2008).
Ajayan, P., Kim, P. & Banerjee, K. Two dimensional van der Waals materials. Phys. Today 69, 38 (2016).
Keyes, R. W. The electrical properties of black phosphorus. Phys. Rev. 92, 580–584 (1953).
Qiao, J., Kong, X., Hu, Z.-X., Yang, F. & Ji Wi High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus. Nat. Commun. 5, 4475 (2014).
Shulenburger, L., Baczewski, A. D., Zhu, Z., Guan, J. & Tomanek, D. The nature of the interlayer interaction in bulk and few-layer phosphorus. Nano Lett. 15, 8170–8175 (2015).
Warschauer, D. Electrical and optical properties of crystalline black phosphorus. J. Appl. Phys. 34, 1853 (1963).
Tao, J. et al. Mechanical and electrical anisotropy of few-layer black phosphorus. ACS Nano 9, 11362 –11370 (2015).
Xiang, D. et al. Surface transfer do ping induced effective modulation on ambipolar characteristics of few-layer black phosphorus. Nat. Commun. 6, 6485 (2015).
Rif ﬂe, J. V. et al. Impact of vacancies on electronic properties of black phosphorus probed by STM. J. Appl. Phys. 123, 044301 (2018).
Xia, F., Wang, H. & Jia, Y. Rediscovering black phosphorus as an anisotropic layered material for optoel ectronics and electronics. Nat. Commun. 5, 4458 (2014).
Liu, H., Du, Y., Deng, Y. & Ye, P. D. Semiconducting black phosphorus: synthesis, transport properties and electronic applications. Chem. Soc. Rev. 44, 2732 –2743 (2015).
Carvalho et al. Phosphorene: from theory to applications. Nat. Rev. Mater. 1, 16061 (2016).
Xiao, J. et al. Theoretical predictions on the electronic structure and charge carrier mobility in 2D Phosphorus sheets. Sci. Rep. 5, 09961 (2015).
Long, G. et al. Achieving ultrahigh carrier mobility in two-dimensional hole gas of black phosphorus. Nano Lett. 16, 7768 –7773 (2016).
Yi, Y., Sun, Z., Li, J., Chu, P. K. & Yu, X.-F. Optical and optoelectronic properties of black phosphorus and recent photonic and optoelectronic applications. Small Methods 3, 1900165 (2019).
Huang, L. et al. Waveguide-integrated black phosphorus photodetector for mid- infrared applications. ACS Nano 13, 913 –921 (2019).
Li, B. et al. Black phosphorus, a rising star 2D nanomaterial in the post-graphene era: synthesis, properties, modi ﬁcations, and photocatalysis applications. Small 15, 1804565 (2019).
Pang, J. et al. Applications of phosphorene and black phosphorus in energy conversion and storage devices. Adv. Energy Mater. 8, 1702093 (2018).
Qiu, M. et al. Current progress in black phosphorus materials and their application in electrochemical energy storage. Nanoscale 9, 13384 –13403 (2017).
Luo, M., Fan, T., Zhou, Y., Zhang, H. & Mei, L. 2D black phosphorus-based bio- medical applications. Adv. Funct. Mater. 29, 1808306 (2019).
Wang, G. et al. Introduction of interfacial charges to black phosphorus for a family of planar devices. Nano Lett. 16, 6870 –6878 (2016).
Yasaei, P. et al. High-quality black phosphorus atomic layers by liquid-phase exfoliation. Adv. Mater. 27, 1887 –1892 (2015).
Liu, H. et al. Phosphorene: an unexplored 2D semiconductor with a high hole mobility. ACS Nano 8, 4033 –4041 (2014).
Kang, J. et al. Solvent exfoliation of electronic-grade two-dimensional black phosphorus. ACS Nano 9, 3596 –3604 (2015).
Kaneta, C., Katayama-Yoshida, H. & Morita, A. Lattice dynamics of black phos- phorus. Solid State Commun. 44, 613 –617 (1982).
Morita, A. Semiconducting black phosphorus. Appl. Phys. A 39, 227 –242 (1986).
Ribeiro, H. B., Pimenta, M. A. & de Matosa, C. J. S. Raman spectroscopy in black phosphorus. J. Raman Spectrosc. 49,7 6 –90 (2018).
Ziletti, A., Carvalho, A., Campbell, D. K., Coker, D. F. & Castro Neto, A. H. Oxygen defects in phosphorene. Phys. Rev. Lett. 114, 046801 (2015). 3 0 .K o e n i g ,S .P . ,D o g a n o v ,R .A . ,S c h m i d t ,H . ,C a s t r oN e t o ,A .H .&Ö z y i l m a z ,B . Electric ﬁeld effect in ultrathin black phosphorus. Appl. Phys. Lett. 104, 103106 (2014).
Rudenko, A. N., Brener, S. & Katsnelson, M. I. Intrinsic charge carrier mobility in single-layer black phosphorus. Phys. Rev. Lett. 116, 246401 (2016).
Yan, X., Wang, H. & Esqueda, I. S. Temperature-dependent transport in ultrathin black phosphorus ﬁeld-effect transistors. Nano Lett. 19, 482 –487 (2019).
Simpkins, B. S., Mastro, M. A., Eddy, C. R. Jr. & Pehrsson, P. E. Surface depletion effects in semiconducting nanowires. J. Appl. Phys. 103, 104313 (2008).
Calarco, R., Stoica, T., Brandt, O. & Geelhaar, L. Surface-induced effects in GaN nanowires. J. Mater. Res. 26, 2157 –2168 (2011).
Snider, G. L., Tan, I.-H. & Hu, E. L. Electron states in mesa-etched one-dimensional quantum well wires. J. Appl. Phys. 68, 2849 (1990).
Tian, T. et al. Electronic polarizability as the fundamental variable in the dielectric properties of two-dimensional materials. Nano Lett. 20, 841 –851 (2020).
Tan, I.-H., Snider, G. L. & Hu, E. L. A self consistent solution of Schroedinger- Poisson equations using a nonuniform mesh. J. Appl. Phys. 68, 4071 (1990). ACKNOWLEDGEMENTS V.K. and A.H. acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) through the project “Elektronische Bauelemente auf Basis des 2D-Materials schwarzer Phosphor — lagenanzahlabhängige Eigenschaften ” (grant No. KR 2264/7- 1 and HU 2827/2-1, respectively). The work was supported by the European Union (ERC Advanced Grant 742145 B-PhosphoChem to A. Hirsch, ERC-2018-StG 804110- 2D-PnictoChem to G.A.). G.A. acknowledges the DFG (FLAG-ERA AB694/2-1) the Spanish MICINN (PID2019-111742GA-I00, “Mari ́ad eM a e z t u” Programme for Units of Excellence in R&D CEX2019-000 919-M) and the Generalitat Valenciana (CIDEGENT/2018/001). AUTHOR CONTRIBUTIONS V.K. designed and supervised the study. M.K.G., L.M., T.U. and M.S. carried out device preparation, electrical characterisation, and AFM measurements. M.S. performed the Raman measurements and G.E. the numerical simulations. T.H., A.H., J.S., M.N., M.Sz. and M.K.G. did AFM measurements and analysis. G.A., S.W., V.L., U.M., M.R., F.H. and A.Hi. carried out the ﬂake preparation. M.K.G. and V.K. analysed the data. All authors contributed equally to the paper write-up. FUNDING Open Access funding enabled and organized by Projekt DEAL. M. Kole śnik-Gray et al. 6 npj 2D Materials and Applications (2023) 21 Published in partnership with FCT NOVA with the support of E-MRS COMPETING INTERESTS The authors declare no competing interests. INCLUSION AND ETHICS This research study has included local researchers throughout the research process with well-de ﬁned responsibilities. The research presented is locally relevant and is due to a profound collaboration of the local partners. ADDITIONAL INFORMATION Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41699-023-00384-2. Correspondence and requests for materials should be addressed to Vojislav Krsti ć. Reprints and permission information is available at http://www.nature.com/ reprints Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional af ﬁliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article ’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http:// creativecommons.org/licenses/by/4.0/. © The Author(s) 2023 M. Kole śnik-Gray et al. 7 Published in partnership with FCT NOVA with the support of E-MRS npj 2D Materials and Applications (2023) 21