Scalable Patterning of Encapsulated Black Phosphorus

Atomically thin black phosphorus (BP) has attracted considerable interest due to its unique properties, such as an infrared band gap that depends on the number of layers and excellent electronic transport characteristics. This material is known to be sensitive to light and oxygen and degrades in air unless protected with an encapsulation barrier, limiting its exploitation in electrical devices. We present a new scalable technique for nanopatterning few layered BP by direct electron beam exposure of encapsulated crystals, achieving a spatial resolution down to 6 nm. By encapsulating the BP with single layer graphene or hexagonal boron nitride (hBN), we show that a focused electron probe can be used to produce controllable local oxidation of BP through nanometre size defects created in the encapsulation layer by the electron impact. We have tested the approach in the scanning transmission electron microscope (STEM) and using industry standard electron beam lithography (EBL). Etched regions of the BP are stabilized by a thin passivation layer and demonstrate typical insulating behavior as measured at 300 and 4.3 K. This new scalable approach to nanopatterning of thin air sensitive crystals has the potential to facilitate their wider use for a variety of sensing and electronics applications.


Analysis of Electron Diffraction Data
: a) Diffraction pattern from Figure 1 in the main text without overlays, and (b), (c) with overlays indicating spot assignments for graphene layers and black phosphorus respectively. The intensity ratio of the (110)/(200) spots is equal to ~0.1, corresponding to a thickness of 5 layers. 1 The scale bar represents 2 nm -1 . Figure S3(a) shows the diffraction pattern from Figure 1(f) in the main text without overlays. Before analysis, the background (tails of directly transmitted beam) was removed using the following method. Firstly, the center of the diffraction pattern was identified using the geometry of the observed spots. The intensity observed in each concentric circle around this central point was then integrated to produce a radial intensity profile. A least squares fit was then used to identify a superposition of two power law functions which described the parts of the radial profile where peaks due to diffraction spots were not present, which was then subtracted from the original image. Individual spot intensities were then obtained by fitting a 2D Gaussian-Lorentzian function to the area immediately around each spot. The obtained intensities for the diffraction pattern from Figure  1(f) in the main text are shown in Figure S3(b) and (c), for the spots associated with the encapsulating graphene and the black phosphorus crystals respectively. It has previously been demonstrated that the intensity ratio between the (110) and (200) reflections can be used to determine the thickness of thin, odd layer numbered black phosphorus flakes 1-2 , assuming that the number of layers is constant over the sampled region. In this case, the right (200) reflection is partially occluded by the beam stop and the intensity could not be extracted. Although the nonsymmetric intensities suggest that there is a small degree of sample tilt the mean of the four (110) / (200) intensity ratios is 0.11, a good fit to the expected ratio for 5 layer thickness (0.10 1 ), and in agreement with the thickness estimate of 4-6 layers obtained from the same flake regions optical contrast on the thin SiO x prior to transfer to the TEM.      Within the holes, the diffraction pattern shows only the diffraction spots of encapsulating graphene sheets (two sets rotated by ~21º). Outside of this region, the black phosphorus is still protected from oxidation and is highly crystalline. Approximate total dose for low resolution EDXS mapping (without damaging graphene) was 1.5 x 10 3 e.nm -2 at a dose rate of 3.2 x 10 2 e.nm -2 .s -1 . Approximate total dose for high resolution imaging in the central area (did damage graphene) was 1.2 x 10 5 e.nm -2 at a rate of 9.9 x 10 3 e.nm -2 .s -1 . Scale bar of (a) and (b) overview is 1 µm, (a) and (b) EDXS elemental map is 500 nm, (c) overview is 500 nm, and (c) SAED patterns is 5 nm -1 .  Figure 2c in the main text. For this specimen the sample was imaged after each cut was etched, leading to damage to the graphene and a subsequent degradation of the unetched regions of the flake in S6b. This has been avoided in Figure 2 by using low dose imaging approaches to minimize unwanted electron beam exposure.

