Uncovering the morphological effects of high-energy Ga⁺ focused ion beam milling on hBN single-photon emitter fabrication

The advancement of quantum information technologies (QITs), such as quantum transducers and sensors, requires solid-state qubit platforms that host optically addressable single-photon or spin states that couple strongly to external degrees of freedom and can be easily engineered for integration into hybrid systems.^1–4 Low dimensional materials are desirable solid-state single photon emitter (SPE) hosts because they exhibit unique electronic, magnetic, and mechanical properties impossible for their bulk counterparts.^5–10 Two-dimensional wide bandgap insulator hexagonal boron nitride (hBN) hosts bright and stable room temperature SPEs that couple strongly to applied electric^6,7 and strain^11,12 fields, as well as optically detectable magnetic resonance (ODMR) with spin state splitting sensitive to external magnetic fields.^13–16 Furthermore, hBN functions well as both a dielectric substrate for other 2D materials, such as graphene^10,17,18 and transition metal dichalcogenides (TMDCs),^19–21 and as a monolithic platform for photonic circuits,^22–24 making it a versatile material to integrate into QITs.

Fabrication challenges, such as patterning hBN SPEs^25–32 and direct-writing nanostructures in hBN,^33–35 have been resolved using focused ion beam (FIB) exposure. However, the high energies required to pattern emitters at precise locations and mill nanostructures cause changes to the morphology of thin hBN crystals that may influence further fabrication steps and impact the material’s ability to host SPEs.^21,35 Moreover, while SPEs found in hBN after FIB exposure have been attributed to either vacancies created deep within the crystal or edges created due to milling, it is not understood how changes to surface morphology that occur due to FIB exposure influence the probability of SPE creation.

To uncover the morphological life-story of hBN crystals to support SPEs while undergoing high energy Ga⁺ FIB exposure, we fabricated samples with a range of exposure doses and characterized changes to the surface roughness, milling depth, swelling, and single photon emission photophysics. Then, we analyzed the correlation between morphological factors and photophysical responses. Furthermore, we performed Stopping Range of Ions in Matter (SRIM) calculations³⁶ to investigate how vacancy concentration influences morphology and SPE creation. This work provides a necessary framework to elucidate how pristine hBN crystals become damaged in controlled ways at the atomic and nano-scales to perform as needed for quantum technologies.

To create samples that contain many visible hBN crystals, we performed an exfoliation transfer of high-pressure high-temperature (HPHT) grown hBN crystals onto Si/SiO₂ substrates using heat-release tape and annealed samples at 500 °C to clean residual tape from the surface. Once we identified crystals of area ∼200 μm² [Fig. 1(a)] and thickness in the range of 10–17 nm (supplementary material, Fig. 1), we performed Ga⁺ FIB patterning tests. Based on our previous work²⁵ that created SPEs at precise locations, we patterned arrays of 500 nm diameter circles with a beam acceleration voltage 20 kV and a beam current of 59 pA. The Ga⁺ FIB used for this work had a minimum exposure time of 1.3 ms for the specified pattern size and current, corresponding to the shortest loop time the FIB can direct-write the circle pattern. We adjusted the ion fluence by changing the exposure time, thereby changing the number of passes performed by the software to pattern the circle and tested 1.3 ms (FIB level 1: 1 pass), 2.5 ms (FIB level 2: 2 passes), 3.8 ms (FIB level 3: 3 passes), and 5.0 ms (FIB level 4: 4 passes) [Fig. 1(b)]. We annealed the exposed samples in a home-built chemical vapor deposition (CVD) system used for vertically aligned carbon nanotube (VACNT) growth with a C₂H₄:H₂:Ar flowing gas mixture at 1000 °C for 5 h in order to create a carbon-rich environment that is expected to activate hBN quantum emitters.^37–39 For more information about the emitter fabrication process, refer to the supplementary material.

