Designing a boron nitride polyethylene composite for shielding neutrons

Neutrons are encountered in many different fields, including condensed matter physics, space exploration, nuclear power, and healthcare. Neutrons interacting with a biological target produce secondary charged particles that are damaging to human health. The most effective way to shield neutrons is to slow them to thermal energies and then capture the thermalized neutrons. These factors lead us to consider potential materials solutions for neutron shields that maximize the protection of humans while minimizing the shield mass, and which adapt well to modern additive manufacturing techniques. Using hexagonal boron nitride (hBN) as a capture medium and high-density polyethylene (HDPE) as a thermalization medium, we aim to design the optimal internal structure of h$^{10}$BN/HDPE composites by minimizing the effective dose, which is a measure of the estimated radiation damage exposure for a human. Through Monte Carlo simulations in Geant4, we find that the optimal structure reduces the effective dose up to a factor of 72x over aluminum (Al) and 4x over HDPE; this is a significant improvement in shielding effectiveness that could dramatically reduce the radiation exposure of occupational workers.


I. INTRODUCTION
Neutron radiation-encountered in healthcare, nuclear reactors, and the space environment-is extremely damaging to human health 1,2,3 .Long-term exposure results in cancer, cognitive decline, and degeneration of the circulatory system 4,5 .Neutrons interacting with a biological target produce secondary charged particles, which have a high linear energy transfer (LET, keV/µm) causing subsequent tissue damage 6 .Shielding against neutron radiation is difficult since neutrons can only be moderated or absorbed.The scattering crosssections of neutron captures are typically orders of magnitude smaller at high energies (> 1 MeV) so the neutrons need to be moderated.In addition, thermal neutrons produce damaging secondary radiation through the neutron capture reaction that typically emits prompt gammas, which are also damaging to health, so ultimately the thermal neutrons should be absorbed to eliminate the danger.
In this work, we consider the physical processes of moderation and absorption to address a basic question: To protect human health, what is the optimal distribution of thermalization and capture elements within a shielding material?Prior experimental studies have considered different internal distributions of materials (e.g.continuous blends or layered composites) for a limited number of configurations, finding an improvement in some configurations 7,8,9,10,11,12 .By addressing this question from a computational perspective using extensive Monte Carlo simulations, we explore a larger configuration space, including structures that are appropriate for modern manufacturing techniques, such as additive manufacturing (AM).The result of this work has clear applications.For a) avira@gatech.eduavira@gatech.edub) first@gatech.edufirst@gatech.eduexample, conventional high-mass radiation shielding is employed for many terrestrial applications but would be a poor choice for space environments such as the Moon, due to associated launch costs.Lightweight polymer materials and composites are increasingly employed for neutron shielding applications and are compatible with a variety of AM techniques 13 .AM capabilities have grown rapidly in recent years and are planned for autonomous lunar construction 14 .The shielding composites, explored in this study, would be of particular interest for lunar applications given the lightweight materials and compatibility with AM.
Because the capture cross-section scales inversely with velocity, reducing the neutron energy improves the capture efficiency.To accomplish these two physical processes, moderation and absorption, we select appropriate materials.Effective moderation requires a high rate of momentum transfer to the material, which is best accomplished by compounds with high hydrogen content due to the near equality of neutron and proton masses.A good moderator is capable of slowing down fast neutrons without undergoing significant radiation damage or capturing the neutrons.Water (ρ ≈ 1 g/cm 3 ), graphite (ρ ≈ 2.26 g/cm 3 ), and polyethylene (ρ ≈ 0.90 -0.97 g/cm 3 ) are commonly-chosen moderator materials.High-density polyethylene (HDPE), (C 2 H 4 ) n for n between 3,500 and 9,000, is ideal for a variety of applications due to its low density and large elastic neutron scattering cross-section.Additionally, unlike water and graphite, HDPE has suitable physical properties to be used as a stand-alone flexible and light-weight radiation shielding material across a broad range of operating temperatures, which is potentially useful for shielding neutrons on the lunar surface 15 .
Among low mass isotopes, 10 B, comprising 19.9% of natural boron, is the clear champion for neutron capture with a cross-section that reaches 3800 b for thermal neutrons 16 .The following capture reactions are relevant for 10 B: Here, 10 B captures the thermal neutrons ( 1 0 n th ), converting them into alpha particles (α) and lithium ions, which have a stopping range of only a few micrometers at the relevant energies-much shorter than the neutron attenuation length.The neutron capture reaction of Equation 1a, with emission of a 0.478 MeV gamma ray, has a ∼94% probability and the reaction of Equation 1b has a ∼6% probability.
Prior experimental studies of neutron transmission through boron-containing HDPE composites reveal an improvement in the shielding effectiveness in comparison to pure HDPE 7,12,7,9,10,11 .A high boron-containing material is more efficient at capturing neutrons.B 4 C has a higher boron content than hexagonal boron nitride (hBN); however, Shin et al.,  2014 showed that hBN can be functionalized to improve the solubility and homogeneity of the composite, which enhanced its shielding effectiveness beyond that of B 4 C composites 8 .Hexagonal boron nitride, a quasi-2D material with a large bandgap of ∼ 5.9 eV, has naturally occurring boron (19.9% of 10 B and 80.1% of 11 B, denoted as h Na BN) and can undergo the capture events shown in Equations 1a-1b.However, hBN can be boron-10 enriched (denoted as h 10 BN), further increasing the probability of undergoing the capture reaction 17 .This work provides a comprehensive characterization that directly builds upon previously published experimental works 7,12,7,9,10,11 by using an open-source Monte Carlo (MC) code, Geant4 (GEometry ANd Tracking), to optimize the shielding composition and structures for h 10 BN/HDPE composites 18,19,20 .Geant4 is a widely used and validated Monte Carlo code that offers freedom in creating the detector geometry and a range of physics models appropriate to different situations.In the context of our own research program, the choice of Geant4 was made because it also has the potential to incorporate detailed solid-state effects.
For human safety, we aim to maximize shielding effectiveness by optimizing the boron distribution within polyethylene composites using the effective dose as a figure of merit.The effective dose is a quantity that provides radiation workers with an assessment of the potential hazard to human health by accounting for the radiation effects on human organs and tissue.Defining the effective dose ratio as the effective dose with a shield divided by the effective dose without a shield, we run Geant4 simulations for three types of composites: (1) homogenous blend of HDPE and hBN ("blended"), (2) pure HDPE interlayered with pure hBN ("ideal layered", inspired by 21 ), and (3) pure HDPE interlayered with an hBN/HDPE blended composite ("manufacturable").For all configurations, the Geant4 simulations assume a homogeneous blend of the specified elements and do not consider the impact of component particle sizes or shapes, which could affect the shielding efficiency 22,23 .The influences of particle size, geometry, and distribution are topics for future research; within our present approximation, we expect that the simulations should remain valid for composites with inhomogeneities on a size scale much smaller than the neutron-capture meanfree path, which in h 10 BN is ∼ 10 µm for thermal neutrons.The manufacturable configuration could be fabricated via multi-material fused filament fabrication, an AM technique where thermoplastic filaments are deposited to build up 3D geometries 24,13 .The simulated incident neutrons have a loguniform energy distribution from 15 meV to 20 MeV.During post-processing, we separate the neutrons by different energy regimes, discussed in Section III C, to isolate the physics occurring within each regime 25 .This also widens the applications of hBN/HDPE shielding composites by providing results that can be applied generally to shielding neutrons for an arbitrary distribution of incident energies.We find that alternating layers of capture and moderation materials provide improved radiation protection compared to blended composites with the same composition.From the manufacturable configurations, we find that we can improve the shielding effectiveness up to a factor of 72 relative to Al and up to a factor of 4 with respect to HDPE.This is a significant improvement that could dramatically reduce the occupational radiation dose for workers in high-risk environments, for example, astronauts or healthcare professionals.
Section II discusses the figure of merit that is used to evaluate the shielding effectiveness of all the simulated shielding composites.The structure of the Geant4 simulations is discussed in Section III.We conclude the paper by discussing the results for the three different geometries for the hBN/HDPE composites (blended, ideal layered, and manufacturable) and the corresponding optimal configurations in Sections IV and V.