Defects and Porosity
Imaging of pristine graphene flakes, mechanically exfoliated from natural graphite, using high resolution imaging techniques such as TEM and scanning tunneling microscopy (STM) indicates that the intrinsic concentration of structural defects is extremely low, 4-5 although it is generally only practical to image relatively small areas at the high resolution necessary to pick up such defects, and it is typically challenging to discern whether observed defective sites were present in the original crystal or introduced by energy transmitted to the sample due to the imaging technique. Scanning nanothermometry of high quality mechanically exfoliated graphene 6 (encapsulated between hBN crystals) reveals very few structural defects over µm scale samples with concentrations of monovacancies or sp 3 bonded adatoms being just 2 x 10 5 cm -2 . In pristine flakes, observed defects typically consist of monovacancies, bond rotations (stone-wales defects 7 ) and occasionally divacancies 8 . While monovacancies are relatively mobile at high T (migration energy barrier ~ 1.3 eV), divacancies (~7 eV), Stone -Wales defects (~10eV) and larger defective structures are not able to migrate thermally, 7, 9 and hence larger structural vacancies are rarely observed in natural graphite. The energy barriers are sufficiently low that such vacancies are relatively mobile under electron illumination at 80 keV (a 80 keV electron can transfer <15.8 eV to a C atom). Clusters of point defects do not necessarily exhibit a larger pore size, due to atomic rearrangement of the surrounding atoms; tetravacancies are usually observed to rearrange into extended defect structures with a maximum pore size no larger than that of a divacancy. 10 The lower knock-on threshold energy of the boron atoms 7 in hBN means that B-vacancies are commonly observed as a result of TEM/STEM imaging 11 However, as for graphene, once a vacancy is formed, the under coordinated atoms surrounding the vacancy have significantly lowered sputtering energy barriers, and are preferentially removed. Consequently, as the point defect density increases under prolonged high intensity electron beam irradiation, point defects will coalesce to form nanopores. In graphene these are often circular while in hBN preferential removal of B tends to produce triangular pores at room temperature. 12 Density Functional Theory (DFT) studies predict the penetration barrier for He atoms through a small defect drops exponentially as the defect size increases, meaning that small vacancy cluster defects such as those which may be present in mechanically exfoliated graphene or hBN are effectively impermeable at room temperature. 13 , 14 , Experimental studies have borne this out, with areas of a few µm 2 of mechanically exfoliated graphene or hBN having been shown to be essentially impermeable to all gases including molecular He, 15 as well as oxygen and water vapor. [16][17][18] To a good approximation, the gas permeability of defective graphene can be explained by comparing the size of the relevant pore to the kinematic diameter of the relevant gas (2.65Å for H 2 O and 3.46Å for O 2 ). 19 As the pore size approaches the molecule size, permeability is modulated by the molecule-pore repulsive interaction and the resistance of the pore structure to bond stretching/ rotation 20 . Gas transport in this regime can therefore be strongly affected by the local structure and chemistry of the pore, and the presence of any adatoms or dangling bonds. Transport of water can be strongly modulated by hydrogen bonding to the terminal groups around the hole 21 , and the permeation energy barrier for water can be reduced in the presence of other water molecules due to additional water-water interactions. 22 The smallest defects observed in graphene structures which may allow partial ingress of air and water consist of a hole created by removing a single ring of 6 carbon atoms, with no terminal groups (eg. hydrogen) around the pore, leading to a pore diameter of approximately 0.5 nm. 23 However, theoretically predicted values for the actual permeance of 0.5 -1 nm diameter pores vary considerably and can be strongly modulated by the presence of surface atoms. 22,[24][25] Experimental verification of permeability at this length scale is also very challenging 26 but the wealth of literature of the subject suggests that the oxygen/phosphorus ingress/egress we observe in our experiments will require at least nanosized pores to be created locally in the graphene/hBN.

