Mechanochemical conversion kinetics of red to black phosphorus and scaling parameters for high volume synthesis

Adopting black phosphorus (BP) as a material in electronic and optoelectronic device manufacturing requires the development and understanding of a large-scale synthesis technique. To that end, high-energy planetary ball milling is demonstrated as a scalable synthesis route, and the mechanisms and conversion kinetics of the BP phase transformation are investigated. During the milling process, media collisions rapidly compress amorphous red phosphorus (RP) into crystalline, orthorhombic BP flakes, resulting in a conversion yield of ≈90% for ≈5 g of bulk BP powder. Milling conversion kinetics, monitored via ex situ x-ray diffraction, manifest a sigmoidal behavior best described by the Avrami rate model with each impact of sufficient energy (>25 mJ) producing BP nuclei; the process appears to be limited by grain growth. Using a kinematic model for ball trajectories and impact energies, the optimum milling condition is determined to be an impact energy near ≈25 mJ and a milling dose near ≈100 kJ/gram. Photoexcitation of exfoliated BP flakes reveals emission in the near-infrared, indicating the formation of few-layer BP, a promising advance for optoelectronic device applications.


INTRODUCTION
The surge of research into two-dimensional materials after the discovery of graphene in 2004 has led to investigating the monolayer limit of layered materials that include transition metal dichalcogenides (TMD), hexagonal boron nitride (hBN), MXenes, layered elemental metals (Te, Sn, Ge, As, Se) 1 , black phosphorus (BP), and beyond. Orthorhombic BP, the high pressure and temperature phase of elemental phosphorus, consists of a puckered hexagonal sheet of covalently bonded atoms. In bulk form, the individual sheets are bound together by van der Waal forces; however, new properties emerge due to quantum size effects when exfoliated into the few-layer regime. A property of note is BP's direct bandgap, which changes from 0.3 eV in the bulk to ≈1.7 eV in the monolayer limit 2,3 . Additionally, BP shows promise in optoelectronic device applications 4-6 from its direct bandgap as BP emissions are between 700 and >4000 nm, which is the near-to mid-infrared wavelength range.
Nevertheless, a standing challenge for BP is its synthesis. Traditionally, generating BP necessitates either high pressure and elevated temperature as in the diamond anvil method first pioneered by Bridgman 7 to phase transform either white or red phosphorus (RP), or the use of a low-pressure mineralizer assisted chemical vapor transport (CVT) 8 . While CVT is the established commercial production process, it suffers from high cost and low volume. In the CVT method, RP, tin iodide, and tin are sealed inside a glass ampule that is evacuated to ≈10 −3 Torr, heated above 825 K, and allowed to cool slowly yielding BP crystals over a 1-3 day period that grow outward from an intermediate Sn-P-I compound 9 .
Alternatively, high-energy planetary ball milling (HEPBM) constrains the pressure induced phase transformation inside a steel vessel under an inert atmosphere, alleviating potential safety hazards associated with the CVT method due to overpressurization of the glass ampules or incorrect cooling rates.
Additionally, no expensive and toxic mineralizers (i.e., tin, iodine) are required. In order to facilitate the scale-up of BP synthesis via HEPBM, details such as the conversion rate as a function of impact energy (i.e., ball density and size) and the cumulative energy input per gram of product formed are required. Several reports [10][11][12][13] of HEPBM on BP synthesis have explored the processing effects on synthesis time such as ball to powder mass ratio and rotation speed of the milling vessel but have not reported the scaling parameters. HEPBM is readily transferred to industrial mills capable of kilogram production capacities, provided that the energetics and kinetics of the phase transformation are well understood.
This study determines the necessary impact energy, the apparent rate constants, and the specific milling dose of the RP to BP phase transformation via HEPBM. By using closed-form analytical models of the kinematic motion inside the mill 14 along with estimates for the phase fraction evolution, an empirical relationship is established showing that higher impact energies lead directly to higher conversion rates while the total energy dose is shown to be an independent constant. Furthermore, the low cost nano-crystalline BP upon exfoliation is shown to be suitable for infrared optoelectronic applications, thereby paving the way for commercial scale-up opportunities.