II. FIGURE OF MERIT: EFFECTIVE DOSE
There are three dose quantities, in radiation physics, that are used to assess radiation exposure: absorbed dose, equivalent dose, and effective dose.The absorbed dose is defined as the mean energy imparted by ionizing radiation in matter per unit of mass (the SI unit of absorbed dose is the Gray, with 1 Gy = 1 J kg −1 ).The equivalent dose (H T ) weights the absorbed dose by the radiation type to determine the effect of the ionizing radiation on a specific organ or tissue.The effective dose (E), provided in units of Sieverts (Sv), takes the equivalent dose and sums over the organ/tissue.Mathematically, the effective dose reads where w R is the radiation weighting factor, w T is the tissue weighting factor, and D T,R is the mean absorbed dose from a radiation type, R, in specific organs or tissues, T .The International Commission on Radiological Protection (ICRP) provides a set of protection quantities based on radiation transport simulations for a monoenergetic particle emitted in various directions (rotational, left lateral, right lateral, isotropic, etc.) incident on a human phantom.ICRP provides FIG. 1. Effective dose per fluence (pSv cm 2 ) for a sex-averaged human versus the incident energy, calculated for various radiation types emitted isotropically.ICRP 116 and 123 are used to obtain the conversion factors from energy to effective dose 26,27 the necessary coefficients to convert incident energy to effective dose for sex-averaged models of human phantoms.These coefficients are used to estimate the effective dose an average human receives for a given radiation type and energy.Here, we use the conversion coefficients for an isotropic input distribution as it more accurately represents the radiation in space.ICRP 103 provides the radiation weighting factors 5 ; ICRP 116 and 123 provide the conversion coefficients for a given fluence (Φ) 26,27 .
Figure 1 shows the effective dose per fluence (E/Φ) as a function of the incident energy for different radiation types.The heavy ions have the largest effective dose and can, therefore, be the most damaging for humans compared to other ionizing particles.Gamma rays, positrons, electrons, and protons have a similar effective dose and are also important to shield against, especially for energies above ∼ 1 MeV.Neutrons are damaging across all incident energies, showing the importance of incorporating 10 B into the shielding components.