EEL SI Processing
In order to obtain compositional maps from multidimensional spectrum images, such as EELS maps, the relative strength of various component signals must be quantified across the spatial dimensions. This can be achieved using least-squares fitting to reference or model spectra, 27 if the number of free fitting parameters is sufficiently low and the signal to noise ratio (SNR) of the individual spectra is sufficiently high. Whilst we were able to use this approach to map the intensity of individual elements (P L-edges, C and O K-edge), we were unable to separate the overlapping BP and PO signals; instead we found that separating the contributory components using non-negative matrix factorization (NMF) 28 provided more distinct EEL component maps.
We can attempt to approximate a spectrum image X (2 spatial, 1 spectral dimensions) using a series of n components x i made up of n 'loadings' L and n 'factors' F such that: Here the factors/eigenvectors F have one spectral dimension and the loadings/eigenvalues L have the two navigation dimensions. NMF imposes the additional constraint that all loadings and factors must be non-negative, which improves the physical interpretability of the separated components. For example, consider an EELS map where the entire area contains species A, and all but a small area contains a small amount of species B. When the non-negative constraint is not imposed, the first extracted factor will tend to be a combined 'A and B' signal (covering the entire area), with the second representing a 'A and negative B' signal localized around the area where only species A is present, such that in this area the A+B and the A-B signals are combined to approximate an 'only A' signal. Assuming no additional variance in the image, NMF decomposition with n=2 would more usefully separate the signal into components representing species A and B separately. In reality, obtaining physically interpretable components can require careful choice of n, depending not only on the dimensionality of the signal (2 in this example, quantity of A, and quantity of B), but also on the SNR and the dimensionality of the noise present in the image (from various sources).
The core loss EEL spectral maps presented in figure 3 in the main text and figure S11 in the supplementary information were captured using a 0.25 eV dispersion. Low loss spectra were captured simultaneously using a lower exposure time. The spectrum images were processed using freeware python modules including HyperSpy 1.3, 29 a multidimensional data-analysis toolbox with significant EEL spectrum image processing capabilities. Firstly, they were divided into EFTEM 'slices' and a median filter (adjacent pixels) was applied to each image before reconstruction, to remove x-ray spikes (anomalously high intensity pixels). The spectra were then aligned with subpixel accuracy by fitting the zero loss peak in the low-loss spectra, to remove any spectral drift. To obtain the elemental intensity maps, the aligned spectrum image was de-noised using PCA filtering, [30][31] and then fit to a model consisting of 3 core loss structures (Hyperspy EELSCLEDGE models for P L, C K and O K edges and a power law for the decaying backgrounds between them). To obtain the relative BP and PO intensities, the aligned spectrum image was first cropped to the energy region containing only the P L edges, and the background was removed by fitting to a power law. NMF decomposition (sklearn.decomposition.NMF 32 packaged in HyperSpy) was then used to separate the remaining phosphorus edge signal into individual components. The choice of the number of components, or the 'output dimension' is nontrivial, and we aimed to choose the highest value of n where output components represent physically interpretable spatial variations in the signal rather than noise effects and/or experimental instabilities. For the spectrum image shown in figure 3 in the main text, an output dimension of 4 was chosen, generating the loadings and factors plotted in Figure S12. The individual components were then assigned as being related to the BP or PO signal by visual inspection. For Figure S12, the 1 st and 4 th components were assigned as BP and the 2 nd and 3 rd as PO. The reconstructed BP and PO signals were then generated by summing the relevant components (loadings x factors). For the spectrum image shown in figure S11, an output dimension of 3 was used, generating the components plotted in Figure S13. Here, the 1 st and 2 nd components were assigned as the BP signal and the 3 rd component was assigned as the PO signal. Figure S13: Extracted NMF loadings and factors used to reconstruct the separate BP and PO signals in Figure S11. Visually, the 1st and 2nd components were assigned as BP and the 3rd was assigned as PO. Now, we examine the effects of the choice of the output dimension n on the NMF decomposition results for the two spectrum images. For n=1-5, we show all n reconstructed components x i integrated over the spatial axes in Figure S14(a). Although this is superficially similar to plotting the factors F n , it quantitatively describes the amount of total signal variance provided by each component relative to each other (here factors are effectively 'weighted' by how much of the overall signal variance they represent). In Figure S14(b), we plot the total residual variance after reconstruction with all n components for many n component distributions where, in this case, n=1,2,3, … ,20. We note that this is not a 'traditional' scree plot which plots the cumulative residual variance after reconstruction with i components in one n component decomposition where i=1,2,3,  … ,n). We also plot the residual signal R in Figure S14(c), defined as the difference between the original signal and the sum of all n components (the 'reconstructed' signal) integrated over the spatial axes for each value of n.
Reconstruction with one component is able to account for 87.09% of the total variance in the spectrum image. The vast majority of the area of the map is 'thick' black phosphorus in this region, so merely fitting each pixel with the average spectra (which is very close to a bulk BP spectra) accounts for most of the variance in the image. As this makes the residuals signal and higher order components much smaller than the 1 st component, we choose to define an integration window immediately around the cut areas, as shown in Figure S15(a). All integrated component spectra, residuals spectra, and fractional variance plots in Figure S14 are calculated over this window rather than the wider image to enhance their visibility. The decomposed factors and loadings used were from the wider image, not this reduced area. However this area represents the vast majority of the spatial variation in the signal, and also that where noise effects are more significant (as the sample is effectively thinner, meaning less electron interaction, in the cut regions). Therefore, restricting our plots to this area makes the effects of changing the decomposition output dimension more apparent. The integrated residuals signal for the full image are dominated by background fitting errors/variations (the fitting region used corresponds to the flat region (E<127 eV) at the lower energy range of the above plots), whereas those for the restricted region are dominated by spatial variations (at low n) and stochastic noise (higher n). Figure S14: Effect of changing output dimension during NMF decomposition for the spectrum image shown in figures S11 and S13 in the supplementary material. In this figure, the residuals and components have been integrated over a reduced region immediately around the cuts (shown in Figure S15a) to increase the size of the residuals signal for higher output dimensions.