After the full emitter fabrication process, we performed tapping-mode atomic force microscopy (AFM) to measure the topography of nanoscale surface features and observed qualitatively different morphological features at each exposure level [Fig. 1(c)]. From the AFM height data, we calculated the following morphological characteristics of 1 × 1 μm² regions containing one exposed circle each: average mill depth $Δz−̄$⁠, maximum mill depth $Δz−max$⁠, mill volume ΔV⁻, average swell height $Δz+̄$⁠, maximum swell height $Δz+max$⁠, swell volume ΔV⁺, and exposed edge surface area SA [Fig. 2(a)]. Local unexposed surface height z₀ and rms roughness RMS₀ for each 1 × 1 μm² region were calculated from points outside of the exposed pattern area and used as the zero height for calculations and as the threshold for identifying FIB induced damage, respectively. Milling values $Δz−̄$⁠, $Δz−max$⁠, and ΔV⁻ were calculated from data points with height below z₀ − RMS₀ and characterized the depth of milled locations and total amount of material removed from the exposed regions. Swelling values $Δz+̄$⁠, $Δz+max$⁠, and ΔV⁺ were calculated from data points above z₀ + RMS₀ and characterized the height and total amount of FIB induced swelling observed in exposed regions. Exposed edge surface area calculations were performed using all data points in the region, with a baseline area of 1 µm² subtracted. For more details about the calculations, refer to the supplementary material.

We report the average values for all calculated morphological characteristics for each exposure level in Table I and show the average values of $Δz+max$ and $Δz−max$ for each exposure level in Fig. 2(b), with average values for ΔV⁺ and ΔV⁻ for each exposure in Fig. 2(c). Plots showing the average values for all calculated morphological characteristics as a function of FIB exposure level can be found in the supplementary material. Based on these calculated values and observations made from the data visualized as AFM scan plots in Fig. 1(c), we found that swelling characteristics and milling characteristics evolve differently with increasing levels of FIB exposure.

Starting at the lowest FIB exposure level (1.3 ms), exposed regions contained small and randomly located areas of milling and swelling within the patterned area with $Δz−̄1$ = −3.26 ± 0.15 nm and $Δz+̄1$ = 1.92 ± 0.09 nm. Regions exposed at the second FIB level (2.5 ms) showed qualitatively similar milling and swelling characteristics and calculated values were similar to the first exposure level, with $Δz−̄2$ = −2.89 ± 0.11 nm (t-test: p = 0.1566) and $Δz+̄2$ = 1.81 ± 0.06 nm (t-test: p = 0.3245). The largest change observed was an 8% decrease in the magnitude of $Δz−max$⁠, with $Δz−max1$ = −8.63 ± 0.33 nm and $Δz−max2$ = −7.92 ± 0.22 nm, but the difference between these values was not strong enough to draw any conclusions about the evolution of morphological characteristics from the first exposure level to the second (t-test: p = 0.0837). Therefore, we determine that there are no significant morphological differences between regions exposed for 1.3 ms and those exposed for 2.5 ms.

At the third exposure level (3.8 ms), we began to observe significant changes to morphological characteristics. While regions exposed at the third exposure level were milled similarly to those at the second with $Δz−̄3$ = −3.25 ± 0.14 nm (t-test $Δz−̄$⁠: p = 0.1742), both the average and maximum swell heights increased from the second to the third FIB level (t-test $Δz+̄$⁠: p = 0.0006, t-test $Δz+max$⁠: p < 0.0001), with the mean $Δz+max$ increasing 44% from $Δz+max2$ = 2.57 ± 0.11 nm to $Δz+max3$ = 3.70 ± 0.23 nm.

Regions exposed at the fourth level (5.0 ms) showed swell heights similar to the values calculated at the third (t-test $Δz+̄$⁠: p = 0.4653, t-test $Δz+max$⁠: p = 0.3754). However, ΔV⁺ increased 23% from $ΔV+3$ = 1.89 ± 0.22 × 10⁵ nm³ to $ΔV+4$ = 2.44 ± 0.23 × 10⁵ nm³ (t-test: p = 0.0444). The areas of swelling for regions exposed at the fourth level appeared more localized to the pattern edges, away from the pattern center, whereas swelling for regions exposed at the third level was observed within the pattern center. While the swelling was pushed out to the pattern edges, we observed deeper milling more localized to the circle pattern at this level, with the magnitude of $Δz−max$ increasing 62% from $Δz−max3$ = −8.19 ± 0.27 nm to Δz⁻ max₄ = −13.23 ± 0.32 nm (t-test: p < 0.0001). Furthermore, we observed similar behavior in the exposed surface area SA, where SA₃ = 1.67 ± 0.07 × 10³ nm² increased 147% to SA₄ = 4.12 ± 0.29 × 10³ nm² (t-test: p < 0.0001) most likely due to the increase in directional milling forming connected sidewalls at the fourth exposure level. We conclude that swelling is the dominant morphological characteristic starting at the third exposure level and is a precursor to the deeper directed milling that began at the fourth level.