III. STRUCTURE OF GEANT4 SIMULATIONS A. Design Space for hBN/HDPE Composites
We restrict our design space to address the use of h 10 BN and HDPE to perform the two-step shielding process, neutron capture and thermalization, respectively.Our goal is to explore the design space (blended, ideal layered, and manufacturable) of h 10 BN/HDPE composites using the effective dose as a figure of merit to obtain the optimal boron distribution in HDPE composites.
Blended Geometry.For the blended geometry, two free design parameters are considered: total composite thickness (t) and global weight percentage of h 10 BN (ω global ).The thickness is varied between 1 mm and 150 mm (specifically, 1,3,5,7,10,20,30,60,100,150 mm) to simulate a shielding material for various applications, such as integration into personal protective equipment or lining for habitat shielding.The global weight percentage, ω global , is varied between 0 and 100%, in increments of 10%.The range of weight percent-age is simulated to cover the entire design space of all possible blended configurations even though blended composites with > 30 wt% of h 10 BN are difficult to manufacture via traditional blending methods due to poor particle dispersion and dramatically increased viscosity of the melt at higher loadings, particularly when hBN nanoparticles are used in the manufacture of blended composites 28,29,24 .Including higher weight percentages of h 10 BN, beyond the range of manufacturable composites, also allows this dataset to be used to evaluate the cost-benefit relationship for the implementation of advanced manufacturing methods which could be used to produce composites with higher weight percentages of h 10 BN.For each thickness, a set of 13 simulations are run to cover the design space of ω global , resulting in 130 unique simulations for the blended geometry.
Ideal Layered Geometry.For the ideal layered geometry, we consider three parameters: total composite thickness (t), number of HDPE and h 10 BN layers within a composite thickness (where n is specifically the number of h 10 BN layers), and global weight percentage of h 10 BN (ω global ).This configuration consists of interlayering pure HDPE and pure h 10 BN layers to explore how spatial separation affects the effective dose compared to a blended structure.We explore this parameter space of interlayering pure h 10 BN to evaluate if there is a benefit to the spatial separation, even though it is difficult to make a pure hBN layer of large size due to the restricted growth conditions required to obtain a pure hBN crystal 30 .From the capture reaction in Equations 1a-1b, it is best to terminate the layered composite with an HDPE layer, which is the same thickness as the other pure HDPE layers, so that the heavy charged particles-which would significantly increase the effective dose (see Figure 1)-can be stopped before exiting the composite shielding.For this reason, there is always one less h 10 BN layer than HDPE layers.The number of h 10 BN layers is varied between 1 (2 layers of HDPE) and 99 layers (100 layers of HDPE) where the number of layers is equivalent to the number of periods.The thickness is varied over the same range of thicknesses used for the blended geometry.The global h 10 BN weight percentage, ω global , changes the relative thicknesses of the h 10 BN layers within the internal structure and varied between 5-95 wt% to cover the entire parameter space (where ω global = 0 is a pure HDPE block and ω global = 100 is a pure h 10 BN block).Based on this design space, there are 1210 unique simulations (10 thicknesses, 11 different ω global , and 11 values of n) for the ideal layered geometry.
Manufacturable Geometry.Due to manufacturability constraints of creating thick continuous hBN layers bonded to HDPE layers in the ideal layered structure and high loadings in some blended structures, it is difficult, or impossible in some cases, to manufacture all the geometries that are explored in the blended and ideal layered geometries 30 .To address this issue, additional simulations were run in which pure HDPE was interlayered with h 10 BN/HPDE blended composite layers with ≤ 20 wt% of h 10 BN.A weight percentage of ≤ 20% is chosen as prior experimental studies have shown little change in mechanical or rheological properties of hBN/HDPE composites at these loadings 24 .Additionally, given the min-imal influence of ≤ 20 wt% on the rheological properties of hBN/HDPE composites, it is reasonable to infer that these manufacturable geometries can be adapted for fused filament fabrication, an additive manufacturing (AM) process 13 .Using common dual extrusion technologies, the blended feedstocks, comprising of hBN nanoplatelets blended into HDPE, are interlayered with pure HDPE feedstocks to build up the 3D geometry where the layer thickness is easily controlled by simple programming.Alternatively, these materials could be manufactured into layered films via common film forming and lamination processes to create large surface area composite sheets suitable for various applications, such as habitat construction for lunar missions.
For the manufacturable geometry, there are four design parameters: total thickness of composite material (t), weight percentage of the blended h 10 BN/HDPE layer (ω global ), local weight percentage of h 10 BN within the blended layer (ω local ), and number of mixture layers (n).Refer to Figure S1, in the supplemental material, for a plot of the interdependence of these design parameters on the total areal density.For a subset of simulations, we explore the decrease in shielding effectiveness by including h Na BN instead of h 10 BN, which is cheaper and more available 17 .The thickness of each mixture h 10 BN/HDPE layer (t mix ) is where ρ total is the total density of the composite material, and ρ mix is the density of the mixture layer.The thickness of each HDPE layer (t HDPE ) is Equations 3 and 4 can be reduced to the ideal layered geometry by assuming ω local is 100% and ω global controls the relative thickness between the pure HDPE and pure h 10 BN layers.For the manufacturable configurations, t is varied between 1 mm and 150 mm using the same thicknesses in the blended configuration, ω global is varied between 10 wt% and 95 wt%, n is varied between 1 and 9 in increments of one, and ω local is varied between 1 wt% and 20 wt%.Based on the large design space, we run 5,670 simulations for the manufacturable configuration alone, which produces a few terabytes of data (see Appendix C for a description of the data reduction).
Reference Geometry.Simulations of pure aluminum, lead, and HDPE block are carried out to provide a direct comparison with commonly used radiation shielding materials.These simulations are run for similar areal densities of the h 10 BN/HDPE composites to allow for a direct comparison.