Black phosphorus characterization
Visually, the powder changes color from a burnt umber red to a dull black with various shades of red/brown/black with increasing BP content, as shown in Fig. 1a. The bulk red and black powders exhibit similar morphology and agglomerations. The X-ray diffraction (XRD) patterns for the conversion of micron sized amorphous RP powder into orthorhombic BP are shown in Fig. 1b as a function of milling duration at the highest experimental milling intensity evaluated in this study. The as-received RP exhibits three broad peaks, with the first defined diffraction peak at ≈17°2θ. The powder diffraction file (PDF 03-065-2491) from the Inorganic Crystal Structures Database (ICSD 15 ) for BP is also shown, with relative peak intensities and positions indicated with the vertical line markers. A gradual shift to higher angles of the first sharp diffraction peak of RP occurs. By contrast, the broader high angle peaks generally decrease in intensity and width as the crystalline BP peaks emerge. XRD scans versus milling rpm using an internal silicon standard and the associated calibration curve are included in the Supplementary Figs 1 The scanning electron microscopy (SEM) images of Fig. 2a show that the RP starting powder exhibits an irregular faceted morphology with particle sizes from~0.1 to 100 μm. The converted BP exhibits similar faceted morphology but has a smaller particle size distribution below 10 μm. The powders appear to consist of soft aggregates of BP crystals embedded into amorphous particles which break down into flakes consisting of crystalline domains surrounded by amorphous matrix under mild ultrasonication. To further verify the crystalline BP domains within the amorphous flakes, transmission electron microscopy (TEM) images and selected area electron diffraction (SAED) 10 Fig. 2 Electron microscopy shows particle size reduction through BP synthesis and exfoliation while retaining polycrystallinity. a Scanning electron microscopy (SEM) images presenting faceted morphology of large particle sizes for RP, smaller particle sizes for BP, and narrow particle size for exfoliated BP flakes that consist of crystalline domains embedded in a lager bulk amorphous particle with plate-like morphology (SEM scale bars from left to right are 100 μm, 10 μm, and 2 μm, respectively). b Transmission electron microscopy (TEM) images and a selected area electron diffraction (SAED) pattern of polycrystalline BP manifesting large lattice domains and converted BP crystallinity. TEM scale bars from left to right are 50 nm, 10 nm, and 2 nm −1 respectively. S.V. Pedersen et al.
patterns are taken and presented in Fig. 2b, revealing the presence of crystalline BP domains embedded within the amorphous flakes and agglomerated powders. The volumetric particle size distributions obtained from solution-based laser diffraction particle size analysis (PSA) shown in Fig. 3a shows nonuniform size distribution for RP. This is an expected result considering the large variation in particle sizes observed in the SEM images (see Fig. 2) and the high sensitivity of volumetric PSA for detection of few but large particles. By contrast, the milling comminution process results in a nearly single mode distribution of the converted BP material, while exfoliation breaks up soft BP aggregates and produces flakes consisting of embedded crystal domains within an amorphous matrix.
To complement the XRD and TEM results on the converted material and verify the crystallinity of the exfoliated BP, Raman spectra are collected on the as-received RP, synthesized BP powder, and exfoliated BP flakes as shown in Fig. 3b. The expected BP vibrational modes 3 A 1 g , B 2g , and A 2 g , respectively, are clearly resolved, while the peak at 520 cm −1 is from the silicon substrate. The slightly raised background between 350 to 475 cm −1 is likely attributed to residual RP 11 .
As aforementioned, the direct bandgap of BP that spans 0.3-1.7 eV makes BP suitable for near-infrared wavelength applications. Few-layer (ca. 5 layers) BP exhibits optically pumped emission 2,5,16 near ≈1575 nm, whereas tri-and bi-layer BP emit at shorter wavelengths. The photoluminescence (PL) spectra for ball mill converted BP after exfoliation is shown in Fig. 3c. The broad emission from 950 to 1600 nm indicates the flakes contain bi-layer, tri-layer, and multilayer BP crystalline domains; however, further efforts to fully separate the various domains by lateral size or layer thickness by such means as density gradient centrifugation must be demonstrated prior to use in photonic applications 2,16,17 .