To determine how the targeted morphological characteristics influenced emitter creation, we identified the fluorescent locations via confocal microscopy with 532 nm excitation [Fig. 3(a)] and characterized quantum emission via Hanbury Brown and Twiss experiment, fitting g²(t) data to the three-level system model for antibunching and using the accepted threshold of g²(0) < 0.5 to denote SPE.²¹ To further confirm the quantum nature of emitters, we measured zero-phonon line (ZPL) spectra [Fig. 3(b)] within the expected range of 1.8–2.2 eV for room-temperature visible range emission in hBN under 532 nm (2.33 eV) excitation.^25,39,40 We calculated the average PL of exposed regions at each level by image analysis on confocal PL scan plots segmented into the regions of 4 μm², with 1 patterned site per 1 μm². We also calculated the average SPE yield of each segmented 4 μm² region by counting the number of emitters with g²(0) < 0.5 and dividing by 4, the number of patterned sites. For more details on the average PL and SPE yield calculations, refer to the supplementary material.

We observed PL emission from successfully patterned sites exposed at all levels with region PL and SPE yield dependent on the exposure level. Patterned sites that bordered or overlapped with an hBN flake edge accounted for ∼1% of all patterned sites and exhibited low PL and SPE yield regardless of the exposure level. We expect that partially patterned circles created structural environments inconsistent with the fully contained patterned sites intended. Therefore, we considered circle pattern bordering or overlapping with hBN edges to be unsuccessfully fabricated and excluded those sites from our analysis. Of the regions that hosted non-zero SPE yields, we found that patterned sites exposed at the third level hosted SPEs with a yield of 0.047 ± 0.016 (mean comparison test to zero: p = 0.0027) and the fourth hosted SPEs with a yield of 0.115 ± 0.025 (mean comparison test to zero: p < 0.0001) [Fig. 3(c)]. The regions exposed at the two lowest levels, which did not host SPEs, exhibited brighter regions on confocal scan plots, with PL₁ = 573 ± 27 cts and PL₂ = 1097 ± 56 cts, while the third level exhibited lower PL on average with PL₃ = 186 ± 25 cts (Tukey test: p < 0.0001) and the fourth even lower on average with PL₄ = 46 ± 3 cts (Tukey test: p < 0.0001) [Fig. 3(d)].

We conjecture that brighter PL indicates that many optically active vacancy defects remained within the crystal at exposed regions, so excited locations within the diffraction limited spot exhibited ensemble emission rather than single photon emission. We further expect that higher SPE yields were possible when many optically active defects were removed by milling, reducing the overall background PL. This conjecture is consistent with the highest SPE yield, lowest region PL, and largest milling observed at the highest FIB exposure level. However, it is inconsistent with the non-zero SPE yield, lower region PL, and lower milling observed at the third exposure level. Therefore, the morphological characteristic of deep and directed milling cannot be fully responsible for the optical characteristics of FIB exposed hBN crystals.

When comparing all calculated morphological characteristics to the measured optical characteristics [Fig. 4(a)], we found that the PL intensity at each FIB exposure level was strongly correlated only to the maximum swell height at each FIB level, with c(PL, $Δz+max$⁠) = −0.969 (p = 0.0309), while the SPE yield at each FIB level was strongly correlated only to the average volume of swell areas at each FIB level with c(SPE, ΔV⁺) = 0.9955 (p = 0.0045) [Fig. 4(b)]. Furthermore, the observed increases in material swell height and volume were the only morphological features that correlated even weakly (0.05 < p < 0.20) to the increased SPE yield and diminished PL simultaneously. We expect that the morphological characteristic of swell is indicative of damage levels within the crystal that decrease PL and increase the likelihood of hosting SPEs prior to deep and directed milling.

To elucidate possible shared causes for both the observed swelling due to FIB damage and the observed optical characteristics on the atomic level, we performed SRIM calculations for vacancy concentrations. We performed the monolayer collision steps/surface sputtering calculation, which was considered optimal for the thickness range of flakes studied in this work.³⁶ From the SRIM calculations, we acquired the vacancy number per ion as a function of depth and lateral position N_v/Ga+(x, z) for a 20 keV Ga⁺ beam entering the material at a single point, where SRIM assumes a beam diameter of 0. To determine the ion exposure profile D_Ga+(x, y) for the full 500 nm circular pattern, we modeled the beam as a 2D Gaussian⁴¹ and summed over all points in the pattern following a circular beam path, with the beam pitch set constant as in the patterning software for all patterning experiments [Fig. 5(a)].