B. Implementation of Simulations
Five million neutrons, isotropic in angle, are uniformly incident on one side of the composite material with energies ran-domly sampled from 15 meV to 20 MeV with a log-uniform distribution.Within each Geant4 simulation, a "sensitive detector" is placed on the entry and exit surface of each composite structure to record the particles entering and exiting the material without saving all the interactions occurring within it (see Appendix A for details of Geant4).The composite material is set to have periodic boundary conditions (using the g4pbc module) in the lateral dimensions to simulate an infinite slab.The lateral dimension for all the simulations is 100 cm, which is ∼10x larger than the thickest sample.We choose such a large lateral dimension to ensure that the periodic boundary conditions do not have a significant effect on the simulation results.

C. Post-processing of Geant4 Results
The incident neutrons are separated into distinct energy bins to isolate the physics driving the neutron interactions and widen potential applications of this work by allowing optimization of hBN/HDPE composite for specific energy distribution (e.g.albedo lunar neutrons or neutron reactors).The neutron energy bins, shown in Table I, are fairly standard bins 25 .Within each energy range, the physics behind the dominant neutron interactions leads to their names.
In thermal equilibrium, neutrons approach a Maxwell-Boltzman distribution with the peak near 0.025 eV.Neutrons are referred to as cold neutrons below the peak, thermal neutrons at the peak, and epithermal neutrons above the peak.At an energy between 0.4 eV and 0.6 eV, neutrons have a large capture cross-section with 113 Cd and, therefore, referred to as cadmium neutrons.Above the cadmium energy range is referred to as epicadmium, which ends at 1 eV since the capture cross-section drops by an order of magnitude.After the epicadmium range, neutrons from 1 eV to 10 eV are referred to as slow neutrons since nuclear resonance typically begins at 10 eV.Above 10 eV, there are typically resonances in the capture cross-sections.The cutoff energy for resonance neutrons is not well-defined but typically stops at 300 eV.Between 300 eV and 1 MeV, there are nuclear reactions that produce other particles, such as protons and alpha particles; this range is called the intermediate energies.Fast neutrons cut off at an energy of 20 MeV since the corresponding velocity is 20% of the speed of light.To reduce the number of bins and increase the signal-to-noise ratio, we combine the following energies: thermal with cold, cadmium with epithermal, and slow with epicadmium.

A. Blended Geometry
We explore the shielding effectiveness of a homogeneous distribution of boron in HDPE by varying ω global from 0% to 100% in 10% increments using the effective dose as a figure of merit (described in Section II).We optimize the effective dose ratio ((E/Φ) out /(E/Φ) in ) by calculating the effective  dose with and without the shielding material.Figure 2 shows the effective dose ratio as a function of areal density (g/cm 2 ) for all h 10 BN/HDPE blended configurations.In both panels, the marker size increases linearly as ω global increases from 0% to 100% (the smallest marker corresponds to a pure HDPE block, largest marker corresponds to a pure h 10 BN block), and each color corresponds to a unique energy range, shown in Table I, based on the neutron energy entering the composite.
In Figure 2a, the shielding efficiencies of the h 10 BN/HDPE composites changes with the incident neutron energy and areal density.For example, in the epithermal range, the optimal structure for areal densities below ∼ 2 g/cm 2 is a pure h 10 BN block (corresponding to the largest marker size).This result indicates the effective dose ratio is not affected by the thermalization process and, instead, the neutron capture with 10 B is the dominant physical process that contributes to the decrease in the effective dose.However, as the areal density increases above ∼ 2 g/cm 2 , a pure HDPE block, the smallest marker, would maintain the lowest effective dose ratio (note that the calculations are for discrete thicknesses), showing that the optimal structure changes with areal density.Similar behavior is observed for cold, epicadmium, and resonance neutrons.
Figure 2b shows the effective dose ratio versus areal density for resonance, intermediate, and fast neutrons.The behavior is different for intermediate and fast neutrons since the op-timal structure is a pure HDPE block (smallest marker size) across all areal densities < 300 g/cm 2 .This indicates that the effective dose ratio is reduced the most when all of the available material is devoted to theramlization because: (1) cold and thermal neutrons contribute less to the overall effective dose than more energetic neutrons (Figure 1) and ( 2) at lower neutron energies the probability of undergoing a capture event increases since the capture cross-section is inversely proportional to velocity.

B. Ideal Layered Geometry
To further investigate the role of internal structure, we consider the next simplest modification, material layering, which separates the two essential processes of neutron thermalization and neutron capture into distinct repeating regions.We first consider a complete separation, where the layers are either HDPE or h 10 BN -a situation we refer to as the ideal layering.Figure 3 shows the effective dose ratio as a function of areal density, within the epithermal energy range, for three example composite thicknesses: 7 mm in (a), 30 mm in (b), and 100 mm in (c).The different colors correspond to a unique number of h 10 BN layers, or equivalently, the spatial period of the repeating structure.As a reference, the blended results are duplicated from Figure 2 (shown in black stars).The individual thicknesses of the pure HDPE and h 10 BN layers can be calculated using Equations 3 and 4 by setting ω local = 100% where t mix would, instead, refer to the thickness of pure h 10 BN.
From Figure 3, we find the ideal layered geometry provides better shielding than a homogeneous blend.For t = 7 mm, the one period trilayer structure (h 1 0BN between two HDPE layers) minimizes the effective dose across all areal densities and results converge to the blended configuration as the number of layers reaches 50+.However, for t = 30 mm, the trilayer configuration minimizes the effective dose ratio only above ∼ 4 g/cm 2 .At t = 100 mm, the optimal structure has six to eight periods versus one period for t = 7 mm, as shown in Figure 3c and 3a.This is the first indication that the spatial scale and organization of thermalization and capture materials within a composite shield can significantly affect the macroscopic shielding properties.