Conversion kinetics
Detailed individual reaction steps during mechanical milling processes are challenging to elucidate, although recent progress with synchrotron X-ray and Raman spectroscopy is now capable of providing in situ monitoring of reaction steps 18 . Nevertheless, empirical models based on ex situ measurements can provide general insight into the overall kinetic behavior, such as determining if the process accelerates to completion, undergoes deceleration, or progresses at a constant rate [19][20][21] . For evaluating the kinetics of a process, typically through a series of experiments, the energy input (such as a constant temperature or constant heating rate) is varied systematically and the evolution of the phase fraction, α, as a function of time is used to extract the rate constant, k 19,22 .
In mechanochemistry, the two kinetic models often used for solid-state phase transformations are the classical Johnson-Mehl-Avrami-Kolmogorov (JMAK) [23][24][25] , and a deterministic collision driven model derived by Delogu and Cocco (DC) 20,26,27 . Both models have been used to obtain the rate constants for single-phase systems undergoing transformations or amorphization and have been applied to multiphase systems undergoing mechanical alloying [23][24][25] . After obtaining the apparent rate constants as a function of milling parameters, the relationship between the milling impact energy, E b , and the rate constants can be determined. If a single impact is sufficient to drive the conversion to completion, then a linear dependence of the rate constants and the impact energies is expected 28 .
Simulations and in situ monitoring 27,29-36 of media interactions in high-energy ball mills have shown agreement with the underlying assumptions and estimated impact energies obtained through kinematic models of high-energy planetary ball mills. While existing models differ in determining the media velocity at impact 14,30,31,34,37 , all models assume the kinetic energy imparted to the media through collisions with the vial walls is transferred into the powder during impact, at least to a first order approximation (some losses due to friction, heat generation, and inelasticity are expected). Nonetheless, the impact energy, E b , is simply related to the media mass, m b , and the impact velocity, v in , as: By extension, the cumulative dose, E dose , is a function of the number of balls, N b , the collision frequency, f b , the milling time, t, and powder mass, m p , as: Using the conversion data obtained from XRD, the method described by Galwey et al. 19 was employed to generate a master plot to determine which subset of solid-state reaction models are applicable, as shown in Supplementary Fig. 6. The JMAK model closely matches the dataset, while the DC model follows the data above 30% conversion. As seen in Fig. 4a for 300 rpm, the rate of conversion slows down, and this transition is visually distinct for the 300 rpm run versus the higher intensities of 400-600 rpm. Similarly, for 200 rpm, the conversion occurs slowly with an induction time of nearly 100 h before BP is observed. Transforming the data using the JMAK model, Fig. 4b, or the DC model, Fig.  4c, shows two or three linear segments. In practice, line segments that show different slopes are treated individually with linear regression lines, as it has been shown that these distinct transitions are often due to changes in growth dimensionality 38 .  Raman spectra confirming the amorphous nature of RP, as well as the crystalline nature of both converted and exfoliated BP. c Photoluminescence (PL) emission of exfoliated BP on SiO 2 with a broad spectrum due to the bilayer to multilayer BP in the exfoliated material.
As a reminder, the DC model is, mathematically, a restricted form of the JMAK model that fixes the exponent n = 1, resulting in an apparent S-shape about the linear regression line in the DC linearization plot. Functionally, this restriction prevents making inferences about the nucleation and growth behavior afforded by the JMAK model, but as shown in Fig. 4d the extracted rate constants are still comparable. In all cases, both the JMAK and DC models appear to be valid, and the linear fits model the data up to 90% conversion well with values for the coefficient of determination (a goodness of fit statistic relating the percentage of the predicted response variation that is modeled by the regression line), R 2 , above >95%. Deviations 19 from the model are likely due to the initial estimation of the phase fraction below 10 wt% and above 70 wt%. In these regions, there appears to be a consistent decrease in conversion rate for all milling intensities, implying a mechanistic change that requires higher mechanical impact energy to progress the conversion. Supplementary Figure 7 (plotted on a log scale) includes the 200 rpm curve; the 200 rpm dataset is not shown in Fig. 4c, as the induction time obscures the higher rpm behaviors. Figure 4d shows the linear dependence of the apparent rate constants as a function of the calculated impact energy. The calculated impact energies and milling intensities using the kinematic model as applied to the HEPBM and the ascertained kinetic parameters for both the JMAK and DC models are presented in the Supplementary Table 1 and Table 2, respectively. Figure 5 shows the phase fraction evolution as a function of the total energy imparted to the powder on a mass basis, revealing three distinct behaviors. First, for higher milling intensities (i.e., mechanical work performed) above 300 rpm, the total dose required for complete conversion is the same. Second, below 300 rpm, there is a transition at low phase fractions (about 30 wt %) requiring progressively longer milling times, as well as a plateau near 75 wt%. This plateau indicates a potentially limiting step which cannot be overcome by the impact energy at 300 rpm. Third and finally, at even lower energy (i.e., 200 rpm), nucleation of BP is evident, but only after a long induction time of 200+ hours.