We simulated the vacancy concentration n_v(x, z) in exposed regions by convolving N_v/Ga+(x, z) at a given point with the simulated exposure profile D_Ga+(x) at that point over a 1 μm position range. Finally, to mimic the material conditions upon which optical data was collected, we further convolved n_v(x, z) with the height data collected by AFM at each point, removing areas of the simulated material at locations where the experimental material was milled a depth −Δz. For more information about simulation parameters and how we accounted for milling at each increased FIB exposure level, refer to the supplementary material.

Based on the results of these simulations, shown in Fig. 5(b), we observed that the n_v(x, z) increased steadily with each increase in the FIB exposure level. Furthermore, while the significant increase in milling at the fourth exposure level removed a large volume of regions with high vacancy concentrations, many regions with high n_v remained due to $Δz−max$ less than the starting thickness of material h₀. This is consistent with our hypothesis that although milling away high n_v areas contributed to an overall lower number of optically active defects at the fourth level, the onset of milling distinct circular holes alone cannot fully justify the PL decrease and onset of isolated SPEs we observed. Furthermore, since the level of swelling correlated strongly to these observed optical characteristics, we expect that there exists some threshold vacancy concentration $nvswell$ that leads to the observed increase in swelling and, therefore, the observed changes to optical characteristics. In addition, by combining $nvswell$ with a corresponding threshold vacancy concentration $nvmill$ that leads to the observed onset of directed and localized milling, we expect that we can determine different regimes of material damage and qualitatively describe how they influence the material’s ability to support SPEs at patterned locations.

To find the value for $nvmill$ from our simulated data, we first calculated the onset exposure time t₀₋ by fitting the normalized remaining thickness $1−Δz−̄(t)/h0$ to a complementary error function with a constant offset, finding the intersection t₀₋ of the maximum amplitude with the negative slope linear regime, and calculated n_v(t₀₋) from a linear fit of n_vmax(t). This process is illustrated in Fig. 5(c), which shows the normalized remaining thickness fit plotted with n_vmax(t). We similarly found $nvswell$ by fitting $Δz+̄(t)$ to an error function with constant offset, finding the intersection t₀₊ of the minimum amplitude with the positive slope linear regime, and calculated n_v(t₀₊) from the linear fit n_vmax(t). From these calculations, we found $nvswell∼150$ nm⁻³ and $nvmill∼243$ nm⁻³. Furthermore, we similarly fit $Δz+max$ to calculate a simulated concentration that corresponded to the maximum observed swelling and obtained $nvmax⋅swell∼186$ nm⁻³. Using these threshold values, we simulated the change in the thickness profile as a function of position Δz(x) for the third and fourth FIB exposure levels by plotting the maximum vacancy concentration as a function of position n_vmax(x) with positive increasing values if the swell threshold was crossed and before the maximum swell, positive decreasing values if the swell threshold was crossed and after the maximum swell, and with negative values if the mill threshold was crossed, as shown in Fig. 5(d). For more details about the simulated thickness profile calculations and fits, refer to the supplementary material.

We observed that the simulated height profile was highest toward the center of the exposed region for the third exposure level, whereas the simulated height profile for the fourth exposure level was highest toward the edges, corresponding to the locations where swelling was observed, in good qualitative agreement with our AFM height data [Fig. 5(d)], as well as results from Glushkov et al.²⁶ that show swelling around the edges of circles patterned with high-energy Xe+ FIB. Additionally, results from Glushkov et al. show that emission is localized to the edges of these patterned circles with swelling, thereby demonstrating that optically active vacancy defects that exhibit photophysical properties of SPEs remain in regions where swelling is observed. Because we only observed swelling in exposed regions with $nvswell<nv<nvmill$⁠, we expect that vacancy concentrations within this range correspond to a distinct type of crystallographic defect, causing higher swelling, reduced region PL, and the onset of non-zero SPE yield.