C. Manufacturable Geometry
Given that the ideal layered geometry reduces the effective dose ratio, we explore a design space that is feasible to manufacture, as described in Section III. Figure 4 shows the effective dose ratio versus areal density for the manufacturable configuration within the epithermal energy range.Figures 4a-4c show the results for t = 7 mm, t = 30 mm, and t = 100 mm, respectively (similar to Figure 3).The color bar in Figure 4 is split into three sections that show the ω local at 5% (shades of brown), 10% (shades of green), and 20% (shades of purple).The different shades of each color correspond to the number of h 10 BN/HDPE mixture layers between pure HDPE layers.For the manufacturable configurations, the optimal configuration occurs at ω local = 20% (shown in the shades of purple) for each thickness shown in Figure 4 (this is also observed for all other thicknesses explored in this study).
A comparison of Figure 3 with Figure 4 shows that the effective dose ratios are similar for all three thicknesses, indicating that we are able to reproduce the shielding efficiency of the ideal layered structure by optimizing the manufacturable composite.Similar to the crossover in the optimal number of layers in Figure 3, the optimal manufacturable configuration for 7 mm consists of a trilayer structure at ω local =20%.As the total thickness of the composite increases to 100 mm, the effective dose ratio is minimized when there are five to eight mixture layers (ω local =20% and areal density of ∼ 10.3 g/cm 2 ).This crossover in the manufacturable results indicates that the physical process driving the change in the number of layers is similar to that occurring in the ideal layered geometry.Figure S2 provides a comprehensive characterization of this crossover feature for the manufacturable configuration.
Since h 10 BN is difficult and expensive to manufacture, we run another set of simulations of HDPE interlayered with  S3, are compared with the performance of the 10 B enriched manufacturable configurations.From this, we find that the effective dose ratio decreases up to a factor of ∼ 1.6 when using h Na BN instead of h 10 BN within the mixture layers.This decrease in performance varies with energy range and areal density (see Figure S3).However, it is important to note that, at specific neutron energies (cold, epithermal, intermediate) and large areal densities (> 5 g/cm 2 ), there is not a significant improvement in the shielding efficiency by using 10 B enriched hBN.