DISCUSSION
The reduction in particle size of RP as it converts to BP due to comminution is typical of ball milling and is desirable as it increases the surface area available for reactions and promotes a single mode particle size distribution which benefits down-stream processes such as exfoliation. The broad BP XRD peaks in Fig. 1b indicate a high degree of nano-crystallinity while the broad base of the peak at 17°2θ indicates some amorphous material is still present. Supplementary Figures 8 and 9 show TEM images of randomly orientated lattice planes surrounded by an amorphous matrix within a single exfoliated flake, indicative of the nature of point impacts that have occurred and that presumably each impact will produce a local formation of BP within the bulk RP  Fig. 4 Mechanochemically induced phase transformation kinetics of RP to BP using the JMAK and DC models. a XRD-based phase fraction versus logarithmic time showing representative S-curves for each rpm condition; higher rpm results in faster conversion while a dramatic reduction occurs at and below 300 rpm. An S-shape curve is typical of nucleation and growth dominated solid-state processes. Error bars represent one standard deviation, lines are visual aids. b JMAK linearization plots for extracting rate constants and exponents, lines are linear regression fits for each straight section; the slopes remain relatively constant indicating the conversion mechanism remains constant, while the kinks at higher phase fraction indicate a possible mechanism change. c DC model linearization plot excluding the 200 rpm data, lines are linear regression fits; the slopes become steeper (higher rate constant) with increasing rpm. d Linear regression fits for rate constants versus calculated impact energies; both the JMAK and DC models show the linear dependence of the rate constant on the impact energy confirming the assumption that a single impact is sufficient for converting the entrapped RP powder between media collisions to BP. particle 10,12 . Amorphous RP possesses some degree of short and medium range order; it has been described as being a 3, 4, and 5 sided cage-like structure linked together by phosphorus atoms in a polymer like manner 39,40 . Presumably, the phosphorus cages are collapsed into the buckled planes due to the high localized compressive forces from the media collisions 39,41,42 . The Raman spectra in Fig. 3b on HEPBM synthesized powder (exfoliated in IPA) shows only modes that correspond to BP, which indicates no significant RP is present. Together with the TEM data, the lack of features corresponding to RP indicates that amorphous BP is likely present [43][44][45][46] . High pressure conversion of RP thin films into BP has also reported the formation of amorphous BP 44 . High-energy ball milling has been used to produce amorphous metals and even cycle between crystalline and amorphous phases of alloys by controlling the milling intensity and duration 47 . In light of this, the coexistence of both crystalline and amorphous, highly disordered, BP within ball mill synthesized BP flakes is unsurprising and has been observed through nuclear magnetic resonance spectroscopy 11 .
While some fraction of BP is present in the XRD patterns at 200 and 300 rpm, such conversion rates are economically impractical due to the power consumption; namely, mill motors constitute the highest recurring electrical costs. The 200-rpm milling intensity shows an entirely different behavior which may be due to the tendency of media to slide along the walls of the vessel with few high-energy collisions at low vial rotation speeds. As a result of the less energetic impacts and low frequency of collision, a significantly higher milling dose is required to producing meaningful quantities of BP at 200 rpm.
From Fig. 4d, the x-axis intercept for a negligible rate constant, k = 0, ideally indicates a minimum impact energy, E b , near 0.5 mJ, but these values do not reflect the apparent change above 70 wt%. As such, this minimum value most likely represents the energy required to form stable nuclei but with no practical grain growth. The utility of Fig. 4d is that extrapolation of an expected rate constant due to higher impact energies is feasible, up to a limit set by the media size, density, vessel diameter, and achievable rpm.
The main scaling parameter for HEPBM is obtained from Fig. 5 and this dictates the optimum milling condition. The product of the milling intensity, E b ·f b ·N b , and time, t, is the same for milling at 600, 500, and 400 rpm. That is, above a certain impact energy (E b = 25 mJ), the conversion takes the same total energy input. The trade-off then is the milling time versus potential media wear and design considerations for the mill motor. An apparent specific milling dose of ≈100 kJ/gram is required to achieve a BP conversion above 90 wt%. While a higher milling dose may improve the yield, it may also introduce more defects into the crystalline domains, as well as produce more media wear and possible contaminations from the milling vessel and media.
The results presented here show that the conversion of RP to BP follows the characteristic S-shape of sigmoidal processes indicative of nucleation and growth dominated phase transformation. The rate constants obtained by applying the JMAK and DC models in Fig. 4 show the expected linear behavior as a function of milling intensity; accordingly, higher impact energies and higher collision frequency should result in more frequent nucleation events and, hence, a faster conversion rate. One benefit of the DC model is that it is less sensitive to induction periods than the JMAK model, such that the extracted rate constants better represent the overall conversion rate. However, the JMAK model allows for inferences into the nucleation and growth behavior.
Using the JMAK model, the values of the exponent can potentially provide insight into the dynamics of nucleation and crystal growth. The calculated exponent value is n ≈ 2.3-2.5, suggesting that RP to BP undergoes continuous nucleation and subsequent interfaced controlled transformation in twodimensions 48,49 . The reduced exponent value below 3 can be attributed to the anisotropic growth and blocking effects, or it could be an outcome of a mixture of 3D and 2D growth with site saturation [48][49][50][51] . In theory, polymorphic phase transformations are categorized as either displacive (rapid and diffusionless) or reconstructive (slow and diffusion limited) [52][53][54] . In practice, assumptions regarding diffusion, preferred nucleation sites such as defects, grain boundaries, or dislocations, and experimentally determined crystal growth habit (i.e., spherical, disc, or rod growth) are required to make reliable inferences from the JMAK exponent.
However, the unique amorphous to crystalline transformation of RP to BP does not fit cleanly into any phase transformation category 53 . We argue that RP to BP is likely a rapid diffusionless, displacive transformation, provided that RP bonds are broken either through thermal or strain input. Generally, high applied pressure inhibits diffusion due to the increase in repulsive forces upon shortening the interatomic bond distance 52 but can enable rapid short displacements via re-orientation 52 . By contrast, grain growth itself is typically a thermally driven process. High pressure studies on the white phosphorus (≈1.2 GPa) or RP (≈8 GPa) to BP transition occurred almost instantly after an incubation time of 10-15 min when exceeding a minimum temperature of 200°C 7,41,43,55-58 . This supports the notion of rapid interface propagation instead of a time-dependent diffusion-controlled grain growth. Comparing our mechanochemical results to Xiang et al. 44 and to those obtained by Bridgman 59 using rotating diamond anvils (≈4 GPa), we contend that the plastic strain facilitated by the direct compression of RP between colliding media breaks molecular bonds within the RP polymer chains, leading to the rapid formation of the orthorhombic crystalline phase. The TEM results here of flakes show clear nanometer-sized lattice domains randomly orientated surrounded by regions of an amorphous matrix, consistent with random nucleation but limited grain growth. Similar behavior was observed with static anvil experiments at a higher pressure (≈6.8 GPa) but a much lower temperature (≈25°C) than typical synthesis conditions. Xiang et al. 44 revealed a slow transformation rate of tens of hours in which the rhombohedral (A7) structure of metallic BP nucleates first from RP and persists indefinitely (up to 71 h). However, the material rapidly converts to orthorhombic BP upon decompression 44 . Collisions in HEPBM are effectively microforges where the RP is entrapped between the media in a millisecond timescale collision. Upon the media rebounding, the BP crystallizes 47 .
In support of our hypothesis, the fact that the rhombohedral crystal structure of BP forms primarily under the applied   Figure 10 shows the in situ temperature measurements obtained during mechanochemical conversion; while the global temperature does not exceed 40°C, it is possible that the local impact energy could heat the powder substantially [60][61][62][63] . Still, the temperature rise is also highly dependent upon the material undergoing the treatment (i.e., if the process results in an exothermic heat release such as in the milling of thermite, metal oxide formation, or sulfide reactions) 20,47,64 . We exclude the possibility of high localized temperatures at the point of impacts, as we do not observe substantial grain sizes. Furthermore, recent work 65 comparing torsional or rotating diamond anvil methods to mechanochemistry shows striking similarities in that the generated plastic strain can reduce the pressure/temperature requirements for phase transformations 66 in Si, Ge, boron nitride 67,68 , and zirconium 69 with similar nanostructures. In practice, knowing the rate constant and specific milling dose enables determination of the mill design (diameter, motor horsepower, media size, and quantity) in order to achieve a given production rate. The ascertained specific milling dose of ≈100 kJ/ gram and corresponding minimum impact energy near 25 mJ/ impact provide meaningful parameters for scale-up either by number of milling stations (i.e., 100 g capable planetary ball mills) or by switching to larger tumbler mills (i.e., 1 kg capacity) and running for the appropriate milling time derived by considering the rate constant, impact energy, and desired amount of product; the difference between the planetary mill versus the tumbler mill being the rate of energy deposition versus the synthesis time 70 .
In conclusion, we have ascertained the rate dependency on impact energy and the scaling parameter (milling dose near 100 kJ/gram) required to convert low cost RP into BP by the highly scalable high-energy ball milling method. While the mechanical impacts readily form BP nuclei, the overall process is grain growth limited; further efforts should seek to enhance the grain growth. The exfoliated BP flakes containing embedded crystalline domains within an amorphous matrix produces PL spectra with emission in the near-infrared range, thereby revealing the potential use of this synthesis route provided improved exfoliation methods are developed to isolate layer-dependent BP nanocrystals. Producing commodity-grade BP is an enabling step to employ the nanomaterial in commercial electronic and optoelectronic devices.