Supported by the morphological characteristics calculated from AFM height data, the optical characteristics calculated from confocal PL and antibunching data, and the vacancy concentration characteristics calculated from SRIM and ion exposure simulations, we expect that at some critical vacancy concentration $nvswell<nv<nvmill$ present after FIB exposure at the third and fourth levels studied in this work, vacancies coalesce into larger volumetric crystal defects, or voids, indicated by the swelling level,^26,36,42 that are not optically active at an excitation wavelength of 532 nm [Fig. 6]. The process of void nucleation and growth in irradiated materials is well studied for materials similar to hBN^43,44 and is consistent with recent investigations into nanopore growth in irradiated monolayer hBN.^45–47 Furthermore, the atomic density n/V_c of multilayer hBN with lattice constants a = 2.5 Å and c = 6.6 Å⁴⁸ with AA′ stacking⁴⁹ is ∼166 nm⁻³, which lies within the simulated range for n_v where swelling is observed. This is consistent with our expectation that voids form when vacancies dominate the underlying crystal structure, which can be described as n_v > n/V_c.

While local layer mixing,⁴² Ga⁺ and recoil ion implantation,³⁶ and carbon impurity doping from the emitter activation process may contribute in small parts to the observed swelling, they cannot fully describe the sudden transition in optical characteristics observed between the second and third FIB exposure levels without considering a critical n_v that causes void nucleation. Because voids act as sinks for vacancies,⁵⁰ the presence of voids in higher quantities starting at the third exposure level further supports the observed decrease in region PL and better isolation of SPEs at patterned sites despite similar milling to the first and second exposure levels. Smaller, optically active vacancy defects would no longer be evenly distributed within the crystal, causing background fluorescence during emitter excitation, and would instead be concentrated at the void locations, coalesced as part of the large volumetric defects that would have different optical characteristics and would not change the overall vacancy number, remaining consistent with our vacancy concentration calculations.

Considering the results, calculations, and analysis presented in this work, we propose the following model for the fabrication of isolated SPEs via high-energy Ga⁺ FIB milling of pristine hBN crystals:

With no FIB exposure, pristine hBN hosts no SPEs. Un-exposed and carbon-annealed hBN hosts randomly located SPEs with a density of 0.0125 μm² (mean comparison test to zero: p = 0.0416) likely due to low concentrations of native vacancy defects forming complexes with carbon impurities.^37,51

At the lowest FIB exposure level (1.3 ms), exposed areas host many optically active vacancy defects [Figs. 3(a), 3(b), and 5(b)] with visible-range emission under 532 nm excitation, and we do not observe isolated SPEs and ZPLs with traditional diffraction limited confocal techniques. This is characterized by the high region PL [Fig. 3(d)] and zero SPE yield [Fig. 3(c)]. Some surface level milling and sputtering occurs and few areas with particularly high vacancy concentrations contain voids. This is characterized by randomly located pockets of milling and nonzero swelling within the exposed area [Fig. 1(c)].

At the next FIB exposure level (2.5 ms), exposed areas host more optically active vacancy defects [Figs. 3(a), 3(b), and 5(c)], and we again did not isolate SPEs with traditional confocal techniques. This is characterized by the similar inconsistent sputtering of the exposed surface compared to the first level (t-test: p = 0.1566) [Figs. 1(c) and 2(b)], the 90% increase in region PL (t-test: p < 0.0001) [Fig. 3(d)], and zero SPE yield [Fig. 3(c)]. Areas with particularly high vacancy concentrations host voids. This is characterized by small areas of randomly located nonzero swelling similar to the first level (t-test: p = 0.3245) [Figs. 1(c) and 2(b)].

At the third FIB exposure level (3.8 ms), exposed areas have even higher concentrations of vacancies [Fig. 5(c)]. Many more regions have reached a critical concentration of vacancies and therefore, contain more voids than optically active vacancies or complexes, but voids have not coalesced into structural defects that cause significant material breakage and directed milling. This is characterized by surface level milling and sputtering similar to the previous level (t-test: p = 0.5058) [Figs. 1(c) and 2(b)], the 44% increase in surface swelling height (t-test: p < 0.0001) [Fig. 2(b)] localized toward the center of the exposed area [Figs. 1(c) and 5(d)], the 83% decrease in region PL (t-test: p < 0.0001) [Figs. 3(a) and 3(d)], and the onset of non-zero SPE yield (mean comparison test to zero: p = 0.0027) [Fig. 3(c)].