V. IMPLICATION OF GEANT4 RESULTS
The optimal distribution of 10 B within HDPE can be determined by comparing the performance of the different geometries described in Section IV.From the MC simulations, we found that separating the thermalization and capture process (ideal layered configuration) provides better shielding than a composite blend, which led us to explore the manufacturable configuration.Figure 5 shows the effective dose ratio versus areal density for each energy range in Table I show the 2σ variation from running the Geant4 simulations with five different random-number seeds.
For cold neutrons (E ≤ 0.025 eV), the composition of the optimal blended and manufacturable configurations changes with areal density.Below ∼ 0.3 g/cm 2 , the optimal blended configurations include small loading of 10 B, demonstrating that it is beneficial to capture the cold neutrons (also see Fig- ure 2).As the areal density increases, the optimal blended configurations fall along the HDPE curve, indicating that the optimal structure is entirely HDPE.Cold neutrons have a very high probability of undergoing a capture event, and most capture events produce prompt gammas, which also contribute to the effective dose (shown in Figure 1).This tradeoff of cold neutrons for gammas seems to be inefficient for the blended configuration for areal densities above ∼ 0.4 g/cm 2 .However, the manufacturable configuration has a better shielding performance than the blended configurations, reducing the effective dose ratio by an additional factor of 1.12-1.67across all areal densities.Similar to the blended configurations, the internal composition changes with areal density for the manufacturable configuration, from trilayer to multilayer composite (also shown in Section IV C).This cross-over transition occurs at an areal density of ∼ 5 g/cm 2 (see Figure S2 in the supplemental material for a comprehensive characterization of the transition for all energy ranges).For both the blended and layered composites, the optimal configurations reduces the effective dose ratio up to a factor of 41x over Al and 2x over HDPE, which is a significant improvement.
As the neutron energy increases to the epithermal range (0.025 eV < E ≤ 0.6 eV), the optimal blended configuration no longer lies along the reference curve for HDPE, indicating that the inclusion of 10 B is important for reducing the effective dose ratio.Similar to cold neutrons, the change in the internal configuration from trilayer to multilayer occurs around an areal density of ∼ 5 g/cm 2 , indicating that the driving physics is likely to be similar for both energy regimes.For areal densities above ∼ 0.1 g/cm 2 , the optimal manufacturable configurations reduce the effective dose more than the blended configurations, resulting in an improvement of up to 1.41x over the blended configurations.At ∼ 0.1 g/cm 2 , the effective dose ratio is lower for the blended composite than a layered composite, indicating that there is no benefit to separating the thermalization and capture mediums.For the composites explored in this study, the optimal configurations reduces the effective dose ratio up to a factor of 21x over Al and 3x over HDPE.
For epicadimum neutrons (0.6 eV < E ≤ 10 eV), the optimal blended configurations include 20 wt% of h 10 BN for areal densities below ∼ 0.7 g/cm 2 and between 5 wt% and 20 wt% for areal densities above ∼ 0.7 g/cm 2 .In contrast to cold and epithermal neutrons, the manufacturable configuration is not the ideal structure for all areal densities -below ∼ 0.7 g/cm 2 the blended configuration reduces the effective dose up to 1.20x over the optimal manufacturable configurations, and above this areal density the optimal layered composite reduces the ratio up to 1.30x over the blended.Above ∼ 0.7 g/cm 2 , the optimal manufacturable configuration changes from a trilayer to multilayer configuration around an ∼ 5 g/cm 2 , which is the same cross-over point as the cold and epithermal neutrons.This interesting result shows the importance of optimizing the design space for h 10 BN/HDPE composites as the arrangement of the moderation and capture medium is not obvious and can result in an improvement of 20x and 3x over Al and HDPE, respectively.
The behaviour of resonance neutrons (10 eV < E ≤ 300 eV) is similar to epicadimum neutrons: (1) the blended configurations improve the effective dose ratio up to a factor of 1.19 over the optimal manufacturable configurations at areal densities below ∼ 2 g/cm 2 and (2) above ∼ 2 g/cm 2 , the structure of the layered composites changes from trilayer to multilayer as the areal density increases.The optimal configurations reduce the effective dose up to a factor of ∼ 20x over Al and ∼ 4x over HDPE for resonance neutrons.All three of these energy ranges (epithermal, epicadmium, resonance) show that there is a clear improvement of including h 10 BN into the HDPE composite.For all the areal densities, the layered composite seems to be the best way to distribute 10 B but there are a few instances, at small areal densities, where it is beneficial to have a homogeneous blend of h 10 BN and HDPE.
For the intermediate range (300 eV < E ≤ 1 MeV), the blended and manufacturable composites have the same performance as pure HDPE until the areal density increases to ∼ 3 g/cm 2 .Above ∼ 3 g/cm 2 , the optimal blended or manufacturable configurations have the similar shielding efficiencies, improving the effective dose by 79x over Al and 4x over HDPE.The improvement over Al is larger than the other energy ranges because neutrons are significantly more damag- Multilayer Blended FIG. 6. Summary of the change in internal structure for the epicadmium energy range as areal density increases.The blue and green circles provide a representation of HDPE and hBN concentrations, respectively, and the solid dark blue is pure HDPE.The optimal structure changes from blended (below ∼ 0.7 g/cm 2 ) to trilayer (between ∼ 0.7 g/cm 2 and ∼ 5 g/cm 2 ) to multilayer (above ∼ 5 g/cm 2 ).
ing to human health within this energy range, as shown in Figure 1, resulting in a larger reduction in the effective dose ratio upon thermalizing and capturing these hazardous neutrons.Thus, above ∼ 3 g/cm 2 , there is a benefit to including 10 B within the composite as the intermediate neutrons are able to undergo elastic interactions with HDPE and increase their probability of 10 B capture events.However, within the range of thicknesses studied, there is no benefit to using a blended or layered composite for this energy range.For fast neutrons (1 MeV < E ≤ 20 MeV), there is no benefit to including 10 B across all areal densities explored in this study due to the large thickness of HDPE (> 250 mm) that is required to thermalize fast neutrons.
Across all energy ranges (except for fast neutrons), we discover a dependence where the optimal manufacturable configurations change from a trilayer to a multilayer structure with areal density.This interesting feature could be due to the trade-off of neutrons for gammas and the attenuation of the gammas through the composite.However, there are specific regions where the blended configurations have a better shielding performance than the manufacturable configurations (epicadmium neutrons below ∼ 0.7 g/cm 2 and resonance below ∼ 2 g/cm 2 ).A summary of the change in internal structure for the epicadimum neutrons is shown in Figure 6 as a function of areal density.This change in internal structure shows the importance of simulating over the design space of h 10 BN/HDPE composites as it is experimentally time-consuming to fabricate and test different configurations.Future work involves understanding the underlying physics driving such a change by developing an analytical model.
We note that there are additional methods to further improve the shielding effectiveness for the materials explored in this study.The effective dose could be further reduced with the addition of high-density material, such as bismuth or barium titanate, at the backend of the composite shielding materials to reduce the gammas that penetrate the shielding material and subsequently contribute to the effective dose.However, the addition of a high-Z material layer to attenuate the photons is outside the scope of this work and could affect the inherent flexibility of the h 10 BN/HDPE composites.Considering real shielding solutions, the longevity of the hBN/HDPE composites could be affected by accumulated energy from secondary gamma production.

VI. CONCLUSION
In this study, we use Monte Carlo simulations to design the boron distribution within HDPE composites for shielding neutrons to address a fundamental question: what is the optimal distribution of moderation and capture elements within a composite?We explore three different configurations: blended, ideal layered, and manufacturable.We use a log-uniform neutron source from 15 meV to 20 MeV, separated into standard neutron energy ranges, to provide the flexibility to optimize the composite structures for specific neutron applications (e.g.reactor with specific proportions of epithermal and intermediate neutrons).We find that, depending on the incident neutron energies and areal density, the optimal design of hBN/HDPE composites reduce the effective dose ratio by a factor of 5-79x in comparison to Al and 1.5-4x in comparison to HDPE.This is a significant improvement in shielding effectiveness that could dramatically reduce the radiation exposure occupational workers receive.
From optimizing the internal structure of hBN/HDPE, we discover that the optimal structure changes with areal density and incident neutron energy.In almost all cases, the layered composite reduces the effective dose ratio more than a homogeneous blend.There are some cases where a homogeneous blend is sufficient -composites with an areal density of < 0.7 g/cm 2 and < 2 g/cm 2 for epicadmium and resonance neutrons, respectively.For all other cases, the layered composite is the most efficient way to reduce the effective dose ratio.The layered structure changes from a trilayer (two HDPE layers around one blended h 10 BN/HDPE layer) to a mulitlayer structure as the areal density increases, which is another interesting feature.By exploring a comprehensive design space for hBN/HDPE composites, we learn that the answer to our original question of the best way to distribute moderation and capture elements is nuanced but it is generally beneficial to incorporate h 10 BN within the HDPE shield.The constraints placed on the manufacturable configuration ensure that the explored design space is compatible with additive manufacturing techniques, a preferred tool for space applications.Moreover, the HDPE and h 10 BN/HDPE blends used in the manufacturable configuration are inherently flexible, allowing for easy incorporation into the linings of personal protective equipment.Future work involves fabricating and irradiating the optimal manufacturable configurations to confirm shielding effectiveness in real-world composites in addition to incorporating realistic particle distributions into Geant4 simulations.