METHODS
In the milling experiment, elemental red phosphorus (1-5 g; Alfa Aesar, 97% purity, −60 mesh) and stainless steel media (10 mm diameter, 100 g) were weighed and loaded into 250 mL stainless steel vessels inside an argon atmosphere glovebox (P O2 < 0.1 ppm) and sealed with internal gas pressures between 90 and 95 kPa at room temperature. The ball to powder mass ratio was 20:1. All milling experiments were conducted with the PM100 planetary ball mill (Retsch, Germany) with pauses in the milling cycles (30 min on, 30 min off) to prevent overheating of the unit and to switch the direction of flywheel rotation. Milling conditions (rpm and duration) were varied between 200 and 600 rpm and 0.15 to 25 h, depending upon the intent of the studies. For the kinetic study, each run was completed three times except the 200 rpm run which was allowed to run once for 220 h. The three runs at 300, 400, 500, and 600 rpm were paused periodically to obtain samples for phase fraction analysis. For each sampling interval, the sealed vessel was first transferred into an argon glovebox. An appropriate amount of powder was extracted from the free powder at the bottom of the milling vessel and sealed inside glass vials; no powder was added back into the milling vessel. The vessel was resealed and transferred out of the glovebox to resume the milling run. An initial milling run at 600 rpm for 2 h with 1 gram of red phosphorus was used to condition the media and vessel surfaces before each milling study. The milling vessel was wiped clean between each study but was not cleaned with isopropyl alcohol, thereby preserving the conditioned state of the vessel.
To obtain few-layer BP suspensions, the as-converted BP powders (0.25 g) are solvent exfoliated in isopropanol (30 mL) using a low-energy mixer mill (MTI, SFM-2) with (1 mm) diameter chrome steel ball media sealed inside a 50 mL plastic container for 24 h at a rotational speed of 150 rpm. The suspension is separated from the media using a sieve (−300 mesh) and stored in glass vials. The suspension is drop cast onto silicon substrates for Raman and PL measurements.
Powder XRD data was obtained in the Bragg-Brentano geometry (CuK α , Rigaku MiniFlex 600, Japan) on 0.165 g of powder. The internal standard method using NIST Silicon 640D powder (ICSD PDF 00-027-1402) was used to generate a calibration curve for estimates of the phase fraction obtained.
Powder Raman spectroscopy (Raman) data was obtained (Horiba LabRam HR Evolution, Japan) on densely packed samples (0.01 grams) prepared inside an argon glovebox that were sealed inside a custom-built enclosure with a sapphire window to prevent oxidation and ignition during measurements (both red and black phosphorus in finely divided powder form can be readily ignited in air even when using neutral density filters and low laser power settings). A 50-mW 532-nm laser source (Nd:YAG) was used to generate spectra collected with a Si CCD. Spectra obtained at ×20 objective with double acquisitions for 60 s at ten random locations throughout the bulk surface using the 1800 lines/mm grating.
PL measurements were obtained using a Princeton Aston spectrometer with a NIRvana InGaAs thermoelectrically cooled array and an excitation wavelength of 532 nm at 10 mW. Data were collected for 5 s using a 100X near-infrared objective.
Solution-based laser diffraction PSA was performed on powders ultrasonically dispersed in isopropanol for 1 min and cycled through the Horiba LA-950 particle size analyzer and solution flow cell; volumetric distributions were obtained using dual 405 and 650 nm laser and light diodes. The complex indices of refractions for RP 71 and BP 72 were taken from literature.
SEM was performed on powders pressed onto carbon tape and imaged at low beam currents and accelerating voltages using a FEI Teneo FESEM. Exfoliated BP powders were drop cast onto cut SiO 2 wafers and imaged under similar conditions using secondary and backscattered electron imaging modes.
TEM was performed using a JEOL-JEM-2100, Japan, with Orius SC1000 CCD camera; exfoliated samples were drop cast onto lacy carbon grids and transferred into the vacuum chamber within 5 min. A low accelerating voltage of 100 keV was used with minimal dwell times roughly 30 s for imaging to prevent the amorphization of the flakes.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
All relevant code is available from the authors on request.