(a) shows all factors in an n-component decomposition of the spectrum image integrated over the window, for n=1-5. (b) Residual signal variance (integrated over signal and navigation dimensions (as a fraction of the total signal variance) after decomposition into n components and subsequent reconstruction with all components, as a function of the number of components n. A 'shoulder' is clearly observed between 3 and 4 components. (c) Residual signal variance integrated over spatial dimensions after decomposition into n total components and recomposition with all n, for n=1-10 components.
Reconstruction with two components results in an additional factor corresponding to a relative increase of the broad central peak (around 155 eV) relative to the L 2,3 and L 1 'edges' at ~130 and 190 eV respectively. This is likely due to (slight) changes in the ratio between the 'electron thickness' of the sample (t/λ) and the elemental phosphorus content due to, for example, the level of surface carbon remaining from the polymer transfer film. We did not use a deconvolution to remove the effects of plural scattering (which would have mostly eliminated any of the signal corresponding to the sample thickness/phosphorus content ratio) -the signal to noise level in the 5-80 eV region of the low loss EELS map was too poor. For such a thin sample, the ZLP peak is far more intense than that of the scattering spectrum, reducing the exposure time we were able to use without damaging the camera. Reconstruction with 3 components adds a third physically interpretable component corresponding to oxidised phosphorus. [33][34] Reconstruction with 4 components introduces a 'mixing' sawtooth-like effect between the 2 nd component and the new 4 th component (and with the 1 st component to a lesser degree), corresponding (when reconstructed) to a series of small peak shifts. As the spectrum image was recorded, the position of the spectra on the camera was slightly shifted due to drift. We were able to remove this signal drift by aligning the spectrum with the zero loss peak at the 0 eV position, for each pixel. However, there is a 'gain noise spectrum' present in the images which is dependent on the position on the camera rather than the electron energy. When we have aligned the spectra using the ZLP, this causes the gain noise spectra to appear to shift (in energy) across the image, which is accounted for in the decomposition by the relative strength of the 2 nd and 4 th components. The spatial distribution of the 4 th component closely follows the initial ZLP peak position prior to realignment. Residuals spectra for a 5 component decomposition are similar to that of the 4 component case with the absence of the sharp peak at the 132 eV onset. This component is distributed 'roughly' inversely across the image compared to the 4 th component, which is correlated to the initial ZLP position. We suspect that this component at least partly represents an error during the phosphorus edge background subtraction process, where the shift of the 'gain noise spectra' causes a smooth change in the best fit power law used to represent the background, which was subtracted from the spectra prior to decomposition. This results in a peak in the residuals spectra near the steepest part of the spectrum, and for most of the image this is the P onset around 130 eV.
New factors generated in the n=5-10 range tended to be splits of existing components (such as that seen between the 4 th and 2 nd for n=4 above) correlated with the ZLP peak shift, and even ZLP peak width, but they are not reproduced here for reasons of brevity. Higher order components typically consisted of spatially and spectrally uncorrelated noise signals.
The 'elbow' of the variance plot in Figure S14(b) obviously occurs at the 3 component decomposition, and this is often considered a good guide to the dimensionality of the signal. However, appropriate choice of output dimension n requires careful examination of the output components, and residuals signal spectra. Clearly, there is an appreciable difference between the residuals spectra for 1, 2, and 3 component decompositions, whereas there is very little change after n=4 (reductions in the residuals spectra after n=4 are distributed evenly across the spectrum, with the exception of the 130eV peak changing slightly for the n=5). We chose a value of n=3 to generate the factors used to plot the separate BP and PO signals in figure S13, as the n=4 decomposition and higher n decompositions produce factors corresponding to experimental instabilities and noise rather than physically meaningful spatial variations. We also verified that the residuals signal obtained after n=3 decomposition and reconstruction did not show significant spatial variation using a second decomposition, examined output loadings and factors up to n=20, and constructed total BP and PO images by visually assigning components (as in S10) for n=4-6 and verified that the final maps and representative spectra are quantitatively similar and qualitatively indistinguishable from those presented in Figure S11.  Figure S16 shows the equivalent plots as those in Figure S14, but for the EEL spectrum image reproduced in Figure 3 in the main text, with its NMF components shown in Figure S12. As in the previous example, the decomposition was performed on the wider image, but the spectra and variances plotted in Figure S16 are integrated over a small window around the cut rather than the entire image, to enhance the visibility of the smaller signals (higher order components, residuals) compared to the larger BP bulk signal, and increase the area where the signal spatial variance is more closely related to real spatial compositional variations (ie. the cut region) rather than background fitting errors (bulk BP area). The integration window is plotted in Figure S15(b).
In contrast to the previous image, a component representing PO appears for n=2, with an additional (split) PO component at n=3. n=4 decomposition results in a 4th component consisting of the onset peak around 130 eV, which appeared as the 5th component in Figure S14(a). The 'mixed' signal related to the movement of the 'gain noise' spectra which appeared between the 2nd and 4th (BP) components in the n=4 panel of Figure S14(a) appears here for n=5, largely between the 2nd and 5th components, although stochastic noise effects are relatively more dominant than the gain noise shift in this case due to a smaller and less linear total energy drift, so the characteristic 'sawtooth' pattern of peak shifts is not clearly visible. Figure S16: Effect of changing output dimension during NMF decomposition for the spectrum image shown in figure 3 in the main text and figure S12 in the supplementary material. In this figure, the residuals and components have been integrated over a reduced region immediately around the cut (shown in Figure S15b) to increase the size of the residuals signal for higher output dimensions. (b) Residual signal variance (integrated over signal and navigation dimensions (as a fraction of the total signal variance) after decomposition into n components and subsequent reconstruction with all components, as a function of the number of components n. A 'shoulder' is clearly observed between 3 and 4 components. (c) Residual signal variance integrated over spatial dimensions after decomposition into n total components and recomposition with all n, for n=1-10 components.
Again, new factors generated in the n=5-8 range tended to be splits of existing components correlated with the ZLP peak shift and change in the ZLP width, and higher order components (n=8-20) typically consisted of spatially and spectrally uncorrelated noise signals. As before, the clear 'elbow' of the scree plot in Figure S16(b) occurs at n=3, however the residuals spectra appears to show a visible change between n=3 and n=4, and none afterwards. As n=4 still generated physically interpretable components rather than noise effects, and accounted for slightly more of the total signal variance than the n=3 decomposition, we chose to use an output dimension of 4 to generate the final PO and BP signals shown in main text figure 3. As before, we also verified that the residuals signal obtained after n=4 decomposition and reconstruction did not show significant spatial variation using a second decomposition, examined the output loadings and factors up to n=20, and constructed total BP and PO images by visually assigning components (as in S9) for n=3-6 and verified that the final maps and representative spectra are quantitatively similar and qualitatively indistinguishable from those presented in Figure 3 in the main text. We note that the total reconstructed BP and PO maps are visually indistinguishable from those generated with n=3.