SUPPLEMENTARY MATERIAL
See the supplementary material for Figure S1, illustrating the interdependence of the design parameters in the manufacturable configuration, Figure S2, illustrating the transi-tion from trilayer to multilayer for all neutron energy regimes (cold, epithermal, epicadmium, etc.), and Figure S3, illustrating the improvement factor of using h 10 BN instead of h Na BN within the mixture layer of the manufacturable configurations.on a single-particle basis.For each hit, the kinetic energy, charge, and position vector can be recorded, providing the user with detailed information of the particle's trajectory through a material.This method allows the user to record all the relevant information, then perform post-processing on the simulation results.Alternatively, the user can define specific variables to tally within the code (e.g., deposited energy, surface flux, number of particles), which is useful when there is a specific variable of interest.For the simulations presented in this study, we record the following information associated with each hit: (1) particle type (e.g.neutron, gamma, alpha, etc.), (2) process type (e.g.inelastic collision, elastic collision, ionization, transport, etc.), (3) particle kinetic energy, and (4) particle position in Cartesian coordinates.From this information, the exact interaction and energy deposited into the material can be determined and tracked through the particle cascade.However, to get valid results, it is critical to use the appropriate physics models.For these simulations, we use the following physics constructors: a modified G4HadronElasticPhysics to include NeutronHPPhysics, G4HadronPhysicsShielding, G4EmStandardPhysics_option4, G4EmExtraPhysics, G4IonPhysics, G4DecayPhysics, G4RadioactiveDecayPhysics, and G4StoppingPhysics.

Specifically,
we modify the original G4HadronElasticPhysics to include the NeutronHPPhysics constructor, which computes the phonon density of states for both the coherent and incoherent part and has been verified and found to have a reasonable agreement with other MC simulation codes along with experimental data 31,32,33,34,35 .The G4HadronElasticPhysics models the elastic physics of particles, as opposed to the G4HadronPhysicsShielding, which models inelastic physics.Both of these constructors are required to obtain accurate descriptions of the interactions occurring for all particles, including thermal neutrons.The other physics constructors are used to account for the following physics: electromagnetic effects (G4EmStandardPhysics_option4, G4EmExtraPhyics), ion interactions (G4IonPhysics), decay channels according to the branching ratios (G4DecayPhysics), radioactive decay of isotopes (G4RadioactiveDecayPhysics), and nuclear capture at rest for negatively charged particles (G4StoppingPhysics).
To ensure that we have loaded the appropriate physics constructors, we compare the Geant4 results to another MC simulation code, Monte Carlo N-Particle code (MCNP), a well-established neutron transport code.
at the center of the sphere and emits neutrons isotropically.For HDPE, we use a density of 0.968 g/cm 3 and a neutron energy of 2.5 MeV so that the thermalization process is observed within the sphere.For h Na BN, a density of 2.1 g/cm 3 and a neutron energy of 100 keV so that the capture events are observed.For the MCNP simulation, the ENDF/B-VIII.0cross-section database and thermal scattering card are used.The surface current is tallied at each spherical shell, which is placed at 5 mm increments starting from the center of the sphere and propagating to the surface.

Calculation of Surface Current Tally
We use the surface current tally, denoted as F1 in MCNP, to compare the Geant4 and MCNP simulation results.This tally counts each neutron that crosses a specified surface and weights that neutron by the angle at which it exits the surface.The surface current tally is mathematically represented as where A is the specific surface area, n indicates its normal direction, r s is its distance to the neutron source, E is the neutron energy, and J is the current vector.By computing the absolute value of n • J, MCNP does not distinguish between forward and backward scattering 36 .The same calculation is performed in Geant4 during the post-processing of the simulation results.The entire sphere is set as the SD and the direction is weighted by | n • J|, replicating the process shown in Equation B1.For simplicity, the F1 tally is calculated for neutrons only in both Geant4 and MCNP.

Comparison Results
Figure 7a (top panel) shows the calculated neutron surface current, F1, at a shell radius of 30 mm, from MCNP6 (red) and Geant4 (black) for an isotropic source at 2.5 MeV located at the center of an HDPE sphere with a radius of 100 mm.Here, MCNP is run for 108 neutrons, and Geant4 is run for 250,000 neutrons.We use a smaller number of neutrons in MCNP as the computational time is longer when recording each individual interaction within the sphere.For direct comparison, the surface current is normalized by the number of incident neutrons, n 0 incident .The difference between MCNP and Geant4 surface currents, at a shell radius of 30 mm, is shown in Figure 7a (bottom panel), where we observe a small discrepancy between the two MC simulation codes.In Figure 7b, the difference between MCNP and Geant4 is shown for each radial shell taken at 5 mm radius increments.The data are plotted with a diverging color bar, where white indicates that the two MC codes give the same results, red indicates that MCNP has a larger surface current than Geant4, and blue indicates that Geant4 has a larger surface current than MCNP.From Figure 7b, we observe that MCNP and Geant4 results are in reasonable agreement with each other for all radial depths within the sphere Exit Distributions (within 1% error between the two codes).After performing the same calculation for an hBN sphere and 100 keV isotropic source, we also find that the MCNP and Geant4 results are in reasonable agreement (within 1% error), indicating that the physics in Geant4 is appropriate.