At the fourth and highest FIB exposure level tested (5.0 ms), exposed areas have the highest concentrations of vacancies, but as voids continue to grow, complete breakage and, therefore, milling of the material occur, leaving behind voids and vacancies along the edges of the milled pattern [Figs. 5(b) and 5(d)]. This is characterized by the 62% increase in mill depth (t-test: p < 0.0001) [Figs. 1(c) and 2(c)], the 147% increase in the exposed edge surface area (t-test: p < 0.0001), the 287% increase in milled volume (t-test: p < 0.0001) [Figs. 1(c) and 2(c)], and the similar level of surface swelling height (t-test: p = 0.8108) but 23% increase in swelling volume (t-test: p = 0.0444) [Figs. 1(c) and 2(c)]. Furthermore, voids act as sinks for many remaining vacancies, causing the growth of voids and decrease in isolated vacancies in nearby regions, further driven by diffusion during the high temperature annealing used to activate emitters. Therefore, the number of optically active defects causing background continues to decrease. This is characterized by the 75% decrease in region PL (t-test: p = 0.0174) [Figs. 3(a) and 3(d)] and the 147% increase in SPE yield (t-test: p = 0.0073) from the third level [Fig. 3(c)].

The spatial distribution and relative concentrations of vacancies, complexes, and voids could be confirmed by high resolution electron microscopy characterization techniques, such as transmission electron microscopy (TEM) or scanning transmission electron microscopy (STEM), performed on cross sections sliced from hBN crystals deterministically patterned with high-energy Ga⁺ FIB. Studies that explore the large range of Ga⁺ ion energies and fluences with fine control could quantitatively measure damage thresholds that correlate with the observed morphological characteristics. Furthermore, combining these characterization techniques with computational studies, such as Molecular Dynamics (MD) and Density Functional Theory, could elucidate void nucleation and growth dynamics while also considering stable quantum emitter defect candidates that include carbon incorporated during the annealing process and allow for the creation of SPEs fabricated with this technique.

In this work, we demonstrated that morphological characteristics due to FIB exposure play an important role in the deterministic placement of SPEs by patterning SPEs in hBN with four levels of FIB exposure and characterizing surface topography via AFM. From our analysis, we found that the evolution of swelling most strongly correlates with the number of SPEs we observed and anticorrelates to PL intensity. We expect that swelling indicates volumetric voids and local amorphization of hBN layers, with these conclusions supported by SRIM calculations and the literature that employ MD modeling.^42,50,52 Therefore, we conjecture that while fluorescent vacancy defects are present at all FIB doses, the increased damage indicated by the swelling onset at the third exposure level and the milling onset at the fourth actually decreases the total number of optically active defects at a given milled site, allowing for the localization of single emitters.

Further experimental and computational studies are necessary to identify the ideal damage threshold to localize emitters in hBN without high-cost and high-precision individual defect placement techniques. For suspended structures, such as photonic waveguides and optomechanical transducers, it may be desirable to identify an ion species that can produce a higher yield of localized emitters with minimal sputtering, demonstrated at our third FIB level.^{22–24,53,54} Another possible route could be to optimize the milled edge profile for SPE creation and design photonic devices that rely on the creation of these edges. By identifying the morphological characteristics that can be leveraged to engineer hBN SPEs for quantum technologies, we have provided insight into the evolution of these defects-by-design that can be applied to the expanding field of solid-state single-photon emitter sources.

See the supplementary material for detailed experimental, analytical, and simulation methods and supporting data for pre-exposure hBN AFM, all studied morphological characteristics as a function of exposure time, distribution of ZPL, and distribution of g²(0) values.

The authors thank Brittany Carter, Uriel Hernandez, Viva Horowitz, and Valerie Brogden for discussion of this work. This work was supported by the University of Oregon, the Reneé James Seed Grant Initiative, and the National Science Foundation (NSF) under Grant No. CMMI-2128671. R.K., J.Z., D.M., K.Z., and B.A. acknowledge facilities and staff at the Center for Advanced Materials Characterization in Oregon and the use of the University of Oregon’s Rapid Materials Prototyping facility, funded by the Murdock Charitable Trust. K.W. and T.T. acknowledge support from the JSPS KAKENHI (Grant Nos. 19H05790, 20H00354, and 21H05233).

The authors have no conflicts to disclose.

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