Appendix C: Data Reduction
For each configuration, a separate executable is created for the user to specify the t, ω local , ω global , and number of layers (not applicable for the blended configuration).The Geant4 code is designed to track the following variables on the input and exit surface of the shielding composite: (1) particle number, type of particle (e.g.neutron, gamma, electron, etc.), type of physics interaction that occurred (i.e.elastic scattering, Compton scattering, etc.), energy of the particle, and the position of a particle (x,y,z).Geant4 outputs these relevant variables to ROOT files, a convenient file format for CERNdeveloped code 38 .The post-processing of the Geant4 results is handled in Python.Due to the large parameter size, described in Section III, we utilize PACE, a high performance computing cluster at Georgia Tech, to run ∼ 20, 000 simulations in parallel.
The ROOT files are read into Python using pyROOT where the ROOT files are converted into feather files, a lightweight binary format.All the necessary information for postprocessing is saved into the feather files, including the energy of each particle at the entry and exit surface along with the particle type and interaction type.To reduce the overall size of each feather file, the data from ∼ 100 ROOT files are saved into a single feather file.The feather files reduce the overall size of the stored data files from ∼ 4 TB to ∼ 300 GB.This is a significant reduction without losing any information from the ROOT files.
The feather files are used to create pickle files that contain the final post-processed data.The neutrons are binned by their energy, shown in Table I.The binned data is saved to a pickle file, a convenient way to store data in Python 39 .Figure 8a shows the log-uniform incident neutron source binned by the neutron energy ranges.The distribution at the output surface of the composite is determined using the particle number for each energy range, as shown in Figure 8b.This final step in the data reduction uses all the data from the feather files to create arrays with the relevant figure of merit.This last step reduces the size from ∼ 300 GB to ∼ 2 MB, which is a manageable size to work with on a local machine.

FIG. 2 .
FIG. 2. Effective dose ratio ((E/Φ) out /(E/Φ) in ) versus areal density (g/cm 2 ) for blended h 10 BN/HDPE composites.In both panels, each color corresponds to a different energy range and the marker size increases with increasing ω global .(a) shows the results for neutron energies below 10 eV and (b) shows the results for neutron energies above 10 eV.The optimal composite changes with areal density for energies below the intermediate range; above the intermediate range, a pure HDPE block is the most effective configuration.

FIG. 3 .
FIG.3.Epithermal Range: Effective dose ratio versus areal density for layered h 10 BN/HDPE composites with a total thickness of 7 mm in (a), 30 mm in (b), and 100 mm in (c).In all panels, the color corresponds to a different number of periods, as shown in the colorbar.The blended results are shown in black stars.The ideal layered structure reduces the effective dose ratio more than the blended composite.In addition, the optimal configuration changes from a trilayer structure to a multilayer structure as the total thickness increases to 100 mm.

FIG. 4 .
FIG. 4. Epithermal Range: Effective dose ratio versus areal density for layered manufacturable composites with a total thickness of 7 mm in (a), 30 mm in (b), and 100 mm in (c).The manufacturable data in each panel are split into three categories, ω local = 5% (shades of brown), ω local = 10% (shades of green), and ω local = 20% (shades of purple).The different shades correspond to the different numbers of mixture layers or the spatial period of the repeating structure.For each thickness, the largest loading of h 10 BN of ω local = 20% reduces the effective dose ratio the most.

FIG. 5 .
FIG.5.Effective dose ratio as a function of areal density for the blended (homogenous blend of HDPE and h 10 BN) and manufacturable configurations (interlayering of HDPE with h 10 BN/HDPE).The colorbar shows the different number of periods within the manufacturable configuration.Each panel shows the dose ratio for each energy range explored in this study (shown in TableI).The reference data for HDPE, Al, and Pb are provided in light gray, coral, and yellow, respectively.Error bars represent the 2σ error determined by the statistical variation of five random-number seeds, shown in the respective colors for the blended and reference geometries and in black for the manufacturable configuration.
FIG. 7. (a) Top panel: MCNP (red) and Geant4 (black) neutron surface current comparison at a shell radius of 30 mm from an isotropic source at 2.5 MeV located at the center of an HDPE sphere of radius 100 mm.Bottom panel: The difference between MCNP and Geant4 in the top panel.(b) Stacked plot of the difference between MCNP and Geant4 as a function of energy.The difference is plotted on a diverging color bar, where white indicates that the MCNP and Geant4 surface current values are equal.The surface current is normalized to the total number of incident neutrons.
FIG. 8.(a) Incident log-uniform distribution separated by neutron energy bins, provided in TableI.The output distribution for each energy bin for all particles is shown in (b).Both of these subplots are for a blended composite with a total thickness of 1 mm and ω global =10%.

TABLE I .
25utron energy bins (determined from25 . The optimal blended configurations, determined by the lowest effective dose ratio for ω global ≤ 20%, are shown as purple dotted lines.The manufacturable configurations are shown for an example of ω global =50% and ω local =20%; we note that the results are similar for other values of ω global .For the manufacturable configuration, the colors ranging from blue to orange correspond to the number of periods within each composite, as shown by the colorbar.Reference geometries are shown for comparison: HDPE in light gray, Al in coral, and Pb in yellow.The error bars (generally smaller than the symbol sizes)