THE DEVELOPMENT OF A NOVEL DIAMOND-BASED NEUTRON DETECTOR AND QUANTUM COLOR CENTER FABRICATION FRAMEWORK By Henry Matthew Thurston A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics — Doctor of Philosophy Computational Mathematics, Science and Engineering — Dual Major 2023 ABSTRACT This work investigates the use of single crystal diamond in sensing technologies. First, homoepitaxial chemical vapor deposition (CVD) grown semiconducting diamond is used to build an ultra-fast neutron detector. Diamond based neutron detectors are extant technol- ogy, but typically limited to neutron energies of less than 14 MeV. This work introduces the Positionally Opposed Schottky Semi-Metal (POSSM) solid state neutron detector. The POSSM device employs two boron-doped P-type semiconducting diamonds in conjunction with a lanthanide foil to create a pair of Schottky junction diodes with a shared cathode. Under reverse bias the diamond-Schottky diodes have undetectable reverse bias leakage cur- rent, resulting in a detector with excellent signal-to-noise properties. Preamplifier circuitry has been designed to exploit the favorable properties of the Schottky architecture. Field test- ing of the device at the Los Alamos Neutron Science Center (LANSCE) yielded successful ultra-fast neutron detection with excellent charge conversion efficiency. Several drawbacks were identified in the performance of the POSSM detector, mainly involving durability of the diodes and speed of the preamplifier. An improved design was developed and is presented in this work, though the improved design has not been built. Secondly, this work presents a novel framework for the simulation of quantum color cen- ter formation in diamond. Diamond color centers have been shown to have unique quantum spin properties making them of interest to quantum computing and field sensing. A meso- scale reaction-diffusion framework is developed and a computational solver is built. The Color Center ANnealing And Reaction-Diffusion (CCANARD) program solves the nonlin- ear reaction-diffusion system by linearization using the Gateaux Derivative and the Crank- Nicolson method. CCANARD is benchmarked for computational efficiency and accuracy. The results are presented and analyzed. An experiment is proposed to further test and develop CCANARD. Copyright by HENRY MATTHEW THURSTON 2023 I dedicate this work to my fellow Conservative students and all others who go against the grain of academic popularism in the pursuit of Truth. Your contributions have, and continue to, advance the understanding of mankind. v PREFACE This dissertation covers two largely distinct projects, both involving the development of diamond-based technology. In order to present this research in the most efficient manner possible, this document is organized into four main sections. Chapter 1 will introduce di- amond as a material, covering the physical properties and natural and artificial diamond synthesis. Chapters 2 through 4 will cover the development of a solid state diamond neutron detector. Chapters 5 through 8 discuss a novel formalism for describing the formation of quantum color centers in diamond and a numerical solver designed to simulate the implant and annealing processes. Finally, chapter 9 summarizes the work and discusses the outlook for the continued development of the technology described herein. vi TABLE OF CONTENTS LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Chapter 1. Background and Introduction . . . . . . . . . . . . . . . . . . . . . 1 1.1 Material Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Type Classification of Diamond . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Formation and Synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.1 Natural Growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.2 Synthetic Growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Diamond as a Semiconductor . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Chapter 2. Detector Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1 Nuclear Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2 Electronic Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.1 Diode Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.2 Measured Diode Performance . . . . . . . . . . . . . . . . . . . . . . 37 2.3 Preamplifier Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.3.1 POSSM Prototype Rev. 0 . . . . . . . . . . . . . . . . . . . . . . . . 40 2.3.2 POSSM Prototype Rev. 1 . . . . . . . . . . . . . . . . . . . . . . . . 42 Chapter 3. Detector Construction . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.1 Diamond Growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2 Diode Assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.3 Device Assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Chapter 4. Detector Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.1 LANL/LANSCE/WNR August 2022 . . . . . . . . . . . . . . . . . . . . . . 55 Chapter 5. Quantum Color Centers in Diamond . . . . . . . . . . . . . . . . . 60 5.1 Point Defects In Diamond . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.2 Properties and Applications of Diamond Color Centers . . . . . . . . . . . . 62 5.2.1 Optical and Spin properties . . . . . . . . . . . . . . . . . . . . . . . 63 5.2.2 Sensitivity to Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.3 Production of Color Centers . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4 Motivation for Color Center Formation Simulation . . . . . . . . . . . . . . . 68 Chapter 6. A Comparison of Monte Carlo Frameworks . . . . . . . . . . . . 69 6.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6.2.1 GEANT4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.2.2 GEANT4/GECO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.2.3 SRIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 vii Chapter 7. Modeling Color Center Formation in Diamond . . . . . . . . . . 87 7.1 Diffusion of Defects in Diamond . . . . . . . . . . . . . . . . . . . . . . . . . 87 7.2 Proposed Reaction Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 7.3 Scheme for Numerical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 97 7.4 Color Center Annealing and Reaction-Diffusion Solver . . . . . . . . . . . . . 100 7.4.1 Specifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 7.4.2 Flowchart and Theory of Operation . . . . . . . . . . . . . . . . . . . 101 7.4.3 Parallelization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7.4.4 Profiling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 7.4.4.1 CPU Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . 111 7.4.4.2 Performance of Jacobi Preconditioner . . . . . . . . . . . . . 114 7.4.4.3 GPU Acceleration . . . . . . . . . . . . . . . . . . . . . . . 115 Chapter 8. Experimental Benchmarking of CCANARD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 8.1 Comparison of CCANARD with Literature . . . . . . . . . . . . . . . . . . . 117 8.2 Benchmarking Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 Chapter 9. Conclusion and Continuing Work . . . . . . . . . . . . . . . . . . . 129 9.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 9.2 Continued Development of the POSSM Detector . . . . . . . . . . . . . . . . 131 9.3 Continued Development of the CCANARD program . . . . . . . . . . . . . . 132 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 APPENDIX A. Design Drawings Rev. 0 . . . . . . . . . . . . . . . . . . . . . . 147 APPENDIX B. Design Drawings Rev. 1 . . . . . . . . . . . . . . . . . . . . . . 148 APPENDIX C. Circuit Schematics . . . . . . . . . . . . . . . . . . . . . . . . . 149 APPENDIX D. Implant and Vacancy Profiles . . . . . . . . . . . . . . . . . . 150 APPENDIX E. CCANARD Quick Start Guide . . . . . . . . . . . . . . . . . 161 viii LIST OF FIGURES Figure 1.1: Two dimensional visualization of carbon atoms covalently bonded to one another resulting in the noble gas configuration with eight electrons in the outer orbital. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Figure 1.2: Two dimensional visualization of tetrahedral-oriented covalent bonds be- tween carbon atoms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Figure 1.3: Ball and stick model of conventional diamond unit cell. Image generated using CrystalMaker software [1] . . . . . . . . . . . . . . . . . . . . . . . 3 Figure 1.4: Common form planes of the diamond cubic. a) (100) b) (110) c) (111) Images generated using CrystalMaker software [1] . . . . . . . . . . . . . 4 Figure 1.5: Flowchart demonstrating the Type classification of diamond based on the presence or absence of nitrogen or boron impurities. A schematic, 2D representation of the diamond lattice shows the aggregation of defects for each type. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Figure 1.6: Schematic diagram of MACVD reactor. Microwaves resonate in a waveg- uide above the reaction chamber and are allowed to enter through a quartz window where they interact with process gases creating a plasma. Free ions in the plasma nucleate crystal growth onto a homoepitaxial substrate. Gas is constantly flowed in and pulled out with a vacuum pump keeping overall pressure around 100 Torr. . . . . . . . . . . . . . . . . . . . . . . 10 Figure 2.1: Diamond detector “sandwich” architecture proposed by Almaviva et al. 2010 [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Figure 2.2: POSSM ultra-fast neutron detector schematic design. . . . . . . . . . . . 16 Figure 2.3: Ball and stick diagram of diamond supercell; (a) (100) (b) (111) [1]. . . . 17 Figure 2.4: Fast (26.1 MeV) Neutron-Diamond interaction pathways detected by Ma- jerle, et al. 2020. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Figure 2.5: Total neutron cross section of gadolinium as a function of neutron energy [3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Figure 2.6: Results of POSSM detector efficiency as function of neutron energy from GEANT4 simulations with (blue) and without (orange) activation layer. 25 ix Figure 2.7: Energy levels in the metal/p-type junction with no applied external po- tential [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Figure 2.8: Energy levels in the metal/p-type junction in the reverse biased configu- ration [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Figure 2.9: Energy levels in the metal/p-type junction in the forward biased configu- ration [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Figure 2.10: Calculated I-V curve of diamond Schottky Diode obtained through finite element solution of the Poisson drift-diffusion equation. The negative current shown in the forward-biased condition is by sign convention choice. 34 Figure 2.11: Measured I-V curves of the POSSM diodes. Both diodes show the charac- teristic exponential increase in current with increasing forward bias. The dashed lines serve a guide for the eye. . . . . . . . . . . . . . . . . . . . . 38 Figure 2.12: Observed preamplifier response using the Rev. 0 circuitry. The input signal is a 2µs square pulse. The time scale is 5µs/Div. . . . . . . . . . 41 Figure 2.13: Calculated preamp circuitry frequency response. . . . . . . . . . . . . . . 43 Figure 2.14: Calculated preamp circuitry noise gain. . . . . . . . . . . . . . . . . . . . 43 Figure 3.1: HTHP Intrinsic diamond (001) face prepared for growth. Images taken at 25x using DICM. a) HMT 001a b) HMT 001b. Water spots are visible in the corners of HMT 001a, redisue from cleaning. . . . . . . . . . . . . . 47 Figure 3.2: HTHP diamond plate undergoing homoepitaxial CVD growth. Note the diamond plate (square object in center of the frame) lies countersunk in a molybdenum sample holder. The intensity of the plasma can be seen to saturate the sensor of the camera near the center top of the image. . . . 48 Figure 3.3: Single crystal growth. Images are taken via DICM at 25x. a) HMT 001a P+ b) HMT 001a P- c) HMT 001b P+ d) HMT 001b P- . . . . . . . . 49 Figure 3.4: The two halves of the diode bracket before final assembly. The Au contacts have been bonded to the plastic bracket pieces, and the diamonds have been fastened, P+ in contact with Au, with a small amount of RTV on the sides of the diamonds. Care has been taken to ensure RTV does not contaminate either the P+ or P- sides of the diamonds. . . . . . . . . . . 51 Figure 3.5: Completed preamp circuitry built on Bakelite perfboard. Through-hole components occupy the top of the board. All traces are connected on the reverse side. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 x Figure 3.6: Completed POSSM detector. Note the G10 panels and perfboard fastened by nylon machine screws. The Diode pack is blue, located near the top at the end of the G10 extension. . . . . . . . . . . . . . . . . . . . . . . . . 54 Figure 4.1: The author’s observation of the POSSM Detector’s energy response during operational testing. No data were recorded during this run. . . . . . . . 56 Figure 4.2: Uncalibrated q-spectrum of (a) POSSM Detector Ch. 2, and (b) 2nd Cividec B8 Diamond Detector . . . . . . . . . . . . . . . . . . . . . . . . 58 Figure 4.3: Uncorrected time of flight spectrum of (a) POSSM Detector Ch. 2, and (b) 2nd Cividec B8 Diamond Detector . . . . . . . . . . . . . . . . . . . 58 Figure 5.1: Point defects in diamond lattice. a) Vacancy. A carbon atom (shown in purple) has been removed from the lattice leaving a void. b) Substitution. A foreign ion (shown in yellow) has replaced a carbon atom at a lattice site. c) Interstitial. The interstitial atom resides in between lattice sites, loosely bound to the surrounding atoms. Images generated using CrystalMaker software [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Figure 5.2: Pseudo-covalent bonding. The dangling bonds left by the vacancy are indicated in yellow. The substitutional nitrogen has five valence electrons. It is covalently bonded with neighboring C atoms, leaving a lone pair, indicated in red, oriented towards the vacancy. This orientation acts like a covalent bond between the substitutional nitrogen and vacancy. . . . . 62 Figure 5.3: Color center structures. a) NV center b) Group III or Group IV color center where X∈{Al, Ga, In, Ti, Si, Ge, Sn, Pb} c) Single vacancy defect. Image taken from [5] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Figure 6.1: Implant depth profiles of 4 He+ into Diamond. . . . . . . . . . . . . . . . 75 Figure 6.2: Implant depth profiles of H+ into Diamond. . . . . . . . . . . . . . . . . 75 Figure 6.3: Implant depth profiles of 14 N+ ions into Diamond. . . . . . . . . . . . . 75 Figure 6.4: Implant depth profiles of 82 Pb2+ ions into Diamond. . . . . . . . . . . . 75 Figure 6.5: Percent difference in average implant depth compared to SRIM. . . . . . 76 Figure 6.6: Vacancy depth profiles created by 4 He+ in Diamond. . . . . . . . . . . . 78 Figure 6.7: Vacancy depth profiles created by H+ in Diamond. . . . . . . . . . . . . 78 Figure 6.8: Vacancy depth profiles created by 14 N+ in Diamond. . . . . . . . . . . . 78 xi Figure 6.9: Vacancy depth profiles created by 82 Pb2+ in Diamond. . . . . . . . . . . 78 Figure 6.10: Percent Difference in average vacancy depth compared to SRIM. . . . . . 79 Figure 6.11: Comparison of SRIM, and GEANT4(/GECO) energy deposition profiles for 1000 keV implants of (a) 4 He+ , (b) H+ , (c) 14 N+ , and (d) 82 Pb2+ . . 80 Figure 7.1: Interstitial and Vacancy mediated diffusion. a) Interstitial (red) mediated diffusion always exchanges with a nearest neighbor (black); b) Vacancy (white) mediated diffusion exchanging with nearest neighbor (black); c) Vacancy mediated diffusion exchanging with second-nearest neighbor . . 92 Figure 7.2: CCANARD program flow chart . . . . . . . . . . . . . . . . . . . . . . . 103 Figure 7.3: The CCANARD “About” informational display reporting that it was com- piled with nvfortran version 22. Note the ASCII ducks, the unofficial logo of CCANARD. (Canard is the French word for duck.) . . . . . . . . . . . 106 Figure 7.4: Multiprocessing performance of CCANARD simulating the annealing of 250 keV nitrogen implants in diamond. The domain is discretized into 51 × 51 × 51 finite difference nodes. Lower is better. CCANARD was compiled with: gfortran -fopenmp -O3 -cpp -s -march=native . . . 112 Figure 7.5: Multiprocessing performance of CCANARD simulating the annealing of 500 keV Proton implants in diamond. The domain is discretized into 51 × 51 × 51 finite difference nodes. Lower is better. CCANARD was compiled with: gfortran -fopenmp -O3 -cpp -s -march=native . . . 113 Figure 7.6: Comparison of computational speed of CG and Jacobi Preconditioned CG with increasing system size. Lower is better. CCANARD was compiled with: gfortran -fopenmp -O3 -cpp -s -march=native . . . . . . . . 115 Figure 7.7: GPU acceleration performance in setting initial conditions. Lower is better. CCANARD was compiled with: nvfortran -mp -Minfo=accel -stdpar=gpu -fast -O3 -cpp . . . . . . . . . . . . . . . . . . . . . . . 116 Figure 8.1: This image, taken from the work of Pezzagna et al. (2010) [6] shows clearly the dependence of color center yield on ion implant dependence in (a) and shows the same relationship in the number of vacancies created as a function of implant ion energy in (b). . . . . . . . . . . . . . . . . . 119 Figure 8.2: Comparison of the energy dependence of color center conversion efficiency modeled in CCANARD with experimental results reported in Lühmann et al. (2018) [7]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 xii Figure 8.3: Color center concentration as a function of annealing temperature at given fluences as simulated by CCANARD. Concentration levels off after 875◦ C, and in some cases even decreases, matching the experimental work of Acosta et al. [8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Figure 8.4: CCANARD simulation of the annealing of 5 keV nitrogen implants for 2 hours at 800◦ C. Conversion efficiency increases with fluence until peaking and then rapidly dropping off. . . . . . . . . . . . . . . . . . . . . . . . . 122 Figure 8.5: Experimentally determined relationship between color center conversion efficiency and fluence using 5 keV nitrogen implants annealed at 800◦ C for 2 hours. The dotted curve is a guide for the eye. Image taken from Pezzagna et al. (2010) [6]. . . . . . . . . . . . . . . . . . . . . . . . . . . 123 Figure 8.6: Final (10 hr) conditions along the sample center line after the 800◦ C (1073 K) anneal. The beam propagates in the x-direction and impinges at the origin. The final color center distribution spans approximately 0.15µm, and is not Gaussian, likely due to a high concentration of vacancies im- peding color center formation in this location [9]. . . . . . . . . . . . . . 126 Figure 8.7: Final (10 hr) conditions along the sample center line after the 1200◦ C (1473 K) anneal. The beam propagates in the x-direction and impinges at the origin. The final color center distribution spans approximately 0.12µm with a tail extending approximately an additional 0.5µm. The distribution is much more Gaussian than in Figure 8.6, likely due to increased mobility in the vacancies. Note that the distribution of interstitials is nearly uniform.127 Figure 8.8: Final (10 hr) conditions along the sample center line after the 1550◦ C (1823 K) anneal. The beam propagates in the x-direction and impinges at the origin. The final color center distribution spans approximately 0.12µm with a tail extending approximately an additional 0.5µm. The distribution is much more Gaussian than in Figure 8.6, likely due to increased mobility in the vacancies. Note that the distribution of vacancies and interstitials is nearly uniform. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 Figure E.1: CCANARD main menu, visible when first entering the program. . . . . 163 Figure E.2: Menu Option 1) Simulation Control. This menu allows the modification of administrative properties of a CCANARD simulation. . . . . . . . . . 164 Figure E.3: Menu Option 2) Physical Domain. This menu allows the modification of the physical and material properties modeled in a CCANARD simulation. 166 Figure E.4: Additional physics options accessed through the More Options... selec- tion of the Physical Domain menu. . . . . . . . . . . . . . . . . . . . . . 168 xiii Figure E.5: Menu Option 3) Discretization and Solver. This menu allows the operator to specify how CCANARD will discretize the domain and solve systems of coupled PDEs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 Figure E.6: Menu Option 4) Annealing options. This menu allows the operator to enter the annealing recipe which CCANARD will simulate. . . . . . . . 172 Figure E.7: Menu Option 5) Output Options. This menu allows the operator to specify and modify CCANARD’s visual output in the form of plots and anima- tions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 Figure E.8: These eight color maps are available in CCANARD. . . . . . . . . . . . 175 Figure E.9: All settings have been verified, indicated by check marks. The operator is able to begin the simulation. . . . . . . . . . . . . . . . . . . . . . . . . 176 Figure E.10: The CCANARD “About” informational display reporting that it was com- piled with nvfortran version 22. Note the ASCII ducks, the unofficial logo of CCANARD. (Canard is the French word for duck.) . . . . . . . . . . 177 Figure E.11: The directory structure of a CCANARD simulation. . . . . . . . . . . . 178 xiv Chapter 1. Background and Introduction Diamond is a solid, structured allotrope of carbon. The earliest written records of diamond date to the third and fourth centuries BC. These historical records describe the magical and robust nature of the material, the first system of assessing the quality of diamond, and the use of diamonds as a method of transfer of wealth [10]. Today, just as 2500 years ago, diamond is regarded as the most precious of gemstones. It is also well known for its hardness and durability, frequently finding a home on the surfaces of polishing and cutting instruments, particularly saw blades, core barrels, and drill bits designed to cut rock or concrete. 1.1 Material Properties Diamond owes these marvelous properties, along with many others, to both its electronic configuration and crystal structure. Carbon has an electron configuration of 1s2 2s2 2p2 , with the four valence electrons spread across the s and p orbitals. In diamond the s orbital mixes with the three p orbitals (px, py, pz) described by the wave function Ψ = 21 (Ψ2s + Ψ2px + Ψ2py + Ψ2pz ). Each carbon atom in the diamond lattice contributes its four valence electrons and receives one valence electron from each of its four neighbors resulting in the so-called noble gas configuration, with eight electrons in the outer shell (Figure 1.1). Each carbon atom is covalently bonded to its four neighbors in what is known as the sp3 bonding which results in a highly directed charge density, tetrahedral in shape. The bond length is 1.545Å[11][12]. This bonding is most easily visualized in two dimensions as shown in Figure 1.2, though the idea generalized to three dimensions gives the full description of bonding in the diamond lattice. 1 Figure 1.1: Two dimensional visualization of carbon atoms covalently bonded to one another resulting in the noble gas configuration with eight electrons in the outer orbital. Figure 1.2: Two dimensional visualization of tetrahedral-oriented covalent bonds between carbon atoms. Diamond crystal belongs to the Fd3̄m(227) space group. The diamond conventional unit cell is the face-centered cubic Bravais lattice with a lattice constant (edge lengths of the unit cell) of a = b = c = 3.567Å. The diamond unit cell is the definitive form of its kind, and the face-centered cubic Bravais lattice is often referred to as the “diamond cubic” even when referring to other elements or minerals. A ball-and-stick model of the conventional diamond unit cell is shown in Figure 1.3. In this figure the full three-dimensional nature of the tetrahedral bonds of each atom to its four nearest neighbors is apparent. 2 Figure 1.3: Ball and stick model of conventional diamond unit cell. Image generated using CrystalMaker software [1] As a naturally occurring mineral, diamond shows adamantine luster and is frequently yellow, brown or colorless, and may, less frequently, contain hints of other hues. It has (111) cleavage perfect in four directions and its common forms include (111), (110), (100), depicted in figure 1.4, frequently with spinel law common (111) twinning [13]. It is the hardest known natural material defining the upper end of both the Vickers and Mohs scales. Diamond has the highest atomic density of any known material, contributing to its hardness and incompressibility. Diamond has also the highest thermal conductivity of any known material. High grade diamond, natural or synthetic, demonstrates wide band optical transparency in the upper x-ray, UV, visible, infrared, far infrared, and into the microwave band. As an electronic material, diamond is a wide band gap semiconductor, with a band gap nearly five times wider than that of silicon. Intrinsic diamond is an insulator. The bulk properties of diamond are presented in Table 1.1. Further discussion of the electronic properties of diamond can be found in Chapter 2. 3 Figure 1.4: Common form planes of the diamond cubic. a) (100) b) (110) c) (111) Images generated using CrystalMaker software [1] 4 Table 1.1: Material properties of diamond Property Value Ref Mohs Hardness 10 Young’s Modulus 1050 GPa [14] Density 3.51 g/cc [13] Atomic Density 1.77x1023 atoms/cc [14][15] Thermal conductivity 1800-2200 W/mK [14] Thermal Expansion Coefficient 10−6 /K [14] Debye Temperature 1860K [14] Band Gap 5.45-5.47 eV [15][14] Dielectric Constant 5.26-5.7 [14][16] Dielectric Breakdown Field 10 MV/cm [15] 1.2 Type Classification of Diamond Pure diamond is made entirely of one element – carbon. However, rarely is diamond, natural or synthetic, completely pure. Often atoms of elements such as nitrogen or boron can replace a small number of carbon atoms in the lattice. While it is possible for other impurities to be found in diamond, classification of diamond is done by the presence of nitrogen or boron. Type I diamonds contain nitrogen impurities and Type II diamonds do not contain nitrogen impurities. Each type is broken into sub-types. The “Type” classification system of diamond is illustrated in Figure 1.5, as originally published in Gems & Gemology [17]. 5 Figure 1.5: Flowchart demonstrating the Type classification of diamond based on the pres- ence or absence of nitrogen or boron impurities. A schematic, 2D representation of the diamond lattice shows the aggregation of defects for each type. Approximately 95% of natural diamonds are Type Ia. These diamonds contain clustered nitrogen impurities and vary in color from near colorless to light yellow. Type Ib diamonds, in contrast, are the rarest diamonds. They also contain nitrogen defects, but as isolated atoms rather than clusters. These diamonds are often bright yellow. Type IIa diamonds contain no significant nitrogen or boron imputirites. They are usually clear, but can also take on a grey-ish hue. Type IIb diamonds contain boron impurities and can be naturally electrically conductive depending on the boron concentration [17]. For electronic and optics applications, generally Type IIa diamond is desired, as this represents the purest diamond available, and the synthesis of Type IIa diamond has been 6 a burgeoning field of research since the mid 1990’s. Prior to the widespread availability of synthetic Type IIa diamonds, diamond application research relied on naturally sourced geologic diamonds. 1.3 Formation and Synthesis 1.3.1 Natural Growth Natural diamonds are found within kimberlite, an igneous rock formed within the Earth’s mantle at depths of 140 to 240km below the surface at pressures in excess of 5 GPa [18]. The environment for kimberlite (and thus diamond) formation is potassic, heavily basic magma, in the presence of anomalous concentrations of exotic minerals. The majority of mined diamonds come from kimberlite pipes, pushed upward into the crust since protozoic times by gas-charged CO2 and H2 0 fluidic systems. Contamination with crustal minerals accounts for the complex chemistry of kimberlite and impurities in diamond. The yield of diamonds from kimberlite ore is approximately 1 part in 8 million [19][20]. For reasons of its immense demand and natural scarcity, synthetic fabrication of diamond was a long-sought process. 1.3.2 Synthetic Growth The first laboratory grown diamonds were manufactured in 1953 and have, since then, been manufactured using the High Temperature/High Pressure (HPHT) process, somewhat anal- ogous to the geologic formation of natural diamond [21]. Surprisingly to researchers of the time, simply applying the equivalent pressure and temperature to graphite did not success- fully produce diamond. Instead, the temperature gradient method was developed, whereby 7 a carbon allotrope (usually graphite) is combined with a small diamond seed and a Group VIII element such as iron or a ferrous alloy is added as a solvent/catalyst. Pressure and temperatures for this process vary by specific method but are generally between 5-6 GPa and 1300-1600◦ C, respectively. In this environment the carbon allotrope is dissolved in the ferrous alloy and free carbon atoms flow against the temperature gradient and nucleate on the diamond seed. Individual methods vary by the geometry of the pressure vessel and manner in which pressure is applied, [22]. The HTHP process is prone to producing diamonds with impurities. Impurities can be easily incorporated through the ferrous catalyst, leech from the HPHT capsule, or be present in the carbon source material. For this reason, the majority of HPHT produced diamonds are of Type Ib with nitrogen concentrations on the order of 100-300 ppm [18]. At this concentration, the diamonds take on a yellow hue. Inconsistency may also be present in HTHP diamonds. A survey of 52 colorless or near colorless (colorless being suggestive of lower nitrogen impurities[18]) HTHP diamonds produced by a major commercial manufacturer showed considerable variance sample-to-sample. 50% of the samples were magnetic, 19.2% were electrically conductive, and the samples showed boron impurity concentrations between 0.0 ppb and 216±30 ppb. Additionally, four samples showed the presence of silicon impurities [22]. This level of inter-sample variability is a hurdle to the design and construction of any device or process which relies on consistent material properties in its engineering. Researchers searched for a method to produce more consistent higher quality diamonds, and in the mid 2000’s Chemical Vapor Deposition (CVD) had come into its own[21]. The CVD process requires the reaction of volatile process gases with a surface to nucleate crystal growth. There are several types of CVD reactors, differing in the method by which the plasma is generated. 8 Hot filament reactors rely on the flow of process gases over a hot filament at less than atmospheric pressure to generate plasma for crystal growth. While this design is simple and inexpensive, it is prone to contaminating the diamond growth with material from the filament, usually tungsten or tantalum. This is undesirable in optical or electronic diamond [23]. Microwave Assisted CVD (MACVD), on the other hand, is a popular method for high purity single crystal diamond growth. In this process microwave power sources are used to generate plasma by coupling microwaves in a resonant chamber through a quartz dielectric window. Process gases are cycled through the reaction chamber which is maintained at less than atmospheric pressure [23][24]. Energy transfer from the microwaves to the process gas molecules through coupling with the gas phase electrons creates a plasma [23]. Free ions and neutral radicals in the plasma are available to nucleate crystal growth on a diamond seed. A schematic of such a system is shown in Figure 1.6. The process gases used in MACVD diamond growth include a carbon source gas, a carrier gas, and a dopant gas. Common carbon source gases include propane (C3 H8 ) or methane (CH4 ). Methane is popular for its availability and purity, though has lower carbon cracking efficiency than propane [25]. Hydrogen acts as a carrier gas and co-reactant, helping high quality diamond rather than graphite nucleate [23]. Additional process gases can be added, such as diborane (B2 H6 ), nitrogen (N2 ), lithium-t-butoxide (LiOC4 H9 ) [26], or triphenylphosphane (P(C6 H5 )3 ) [27], among others [28], to dope the diamond growth with P or N type carriers to produce semiconducting diamond films. 9 Figure 1.6: Schematic diagram of MACVD reactor. Microwaves resonate in a waveguide above the reaction chamber and are allowed to enter through a quartz window where they interact with process gases creating a plasma. Free ions in the plasma nucleate crystal growth onto a homoepitaxial substrate. Gas is constantly flowed in and pulled out with a vacuum pump keeping overall pressure around 100 Torr. 1.4 Diamond as a Semiconductor Diamond is a Group IV semiconductor, joining other venerable semiconductors such as silicon (Si) and germanium (Ge). However, diamond demonstrates properties that make it more desirable for use in highly demanding applications and austere environments. Table 1.2 compares the semiconducting properties of diamond with other common semiconductors and lists the application advantages of diamond. This work focuses on applying semiconducting diamond to the task of neutron detection. As such, application specific detail of semiconducting diamond will be provided throughout, primarily in Chapters 2 and 3, covering the design and construction of the device. 10 Table 1.2: Comparison of semiconductor properties of common semiconductors Diamond [29] Si [29] Ge GaN [29] 4H-SiC [29] Application benefit Band Gap (eV) 5.47 1.1 0.67 [30] 3.44 3.2 High temperature Breakdown Field (MV/cm) 10 0.3 0.1 [31] 5 3 High voltage e- saturation velocity (x107 cm/s) 2 0.86 1.3 [31] 2.5 3 High frequency Hole saturation velocity (x107 cm/s) 0.8 n/a n/a n/a n/a e- mobility (cm2 /Vs) 2000-4000 1450 3800 [30] 440 900 Hole mobility (cm2 /Vs) 1800-3800 480 1800 [30] 200 120 Thermal Conductivity (W/mK) 2200 150 55 [30] 130 500 High power 11 Chapter 2. Detector Design The Positionally Opposed Schottky Semi-Metal (POSSM) Ultra-fast Neutron Detector is a solid state particle detector built from semiconducting diamond optimized for neutron detection with non-negligible efficiency in the ultra-fast regime above 20MeV. The architecture of the POSSM detector is inspired by and borrows heavily from extant semiconducting diamond thermal neutron detectors. The POSSM detector pushes detection sensitivity beyond the current body of work and into higher energy regimes through the use of Semi- and Inner Transition Metals to form a rectifying metal-semiconductor interface. Using this interface, two Schottky diodes are built around a shared cathode. Electronically speaking, the detector resembles two mutually opposed Schottky diodes. The advantages of the metal-semiconductor junction over the more common P-N junc- tion include faster response and lower bias voltages [32][33]. Schottky barrier devices have been used with success for particle detection even with zero bias voltage [33]. The choice of diamond as a semiconductor material was driven by electrical, atomic, and nuclear consid- erations. The nuclear and atomistic/electronic properties of diamond are discussed in detail in Sections 2.1 and 2.2, respectively, but the wide band gap and large carbon cross section make diamond a desirable material for neutron detection, the former reducing signal noise and the latter enhancing detector efficiency. Other material properties of note in diamond include the intrinsic breakdown potential, of which diamond has the highest of any semiconductor at 10MV/cm. Similarly, diamond has excellent thermal properties with room temperature conductivity of 2200 W/mK [29]. Diamond truly is one of the most remarkable naturally occurring materials and perfectly suited for high power applications. However, bias voltages in diamond detectors are typically 12 an order of magnitude lower those than usually required in a scintillator and photo multiplier tube detection chain. Research cites bias voltages between 2 and 5 V/µm for a total bias on the order of hundreds of volts, depending on detector thickness [15][34][35]. The design of the POSSM detector was governed by two design philosophies: 1. Develop a device that achieves appreciable neutron detection efficiency at and above 50MeV. 2. Keep the device as small as possible to allow easy incorporation into existing nuclear physics experiments. The POSSM detector should provide neutron detection capabil- ities to experiments that otherwise might not have such capability. These design philosophies are somewhat contradictory, as efficient neutron detection, espe- cially in the ultra-fast regime, requires significant detector volume. The leading ultra-fast neutron scintillator arrays have volumes on the order of cubic meters. After a rigorous design process, the architecture of this detector reduces the form factor from contemporary high energy neutron detectors by nine orders of magnitude (from cubic meters to cubic millimeters) with only a two order of magnitude reduction in gross efficiency. The small form factor of the detector will provide desirable capability to nuclear scientists, including, but not limited to, unparalleled spatial resolution and self-vetoing, allowing neu- tron and charged particle discrimination without the need of an auxiliary detector. The device will provide neutron detection capability to nuclear experiments which may otherwise be unable to take advantage of larger scintillator detectors, radiation-hard beamline mon- itoring capabilities across a wider range of neutron (or charged particle) energies, among other applications such as nuclear medicine and port screening. This design builds on the work of Almaviva et al. (2010) [2], who improved existing diamond neutron detector designs by incorporating a “sandwich” design. The design, shown 13 in Figure 2.1, consists of CVD diamond grown homoepitaxially on an HTHP seed. The first growth grew 25µm of heavily boron doped P-type diamond, on top of which was grown 55µm of intrinsic diamond. This formed the active region of the detector. Almaviva et al. used a 1µm thick 6 Li activation layer. The authors’ reason for developing the “sandwich” architec- ture was to electrically isolate the HTHP seeds from the inner workings of the detector. The HTHP seeds were assessed to have a higher defect rate and, the authors wanted to prevent any defects in the HTHP seed from any deleterious effects in the workings of the detector [2]. Figure 2.1: Diamond detector “sandwich” architecture proposed by Almaviva et al. 2010 [2]. The current work began with the sandwich design. The “sandwich” was retained through- out the design process and into production, but for different reasons that those cited by Al- maviva et al. The POSSM detector would also use homoepitaxial growth on an HTHP seed, but the HTHP seeds acquired through New Diamond Technology (Saint Petersburg, Russia) were of electronic grade with low ppb nitrogen impurities. This level of quality is assessed to not impart any deleterious effects in the performance of the detector. In fact, the quality of 14 this material allows it to be used for the active region of the POSSM device. The sandwich design means the POSSM detector is basically two detectors in one, as the intrinsic diamond in each half of the “sandwich” is its own active detector region. This means the detector has two channels, both of which can independently detect radiation events. By convention, the channel which faces the incident beam is channel 1 and the back layer forms channel 2. The active area of the detector, the intrinsic diamond, is oriented such that the preferred direction of incident beam (See Figure 2.2) will be nominally oriented along the ⟨1̄00⟩ axis, normal to the (100) plane. There are three reasons for this: 1. HTHP seeds with this orientation are the easiest to acquire. 2. ⟨001⟩ (by symmetry, analogous to ⟨100⟩) oriented growth results in shorter growth times and fewer defects such as crystal twinning [25]. 3. The (100) face provides the greatest “surface-projected atomic density” (SPAD), the two-dimensional atomic density if all the atoms in the bulk material were projected onto a single plane. This is in contrast to surface density, the two-dimensional atomic density of only atoms in a physical plane. A greater surface-projected atomic density should result in better neutron detection efficiency, as it simply gives an incident neutron more possible carbon atoms from which to scatter. SPAD was calculated for the (111) face (the face with the highest surface density) and for the (100) face via conventional technical drawing take-off methods using scale models and a Keuffel & Esser model 4296 compensating polar planimeter. The (100) face was determined to have an SPAD of 0.7228 ± 0.0006 atoms/Å2 and the (111) face has an SPAD of 0.6359 ± 0.0018 atoms/Å2 . Figure 2.3 shows scale diagrams of each crystal face. The prototype POSSM design calls for 1µm of heavily boron doped P-Type diamond. This material is to be so heavily doped as to have near metallic conductivity for use as a 15 readout connector. This requires a concentration of ∼ 1020 atoms/cc. On top of this layer would be 500µm of intrinsic diamond. The intrinsic diamond will be the electronic grade HTHP seed. On the other side of the intrinsic diamond, 3µm of lightly boron doped P-type diamond will be grown via CVD. Two such diamonds are grown, and a 10µm gadolinium foil is sandwiched between the lightly doped P- layers. Readout contacts are gold. This is shown schematically in Figure 2.2. Please see Appendix A for detailed design documents. Figure 2.2: POSSM ultra-fast neutron detector schematic design. 16 Figure 2.3: Ball and stick diagram of diamond supercell; (a) (100) (b) (111) [1]. The diamond “sandwich” is mechanically fastened via a 3D printed bracket. The bracket is a two-piece clam shell design essentially squeezing the diode assembly together and en- suring physical and electrical contact of components. The two-piece bracket is held together by 2x #2-56 nylon machine screws torqued to a stress of 16 psi. The bracket is printed of polyetheretherketone (PEEK) plastic, an industrial plastic known for its strength, excellent electrical properties, and most importantly in this application, radiation hardness. 1/16” G10 fiberglass board forms the body of the detector, holding the diode bracket assembly, onboard electronics, and input/output connectors while also offering protection to the sensitive components of the detector. In the prototype design, all onboard electronics are built on standard single-sided perf- board with conventional 0.1” pad spacing. High voltage is connected via an MHV terminal mounted to the G10 body. Common ground, ±12V, and detector signal are interfaced through female LEMO 00 series 50Ω coaxial connectors. Raw signal taken directly from the diodes without preamplification and the preamplified signal for both channels are output 17 using the same LEMO 00 series connectors. 2.1 Nuclear Design The potential of diamond as a material for particle detection, specifically for the detection of alpha particles was first noted in the late 1940’s [36]. In the 1970’s, work performed largely by notable Soviet scientists demonstrated the capability of type-II natural diamonds in alpha radiation detection in austere environments [37]. Synthetic diamond fabrication techniques, both CVD and HTHP, were improved sufficiently such that by the 1990’s electronic-grade artificial diamond could be made with reliability sufficient to fuel further research into dia- mond based radiation detection [38]. Modern diamond detectors have demonstrated 100% charge-collection efficiency in the detection of charged particles [39]. In recent years, interest has grown in using diamond-based detectors for the detection of neutrons. Most of this work has focused on the detection of thermal neutrons [15][34][2]. Diamond is particularly well suited for neutron detection. It is among the most atomically dense of naturally occurring or man-made materials, and has a fast neutron cross section of 2.31b (approximately 77% that of leading engineered scintillators) resulting in excellent neutron attenuation [15]. Because of the electronic properties of diamond (see Section 2.2), diamond detectors have shown remarkable energy resolution in the detection of thermal neutrons, as low as 1.5% FWHM [40]. This value represents the ratio of the FWHM meaure of the detected signal to the energy of the incident particle. The lower this value, the better the detector is able to distinguish between particles of simlar energies. Further, because diamond demonstrates excellent radiation hardness [39], diamond detectors are frequently employed in low-energy/high-flux environments, such as reactor radiation monitoring, spallation source 18 beam monitoring, and radiation-surveillance [2][38][40][41]. Whereas the neutron, being a neutral particle, cannot interact coulombically, and thus cannot create electron-hole pairs in a semiconducting medium, any solid state neutron detec- tor, including those of diamond construction, must rely on nuclear interaction to transitively generate electron-hole pairs which can then be read out as a signal. (Please see Section 2.2 for a detailed treatment of the solid state physics of particle detection.) Three events can occur when a neutron impinges on a diamond detector [38]: 1. The neutron passes through the diamond crystal without interaction. The neutron is not detected. 2. The neutron undergoes nuclear elastic scattering, displacing a carbon atom from the lattice. If the displaced carbon atom ionizes or induces additional ionization in the lattice, the neutron is transitively detected. 3. The neutron is absorbed by a carbon atom through nuclear interaction. Daughter particles are subsequently ejected and create electron-hole pairs through ionization of lattice atoms. The neutron is transitively detected. Existing literature cites option 3 as the predominant vector of transitive neutron detection [2][15][34][40]. Through an examination of the response 17-34MeV neutrons in a diamond detector, Majerle, et al. (2020) [41] identified four neutron-carbon reaction pathways. The data they present suggest the presence of two additional pathways, which the authors do not identify. Their results are shown below in figure 2.4. The (n, d) pathway dominates the other identified pathways in probability of occurrence at an incident neutron energy of 26.1 MeV [41]. 19 Figure 2.4: Fast (26.1 MeV) Neutron-Diamond interaction pathways detected by Majerle, et al. 2020. Literature is sparse on the subject of ultra-fast neutron-carbon reactions. A search of the EXFOR database for neutrons against a diamond target at energies above 50 MeV yields no results [42]. It is well known that neutron cross sections can vary drastically with neutron energy so additional efforts were made to understand the expected performance of diamond subjected to neutrons in excess of 50 MeV. Geometry and Tracking 4 (GEANT4) simulations[43] were used to analyze neutron- diamond reactions between 50 MeV and 150 MeV. The simulations consisted of a custom diamond target material invoked through an instance of the G4Material class, using G4 C from the NIST Materials Manager set to a density of 3.51g/cc[13]. Crystalline structure was incorporated through the GEANT4Crystal Object (GECO)[44] code using an extension of the G4Material class. Results are tabulated in Table 2.1. The (n, α) and (n, p) pathways are dominant throughout this energy regime, with the (n, α + p) and (n, α + d + p) pathways becoming significant as energy increases. The presences of one or more of these peaks in 20 detector output will indicate the successful detection of neutrons. Table 2.1: Expected inelastic reaction pathways for ultra-fast Neutron-Diamond interactions Reaction Products Q-Value (MeV) [45] 50 MeV 100 MeV 150 MeV 9 Be + α -5.7 56.92% 36.40% 28.06% 12 B +p -12.59 32.38% 36.76% 38.74% 11 B +d -13.73 3.13% 0.00% 0.40% 10 B +t -18.93 0.00% 0.00% 0.00% 7 Li +α+d -22.4 0.26% 1.10% 0.79% 8 Li +α+p -22.59 2.35% 11.76% 17.00% 8 Li +α+p -22.59 2.35% 11.76% 17.00% t + d + 2α -24.86 1.04% 0.37% 0.00% 8 Be +d+t -24.96 0.00% 0.00% 0.40% 10 Be + p + d -24.96 1.57% 0.37% 1.18% 5 He +α+t+p -27.82 1.83% 4.78% 2.37% 6 He +α+d+p -32.37 0.52% 8.46% 11.06% Extant diamond-based fast neutron detectors frequently employ an “activation layer,” a thin layer of, usually lithium, to enhance the neutron cross section of the device as a whole and help tune the device for improved efficiency in a particular energy regime. The idea behind the activation layer is that a thin layer of a material with a large neutron cross section will result in more neutron capture events subsequently producing more ionizing radiation which can be detected through electron-hole pair generation in the semiconducting 21 diamond [2][34]. The POSSM Ultra-fast Neutron Detector also employs an activation layer. The activation layer material chosen for this detector is gadolinium (Gd). In the early design process, Cadmium (Cd) was considered. Cadmium is frequently used in neutron detectors as an activation layer or moderator. This is most often done in conjunction with a scintillator detector [46][47]. Cd would serve well as an activation layer, with a large neutron cross section and excellent neutron attenuation properties. It also has the electronic properties required to serve as the shared cathode in the Positionally-Opposed Schottky barrier architecture. However, concerns with the toxicity and possible persistent activation after neutron exposure led to Gd being chosen for the construction of the POSSM detector instead. gadolinium has the largest thermal neutron cross section of any element, of 49,000b.[48] Gd’s neutron cross section drops off quickly with increasing neutron energy[49], but retains a nearly constant value of approximately 5.5b through 20MeV [50] [51]. The neutron cross section of gadolinium up to 20MeV is shown in Figure 2.5. Reliable data for higher energies was not found, though the constant value through 20MeV suggests Gd may perform well even into the ultra-fast regime. Gadolinium is found in seven naturally occurring stable isotopes in abundances as presented in Table 2.2. 22 Table 2.2: Natural abundances of stable gadolinium isotopes Isotope Abundance 152 Gd* 0.20% 154 Gd 2.15% 155 Gd 14.73% 156 Gd 20.47% 157 Gd 15.68% 158 Gd 24.87% 160 Gd 21.90% *semi stable. t1/2 = 1.08 × 1014 yr Figure 2.5: Total neutron cross section of gadolinium as a function of neutron energy [3]. The thickness of the activation layer must be optimized. In general, a thicker activation layer will result in more neutron interactions, but an activation layer which is too thick 23 may stop charged daughter particles produced from n-Gd reactions from entering the semi- conducting diamond. All the predominant n-Gd reactions produce α-particles, and many produce only α-particles plus a heavy ion [45]. In all of these cases, the α-particles will be the shortest traveled (ignoring heavy ions which will have only negligible energy), and will therefore govern the activation layer thickness. A 1.85 MeV α-particle (lower end of typical n-Gd reactions)[45] has a projected range of 4.758µm in Gd [52]. Again, GEANT4 simulations were used to aid in design. Simulating 50MeV neutrons incident on thin layer of Gd showed an optimized thickness between 5µm and 10µm. Origi- nally the design intent had been to sputter or otherwise electronically deposit the Gd layer onto the semiconducting diamond. Constructability concerns about the adhesion and fea- sibility of sputtering a layer of the specified thickness led to a modification of the design to use Gd foil. 10µm thick foil was chosen based on availability. The foil was purchased from Goodfellow, PN: GD00-FL-000110. Design drawing are included in their entirety in Appendix A. The entirety of the detector design was simulated in GEANT4. The detector was mod- eled as a stack of absorbers consisting of amorphous G4 C or G4 Gd. Note that in these simulations, amorphous carbon was used for the diamond absorbers rather than the custom diamond material described above. This will cause a systematic underestimate of detector efficiency. Incident neutron energy was varied from 10 MeV to 150 MeV in increments of 10 MeV. After each simulation, data was analyzed to ascertain detection efficiency. Because of the excellent charge collection efficiency of diamond diode devices, any ionizing radiation produced by a neutron interaction in or otherwise entering either of the active regions of the detector was scored as a successful detection. Total detector efficiency was then calculated as the ratio of successful detections to the number of incident neutrons. This was done 24 both with and without the activation layer. The results of these simulations are presented in Figure 2.6. The detector with the activation layer consistently performs by a margin of between 2 - 5% beyond what the semiconducting diamond can do alone. A notable boost in performance occurs at 35MeV. The reason for this is unknown. The 35 MeV instance was simulated multiple times and the significant Gd enhanced performance is consistent – it appears not to be an artifact of psuedo-random number generation in the Monte Carlo simulation. Figure 2.6: Results of POSSM detector efficiency as function of neutron energy from GEANT4 simulations with (blue) and without (orange) activation layer. The POSSM Ultra-fast Neutron Detector can be configured in self-vetoing mode to aid the nuclear experimentalist in particle discrimination. The “sandwich” design allows the signal from the first active region of the detector to serve as a veto. Below a certain energy level, the energy level at which energetic protons would be expected to pass through the first active region, charged particles will be detected and stopped only in the first detector region and only neutrons will be detected in the activation layer or second region of the detector. Channel 2 can be vetoed or triggered by channel 1. This removes the need for a dedicated 25 veto scintillator, simplifying experimental setup. Stopping and Range of Ions in Matter (SRIM) simulations conducted by Ryc, et al. (2016) indicate the limiting particle energy for the self-vetoing configuration is the 10 MeV proton [13]. Reactions producing protons of a higher energy will not allow for reliable use of self-vetoing configuration because charged particles are expected to penetrate into the activation layer or second region of the detector. Future design revisions will include thicker active regions (thicker diamonds), limited only by the semi-conducting physics of the diamond material. This will improve overall detector efficiency and allow the self-veto configuration to operate at higher particle energies. Future designs could incorporate asymmetry in the thickness of the first and second regions of the detector, permitting tuned efficiency in each channel, independently. 2.2 Electronic Design 2.2.1 Diode Design The fundamental function of this device is that of the Schottky diode, the ubiquitous Metal- Semiconductor diode named for the work of the German physicist Walter Schottky. Schottky barrier particle detectors represent an extant group of detector archetype[53]. While not as widespread as their p-n and p-i-n counterparts, Schottky detectors are employed in the detection and monitoring of charged particles, and have been adapted to the detection of fast and thermal neutrons [54]. Schottky barrier detectors are particularly well suited for fast response time and low signal applications, given a Schottky barrier height in excess of 0.5 eV [32][33]. The operating principle of the Schottky barrier detector is like that of other semicon- 26 ductor particle detectors, relying on incident radiation to create electron-hole pairs within the depleted region of the biased device. The electrons and holes will drift against and with the electric field, respectively, creating a signal which can be read out [12]. By knowing the electron-hole formation energy of the semiconductor and the magnitude of the measured signal, one can calculate the energy of the incident particle. In the case of neutron detection, this is two-step process, first relying on a neutron-nuclear interaction to produce ionizing radiation which can then be detected through the afore-described electron-hole mechanism. The subject device consists of boron-doped diamond, a p-type semiconductor. The ac- tivation layer, a 10µm layer of gadolinium, forms the rectifying contact (Figure 2.7). The device is designed to be used in the reverse biased configuration. In this discussion, “reverse biased” shall describe a reference voltage applied to the p-type semiconductor and a corre- sponding positive voltage applied to the metal side of the rectifying contact (Figure 2.8). “Forward biased” shall be taken to mean the opposite (Figure 2.9). In the reverse biased configuration, a perfect Schottky device would not allow the flow of any current. In the absence of an applied external voltage, there will still exist a built-in potential between the metal and semiconductor junction due to a mismatch in energy levels in the band structure of each material. The electric field generated from this built-in potential exists entirely within the semiconductor. Within this space-charge region the acceptors are nearly completely full, and consequently current carrying holes within the valence band of the semiconductor are nearly all excluded from this region. Thus, holes cannot traverse from the semiconductor to the metal without energy sufficient to overcome the barrier, approximately the difference between the work functions of the semiconductor and metal, respectively [4]. This is presented in detail below. As reverse bias is applied, electrons in the semiconductor are attracted toward the junc- 27 tion through coulomb attraction and fill additional holes in the valence band of the semi- conductor. The occupancy of these additional states results in upward band bending in the semiconductor in realtion to the metal thus expanding the extent of the depletion region and effectively raising the height of the Schottky barrier [4]. Figure 2.7: Energy levels in the metal/p-type junction with no applied external potential [4]. Figure 2.8: Energy levels in the metal/p-type junction in the reverse biased configuration [4]. The reverse is also true. In the forward biased configuration, holes are injected into the semiconductor, raising the energy of the valence band. The Schottky barrier is effectively 28 lowered, and when external voltage is sufficient to overcome the Schottky barrier, current flows. In electronics design, the potential necessary to bring the diode into this conducting state is referred to as “turn-on voltage.” Once in a conducting state, the current through the diode is impeded only by the forward resistance of the device. For an ideal diode this value is zero. Figure 2.9: Energy levels in the metal/p-type junction in the forward biased configuration [4]. The easiest way to derive (and visualize; see above figures) the current properties of a metal-semiconductor junction is to analyze total current flow as the sum of four components: 1. Current I1 is comprised the flow of electrons from the metal into the semiconductor. 2. Current I2 is comprised the flow of electrons from the semiconductor into the metal. In equilibrium this is equal and opposite to I1 . 3. Current I3 describes the flow of holes from the lower portion of the metallic valence band into the lower band of the semiconductor. In the metal/p-type junction this is governed by a relationship similar to Richardson’s equation of thermionic emission [4]. 29 4. Current I4 describes the flow of electrons from the valence band of the semiconductor into the metal. The values of the currents are given, in units of amperes\cm2 , as [4]: −qDn −qθ I1 = √ · 2U T 3/2 e kT (2.1) Dn τn qDn −qθ I2 = √ · 2U T 3/2 e kT (2.2) Dn τn −q(φs −φm ) I3 = QT 3/2 e kT (2.3) −q(φs −φm ) I4 = −QT 3/2 e kT (2.4) Combining these current components, we can calculate saturation current, Isat = I2 + I3 . (2.5) And finally we can calculate total current in the presence of an external potential, Va , qVa   Itot = Isat 1 − e kT , (2.6) where Dn , the diffusion constant, is given by the Einstein relationship 30 kT Dn = µn (2.7) q cm2 µn = 2000 (Carrier mobility [55]) (2.8) V ·s τn = 1.2 ∗ 10−11 s (Carrier lifetime [29], [56]) (2.9) q = 1.602176634 ∗ 10−19 C (Elementary charge) (2.10) amperes Q ≈ 60 (Richardson constant [57]) (2.11) cm2 K 2 1 U= (Constant relating energy levels in upper to lower bands [12]) (2.12) T 3/2 T = 300K (Temperature; calculation done at room temperature 300K) (2.13) k = 1.380649 ∗ 10−23 (Boltzmann constant) (2.14) φm,s ≡ work function; subscript denotes m = metal, s = semiconductor (2.15) θ ≡ distance between bottom of upper band and Fermi level; see Figure 2.7 (2.16) While this is, strictly speaking, all that is necessary to calculate the I-V performance of the presently considered detector, it relies on values which must be obtained through additional assumptions or calculations. A value of φs = 3.77 eV was determined for the work function of the lightly boron doped diamond. This value is based on the work of Diederich, et al. (1998) [58], who found a relation between work function with crystal orientation, annealing temperature, and growth termination method. Using their tabulated results for the (100) face with low boron doping and hydrogen termination, we perform a linear regression to determine the relation between work function and annealing temperature in degree centigrade, φs = 0.0014(T ) + 3.7363 [eV] with a mediocre goodness of fit r2 = 0.8181. Using this formula we arrive at φs = 3.77 eV for the (100) face hydrogen-terminated 31 unannealed CVD-grown single-crystal diamonds used in the device manufacture. Despite the lackluster goodness of fit in the linear regression, the calculated value fits nicely with the data of Diederich, et al. (1998) [58] which included measurements of the diamond work function following anneals at a variety of temperatures, reporting a value of φs = 3.9 eV for the hydrogen-terminated unannealed (100) face. As acceptors are added to the semiconductor through p-type doping, one may intuitively suspect that the thermodynamic energy required to add an additional electron to the semi- conductor decreases. This requisite amount of energy is known as the Fermi level and for the purposes of this report will be measured and/or reported in terms of energy above the top of the valence band and denoted by ψ. Literature identifies the Fermi level in diamond between 1.7 → 2.0 eV with strong Fermi level pinning [59] [60]. This suggests that, in contrast to the naive picture, doping concen- tration and the choice of metal at the metal-semiconductor have little to no effect on the Fermi level. More recent research, however, presents calculations suggesting that Fermi level pinning in diamond is not as strong as previously reported, and cites additional research that there is a correlation between doping concentration and Fermi level [58]. In the face of conflicting literature, it was deemed judicious to perform our own cal- culations for the basis of our detector diode performance specifications. This was initially planned to be done through a series of density functional theory simulations, but given time constraints of manufacture, a simple model relating the number of acceptors and Fermi-Level was adopted [4], simply to sanity-check design specifications: Fermi Level, ψ = .6221∗Na−.02 , where Na denotes the acceptor concentration per cc of semiconductor material. Gadolinium (Gd, atomic number 64), an Inner Transition Metal, was chosen for the metallic portion of the metal-semiconductor junction. The rationale for choosing this ma- 32 terial was driven by nuclear considerations (discussed in Section 2.1), though it performs suitably in the diamond Schottky diodes. Gadolinium has a work function of φm = 2.90(6) eV [61] and electrical resistivity of ρ = 1310nΩ · m [62]. This high resistivity compared to a True Metal suggests Gd will not perform as well as a True Metal in a Schottky diode, increasing forward biased resistance. However, the nuclear property advantages makes this trade-off desirable in the pursuit of detector sensitivity to ultra-fast neutrons. As designed, the diamond Schottky diodes are calculated to have excellent current-voltage properties, shown here in Figure 2.10. Calculations were performed both through direct cal- culation of performance characteristics and through finite element analysis to solve the semi- conductor drift-diffusion-Poisson equation. This was accomplished using the NIST Sesame Optoelectric modeling software [63][64]. Table 2.3 shows the calculated properties of the diamond diodes as designed. 33 Figure 2.10: Calculated I-V curve of diamond Schottky Diode obtained through finite element solution of the Poisson drift-diffusion equation. The negative current shown in the forward- biased condition is by sign convention choice. 34 Table 2.3: Semiconducting and electronic design values of POSSM detector diamond diodes. Property (Unit) Value Schottky Barrier Height (eV) 0.8521 Fermi Level (eV) 0.3118 Built-in Potential (V) 1.47 Depletion Depth at 150V reverse bias (µm) 95.96 Reverse Bias Leakage Current (A/cm2) 7.19E-10 Forward Bias Turn-On Voltage (V) 3.5 Ideality Factor (unitless) 6.98 M-S junction capacitance (pF) 2.83 Electronic grade HPHT single-crystal diamond plates were purchased from New Dia- mond Technology. A single plate with nominal dimensions of 3.5 x 3.5 x 1 mm was chosen for construction of the diamond diodes. The plate was laser cut resulting in two plates each of approximately 400 µm thick. The plates were mechanically polished and CVD growth of boron doped diamond was conducted on each side of each plate, a lightly doped semi- conducting P- layer on one side, and, to guaratee uniformity in the electric field induced in the intrinsic layer under bias, a heavily doped near-metallically conductive P+ layer on the other side. All sample preparation and growth was conducted by staff of the Engineering Research Center, Microwave and Plasma Processing Laboratory, located aboard Michigan State University. The diamond growth results are tabulated in Chapter 3, Table 3.1. See Chapter 3 for full diamond growth and device construction details. 35 The diamond diodes were then assembled using a piece of 10 µm thick Gd foil as a shared cathode. (See Appendix A for as-built drawings). The foil was mechanically fastened between the P- diamond material via a plastic bracket to a pressure of 16 psi, limited by the yield strength of the #2-56 nylon machine screws used to fasten it. A tail of Gd foil was mechanically clamped to Au foil which served as a solder pad for connection. Similar Au foil pins were fashioned as cathodes. These were, likewise, mechanically clamped to the P+ sides of the diamond plates ensuring an ohmic connection. The final product is the active portion of the POSSM Neutron Detector, electronically speaking, two mutually opposed Schottky diodes. The diodes were then characterized through I-V measurements and Capacitance mea- surements. Diode ideality factor was calculated. Diode ideality factor is a measure of how well a diode follows the ideal diode equation proposed by William Shockley:  Vq  I = Isat e kT − 1 (2.17) where I ≡ Current through diode (2.18) Isat ≡ Reverse bias saturation current (2.19) V ≡ Voltage across diode (2.20) q ≡ Elementary charge (2.21) k ≡ Boltzmann’s constant (2.22) T ≡ Temperature. (2.23) 36 At normal temperatures eV q/kT −1 ≈ eV q/kT so the diode equation can be approximated as a linear relation ln(I) = ln(Isat ) + VkTq for a perfect diode. For a real diode this equation takes the form of ln(I) = ln(Isat ) + n1 VkTq with ideality factor n. This can be obtained from a linear regression of ln(I) vs V . Calculated in this manner, ideality factor assumes a diode of a single junction, either P-N or M-S. In a more complicated device, such as the POSSM dia- mond diodes, the overall ideality factor is a sum of the ideality factors of each junction. The diodes discussed in this work contain two homojunctions, one heterojunction, and one ohmic contact. The ohmic contact does not contribute to diode behavior, so overall ideality factor is the sum of the ideality factors of the homo- and heterojunctions in the device [65]. For a single diode, ideality is pinned between values of 1 and 2, caused by competition between carrier drift diffusion and the Sah-Noyce-Shockley generation-recombination process. The theoretically calculated ideality factor of 6.98 suggests ideality factors of slightly over 2 for each of the homo- or heterojunctions. This would suggest recombination manifests through band to band and Shockley-Read-Hall processes, governed by both carrier types [65]. This is in line with other studies of M-S junction diodes [66], though deviates considerably from the measured properties of the POSSM diodes. 2.2.2 Measured Diode Performance The results of experimental diode characterization of the POSSM diodes are tabulated in Table 2.4. Measured I-V curves are shown in Figure 2.11. The measurements were recorded at 19◦ C. 37 Table 2.4: Characterization measurements of mutually opposed diamond Schottky diodes. HMT 001a (Ch. 1) HMT 001b (Ch. 2) Forward Bias Turn-On 5.1 V 5.1 V Ideality Factor 193 193 Total Capacitance (0 bias) 94 pF 93 pF Reverse Bias Leakage Current (-150V) Undetectable (< 0.1mA) Undetectable (< 0.1mA) Figure 2.11: Measured I-V curves of the POSSM diodes. Both diodes show the characteristic exponential increase in current with increasing forward bias. The dashed lines serve a guide for the eye. These diodes show considerably worse performance than predicted in the design calcu- lations in forward bias. The high ideality factors are indicative of partially ohmic behavior. There are two reasons for this: First, the Gd cathode resists current in the forward bias con- figuration compared to a True Metal. Secondly, the mechanical fastening method does not 38 guarantee a uniform, optimal electrical connection. Applying additional pressure to the me- chanically clamped bracket improved the ideality factor by nearly 50%. In future iterations the clamping bracket should be redesigned to apply additional pressure to the diode stack. Further it is likely that the Gd foil oxidized resulting in a (semi)metal oxide-semiconductor junction rather than a (semi)metal-semiconductor junction. Recent advances suggest that this may actually be beneficial, resulting in a higher Schottky barrier [67]. This may simul- taneously increase junction capacitance because GdO3 will become more strongly polarized under bias [68]. However, the electrical effects of GdO3 specifically (or transition-metal oxides in general) are not fully understood at this time. The existing body of work has demonstrated p-type diamond Schottky diodes with ex- cellent forward biased performance and report ideality factors, n ≈ 3. However, these works suggest forward biased performance is heavily temperature dependent, improving with tem- perature increase [23][60]. Notably, even these high performing diodes passed significantly less forward current in room temperature conditions, in agreement with observations of the POSSM diodes. Ultimately, it is assessed that a different choice of metal and construction process would dramatically improve forward biased performance, if that were the objective of the work. The POSSM diodes perform acceptably well in reverse bias, which is the configuration in which they will be used as detectors. Throughout the range of test voltage (up to 150 V reverse bias) the reverse bias leakage current was undetectable in each diode via direct measurement. Small changes were made to the diode package design in the POSSM Prototype Design Rev. 1 package, to improve reliability and durability. These changes consisted of replacing the gold contact interfaces at the P+ material with conductive paint and relocating the nylon 39 machine screws to provide more even pressure on the diode package. See Appendix B. 2.3 Preamplifier Design Onboard electronics were designed to optimize pre-amplification of the diamond diode signal. The electronics designed and built for the POSSM Prototype Design Rev. 0 consisted of two (one for each detector channel) single stage operational amplifier circuits. See Appendix C for circuit diagrams. 2.3.1 POSSM Prototype Rev. 0 The preamplifier electronics for the Rev. 0 prototype were built around the Analog Devices AD823 rail-to-rail operational amplifier following the architecture of a typical photodiode trans-impedance amplifier. The Rev. 0 circuitry design was obtained largely through an iterative design process achieved through breadboard testing and revision. It became apparent shortly into the design process that the circuit would have to be a compromise of speed, gain, and noise. While this may be true of any amplifier circuitry, it was particularly salient in this design because the AD823 is not situated for optimal performance in any of these metrics in the regimes required for a particle detector preamp. The completed Rev. 0 circuitry was tested using off-the-shelf Infineon Technologies Schot- tky diodes (PN# IDH03SG60CXKSA2) in place of the POSSM diodes which had not yet been prepared. The design provided a maximum gain of 5 with a rise time of approximately 60nS. In order to achieve this rise time with the AD823, a feedback capacitance of 1pF was required. This increased peaking, but it also resulted in significant backswing. See figure 2.12 for the observed preamp response. 40 Figure 2.12: Observed preamplifier response using the Rev. 0 circuitry. The input signal is a 2µs square pulse. The time scale is 5µs/Div. Once the POSSM diodes had been prepared and inserted into the detector, informal tests were conducted using a 5µCi source. Sufficient testing was conducted to show that the detector would produce a signal on an oscilloscope. No further effort was made to match the preamplifier circuitry with the POSSM diodes junction capacitance. Neutron beam testing of the POSSM device was conducted at the Weapons Neutron Research Facility (WNR) located at the Los Alamos Neutron Science Center (LANSCE) aboard Los Alamos National Laboratory and demonstrated overall success of the detector and preamplifier. The testing did however, accentuate shortcomings in the electronics design, namely in the response speed of the amplifier, necessitating a Design Revision 1. Please see 41 Chapter 4 for a full discussion of detector testing and results. 2.3.2 POSSM Prototype Rev. 1 The Design Rev. 1 preamplifier consists of multi-stage trans-impedance amplifier circuitry based around the Analog Devices AD8056 two-channel high speed video operational ampli- fier. The previous Design Revision was based on the Analog Devices AD823, but preliminary field testing showed this opamp to be too slow for the requirements of the POSSM detector, necessitating the Design Rev. 1. The Design Rev. 1 preamplifier circuitry benefited from a more concerted and rigorous design process. Field testing did not simply reveal shortcomings in the Design Rev. 0, but also provided the basis for a set of design specifications dictating the required performance of the device. Specifically, the Design Rev. 1 preamplifier was required to produce a pulse with a rise time of less than 20ns and pulse width of 45ns (a more than 3x bandwidth improvement over the AD823). Over-volt and over-current protection were required for both the high and low voltage buses. Low voltage input was required to be single-sided supply rather than requiring a dual supply of both positive and negative voltages. The design process incorporated software aided design, significantly reducing cost and turn-around time compared to the time intensive method of, what was essentially, trial and error on a breadboard employed in the Design Rev. 0. The Design Rev. 1 was only breadboarded once a suitable design had been found in software. The Design Rev. 1 preamplifier circuitry produces an inverted wave form with a 12nS rise time at a bandwidth of 52.6MHz, a signal-to-noise ratio of 52.7, and a 66.9mV output offset. Total noise was calculated to be 615µV rms. Trans-impedance gain was calculated to be 5.01kV/A at the design operating frequency of 52.6MHz [69][70][71]. These parameters 42 meet or exceed the performance requirements. See the drawings in Appendix C for circuit schematics and Design Rev. 1 drawings in Appendix B the perfboard layout. Bode plots are presented below in Figures 2.13 and 2.14. Figure 2.13: Calculated preamp circuitry frequency response. Figure 2.14: Calculated preamp circuitry noise gain. Both MHV and SHV terminals have been provided for flexibility in experimental set-up. Additionally, high voltage diodes were incorporated into the high voltage circuit, in series with each HV bulkhead connector. This is a safety measure to prevent the energization of 43 whichever connector is not in use. The high voltage input is protected via a 360V crowbar circuit triggered by a 600V Silicon Controlled Rectifier (SCR) and zener diode array. Trig- gering the crowbar draws excessive current through the voltage input firing a quick-blow fuse rendering the circuit safe. A snubber capacitor prevents spurious triggering of the circuit, such as when switching on the power supply. Over-current protection against a possible short failure in the diamond diodes is provided through the use of a 1/2 Watt, 50MΩ resistor. This resistor is capable of reducing the current through a short circuit to the µA level and dissipating the corresponding power. Over-voltage protection in the low-voltage input is provided through similar means. A low voltage crowbar circuit is designed to trigger at 30 volts using the same mechanism as the high-voltage crowbar. The single-side supply feeds two VX78-500 DC switching regulators. One is configured to provide +5V in reference to common ground and the other is configured to provide -5V. These outputs provide the dual-sided power requirements of the AD8056 opamp packages. Low voltage over-current protection is provided through a quick-blow fuse on the input. See Appendix C. Board layout was done using best practices for signal fidelity, minimizing interference, cross talk, and parasitic capacitance [72]. 44 Chapter 3. Detector Construction 3.1 Diamond Growth Detector construction began with the purchase of HPHT diamond seeds. The diamonds selected were single crystal electronic grade intrinsic diamond plates from New Diamond Technology. These samples represent some of the highest quality intrinsic diamond available and show minimal defects and low ppb nitrogen impurity. The diamond plate purchased was nominally 3.5 x 3.5mm in cross section and 1 mm thick. The plates were oriented with the faces aligned to the (001) crystal plane ±3◦ , and the edges aligned along the (110) and (011) planes. The plate was then cut in half along the (001) plane using a Nd:YAG Bettonville Ultra- shape 5xs IR cutting laser, resulting in two plates each of approximately 3.5 mm x 3.5mm x 0.45mm. The laser cutting kerf accounts for roughly 50µm of material lost. Additionally, the laser cut results in some amount of surface damage which must be polished before CVD growth can be conducted. The plates were designated HMT 001a and HMT 001b. Mechanical polishing was conducted using a Coborn PS1 polishing bench with a cast iron scaif. The scaif is seeded with 5-10µm natural diamond powder and a few drops of polishing oil, and then spun at 3450 RPM, resulting in contact speeds of between 28.5 - 32.5 m/s, depending on radial distance. Planetary motion of the polishing wheel slowly adjusts radial distance to reduce polishing grooves and preserve scaif lifetime. The sample was loaded to 500 kPa pressure against the polishing wheel [73]. Both of the laser cut faces (one on each sample) were polished. The untouched faces were determined to be sufficiently smooth as provided from New Diamond Technology. Following 45 the polishing process, each plate was approximately 400µm thick, the polishing process removing the remaining damage from laser cutting. The diamonds were then cleaned in a first wash of hydrochloric acid and a second wash of combined sulfuric/nitric acid. All cutting, polishing, and cleaning activities were performed by staff engineers at Fraunhofer USA - Center for Coatings and Diamond Technologies, co-located aboard the Michigan State University campus in East Lansing, MI. After cutting and growth preparation (see Figure 3.1), the P+ (heavily boron doped) layers were grown in the ⟨001⟩ direction. The homoepitaxial growths were conducted in a completely bespoke MSU-built tunable cavity type 2.45 GHz microwave plasma assisted chemical vapor deposition system (MPACVD). This reactor, designated DS1, is a dedicated heavy boron doping system using 1000 ppm diborane in hydrogen gas and methane as process gases. A picture of such a diamond in the growth process is shown in Figure 3.2. Each diamond had 1.4 ± 0.3µm of ∼1020 B/cc P+ single crystal diamond grown. The thickness of the material added was determined through high precision mass measurements [73]. Fraunhofer USA personnel conducted microscopic photography analysis using Differential Image Contrast Method (DICM) and performed 4-point sheet resistance measurements of the P+ growth. 46 Figure 3.1: HTHP Intrinsic diamond (001) face prepared for growth. Images taken at 25x using DICM. a) HMT 001a b) HMT 001b. Water spots are visible in the corners of HMT 001a, redisue from cleaning. Following completion of the P+ growth, the diamonds were again cleaned and prepared for P- growth. The intrinsic side of each diamond (the face opposite the P+ growth) received approximately 3µm of lightly doped, ∼ 1015 B/cc, P- single crystal semiconducting diamond. The P- growths were conducted in a heavily modified Lambda Technologies MPACVD system [73][24]. 47 Figure 3.2: HTHP diamond plate undergoing homoepitaxial CVD growth. Note the diamond plate (square object in center of the frame) lies countersunk in a molybdenum sample holder. The intensity of the plasma can be seen to saturate the sensor of the camera near the center top of the image. Final growth results are tabulated below in Table 3.1, and DICM results are shown in Figure 3.3. Table 3.1: Measured electric properties of semiconducting diamonds. HMT 001a (Ch. 1) HMT 001b (Ch. 2) P+ thickness 1.4(.3)µm 1.4(.3)µm P+ concentration 1020 cm-3 1020 cm-3 P+ sheet resistivity 908 Ω/sq 994 Ω/sq P+ resistivity 0.13 Ω·cm 0.14 Ω·cm P- thickness 3.49 µm 2.79 µm P- concentration 1015 cm-3 1015 cm-3 48 Figure 3.3: Single crystal growth. Images are taken via DICM at 25x. a) HMT 001a P+ b) HMT 001a P- c) HMT 001b P+ d) HMT 001b P- From the figure above, it is apparent that the P- growth was not as smooth and uniform as the P+ growth. There are a variety of explanations for this, though an exhaustive analysis of the growth process is outside of the scope of this work. Due to time constraints, no further action was taken into characterization of the diamond growths, exempli gratia, X-ray Diffraction or Birefringence measurements. See references [25][74][75]for more discussion on these topics. 3.2 Diode Assembly The POSSM diode assembly was constructed using the semiconducting single crystal dia- mond plates described in Section 3.1 and a two-piece 3D printed bracket. Because of material availability, polyactide palastic (PLA) was used in the prototype. This deviates from the 49 design specification of PEEK plastic. In a producction environment this substitution is ex- pected to reduce longevity of the device by nearly 50%, however it is not assessed to have a negative impact in the performance of the prototype[76][77]. The diode bracket used in detector construction was an older design and included a key-way cut into the top bracket piece and fastening screws located 1mm off center. This necessitated a small shim to ensure proper pressure was applied to the Gd-Au cathode junction. This is shown in the Rev. 0 drawings in Appendix A. Refer to Appendix B for design drawings which show the proper diode bracket which does not necessitate the shim. The Au anode connection pins, the contacts contacting the P+ diamond, were fashioned from a 0.01” thick gold foil (McMaster-Carr product# 6375N113). The foil was cut by hand using conventional kitchen shears, into a “Γ” shape of roughly 1mm in width. The gold contacts were then bonded to the plastic bracket using cyanoacrylate adhesive. The diamonds were placed with the P+ side on the Au contacts and bonded in place by a small amount of Room Temperature Vulcanizing (RTV) compound, an electrically insulating rubber, along the edges of the diamond. Care was taken to ensure that RTV did not seep between the diamond and the Au contact surface, nor get on the exposed P- surface. This is seen in Figure 3.4. Note that the legs of the Au anode pins have been bent outboard for increased physical separation to prevent inadvertent shorts during the assembly procedure. These can be bent back straight after assembly. Also note that solder has been adhered to the Au leads to facilitate future connections. Sn-free solder should be used in contact with Au because the gold is soluble in liquid tin and can be dissolved. A solder temperature of 300◦ C was used for these contacts. At this point, a nominally 5mm x 10mm piece of Au foil was cut and bonded to the flat piece of the two-part bracket using cyanoacrylate adhesive, electrically isolated from the 50 anode pins. This piece forms the basis of the shared cathode pin. Figure 3.4: The two halves of the diode bracket before final assembly. The Au contacts have been bonded to the plastic bracket pieces, and the diamonds have been fastened, P+ in contact with Au, with a small amount of RTV on the sides of the diamonds. Care has been taken to ensure RTV does not contaminate either the P+ or P- sides of the diamonds. A piece of 10µm thick gadolinium foil, 3.5mm x 10mm, was cut using conventional kitchen shears. The two halves of the diode bracket, each with one diamond were then joined using the specified nylon machine screws. The piece of Gd foil was placed between the diamonds before the diode pack was fully fastened. Tightening the nylon machine screws to the specified torque, sufficient to exert 16psi of pressure on the diamond sandwich, squeezes the Gd foil between the P- layers of the diamond forming rectifying contacts. The raised edge of the smaller piece of the detector bracket clamps the tail end of the Gd foil against the Au Cathode pad forming an Ohmic connection. Once assembled, the excess Gd tail can be trimmed using a sharp razor knife. See Appendix A for design drawings. When not actively being handled, the Gd foil was stored in 99% isopropyl alcohol to prevent oxidation of the material. 51 Once the POSSM diode pack had been physically assembled, remaining voids inside of the detector bracket were filled with RTV to provide mechanical strength to the internal components, provide mechanical shock resistance, hermetically seal the internals, and prevent shifting of internal components degrading performance or causing short circuits. The RTV was given 48 hours to fully vulcanize before the diodes were tested. The diode pack was assembled and tested at Hillsdale College, located in Hillsdale, MI. 3.3 Device Assembly The frame for the POSSM detector was cut from 1/16” G10 fiberglass board. The two frame pieces provide physical strength for mounting the detector, protect the detector electronics, and provide a surface for clamping or otherwise physically mounting the device. All onboard electronics were assembled on a single piece of 7cm x 9cm single sided Bakelite perfboard using conventional 0.1” pad spacing. Where available, through-hole components were used. DIP-8 sockets were installed where necessary for easy insertion and replacement of PDIP-8 operational amplifiers or DIP-8 adapters as appropriate. The completed preamp circuitry is shown below in Figure 3.5. Connections between components were made using 22AWG enamel coated copper magnet wire. Junctions were soldered using rosin core 60/40 Si,Pb electronics solder at a temperature of 350◦ C. The diamond diode pack was pressure fit into a rectangular hole in an extended protrusion in the G10 base. The diode pack was fixed in place with a small amount of cyanoacrylate adhesive. The diode leads were connected to the onboard preamp circuitry via 22AWG enamel coated magnet wire. RoHS compliant components were used throughout, but assembly using leaded solder means the detector is not RoHS compliant. However, it is 52 exempt from RoHS regulation under Article 2, paragraph 4(j) [78]. Figure 3.5: Completed preamp circuitry built on Bakelite perfboard. Through-hole compo- nents occupy the top of the board. All traces are connected on the reverse side. Two versions of the on-board circuitry were designed. The first was based around the Analog Devices AD823 operational amplifier and built as described above. Preliminary field testing showed this preamp was too slow for the application. The preamp circuitry was redesigned to the parameters described in Section 2.3.2, though this updated version has not yet been built. The G10 frame and perfboard were joined according to design drawings using #4-40 nylon machine screws and 1/4” spacers, completing final assembly of the detector. The 53 finished detector is pictured below in Figure 3.6. All detector and electronic assembly was conducted at Hillsdale College. Rev. 0 draw- ings (as-built) are presented in their entirety in Appendix A and electronic circuit design schematics are presented in their entirety in Appendix C. Figure 3.6: Completed POSSM detector. Note the G10 panels and perfboard fastened by nylon machine screws. The Diode pack is blue, located near the top at the end of the G10 extension. 54 Chapter 4. Detector Testing 4.1 LANL/LANSCE/WNR August 2022 In August 2022, The POSSM detector was tested at Weapons Neutron Research Facility located at the Los Alamos Neutron Science Center, part of Los Alamos National Laboratory. The testing was done in conjunction with the MoNA collaboration. The POSSM detector was placed in parasitic mode with the MoNA experiment, itself running in parasitic mode behind an unrelated experiment. The detector was set up in a target shed located 90 meters from the spallation neutron source and approximately 15◦ off-axis, exposed to a neutron beam with a nominal energy spectrum of 20 to 800 MeV. (The small off axis angle of the target shed results in a maximum neutron energy of approx. 790 MeV at the detector.) Due to schedule constraints, no local testing of the fully-assembled POSSM detector had been conducted. Consequently, the first goal at LANSCE was simply to see if there was a signal from the detector. To get the purest beam possible without the interference of any charged particle production from neutron interactions with the MoNA scintillator detectors, the POSSM detector was placed ahead of the MoNA detectors, approximately 6 inches from the neutron collimator. In this configuration signals were observed from both detector channels on an oscilloscope, and energy spectra were observed in real-time on the DAQ computer. The energy spectra showed approximately the expected shape, and, as expected, a charged particle peak was observed in channel two. At this point, no data were recorded electronically though quick sketches were recorded by the author. The sketches of the DAQ histograms recorded in the author’s notebook are shown in Figure 4.1. 55 Figure 4.1: The author’s observation of the POSSM Detector’s energy response during op- erational testing. No data were recorded during this run. Following the successful operational test of the POSSM detector, the device was attached to the rear of the MoNA scintillator stack, so as to avoid interference in their further data collection. The experiment was facilitated by equipment which was sourced locally (from in and around LANSCE Bldg. 17). All of the available high voltage supplies and cables provided SHV connectors, and the POSSM Prototype Design Rev. 0 included only MHV high voltage input. Luckily an MHV-SHV adapter was sourced to allow the requisite HV hook up. Two laboratory bench type DC power supplies (of different models) were acquired to provide low voltage power to the preamplifer. Due to the author’s lack of familiarity with one of the power supplies used, during the reconfiguration of the detector position, he inadvertently put 48 VDC across the +12VDC preamplifier input. The power state was immediately remedied, however the excessive voltage had permanently damaged one of the AD823 opamps at the heart of the preamplifier. During further data collection, no signal could be obtained from channel 1, though channel 2 functioned without error. The MoNA Collaboration experiment involved using two Cividec Instrumentation B8 open-design diamond detectors, commercially available single crystal diamond detectors de- signed for alpha and heavy ion spectroscopy [79] placed in the neutron beamline one approx- imately one half inch behind the other. These detectors provide a good frame of reference 56 with which to compare the performance of the POSSM device. Though not designed for neutron detection, the dimensions of the Cividec B8 diamond is nominally 4mm × 4mm × 500µm, very comparable to the size of each POSSM channel. Only data from the second of the two MoNA diamond detectors will be compared, as this is the most equivalent to the channel 2 POSSM diode. The experimental data were analyzed using the ROOT Data Analysis Framework, version 6.26/04 [80]. Figure 4.2 shows the uncalibrated q-spectrum, and Figure 4.3 shows the uncorrected time of flight for a complete 1800 ns pulse of the WNR beamline. Presented here is the raw data, uncalibrated and uncorrected, so only generalized com- parisons can be made. In the case of the q-spectrum, the POSSM device shows a peak at an arbitrary value of 10 as opposed to the commercial diamond detector’s peak at or near zero. This suggests that the POSSM detector is more sensitive to higher energy neutrons. Additional evidence of higher energy sensitivity is found in the integral of the histograms, shown in Figure 4.2. In the case of the POSSM detector, 20,057 events result in an area under the curve of 15,541, contrasted with the commercial diamond detector’s 52,155 events and 2,338 integral. These numbers demonstrate the POSSM’s superior charge conversion and collection capabilities, attributable to the presence of the activation layer and Schottky barrier architecture. While it is not as pronounced as observed in the first round of operational testing, the bimodal nature of the Channel 2 q-spectrum as shown in Figure 4.1 is possibly visible in Figure 4.2(a) and this figure shows a clear broadening on the right shoulder of the curve. It is believed this is due to the detection of charged particles generated in the activation layer. 57 Figure 4.2: Uncalibrated q-spectrum of Figure 4.3: Uncorrected time of flight (a) POSSM Detector Ch. 2, and (b) 2nd spectrum of (a) POSSM Detector Ch. 2, Cividec B8 Diamond Detector and (b) 2nd Cividec B8 Diamond Detec- tor The POSSM’s time of flight spectra suggest peak neutron detection between 1000-1200 ns. Based on published time of flight data for the WNR beamline [81], this corresponds to neutron energy 20 MeV < En <100 MeV . Gamma flash is expected at 300 ns [81]. This is not distinctly evident in the POSSM’s time of flight spectrum, however there is a broad peak centered approximately at 200 ns. It is likely that this peak is attributable to gamma flash and “wrap-around” neutrons and the broadening is a result of the preamplifier being too slow to cleanly resolve events. The time of flight spectra generated from the commercial diamond detector’s data does not show any clear structure, though a small peak is seen at 300 ns. 58 The testing that has been conducted with the POSSM detector at the time of writing is inconclusive, and more testing with a fully functioning device is required before any definite claims of performance can be made. Nevertheless, the limited data available suggest that the POSSM detector performs at or near its design specifications. It appears to offer increased efficiency in the detection of ultra-fast neutrons over presently available diamond detector technology. The activation layer appears to function as designed, increasing neutron inter- action, manifesting in the greater charge conversion observed in the POSSM’s q-spectrum. Proof of concept of the POSSM architecture has been demonstrated. Additional work in the development of this technology is discussed in Chapter 9. 59 Chapter 5. Quantum Color Centers in Diamond Quantum color centers (hitherto referred to simply as color centers) observed in diamond are a point-defect similar to F-centers observed in alkali metals whereby one or more un- paired electrons occupy a vacancy in the crystal lattice. Color centers demonstrate a variety of interesting, unique, and useful properties of interest for quantum technologies such as computing and encryption, or hypersensitivity to changes in electric or magnetic fields. Color centers absorb light at frequencies at which the crystal lattice is transparent and exhibit photo-luminescence (PL) at a different frequency, earning the “color center” name. 5.1 Point Defects In Diamond There are three unique single-atom point defects which can exist in the crystal lattice. They are the vacancy, substitution, and interstitial. These defects are represented in Figure 5.1. The vacancy is, as its name suggests, a vacancy in the lattice; one of the host material’s compositional atoms has been removed from its lattice site leaving a void. The substitution or substitutional ion, is a foreign atom which has replaced one of the host material’s com- positional ions in the lattice. Exempli gratia, a nitrogen atom which has replaced a carbon atom in diamond is said to be substitutional. Finally, an atom can reside in the bulk of the material, but not at a lattice site. Such an atom residing in between lattice sites is an interstitial. 60 Figure 5.1: Point defects in diamond lattice. a) Vacancy. A carbon atom (shown in purple) has been removed from the lattice leaving a void. b) Substitution. A foreign ion (shown in yellow) has replaced a carbon atom at a lattice site. c) Interstitial. The interstitial atom resides in between lattice sites, loosely bound to the surrounding atoms. Images generated using CrystalMaker software [1]. When concentrations of specific substitutional ions and vacancies coexist in the lattice, unpaired electrons belonging to the substitutional ions(s) can occupy the vacancy. This results in a “pseudo-covalent” bonding between the vacancy and the electron’s donor ion. This has been described as a “solid-state artificial atom” with ionic charge and a hierarchy of discrete internal energy states [82][83]. A substitutional N atom (five valence electrons) located adjacent to a vacancy will form three covalent bonds with its three nearest neighbor C atoms. Its remaing two valence electrons are a lone pair configured towards the vacancy and overlapping with the three dangling bonds of the vacancy-adjacent C atoms, as shown in figure 5.2. This, in effect, appears as a covalent bond between the substitutional N and the vacancy defect [84]. More complex color centers have been documented, for instance, the nickel-nitrogen-vacancy center (NE8) which consists of a substitutional nickel ion and four substitutional nitrogen ions bonded to each other and two vacancies [85]. 61 Figure 5.2: Pseudo-covalent bonding. The dangling bonds left by the vacancy are indicated in yellow. The substitutional nitrogen has five valence electrons. It is covalently bonded with neighboring C atoms, leaving a lone pair, indicated in red, oriented towards the vacancy. This orientation acts like a covalent bond between the substitutional nitrogen and vacancy. 5.2 Properties and Applications of Diamond Color Cen- ters Color centers demonstrate a variety of physical properties which make their research of interest to applications in quantum technologies for industry, military, and research purposes. Some color centers, such as the negative nitrogen-vacancy (NV-) center and the silicon- vacancy (SiV-) center have been well characterized experimentally and theoretically [82] and show promise as high temperature qubits for use in quantum information processing. These and other color centers also have extreme sensitivity to changes in magnetic, electric, temperature, and strain fields providing a tool for precise quantum sensing in applications 62 such as single-molecule Nuclear Magnetic Resonance (NMR) measurements [85][86]. 5.2.1 Optical and Spin properties The NV- center is the most studied solid state single photon emitter and has been intensely studied for its use in quantum computing, favored primarily for its excellent single electron spin coherence time (T2 ) at high temperature [5]. Experiments have demonstrated values of T2 ≈ 0.6s at 77K [87], and the relaxation time of the population of spins towards equilibrium, T1 ≈ 5ms in natural diamond at room temperature [85]. This is sufficient time to manipulate the spins and entangle two NV- centers [88]. Entanglement of two NV- centers in two separate diamonds has been demonstrated at distances of up to 3m [86]. The spin state of the NV- center can then be read via light or RF fields. The ability to prepare, manipulate, and read out the spin state is necessary in the development of a functioning qubit. While the NV- center is undeniably a great qubit candidate, it is not without its chal- lenges for the application. The excellent spin-coherence time of the NV- center can be reduced by coupling to the nuclear spins of the 13 C present in diamond [5]. Further, only about 4% of the fluorescence of the NV- center occurs in the zero phonon line (ZPL) and other color centers may offer order-of-magnitude greater branching into the ZPL resulting in higher signal-to-noise ratio on readout [5][82]. Further, the NV- center has been shown to sometimes lose its charge state under optical excitation and consequently lose some of its desirable characteristics documented above [83][89]. As optical excitation is necessary to read out the spin state of the color center, the reading process risks destroying the qubit. Re- search suggests that post-annealing treatments may mitigate this issue, nevertheless this has prompted research into the suitability of other diamond color centers in quantum information [83]. 63 The Group IV vacancy centers (SiV-, GeV-, SnV-, PbV-) have been extensively charac- terized and show promise for use as qubits. Enhanced spin orbit coupling for the heavier elements manifests in a larger energy level splitting in the ground state for the singly neg- atively charged color center which increases with the mass of the substitutional ion. This results in larger spin coherence times at higher temperatures [5]. Because the Group IV color centers exhibit inversion symmetry, they have no permanent electric dipole moment, unlike the NV- center. This results in optical transitions less sensitive to electric field noise near the surface of the host material [82]. Resistance to electrical noise is essential for a micro- or nano-scale quantum device. The Group IV color centers suffer slightly lower PL quantum efficiency than NV-, but more than make up for this loss with order-of-magnitude greater branching ratio into the ZPL [82], exempli gratia, the SiV- color center has 60% of its luminescence in the ZPL [5]. Acoustic phonons can be detrimental to the performance of the Group IV color centers [5]. Acoustic phonons can serve as a dephasing mechanism effectively limiting the spin- coherence time of the color center. Several techniques are showing promise for reducing the phonon density of states and mitigating this problem which will be discussed in Section 5.3. The Group III vacancy centers (AlV-, GaV-, InV-, TiV-) have begun to be studied. A theoretical work published in 2019 relies on computer simulation and DFT to analyze the traits of these color centers [82]. This work determines that the Group III color centers share the split-vacancy stability of the Group IV color centers. In many respects they are similar to their Group IV neutral counterparts, but the Group III color centers have s spin-1 ground state. The Group III color centers show different symmetry depending on charge: D3d symmetry for the ±1 charge state and C2h symmetry for the −2 charge state. Both of these configurations have inversion symmetry precluding the existence of a net dipole 64 moment resulting in optical transitions with the same resistance to electric field noise as seen in the Group IV centers [5][82]. See Figure 5.3 for the atomic configuration of Group III and Group IV color centers. Figure 5.3: Color center structures. a) NV center b) Group III or Group IV color center where X∈{Al, Ga, In, Ti, Si, Ge, Sn, Pb} c) Single vacancy defect. Image taken from [5] 5.2.2 Sensitivity to Fields In addition to their quantum spin properties making them good candidates for high temper- ature qubits, diamond color centers also demonstrate potential as quantum sensing elements taking advantage of excellent sensitivity to changes in a variety of fields [85]. For example, the nitrogen-vacancy center has been shown to be able to detect nanoTesla scale magnetic fields with 10nm resolution [90]. This sensitivity and resolution is suitable for single particle observations [90]. Similar results have been shown in electric field sensitivity, again using the nitrogen-vacancy center to detect an electric field of 202 Vcm−1 Hz−1/2 . This is sensitive enough to detect the electric field induced by an elementary charge at a distance of 150nm [91]. Temperature nanoprobes have also been built around the NV center with a temperature resolution of 1mK and spatial resolution of a few tens of nanometers.[92] The specifics of the methods vary, but all these measurements take advantage of the unique spin properties 65 of the NV color center and use some kind of optical readout. For this reason, the same weaknesses of the NV center as discussed in Subsection 5.2.1 apply. This has motivated re- search into more exotic color centers to enhance desired sensitivities and reduce undesirable characteristics. Some exotic color centers of interest include the oxygen, calcium, fluorine, magnesium, phosphorus, osmium, and scandium color centers, among others [7]. 5.3 Production of Color Centers The formation of color centers in diamond require three things: vacancy defects, substitu- tional defects, and thermodynamic energy sufficient to form a bond between these defects. This can happen in nature, and indeed NV- is regularly found in naturally occurring Type Ib diamond, though the concentration of naturally occurring color centers is generally too low to be of use in quantum devices [83]. To create artificial color centers in diamond, one must first ensure a sufficient quantity of vacancies and substitutional ions. Substitutional ions can be placed through implanting via an energetic ion beam or included through doping during the diamond growth process. Vacancies are naturally created through the implant process as heavy implant ions knock carbon atoms out of the lattice. Other methods have been used to create vacancies such as irradiation with neutron beams [9][93]. Following irradiation/implantation the diamond samples are annealed at temperatures typically ranging from 600◦ C to 1500◦ C for a period of minutes up to days [58][94]. The purpose of annealing is twofold. First, annealing provides thermal energy to the system to promote pseudo-covalent bonding between the vacancies and substitional ions. Secondly, annealing at temperature enhances diffusion of defects in the lattice promoting 66 thorough mixing. As defects diffuse through the crystal lattice vacancies and interstitials can recombine, and vacancies and substitiutional ions can form color centers when diffusing in appropriate vicinity of one another. A detailed discussion of the physics of self-diffusion in diamond will follow in Chapter 7. Following the annealing process the diamond sample can be checked for the presence of color centers by observing PL induced through optical excitement. Color center concentra- tion can be mapped through this technique. It has been shown that reducing the phonon density of states will mitigate signal-noise and interference issues caused by acoustic phonons on the performance of Group IV color centers. One method to accomplish this is to reduce the concentration of defects other than color centers in the crystal lattice, effectively increasing the volume of the conventional unit cell relative to the size of the crystal sample. If the crystal were to contain zero defects, the unit cell would reduce to the primitive unit cell of the diamond cubic. If the crystal were to contain a single defect, then the unit cell must, by definition, be the entire crystal. Roughly speaking, the volume of the unit cell is the volume of the crystal divided by the number of defects present. This assumes a small, but non-zero number of uniformly spaced defects [5]. An alternative solution to reduce the phonon density of states is to use heavier implant ions, thus changing the properties of the resulting color centers [85]. The challenge of applying present color center technology to practical quantum devices lies in controlling the concentration gradient and orientation of the color centers. Some efforts have been made to improve the precision of manufacture, such as localized annealing through laser impulse heating. While this technique can, in principle, create a single NV- center at a desired location within about 200nm, literature reports success rates of less than 50% [95]. Little research has hitherto been done into laser annealing of the more exotic color 67 centers. 5.4 Motivation for Color Center Formation Simulation In an effort to start bridging the gap between simply producing and truly engineering color centers, a reliable process of developing an implant and annealing “recipe” to obtain a specified color center outcome must be developed. Such a process must be suitably fast and inexpensive to allow an iterative design process. To these ends a computer simulation is developed to simulate ion implantation, annealing, and color center formation in silico. Such a simulation will take advantage of existing Monte Carlo transport codes such as SRIM and GEANT4 to model ion implantation and then bespoke code to simulate the diffusion of defects, recombination of vacancies and interstitials, and the formation of color centers. Any such simulation must be able to handle fine meshes over a large domain, scale across multi-processing architectures, and provide an intuitive and easy-to-use user interface. A full description of the simulation and the newly proposed mathematical framework for meso-scale color center formation is provided in Chapter 7. 68 Chapter 6. A Comparison of Monte Carlo Frame- works 6.1 Motivation A color center formation simulation must begin with the implant of the desired ion(s) at the desired concentration(s). The post-implant condition of the crystal provides the initial conditions for the annealing process. An accurate simulation necessitates accurate initial conditions, so the question of accuracy must be posed to the available Monte Carlo frame- works. The implementation of stopping power calculations varies across commonly used general- purpose Monte Carlo codes. This results in sometimes widely varying results, particularly when approaching the lower (E < 1 MeV) energy regime [96]. As this energy regime is particularly important to the study of materials science, space science applications, radiation damage assessments, semiconductor device development, and sample preparation, accurate simulation is essential to a variety of problems. Previous work has investigated and compared different Monte Carlo frameworks [96][97]. Much of this work has focused on an energy regime (E > 1 MeV) [96] that is too high to be of application to many segments of the materials science, space science, and device manufacturing communities. Often, when the results of multiple Monte Carlo frameworks at lower energies have been compared, it has been necessitated by specific work, such as in the case of investing the structural integrity of graphite modified nuclear irradiation in high power graphite-moderated nuclear reactors [97]. Thus, these works offer individual cases rather than a comprehensive analysis of the performance of various Monte Carlo frameworks. 69 A more comprehensive survey can be of interest to researchers tasked with implementing Monte Carlo particle transfer methods in their works. The current chapter of this work compares the low energy results of two of the most commonly used general-purpose Monte Carlo frameworks. The authors analyze low energy ion implants into diamond, silicon, and germanium using the Stopping and Range of Ions in Matter (SRIM) 2013 code [98] and Geometry and Tracking 4 (GEANT4) version 4.10.07.p03 [43]. Additionally, whereas many materials of interest, including all of the Group IV semi- conductors, are defined by a crystalline structure, comparisons are also made using the GEANT4 Crystal Object (GECO) code, included in the GEANT4 kernel beginning with version 4.10.01 [44]. SRIM and GEANT4 do not model channeling effects which manifest in the crystal lattice, and the addition of GECO to the GEANT4 kernel was to “[allow] the manipulation of particle trajectories by means of straight and bent crystals and the scal- ing of the cross sections of hadronic and electromagnetic processes for channeled particles” [44][99][100]. 6.2 Methods In this work, four implant ions (H+ , 4 He+ , 14 N+ , 82 Pb2+ ) are evaluated at five incident energies ranging from 100 keV to 1000 keV. These implant schemes are repeated with targets of diamond, silicon, and germanium, using each of SRIM, GEANT4, and GEANT4/GECO. Implant profiles, vacancy profiles, and energy deposition along the beam axis are presented and compared among each other and against experimental data found in literature where appropriate and available. All simulations consisted of 10,000 incident ions. GEANT4 and GEANT4/GECO sim- 70 ulations and data analysis were performed on a Ryzen 3950x with 128 GB of DDR4 RAM running a highly modified Calculate Linux 64-bit operating system. SRIM simulations were performed on an Intel 9th generation Core i-7 running 64-bit Windows 10. A random sample of 45 simulations (25% of the total 180 simulations) was chosen, and these simulations were repeated two times additionally to infer inter-simulation statistics. In total, 1.8 TB of output data was produced and processed. 6.2.1 GEANT4 The GEANT4 simulations consisted of a rectangular-prism absorber of 9µm × 3µm × 3µm in the x, y, and z-directions, respectively. The beam propagates along the x-axis. In the case of 1000 keV protons, the targets were increased in size to 20µm × 20µm × 20µm to ensure the implants stopped within the target material without transmission. Amorphous absorber materials were sourced from the G4 NIST Materials Manager. The G4 Ge and G4 Si standard materials were used for germanium and silicon, respectively. A custom amorphous diamond material was built via an instance of the G4Material class, using carbon given a density of 3.51 g/cc [13]. The G4ParticleGun class is used to generate the desired ions, including specifying the charge state of the ion. A +1 charge state was specified for protons, alphas, and nitrogen. A +2 charge state was used for lead. Modular physics were employed, using the FTFP BERT HP physics list for hadronic and neutron physics. The GEANT4 EMStandardSS physics list was used for coulomb scattering and other electromagnetic interactions. The single scattering model was chosen for its ability to employ theory-based cross sections [101], rather than relying on empirical data, which is scarce for the particle/target/energy combinations of interest. Using the single scattering model, each elastic scatter is simulated, increasing the number of steps of charged particles, 71 negatively impacting CPU performance [101]. The discrete charged particle fidelity offered by the complete simulation of elastic collisions, however, is necessary for reliable simulation of vacancy creation. 6.2.2 GEANT4/GECO The same absorber geometry was applied to GEANT4/GECO as in the amorphous GEANT4 case. In implementing the crystalline structure, intra-lattice fields were calculated using DYNECHARM/ECHARM [102] software, or taken from the ECHARM author’s github [103] repository if available. Crystal structure was added to the absorber materials used in the amorphous case via the extensions to the G4Material class included in the GECO package. All crystals were oriented with the beam impinging on the (110) face. The same physics lists and G4ParticleGun instance were used as in the amorphous GEANT4 case. 6.2.3 SRIM SRIM 2013.00 was used in this work. SRIM is a collection of software packages which eval- uate the physics of ion penetration and transport in matter. It includes a Monte Carlo simulation of ion implants into layers of elemental or compound material. The material density, atomic dislocation energy, lattice binding energy, and surface binding energy are directly set by the user [98]. Crystal structure is not explicitly implemented in SRIM. Table 6.1 shows the material properties used in this work. SRIM does not allow for the explicit setting of ion charge state. 72 Table 6.1: Material properties used for target materials in SRIM simulations Material Element Density Displacement Lattice Binding Surface Binding (g/cc) Energy (eV) Energy (eV) Energy (ev) Diamond C 3.51 45.00 5.00 7.50 Silicon Si 2.32 15.00 2.00 4.70 Germanium Ge 5.35 15.00 2.00 3.88 Target depth was adjusted based on ion species stopping range at 1 MeV. Subsequent simulations at lower energies for each ion species were kept fixed at the selected 1 MeV implant depth. Full cascade damage modeling was used with 50 eV Energy Interval. 6.3 Results Raw output from all three simulation frameworks was post-processed to obtain x, y, and z profiles for implant ions and created vacancies. SRIM directly calculates vacancies and recombinations. Raw output from GEANT4 and GEANT4/GECO was parsed and lat- tice atoms which had been given energy following elastic scattering beyond the dislocation energy were assumed to form vacancies. Their initial positions (or vertex) positions were recorded. Ion positions were directly recorded from GEANT4 and GEANT4/GECO output and from the EXYZ.txt file produced by SRIM. Energy deposits were recorded from the histomanager.cc class in GEANT4 and GEANT4/GECO and recorded from EXYZ.txt in the case of SRIM. 73 Using the random sample of 45 implant scenarios, each with a total of three simulations, statistics were calculated to verify statistically significant reliability of the results. In other words, the author sought to confirm that the results reported herein are representative of their respective Monte-Carlo frameworks and in no case a quirk of the random number generation. In each case, the results of each simulation were within statistical agreement to a 95% confidence interval. The implant and vacancy depth profiles for all ion and target combinations were statis- tically fit. The results of these fits are recorded in Appendix D. Fits were done using the Python Fitter package [104] and the scipy.stats package [105]. The Scipy documentation should be consulted for the meaning of the fit parameters. The implant depth profiles for the diamond targets are plotted in Figures 6.1 to 6.4. 74 Figure 6.1: Implant depth profiles of Figure 6.2: Implant depth profiles of H+ 4 He+ into Diamond. into Diamond. Figure 6.3: Implant depth profiles of Figure 6.4: Implant depth profiles of 14 N+ ions into Diamond. 82 Pb2+ ions into Diamond. As can be observed in figures 6.1 to 6.4, there are no significant differences between the results obtained using GEANT4 with or without the GECO crystal structure classes; the average implant depths vary by less than 1% between these two models. In literature, channeling is observed when the implant beam is aligned with particular axes of the crystal target, allowing implant ions to travel along rows of planes of target atoms. This can result in underestimation of ion range if channeling effects are not accounted for in simulation [99]. No effort was made in this work to specifically observe crystal channeling effects. However, 75 all particle beams were given an angle of incidence ϕ = 0.0◦ , which should result in maximum implant depths due to channeling effects [106][107][108]. Figure 6.5 shows the percent difference in implant depth of GEANT4 and GEANT4/GECO compared to SRIM. This is calculated using the mean result of each framework as percent GEANT − SRIM difference, ϵ% = 100 . SRIM Figure 6.5: Percent difference in average implant depth compared to SRIM. In the case of the three lightest ion implant simulations, SRIM systematically shows slightly deeper average penetration than GEANT4 or GEANT4/GECO, though the extent of that difference varies greatly by case. It is only in the case of lead implants that GEANT4 and GEANT4/GECO estimate deeper penetration, and this is only the case for diamond and germanium target materials. In the case lead implants into a silicon target, SRIM demonstrates significantly deeper penetration. Differences in mean penetration depth are tabulated for 1000 keV implants in Table 6.2. For purposes of comparison the data for 76 implants into Si and Ge are also presented. Vacancy depth profile histograms are plotted for each of the diamond simulations in Fig- ures 6.6 to 6.9. Figure 6.10 displays the percent difference in vacancy depth of GEANT4and GEANT4/GECO compared to SRIM. Table 6.2: 1000 keV implant mean depth. Percent difference from SRIM are shown for GEANT4(/GECO). Negative percent difference indicates shallower mean penetration. 1000 keV mean SRIM GEANT4 % Difference GEANT4/ % Difference implant depth of from SRIM GECO from SRIM 4 He+ into Diamond (µm) 1.72 1.61 -6.4 1.61 -6.4 H+ into Diamond (µm) 8.04 7.99 -0.62 7.99 -0.62 14 N+ into Diamond (µm) 0.71 0.61 -14.08 0.61 -14.08 82 Pb2+ into Diamond (µm) 0.13 0.2 53.85 0.2 53.85 4 He+ into Si (µm) 3.5 3.54 1.14 3.54 1.14 H+ into Si (µm) 16.56 8.98 -45.77 8.98 -45.77 14 N+ into Si (µm) 1.45 1.38 -4.83 1.38 -4.83 82 Pb2+ into Si (µm) 0.25 0.73 192.0 0.72 188.0 4 He+ into Ge (µm) 3.06 3.09 0.98 3.09 0.98 H+ into Ge (µm) 11.55 11.41 -1.21 11.41 -1.21 14 N+ into Ge (µm) 1.15 1.14 -0.87 1.14 -0.87 82 Pb2+ into Ge (µm) 0.16 0.40 150.00 0.40 150.00 Vacancy depth profiles generally show the same trends as seen in the implant profiles, 77 however GEANT4 and GEANT4/GECO consistently estimate significant surface damage that is not seen in the SRIM vacancy results. The vacancy profiles produced by GEANT4 and GEANT4/GECO for the three lightest implants show peak vacancy concentrations in the first 100Åof sample depth. Figure 6.6: Vacancy depth profiles cre- Figure 6.7: Vacancy depth profiles cre- ated by 4 He+ in Diamond. ated by H+ in Diamond. Figure 6.8: Vacancy depth profiles cre- Figure 6.9: Vacancy depth profiles cre- ated by 14 N+ in Diamond. ated by 82 Pb2+ in Diamond. 78 Figure 6.10: Percent Difference in average vacancy depth compared to SRIM. Energy deposition, presented below in Figure 6.11, shows good agreement across all three frameworks in the H+ and 4 He+ implant cases, though with a slightly higher Bragg peak appearing in GEANT4(/GECO) in the 4 He+ implant case. The 14 N+ energy depo- sition curves start to show significant deviations between GEANT4(/GECO) and SRIM. GEANT4(/GECO) show more stopping power of 14 N+ implants than SRIM near the sur- face of the material, manifesting in the increased concentration of implant ions near the surface and the surface damage seen in the vacancy profiles. In the case of H+ implants, this high concentration of short ranged ions collecting near the surface of the target material is observed experimentally in the work of Kaminsky et al. [109]. 79 Figure 6.11: Comparison of SRIM, and GEANT4(/GECO) energy deposition profiles for 1000 keV implants of (a) 4 He+ , (b) H+ , (c) 14 N+ , and (d) 82 Pb2+ . 82 Pb2+ generally shows the largest difference between GEANT4(/GECO) and SRIM. In the case of implants, 82 Pb2+ straggling decreases with increasing energy in the GEANT4 simulations, both amorphous and GECO (Figure 6.4. This is in contrast to SRIM which predicts the expected behavior of straggling increasing with implant energy. The GEANT4 and GEANT4/GECO 82 Pb+2 induced vacancy profiles show righ-skewness at all energy lev- els. The surface damage increases with increasing energy. The SRIM produced 82 Pb+2 vacancy profiles show right-skewness below 1000 keV; The 1000 keV profile shows left- skewness. Surface damage decreases with increasing energy. The SRIM vacancy profile peaks are 2× deeper than the GEANT4(/GECO) vacancy peaks. The 82 Pb2+ energy depo- sition profiles show nearly twice the stopping power in SRIM, however the the GEANT4 and 80 GEANT4/GECO curves show a more gradual tail, explaining the significant right-skewness observed in the Figure 6.9 and the greater straggling observed in Figure 6.4. These differences are caused by differences in implementation of nuclear stopping power calculations. SRIM relies on the derivation of parameterized “universal interatomic poten- tial” and corresponding “universal nuclear stopping power” formulae building on the work of Lindhard, Moliere, and others[98]: Universal nuclear stopping power, 8.462×10−15 Z1 Z2 M1 Sn (ϵ) Sn (E0 ) = [98] (M1 +M2 )(Z10.23 +Z20.23 ) Where ϵ, the reduced energy, is given by: 32.53M2 E0 ϵ= [98] Z1 Z2 (M1 +M2 )(Z10.23 +Z20.23 ) and the reduced stopping power, Sn (ϵ) is piecewise defined [98]:  ln(1+1.383ϵ)  ϵ ≤ 30   2(ϵ+0.01321ϵ0.21226 +0.19593ϵ0.5 )  Sn (ϵ) =  ln(ϵ)    2ϵ ϵ > 30 GEANT4(/GECO) uses the Bethe-Bloch formula with a variety of corrections, particu- larly in the low energy regime (E < 2 MeV) where the accuracy of the Bethe-Bloch formula 81 degrades. This is further parameterized by data provided by the ICRU’49 report [110] or NIST data. In the current work, NIST data are used for the low energy parameterization invoked through the G4NIST Materials Manager [101]. The parameterized Bethe-Bloch for- mula is given as: 2mc2 β 2 γ 2 Tup     2  Tup  2Ce dE = 2πre2 mc2 nel z 2 ln − β2 1+ −δ− + S + F [101] dx β I2 Tmax Z The S and F terms are higher order corrections, including spin, shell correction, density correction, Mott correction, and finite size correction. Of note in this analysis is Z, the atomic number of the target material. The higher order corrections result in a additive term of F ∼ Z 3/2 . These competing nuclear stopping models place different weights on the charge and atomic mass of the implant ion and target material manifesting in larger differences as atomic weight increases, explaining the trend of growing differences in stopping power observed in the current study in 14 N+ and 82 Pb2+ . Finally, the most striking difference observed in this work is the vacancy profile created by 1000 keV proton implants in the GEANT4(/GECO) cases. Although the carbon vacancy in diamond (GR1 center) can be observed in PL, the correlation between PL intensity and vacancy concentration requires careful calibration and can be masked by the PL of other defects [7]. The author of this dissertation and his collaborators believe the drastic differ- ences seen in the 1000 keV H+ vacancy profiles as given by GEANT4 and GEANT4/GECO compared to SRIM is erroneous. Additional testing revealed that this profile develops only in proton implants above 500 keV. It is assessed that the differences in vacancy profiles are due to differing knock-on cascade models. When a host material atom is given sufficient energy due to a scattering event, 82 it is removed from its lattice site and any energy greater than the lattice binding energy becomes kinetic energy. This energetic host atom is then able to scatter with additional host atoms, each of which in turn gains kinetic energy and a knock-on cascade is formed. When simulating 1000 keV H+ implants, SRIM typically calculates on the order of ∼300 vacancies created for each implant ion. The majority of these vacancies are created through knock-on atom cascades. In the corresponding GEANT4 and GEANT4/GECO simulations, this value is an order of magnitude less, typically only ∼30 vacancies per implant ion, and the vast majority of vacancies (∼98%) being created by the implant ions directly. However, as seen in Figure 6.11, the stopping power curves for H+ ions are in good agreement across all three frameworks. Thus the energy deposition in the GEANT4(/GECO) simulations must be going into other physical processes than kinetic energy of host material atoms. Attempts to isolate the cause of this included swapping the hadronic physics to use the Binary Cascade model in place of the Bertini Cascade model used in the FTFP BERT HP physics list. The GEANT4 literature suggests the Binary Cascade model performs better at low energies than the Bertini Cascade model [101], however no significant change in vacancy profile was observed. Scintillation flourescense, internal conversion, nuclear deexcitation, δ-ray production, and particle-induced xray emissions were each toggled on and off with no significant change in the GEANT4(/GECO) vacancy profiles. Previous work by Vagena, et al. (2020) [96] finds that as implant energies decrease, the relative spread in results from common Monte Carlo simulation frameworks increases [96]. Their work was limited to energies above 1 MeV/u, and the current work finds that this trend continues into the sub-MeV range and deviates greatly from experimental work, though experimental comparisons could only be made in a handful of cases. Comparing the simulation results of this work with a survey of reported experimental 83 data found in literature one finds that all of the Monte Carlo simulation frameworks used in this work systematically underestimate particle range in the cases available for comparison. Comparisons to experimental data are presented in Table 6.3. In order to make comparisons across a greater range of experimental data, the results of the simulations have been inter- polated for intermediate values. Table 6.3: Comparison of Monte Carlo simulations with published experimental data. Neg- ative percent difference means an underestimation of experimentally observed ion range. Case Average Implant Depth Method GEANT4 % GEANT4/GECO SRIM % – Experimental (µm) Difference % Difference Difference 200keV H 1.10 [111] SIMS -82.16 -82.16 -81.05 into Diamond* 120keV H 0.81 [109] SIMS -30.64 -30.61 -35.24 into Diamond* 150keV N 0.20 [111] SIMS -21.79 -22.01 -17.94 into Diamond 200keV N 0.22 [111] SIMS -10.78 -10.82 -5.27 into Diamond* 640keV N 0.60 [112] EPR -24.13 -24.11 -8.45 into Diamond* 100keV H 0.80 [113] SIMS -39.94 -39.85 -44.63 into Silicon *Intermediate value interpolated from simulation results GEANT4, GEANT4 with GECO crystal object classes, and SRIM are all Monte Carlo 84 simulation frameworks used for the study of interaction of ions with, and passage through, matter widely applied in high energy physics, nuclear physics, medical and space physics, and materials science. These frameworks have shown to lose accuracy at lower energies and this has prompted the development of additional physics lists in GEANT4 [114][115]. However, the single scattering models necessary for analyzing lattice damage and vacancies appear to be less well suited to the sub-MeV regime. Both SRIM and GEANT4(/GECO) results differ greatly from the limited experimental data available, all underestimating the range of the implanted particles. It should be noted that a much more extensive work by Ziegler, et al. (2010) [116] comparing thousands of available experimental measurements has found SRIM implant depth profiles to be typically within 5% of experimentally reported values, though this claim is made simulating target materials of Ag, Al, and H2 O. The deviations seen in the current work may be attributable to ion channeling in the diamond cubic not being faithfully modeled in silico. In agreement with the results presented here, Ziegler, et al. (2010) [116] also find that there is significantly more spread in the experimental results for heavier ions and that SRIMs accuracy decreases at energies below 100 keV/u. Ziegler, et al. attribute this to the common use of Inverted Doppler Shift Attenuation (IDSA) measurements in experimental work. Bet- ter agreement is found with the few available experiments using different methods. IDSAF functions, fundamentally, by letting an implant ion undergo nuclear decay in the target mate- rial and observing emitted gammas. This method relies heavily on knowledge of the nuclear lifetime and kinematics of the nuclear decay and recoiling particle(s) [110]. It is interest- ing that [116] attributes the largest divergence between SRIM simulation and experimental values to nuclear effects, giving credence to the explanation above imputing the difference 85 observed in the simulated stopping power of lead to the calculation of nuclear effects. GEANT4(/GECO) and SRIM calculate scattering in similar ways, though with manifest differences in nucleus-nucleus interactions. However, the building block structure of physics in GEANT4 allows the end user to tailor the simulation to a specific situation, allowing, in principle, results more representative of the real-world situation under consideration. This also requires the end user to have thorough knowledge not only of the physics applicable to the situation but also the inner workings of GEANT4 to ensure the simulation is designed with fidelity to the real world. Despite this additional power, GEANT4(/GECO) does not appear to offer additional accuracy over SRIM in the sub-MeV regime in the applications herein investigated. Based on the results of the comparative study described in this chapter, SRIM will be used in the remainder of this work, though development on the GEANT4/GECO simulations will continue to hopefully improve future results. 86 Chapter 7. Modeling Color Center Formation in Diamond When modeling the formation of color centers there have historically been two approaches: the micro-scale, analyzing the interaction of individual defects with one another and neigh- boring lattice atoms using DFT, and the macro-scale, treating the diffusion of defects en masse using the framework of Fick’s 2nd Law. While both of these approaches offer physical insight into the machinery of color center formation, neither offers a complete picture. The present work develops a meso-scale approach which starts from the macro-scale and works down to near-micro scale resolution. This is accomplished through casting the color cen- ter formation problem in a reaction-diffusion formalism wherein the macroscopic diffusion of defects through the lattice is solved simultaneously with a statistically driven reaction term describing the formation of color centers in many small tributary volumes within the discretized sample domain. A meso-scale simulation and numerical solver, the “Color Center ANnealing and Reaction- Diffusion” simulation (CCANARD) has been developed to simulate the annealing process and color center formation using the reaction-diffusion formalism and is presented in Section 7.4. 7.1 Diffusion of Defects in Diamond The standard mathematical tools for the analysis of a diffusive profile are the equations of diffusion proposed by Adolf Fick in the mid-nineteenth century [117][118]. Fick’s Laws as they’ve come to be known, are partial differential equations which relate concentration, 87 diffusivity, and flux. Most generally, they take the form of: ⃗ (⃗r) J = −D(⃗r, T )∇U (7.1) ∂U ⃗  ⃗  = ∇ D(⃗r, T )∇U (⃗r) (7.2) ∂t where J is the chemical or particle flux, D is diffusivity as a function of spatial coordinate ⃗r and temperature T , and U is the chemical or particle concentration. Equation 7.1 is known as Fick’s First Law, and Equation 7.2 is known as Fick’s Second Law and will be referred to as such in the remainder of this document. Both of these equations can be simplified by the physics of the problem and through reasonable simplifying assumptions. In the case of the diffusion of defects in diamond, this work assumes isotropic diffusivity with temperature dependence only. This assumptions ignores the effects of stresses and the associated strains on diffusivity, however the existing body of literature suggests these effects are second order to temperature [119]. The work of Daw et al. (2001) [120] conducts a full treatment of diffusion of dopant atoms in the diamond cubic crystal structure using silicon media. Their work derives the diffusion equation in the presence of defect induced stresses using a tensor approach and compares with experimental data. They conclude that while the induced stresses do manifest in corrections to the diffusive model, it is highly situational, finding that when the crystallographic direction of growth is oriented along the ⟨100⟩ axis, that diffusion of defects is isotropic, and that growth along the ⟨111⟩ axis results in strong anisotropies in the diffusion tensor [120]. It is well understood that growth along the ⟨111⟩ axis is prone to the inclusion of additional defects such as crystal twinning due to reflective symmetries in the diamond cubic [121]. These defects introduce 88 anisotropic stress, so it is reasonable to find greater stress induced anisotropy in diffusion in a sample with this growth orientation. The same study also finds that the stress induced effects on diffusion are highly dependent on the defect or defects being diffused, noting higher stress dependence in the diffusion of vacancies than in the diffusion of implant structures, where the correction factor drops to less than 3%. The authors note that the anisotropies observed are most prominent at low temperatures and high external pressures. However, despite sharing the same crystal structure, diffusion mechanisms in diamond are different from those in silicon due to the ability of carbon to easily form multiple bonds (allotropes), a much greater electronic density in diamond, and a much larger bulk modulus in diamond which reduces strain from external stress in diamond [122]. Under the guidance of the Air Force Research Laboratory (AFRL) (Wright Patterson AFB, OH) the present work is studying color center formation in diamond crystals grown along the ⟨100⟩ axis and at intermediate-to-high temperatures at negligible gauge pressure. Additionally the material properties of diamond reduce the effects of stress induced correc- tions to diffusivity compared to the more thoroughly studied silicon. For these reasons this author and his collaborators believe it is reasonable to disregard higher order stress effects on diffusion at this time. Given these assumptions, Fick’s second law reduces to: ∂U ⃗2 U (⃗r) = D(T )∇ (7.3) ∂t Further, we assume zero-flux boundaries along all edges of the physical sample. This is easily understood as the absence off-gassing of particles during the annealing process. In other words, the ions that are implanted in the diamond sample will remain in the 89 diamond sample during annealing due to the low diffusivity observed in the implant ions. This simplifies Fick’s first law to: J =0 (7.4) The much stronger dependence of diffusivity on temperature, however, cannot be dis- regarded. Diffusivity is strongly correlated to temperature through the Arrhenius equation [94][123]: −E D = D0 e kB T (7.5) where D is the diffusivity coefficient, D0 is the maximal diffusivity coefficient (id est, D as T → ∞), E is the diffusion activation energy [123], kB is Boltzmann’s constant, and T is temperature. The diffusivity coefficient is a macroscopic property of the material being diffused. Hence we require a different value for each defect undergoing diffusion, including different values for every implant ion we wish to evaluate. As this is a macroscopic material property, ascer- taining values becomes an experimental problem relying on results of experiments reported in literature. Unfortunately, only values of diffusivity for the most studied implants and defects, namely nitrogen and vacancies, are present in the literature [124]. Thus it becomes necessary to calculate diffusivity from Equation 7.5. This requires reliable values of D0 and E. Values for D0 are widely reported in literature and vary widely from source to source. Reported values span a range of five orders of magnitude, though most reported values are near D0 ∼ 10−6 cm2 /s [94][122][123]. 90 Activation energy values are also available in literature for vacancies and nitrogen, but can be calculated for any implant ion using DFT [83][94][123][125]. This effectively removes the burden of determining reliable diffusivity coefficients from the experimental realm and transfers it to the theoretical realm. Consequently, the diffusion model of defects in diamond becomes much more widely applicable to any combination of defects, provided one has access to reliable DFT simulations. In this work the author and his collaborators use the following default values in CCA- NARD: Maximal theoretical diffusivity: D0 = 3.6 × 10−6 cm2 /s [122][126] Activation energy of implanted nitrogen: Eion = 5.6eV [94] Activation energy of vacancies: Evacancy = 2.6eV [83] Activation energy of carbon interstitials: Einterstitial = 1.839eV [125] The diffusion activation energy of color centers is assigned an arbitrarily high value of 10eV because it is believed that once formed color centers will not diffuse at annealing temperatures. These values provide results that are consistent with expectations of color center conver- sion rates. A full discussion of CCANARD’s results and planned bench marking experiments occurs in Chapter 8. 7.2 Proposed Reaction Terms Diffusion of defects in diamond can proceed either through interstitial or vacancy mediated pathways. In the case of interstitial mediated diffusion this exchange always takes place 91 between nearest neighbors. Vacancy mediated diffusion can take place with either nearest or second-nearest neighbor exchange [94][117]. This can be easily visualized in two dimensions on a square lattice as shown in Figure 7.1. An interstitial will always have a nearest neighbor on a diagonal and a vacancy will always have a nearest neighbor on an axis and a second- nearest neighbor on a diagonal: Figure 7.1: Interstitial and Vacancy mediated diffusion. a) Interstitial (red) mediated diffu- sion always exchanges with a nearest neighbor (black); b) Vacancy (white) mediated diffusion exchanging with nearest neighbor (black); c) Vacancy mediated diffusion exchanging with second-nearest neighbor For a color center to form, an implant ion and vacancy must be nearest neighbors. Given that a substitutional ion occupies a lattice site, knowing the number of available lattice sites (the atomic density of the material multiplied by the volume of interest) and number of available vacancies in the volume of interest, we can calculate the probability that a vacancy will reside in a nearest neighbor position to the substitutional ion. Assuming Nion number of substitutional ions are uniformly distributed in the volume of interest containing N lattice sites, the probability that any given lattice site contains a 92 substitutional ion is given by Nion Pion = . (7.6) N Given that a lattice site contains a substitutional ion, the probability that a vacancy occupies one if its nearest neighbor sites is given as Nvac Pvac = b , (7.7) N −1 assuming that the vacancies are also uniformly distributed among the N lattice sites. The denominator takes the value of N − 1 instead of N because the initial substitutional ion occupies a lattice site. The prefactor b gives the number of nearest neighbors. In three dimensions, b takes on the following values.   6; interior lattice sites           5; face lattice sites  b= (7.8)       4; edge lattice sites      3; corner lattice sites  The probability that any given lattice site contains a substitutional ion and that a neigh- boring site contains a vacancy is the probability of the existence of substitutional-vacancy pair. Performing this calculation for every available lattice site should, then, yield an ex- pected number of pairs in the domain of interest. This expression is given by Nion Nvac b̄ , (7.9) N − Nvac N − Nion 93 where b̄ represents the average of nearest neighbors per lattice site. However, as the maximum possible number of pairs is limited by the smaller of either Nion or Nvac the expression must be modified as follows:   Nion Nvac Nion Nvac Nion + Nvac − |Nion − Nvac | b̄ min(Nvac , Nion ) = b̄ N − Nvac N − Nion N − Nvac N − Nion 2 (7.10) Finally, to prevent double counting, b̄ is divided from the expression yielding an expression for the expected number of substitutional-vacancy pairs as a function of the number of substitutional ions, vacancies, and total number of lattice sites:   Nion Nvac Nion + Nvac − |Nion − Nvac | ⟨Psv ⟩ = (7.11) N − Nvac N − Nion 2 Two key facts allow the same derivation to also apply to interstitial-vacancy pairs. As seen in Figure 7.1, the nearest neighbor atoms to interstitials are different than for sub- stitutional ions. However, as b̄ is ultimately divided out of the expression, this does not matter. Secondly, while the derivation of Equation 7.11 relies on the occupancy of specific lattice sites, this is not strictly necessary to the derivation. Interstitials, by definition, do not reside at lattice sites but rather in between them. If we consider Nint to represent the number of inter-lattice point interstitial sites, then we can express the expected number of interstitial-vacancy pairs as follows.   Nint Nvac Nint + Nvac − |Nint − Nvac | ⟨Piv ⟩ = (7.12) N − Nvac N − Nint 2 Assuming sufficient thermodynamic energy at annealing temperatures to facilitate color 94 center formation, ⟨Psv ⟩ Ucc = , (7.13) V– where Ucc is the concentration of color centers and V– is the volume of the sample domain of interest. Thus the reaction terms describing the change in concentration of interstitials, substitu- tionals, vacancies, and color centers are as follows.   Uion Uvac Uion + Uvac − |Uion − Uvac | ∆Ucc = k1 (7.14) Na − Uvac Na − Uion 2   Uion Uvac Uion + Uvac − |Uion − Uvac | ∆Uion = −k1 (7.15) Na − Uvac Na − Uion 2   Uion Uvac Uion + Uvac − |Uion − Uvac | ∆Uvac = −k1 Na − Uvac Na − Uion 2   Uint Uvac Uint + Uvac − |Uint − Uvac | − k2 (7.16) Na − Uvac Na − Uint 2   Uint Uvac Uint + Uvac − |Uint − Uvac | ∆Uint = −k2 (7.17) Na − Uvac Na − Uint 2 Here, Na is the atomic density of the host material, and k1 and k2 represent the number of neighbors in a given number of coordination spheres. These terms can be adjusted to modify interaction range. Id est, if one wants to assume that defects can bond over ranges longer than nearest neighbor, these k terms can be adjusted. The k terms are derived from geometry assuming a cubic lattice and are simply the ratio of the number of nearest neighbors in a given coordination sphere to the number of nearest neighbors in the first coordination 95 sphere. (2n + 1)3 − 1 k1 = (7.18) 26 (2n)3 k2 = (7.19) 8 With a reasonable assumption that the concentration of defects, Ux << Na , then Equa- tions 7.14 - 7.17 simplify to   Uion Uvac Uion + Uvac − |Uion − Uvac | ∆Ucc = k1 (7.20) Na Na 2   Uion Uvac Uion + Uvac − |Uion − Uvac | ∆Uion = −k1 (7.21) Na Na 2   Uion Uvac Uion + Uvac − |Uion − Uvac | ∆Uvac = −k1 Na Na 2   Uint Uvac Uint + Uvac − |Uint − Uvac | − k2 (7.22) Na Na 2   Uint Uvac Uint + Uvac − |Uint − Uvac | ∆Uint = −k2 (7.23) Na Na 2 96 Having derived the diffusive terms in Section 7.1 and having derived the reaction terms in Section 7.2, the complete reaction-diffusion formalism can be expressed as a system of coupled non-linear PDEs:          Uion  Dion   Uion   ∆Uion                  ∂ Uvac      Dvac  ⃗2 Uvac  ∆Uvac       = ∇  + (7.24) ∂t            Ucc   Dcc   Ucc   ∆Ucc                  Uint Dint Uint ∆Uint 7.3 Scheme for Numerical Solution For notational simplicity Equation 7.24 is rewritten as ∂ ⃗2 U + kS U = D∇ (7.25) ∂t This system will be solved with a one-step in time Crank-Nicolson (C-N) scheme. C-N was chosen for its unconditional numerical stability [127], its second order accuracy in space and time, O(∆t2 + h2 ), (when used with second order spatial derivative approximations), and good performance in reaction-diffusion problems [128][129]. C-N has been shown to outperform other commonly used methods in solving stiff systems such as those found in diffusive problems [130]. The reaction-diffusion C-N scheme with 2nd order central difference approximation is given in [128] as un+1 − un D h n+1 i kh i = Au + Aun + S(un+1 ) + S(un ) (7.26) ∆t 2 2 97 where A is the three-dimensional discrete Laplacian operator matrix and superscripts denote time steps. Every solution un+1 is some incremental change from un , given by un+1 = un + ∆u. (7.27) Substituting Equation 7.27 into the linear part of Equation 7.26 yields ∆u + un − un D kh i = [A(∆u + un ) + Aun ] + S(un+1 ) + S(un ) , (7.28) ∆t 2 2 and with algebraic rearrangement, D∆t n k∆t h n+1 n i ∆u = [A∆u + 2Au ] + S(u ) + S(u ) . (7.29) 2 2 The implicit, nonlinear S(un+1 ) term must now be linearized through Taylor expansion and substitution of the Gateaux derivative, S(u + λ∆u) − S(u) d DS(u; ∆u) = lim = S(u + λ∆u)|λ=0 . (7.30) λ→0 λ dλ The leading terms of the Taylor expansion of S(un+1 ) are ∂S(un ) S(un+1 ) ≈ S(un ) + ∆u. (7.31) ∂un 98 Equation 7.30 is equivalent to the second term of Equation 7.31. Making this substitution yields S(un+1 ) ≈ S(un ) + DS(un ; ∆u) (7.32) Substituting Equation 7.32 into 7.29 gives D∆t k∆t ∆u = [A∆u + 2Aun ] + [DS(un ; ∆u) + 2S(un )] . (7.33) 2 2 Factoring a factor of ∆u out of Equation 7.30 allows it to be rewritten as DS(u; ∆u) = DŜ(u)∆u. (7.34) Substituting Equation 7.34 into Equation 7.33 and combining like terms yields D∆t k∆t ∆u − A∆u − DŜ(u)∆u = D∆tAun + k∆tS(un ). (7.35) 2 2 Factoring out ∆u gives the form   D∆t k∆t I− A− DŜ(u ) ∆u = D∆tAun + k∆tS(un ). n (7.36) 2 2 Equation 7.36 is a linear system which can be solved for ∆u. Then the solution is updated: un+1 = un + ∆u. (7.37) This process is repeated to step through the time domain by time increment ∆t until a 99 a final solution is reached. 7.4 Color Center Annealing and Reaction-Diffusion Solver The Color Center ANnealing And Reaction Diffusion (CCANARD) solver is a program which has been developed to efficiently solve Equation 7.24 using the numerical scheme derived in Section 7.3. This section will discuss the properties and design of CCANARD and analyze its com- putational performance. Instructions in the use of CCANARD are included in Appendix E. Commentary on the accuracy of CCANARD and experimental benchmarking will be discussed in Chapter 8. 7.4.1 Specifications CCANARD has been developed under the guidance of the Air Force Research Lab, Materials & Manufacturing Directorate. The guidance provided by AFRL required that CCANARD be able to simulate the formation of color centers in diamond during the annealing process in the full three-dimensions to a resolution on the order of nanometers using an implant profile provided by either GEANT4 or SRIM as the initial conditions. The output must show the color center concentration at any point in the sample. These requirements mean that CCANARD must have a focus on computational efficiency, as the system sizes necessary to meet the provided requirements are in the hundreds of millions. For example, a 5µm × 5µm × 5µm sample discretized to 10nm grid spacing 100 would result in 125,751,501 finite difference nodes. Since each of the four equations shown in Equation 7.24 must be solved at each node, the overall system size is slightly over half a billion! CCANARD employs a custom sparse-matrix implementation to minimize memory use, efficient linear solvers, and extensive use of parallelization, both CPU and GPU, to take maximum advantage of concurrency and reduce wall time. Additionally, in an effort to make CCANARD as general and widely applicable as pos- sible, many user-adjustable settings are present, including all physical parameters to allow CCANARD to be used to model any desired annealing situation. 7.4.2 Flowchart and Theory of Operation In keeping with the tenets of good software design [131], CCANARD follows a modular architecture with five program segments consisting of seven modules of modern Fortran written in compliance with the 2008 standard [132][133], and a series of Python and C functions. Figure 7.2 shows the CCANARD program flowchart. The main program module, encapsulated in file main.f08 controls program flow and calls all other routines. Upon starting execution the main program module sets all the default variable values for the program through a subroutine call to set defaults, a member subroutine of the setup module. Program control is then handed over to the user interface. Once user input is complete and recorded, the main program resumes control and through a number of subroutine calls reads input, allocates and initializes variables, runs the solver, records results, and elegantly ends program execution. The user interface in an interactive command-line interface (CLI). User options are orga- nized into groups accessible through user input, specifically by the user entering an integer 101 number corresponding to the menu option desired. The file sys keyin.c contains several void functions which are used to interface with the terminal to allow for reading single key strokes. This means the user does not have to press [Enter] to tell the program to read the stdin buffer, resulting in a more seamless, intuitive, and faster user interface experience. The user is required to view and confirm or change all program settings. This is designed to prevent simulations form being run with erroneous input or settings. Only once the all settings have been confirmed does the user have the option to begin the simulation. All user input is type checked, and the program will gracefully handle mismatched data types by alerting the user to the error and requesting it be fixed. Settings are not modified until the correct data type has been received. The user can also access the program details in an “About” menu. Notably, this includes the compiler used to compile the running version of CCANARD. This feature requires that the Fortran code be compiled with the use of the C Preprocessor. By automatically, at compile time, detecting the compiler used, certain program features can be locked or made accessible in line with the capabilities of each compiler. Figure 7.3 displays the “About” window showing that the running instance of CCANARD was compiled using the Nvidia Fortran compiler, nvfortran, version 22. 102 Figure 7.2: CCANARD program flow chart Once user input is complete, the main program makes a system call to execute the Python script which reads the SRIM data and translates it for parsing into CCANARD. The main program then makes several subroutine calls to the setup module (set up.f08) which then completes setting the initial conditions. The translated SRIM files containing the 103 (x, y, z) positions of each ion or vacancy are read into Fortran arrays. The arrays are then parsed and each particle is binned to a finite difference node in a sample domain established from the user’s input. Each node is assigned a tributary volume so that the entire sample domain is covered, but no two tributary volumes overlap. The tributary volume for any given node extends half the distance to each neighboring node. Any particle (or vacancy) whose coordinates fall into a node’s tributary volume will be represented in the value at that node. Id est, each node represents the concentration of defects in its tributary volume. SRIM assumes a point source particle gun. This results in an implant (and vacancy) profile generated from an infinitely thin particle beam. Once the SRIM data has been binned a convolution is applied to extrapolate the effects of beam width on the implant ion (vacancy) profile. The beam width can be set by the user to adapt the simulation to a desired experimental parameter. Finally, the initial implant and vacancy concentrations are scaled according to either the user’s desired implant fluence or the user’s desired maximum concentration. Whichever value is not specified is calculated and reported to the user to aid in designing future experiments. Maximum concentration is a function of fluence and implant energy. If the distribution of implant ions is assumed to be Gaussian or near-Gaussian (see Chapter 6) the maximum concentration is equivalent to the maximum height of the distribution. The height of the normal distribution is a function of standard deviation given as: 1 Pmax = √ . (7.38) σ 2π In ion implant physics, the second moment about the mean implant depth is referred to as straggle and is heavily energy dependent. For a given fluence, higher energy implants 104 have more straggling meaning a wider distribution (larger standard deviation) resulting in a lower maximum ion concentration compared to those implanted at a lower energy. (See the results presented in Chapter 6.) For a given implant energy, concentration and fluence can be related [98] by M C=F , (7.39) (1L)N where C is the concentration, F is the fluence, M is the number of ions in volume L3 of host materials, N is the number of implant ions, and L is an arbitrary unit of length. The factor of (1L) is necessary for dimensional compliance, but will be omitted for sim- plicity in the remainder of this work. In order to make this formula more general an energy dependent “straggle factor” term g(E) is incorporated: M C=F g(E). (7.40) N Multiplying both the numerator and the denominator by volume−1 puts the numerator in terms of atomic density, ρ, and the denominator in terms of ion concentration. M V–1 ρ C=F g(E) = F g(E) (7.41) N V–1 C Finally, p C2 C= F ρg(E) ⇐⇒ F = . (7.42) ρg(E) To avoid divergence as σ → 0, a reference standard deviation is set. In CCANARD 105 σref = 0.03249372µm, obtained from 1 keV implants into diamond. This value is the same for both 100 keV proton implants and 100 keV nitrogen implants, and is therefore assessed as a reasonable value representative of the lightest half of the periodic table. For any given set of implant coordinates passed into CCANARD, the standard deviation is calculated and therefor the straggle factor can be determined: 1 √ σ 2π σ g(E) = = ref (7.43) 1 σ √ σref 2π Thus, the initial conditions are set. Figure 7.3: The CCANARD “About” informational display reporting that it was compiled with nvfortran version 22. Note the ASCII ducks, the unofficial logo of CCANARD. (Canard is the French word for duck.) The main program then executes the time-step iteration, calling a variety of mathemat- ical subroutines to solve the non-linear system given in Equation 7.24 at each step. These mathematical subroutine are contained in the math funcs module (math functions.f08). This module relies on a bespoke sparse matrix implementation using Fortran Derived Types, Type Extensions, and Type Bound Procedures. Limited polymorphism is employed to gen- 106 eralize the use of the sparse matrix derived types, allowing for storage in both Compressed Sparse Row (CSR) format and Coordinate List (COO) format [134]. CSR storage of the 3-dimensional discrete Laplacian matrix scales as O(15N ) where N denotes the number of equations in the system. This outperforms COO storage, O(21N ), and greatly outperforms dense storage, O(N 2 ). Matrix storage size is a particular concern in CCANARD because the 3D discrete Lapla- cian matrix scales as n6 where n denotes the number of finite difference nodes along one axis of the domain mesh. The CSR format allows storage of a system with hundreds of millions of nodes. A system discretized into 501 nodes in each direction, such as that described in Section 7.4.1, can be run on a high-end workstation without having to resort to a cluster or super computer. Memory and computational scaling will be discussed in Section 7.4.4. Multiple steps are needed to solve the non-linear system at each time step. First, the annealing temperature must be calculated from the annealing recipe and the diffusivity constant, D, must be calculated for each species. The system is then linearized using the Gateaux derivative and the left-hand-side matrix can be built in CSR format. The right- hand-side can be formed and stored as a 1-dimensional array. The linearized system is then solved using an iterative Krylov subspace method. Based on the symmetric and definite structure of the discrete Laplacian matrix, the Conjugate Gra- dient (CG) method offers the most efficient iterative algorithm [135]. In CCANARD, Jacobi preconditioning is applied to the CG algorithm by default. A library of bespoke sparse ma- trix mathematical routines is contained in the sparse tools module (sparse tools.f08). Many of these were written from scratch, and others were adapted from the open source spmat and SPARSEKIT libraries [134][136]. The sparse matrix routines were customized to the sparse matrix Derived Types and Type Extensions developed for CCANARD with a fo- 107 cus on computational efficiency: reducing nested loops to the extent possible, ordering loop indices to take advantage of column-major storage in Fortran, and maximizing concurrency. Fortran’s pass-by-reference nature eliminates duplication of memory objects in subroutine calls through assignment of the ALLOCATABLE attribute to arrays rather than requiring C-style pointers. Aside from preventing memory object duplication, declaring arrays ALLOCATABLE assists in memory management, places arrays on the heap rather than stack (necessary for storing large arrays), and offers better compiler optimization than C-style pointers [132]. All this results in highly efficient mathematical functions allowing larger systems to be solved within tolerable computational time. Once the system has been solved at the current time step, the main program takes control to increment the time and record results and update the user as desired. The system is then solved for the next time increment, and the process repeats until the end of the specified annealing recipe. Once the simulation is complete, the main program makes several subroutine calls to record the results, plot graphs, and animate graphs as requested. Plotting is handled natively in Fortran using the ogpf library, a fully object-oriented Fortran wrapper for GNUPlot [137]. ogpf is written for a Microsoft Windows environment, so several small changes had to be made, namely the handling of escape characters in plot titles. An additional function was written to allow dynamically changing the GNUPlot terminal so plots could be both saved as images and displayed on screen. Animations are generated from a series of plots saved as .png images using a system call to ffmpeg using high quality mpeg4 encoding. Eight user-selectable color maps are available for producing output visualization. These have been chosen to provide diverging, perceptually uniform color maps with linearly in- creasing luminosity for the highest perceptual resolution. These maps have been chosen to 108 provide the end user aesthetically pleasing color schemes that maximize defect distribution fidelity on screen and on paper[138]. Examples of the default “Bent Cool Warm” color map can be seen in the CCANARD simulation results presented in Section 8.2. 7.4.3 Parallelization Parallelization in CCANARD is achieved through three methods: single instruction multiple data (SIMD) vectorization, multithreading using the OpenMP Application Programming Interface (API), and compiler driven parallelization of native language constructs, namely the Fortran DO CONCURRENT construct in the case of CCANARD. Automatic vectorization realizes data-level parallelization through compiler analysis and optimization of program characteristics. Most modern processors support SIMD instruction set extensions and most modern compilers support some level of automatic vectorization optimization [139]. CCANARD assumes the presence of these capabilities and uses auto- vectorization to great extent for the mass updating of arrays. Additional auto-vectorization is left to the whims of the compiler and optimization level specified at compile time[140]. OpenMP is a parallel computing API developed continuously since the early 1990’s de- signed to bring parallel execution to serial C/C++ and Fortran applications through the use of compiler directives in a shared memory multiprocessing environment. OpenMP uses a multithreading fork-join model allowing the program to fork (split) into regions of parallel execution and join (merge) back into regions of serial execution [139]. OpenMP worksharing-loop constructs are used frequently in CCANARD, both in the form of simple parallel execution and in reduction operations. OpenMP allows for nested multithreading, and this is used effectively in the sparse matrix functions where it is necessary to iterate a nested loop over the entries of a CSR matrix inside a loop over the the beginning 109 and ending indices of each row. Setting the initial conditions is one of the most computationally intensive parts of CCA- NARD. This is due to nested loops iterating over all particles and all finite difference nodes. The computational burden is lessened through short-circuit analysis and exiting the inner- most loop when a particle is placed in a bin, but the process remains the longest running single process in CCANARD. Once binned, the data must be convolved in three dimensions as described in Section 7.4.2, adding to the wall time required to set the initial condi- tions. Fortunately, this region is embarrassingly parallel in the sense that each of these processes must occur for each of implant ions, vacancies, and interstitials. As this requires slightly different instructions to be run for each species, it is handled through the use of the OMP SECTIONS construct allowing for multiple different instruction to be run on multiple data simultaneously. Limited GPU acceleration is available in CCANARD, though presently limited to modern (2008 and later) Nvidia GPU hardware and the Nvidia HPC SDK suite of compilers of version 20.11 or greater. This GPU acceleration is provided as a compiler optimization to standard Fortran 2008 or later [141]. The Fortran DO CONCURRENT statement is used to assert to the compiler that the loop iterations have no interdependence [132]. This allows compilers to perform automatic parallelization. If instructed to do so, the Nvidia Fortran compiler, nvfortran, can offload DO CONCURRENT code to the GPU. Host-device data transfer limitations still apply, so this optimization must be used carefully, and in a way to maximize computation on data once it has been transferred to the GPU. At the time of writing only the setting of the initial conditions, the binning of particles and convolution is offloaded for GPU acceleration. Given that the complete array of particle coordinates and the finite difference nodes can be transferred to the GPU and the entire binning and convolution process can be 110 computed, this is a process ideally suited to GPU offload. Parallelization is controllable by the end user through the CCANARD CLI. User input is checked against system capabilities before acceptance. Exempli gratia, the user is not able to request more processors than the system has available. By default GPU acceleration is OFF, and the user is only able to turn it ON if the option is available, id est, if CCANARD was compiled with nvfortran. 7.4.4 Profiling CCANARD has been profiled for computational efficiency and multiprocessing scaling. Through- out the development process, CCANARD was profiled using gprof, the GNU profiler. This was used to find bottlenecks and focus development efforts. The results of these efforts are presented below. All of the profiling discussed in this subsection was conducted on a com- puter using an AMD “Ryzen 9” 3950x 16-core (32-thread) overclocked to 4250 MHz on all cores, 128GB of DDR4 RAM, an Nvidia GeForce RTX2060 Super graphics card (2176 stream processors/272 Tensor Cores), running a highly customized and stripped down version of the Calculate Linux operating system [142]. 7.4.4.1 CPU Scaling The multiprocessing efficiency of CCANARD was analyzed using two test cases each for the CG and PCCG algorithms. The first test case simulated the annealing of a diamond after implantation of 250 keV nitrogen ions at a fluence of 1015 N/cm2 . The results are presented in Figure 7.4. The annealing recipe used was a 10 hour anneal, ramping from 0◦ C to 800◦ C in one hour and holding at 800◦ C for 9 hours. The system was discretized into a 51 × 51 × 51 node domain. The iterative solver convergence tolerance was set to 10−6 . 111 Up to 32 processors were used. The hardware used in testing only has 16 physical pro- cessor cores, so the test using 32 processors takes full advantage of AMD’s “simultaneous multi-threading” technology. In all test cases, the 32 processor trial obtained the same or worse results than the 16 processor case. This is attributed to the fact that this requires two threads to be simultaneously running on each processor core. This may be advanta- geous when each thread is only demanding a portion of the processor’s capacity. However, the profiling results presented herein suggest that CCANARD is efficiently using nearly all available processing power in each physical processor core leaving little available to run a second thread. Figure 7.4: Multiprocessing performance of CCANARD simulating the annealing of 250 keV nitrogen implants in diamond. The domain is discretized into 51 × 51 × 51 finite difference nodes. Lower is better. CCANARD was compiled with: gfortran -fopenmp -O3 -cpp -s -march=native The second case profiled for multiprocessor scaling was 500 keV protons implanted into diamond at a fluence of 1015 H/cm2 . The same annealing recipe was used. The results are shown below in Figure 7.5. CCANARD’s implementation of PCCG scales somewhat better than the CG implemen- 112 tation, however in the test cases herein presented, CG outperforms its preconditioned coun- terpart. This phenomenon will be explored in the next subsection. Table 7.1 tabulates the results shown in Figures 7.4 and 7.5. Table 7.1: CPU scaling results Implant Algorithm CPU scaling 250 keV nitrogen CG x−0.661 250 keV nitrogen PCCG x−0.7 500keV Proton CG x−0.571 500keV Proton PCCG x−0.376 Figure 7.5: Multiprocessing performance of CCANARD simulating the annealing of 500 keV Proton implants in diamond. The domain is discretized into 51 × 51 × 51 finite difference nodes. Lower is better. CCANARD was compiled with: gfortran -fopenmp -O3 -cpp -s -march=native 113 7.4.4.2 Performance of Jacobi Preconditioner In the test case of 250 keV nitrogen, the Jacobi preconditioner appears to slow down conver- gence, resulting in slower computation time. There is additional overhead in implementation of PCCG, such as extracting and inverting the diagonal of the LHS matrix to form the pre- conditioner. Each iteration of the preconditioned solver requires one additional matrix-vector multiply. However, the performance difference cannot be attributed to this small amount of computational overhead. On average, PCCG requires 10 iterations to solve the test system in each time step, in contrast to the average 8 iterations required by the un-preconditioned CG implementation. Clearly, in this test case, Jacobi preconditioning does not sufficiently reduce the condition number of the LHS matrix. Perhaps a more complex preconditioner would yield better results, however with only 8 iterations required by CG, performance is acceptably good without any preconditioning. The second test case of 500 keV Proton implants is more interesting. As in the first case, PCCG shows better scaling than CG and worse performance for a number of processors less than 16. In the 16 processor case PCCG slightly outperforms CG. This is not seen at 32 processors, though this is attributed to running multiple threads on each core as described in the previous subsection. In the process of running CCANARD for testing and experimental planning it appears that the performance of the Jacobi preconditioner is very situationally dependent. CCA- NARD is designed to be highly versatile, readily applied to any number of problems involving diffusion and/or annealing of defects in crystal. The condition number of these possible sys- tems could vary greatly, depending upon the distribution of initial defects, the discretization and domain sizes chosen by the user, the interaction ranges selected for color center forma- tion, the annealing recipe, the activation energies chosen for defect mobility, and many other 114 user-definable physical properties. The inclusion of the Jacobi preconditioner in CANNARD was to provide additional versatility of the tool. In general, from the author’s qualitative analysis during testing, it seems PCCG compares better to CG in terms of performance as system size increases and convergence tolerance is decreased, though this is also highly situational and does not manifest in the test case herein presented. Scaling of the iterative solvers with system size is presented below in Figure 7.6 for the 500 keV proton implant case described above. This particular testing indicates that the performance of PCCG and CG is nearly equivalent for small and intermediate sized systems. Only as system size increases does a performance difference emerge with CG outperforming PCCG. Figure 7.6: Comparison of computational speed of CG and Jacobi Preconditioned CG with increasing system size. Lower is better. CCANARD was compiled with: gfortran -fopenmp -O3 -cpp -s -march=native 7.4.4.3 GPU Acceleration At the time of writing, GPU acceleration in CCANARD is limited to the embarrassingly parallel problem of binning the implant ions and performing the beam-extent convolution. The concurrent nature of this problem shows better CPU scaling than the time step itera- 115 tion, approaching the theoretical linear reduction in wall time as processor count increases. Despite the excellent CPU parallelization, this problem still benefits from GPU offload. This benefit is presented in Figure 7.7. Figure 7.7: GPU acceleration performance in setting initial conditions. Lower is better. CCANARD was compiled with: nvfortran -mp -Minfo=accel -stdpar=gpu -fast -O3 -cpp The GPU performance scaling reduces with system size, ultimately paralleling CPU scaling. The process of binning and convolving the initial coordinates is performed in an OMP SECTIONS work sharing construct. This means that in the current version of CCANARD only one CPU core will access the GPU at a time, effectively making the process CPU bound, though still taking advantage GPU speed-up as shown above. 116 Chapter 8. Experimental Benchmarking of CCANARD Before it can be used in scientific endeavors, CCANARD must be tested for accuracy. This is accomplished in two parts. The first is comparing CCANARD results to experimental results found in literature, and the second is through a dedicated experiment designed to test CCANARD. 8.1 Comparison of CCANARD with Literature The first exercise of testing the accuracy of CCANARD was studying the effect of implant en- ergy with color center conversion efficiency. The work of Lühmann et al. (2018) [7] and Pezza- gna et al. (2010) [6] indicate that NV color center conversion efficiency, or yield, is dependent total color centers formed on implant ion energy, where conversion efficiency is defined, η = . total implantations Using an annealing scheme of 2 hours at 800◦ C, the conversion efficiency of a variety of nitrogen implant energies are reported, ranging from less than 1% at 5 keV up to nearly 50% for 18 MeV implants [6]. It is assessed the higher conversion efficiemcy at higher energies is due to the increased yield of vacancies at higher energy. Pezzagna et al. (2010) [6] point out that the curve of vacancies created per ion (according to SRIM simulation) is nearly coincident with the curve of color center conversion efficiency as a function of implant energy. This is shown in Figure 8.1, taken from Pezzagna et al. (2010) [6]. The work of Lühmann et al. (2018) [7] reproduced these results. The same experiment was repeated in CCANARD. Results are presented in Table 8.1 and Figure 8.2. CCANARD matches the literature values perfectly at 100 keV, but quickly diverges from the results 117 presented in [7] and [6]. At an implant energy of 500 keV, CCANARD underestimates the yield by a factor of nearly 3, and at 1000 keV, underestimates the yield by a factor of more than 5. These simulations were constrained to a domain of spatial extent 1µm × 1µm × 1µm, discretized into a 51 × 51 × 51 node grid resulting in 200 Å resolution. The interaction range between nitrogen ions and vacancies was set to 17 nm. A fluence of 1014 N/cm2 and beam spot size diameter of 10 nm were used. Other physical parameters were left at the default values specified in Section 7.1. Despite the differences in color center conversion efficiency, the average number of va- cancies produced per implant ion in the CCANARD simulations agrees well with the values shown in [6]. This is reasonable as SRIM is used in both, and this bears witness that the discrepancy observed in the results is not attributable to the simulation input. 118 Figure 8.1: This image, taken from the work of Pezzagna et al. (2010) [6] shows clearly the dependence of color center yield on ion implant dependence in (a) and shows the same relationship in the number of vacancies created as a function of implant ion energy in (b). Table 8.1: Tabulated values of color center conversion efficiency at given implant energies as measured in Lühmann et al. (2018) [7] and simulated in CCANARD. Vacancies produced per implant ion as determined by SRIM and used in CCANARD are presented and in agreement with Pezzagna et al. (2010) [6], shown in Figure 8.1 keV % Conversion [7] % Conversion CCANARD Vacancies/ion 100 10.45 10.11 186 500 11.5 4 311 1000 11.75 1.6 341 119 Though the reaction term (see Section 7.2) is designed to be energy agnostic relying only on the local concentration (itself an implicitly energy dependent quantity as shown in the Chapter 6 results) of defects, these results, especially in light of agreement in the number of vacancies created, suggest the necessity of incorporating an explicitly energy-dependent term or otherwise modifying the reaction term. Figure 8.2: Comparison of the energy dependence of color center conversion efficiency mod- eled in CCANARD with experimental results reported in Lühmann et al. (2018) [7]. The experimental work of Acosta et al. (2009) [8] shows that the maximum concentration of color centers saturates at an annealing temperature around 875◦ C. Annealing at tempera- tures beyond this does not result in a significant increase in color center concentration. This was tested in CCANARD using 150 keV nitrogen implants at eight fluences between 1012 and 5 × 1015 N/cm2 annealed at three temperatures: 700◦ C, 875◦ C, and 1050◦ C. All other program settings were as previously described. The results are shown below in Figure 8.3 where a plateau clearly manifests after 875◦ C for all implant fluences, in excellent agreement with Acosta et al. (2009) [8]. 120 Figure 8.3: Color center concentration as a function of annealing temperature at given fluences as simulated by CCANARD. Concentration levels off after 875◦ C, and in some cases even decreases, matching the experimental work of Acosta et al. [8]. The work of Pezzagna et al. (2010) [6] also determined a relation between color center conversion efficiency and implant fluence. Up to a certain fluence color center yield increases, but after an inflection point drops off. The conversion efficiencies were recorded for all the above simulations and generally agreed with Pezzagna et al. (2010) [6], increasing to a maximimum conversion rate and rapidly dropping off with increasing fluence. This was hard to see in the cases of the 875◦ C and 1050◦ C anneals, however both final points (fluence = 5×1015 ) plot on the order of 1% less than the penultimate points at fluence = 1015 . These data suggest that the fluence at which peak color center formation occurs increases with annealing temperature. To make a direct comparison with Pezzagna et al.’s (2010) [6] work, a series of CCANARD simulations were run using 5keV nitrogen implants. The results are shown in Figure 8.4. The corresponding experimental results from Pezzagna et al. (2010) [6] are presented for context in Figure 8.5. 121 Figure 8.4: CCANARD simulation of the annealing of 5 keV nitrogen implants for 2 hours at 800◦ C. Conversion efficiency increases with fluence until peaking and then rapidly dropping off. The general shape of Figure 8.4 matches that of Figure 8.5, however the conversion efficiency differs by up to an order of magnitude, in this case with CCANARD overestimat- ing conversion efficiency. And while both plots show a peak followed by a rapid drop-off, CCANARD estimates the peak at a fluence of 3.5 × 1015 , compared to the experimentally determined peak at a fluence of ∼ 2 × 1013 . 122 Figure 8.5: Experimentally determined relationship between color center conversion effi- ciency and fluence using 5 keV nitrogen implants annealed at 800◦ C for 2 hours. The dotted curve is a guide for the eye. Image taken from Pezzagna et al. (2010) [6]. The similarities observed between the simulated and experimental conversion efficiency dependence on fluence suggests that the formula used in CCANARD to convert between fluence and concentration (see Section 7.4.2) is accurate. However the difference in overall conversion efficiency is further evidence that CCANARD does not fully capture the implant energy dependence observed in [6] and [7], which would suggest conversion efficiency of ∼ 1% from 5 keV nitrogen implants. 8.2 Benchmarking Experiment While the comparison of CCANARD with literature has been elucidating into the strengths and weaknesses of its use as a prediction engine, a dedicated experiment, or series thereof, 123 would reveal much more about about CCANARD’s capabilities, and help gather data for implementing empirically-driven corrections to CCANARD’s reaction-diffusion framework. To these ends, an experiment is proposed, and at the time of writing, being planned in conjunction with AFRL. The experiment will consist of three high-purity diamond sam- ples implanted using the nanobeam facility at the Center for Integrated Nanotechnologies (CINT), part of Los Alamos National Laboratory. Implants shall be 150 keV nitrogen at a fluence of 1014 N/cm2 . Multiple sites will be implanted on each sample for good experimen- tal statistics. Each sample will be annealed at a different temperature: 800◦ C, 1200◦ C, and 1550◦ C to test CCANARD across the entire range of typical annealing temperatures. Each sample will be subjected to a 2 hour ramp in temperature, taking the sample from room temperature up to its full annealing temperature. The samples will be measured using PL and confocal microscopy at four hours of total annealing time and every two hours thereafter. Color center concentration and gradient will be recorded. To the extent possible, concen- tration and gradient will be recorded for nitrogen ions and vacancies [7][143]. If possible, SIMS measurements will be conducted on the samples following all other measurements to ascertain an accurate profile of the nitrogen implants. These experiments have been simulated in CCANARD using a domain of spatial extent .5µm × .5µm × .5µm, discretized into a 51 × 51 × 51 node grid resulting in 100 Å resolution. The interaction range between nitrogen ions and vacancies was set to 17 nm. A beam spot size diameter of 10 nm was used. Other physical parameters were left at the default values specified in Section 7.1. The simulated results are tabulated in Table 8.2. 124 Table 8.2: Simulated results of planned experiment. Time (hr) Temp (K) Max. Ion con. Max. Vac con. Max. Int con. Max. CC con. Conversion % initial – 4.500E+18 3.500E+21 3.500E+21 0 – 4 1073 4.000E+18 3.500E+21 3.000E+21 4.000E+17 6 1073 4.000E+18 3.000E+21 1.200E+18 5.000E+17 8 1073 4.000E+18 3.000E+21 8.000E+17 6.000E+17 10 1073 4.000E+18 3.000E+21 6.000E+17 7.000E+17 8.9 4 1473 4.000E+18 7.000E+18 3.481E+17 5.000E+17 6 1473 4.000E+18 3.000E+18 3.475E+17 6.000E+17 8 1473 4.000E+18 1.800E+18 3.469E+17 6.000E+17 10 1473 4.000E+18 1.200E+18 3.469E+17 6.000E+17 7.18 4 1823 4.000E+18 3.450E+17 3.559E+17 3.500E+17 6 1823 4.000E+18 3.435E+17 3.553E+17 3.500E+17 8 1823 4.000E+18 3.428E+17 3.546E+17 3.500E+17 10 1823 4.000E+18 3.422E+17 3.540E+17 3.500E+17 4.14 Based on the simulation, the experiment should reveal rapid saturation of color center conversion at the two higher temperatures. Longer annealing times will not yield significantly more color centers. In the 800◦ C (1073 K) anneal sample, color center concentration will increase throughout the entire annealing period, surpassing the color center concentrations in the higher temperature anneals. The simulated final color center distributions are presented below in Figures 8.6 - 8.8. 125 Figure 8.6: Final (10 hr) conditions along the sample center line after the 800◦ C (1073 K) anneal. The beam propagates in the x-direction and impinges at the origin. The final color center distribution spans approximately 0.15µm, and is not Gaussian, likely due to a high concentration of vacancies impeding color center formation in this location [9]. 126 Figure 8.7: Final (10 hr) conditions along the sample center line after the 1200◦ C (1473 K) anneal. The beam propagates in the x-direction and impinges at the origin. The final color center distribution spans approximately 0.12µm with a tail extending approximately an additional 0.5µm. The distribution is much more Gaussian than in Figure 8.6, likely due to increased mobility in the vacancies. Note that the distribution of interstitials is nearly uniform. 127 Figure 8.8: Final (10 hr) conditions along the sample center line after the 1550◦ C (1823 K) anneal. The beam propagates in the x-direction and impinges at the origin. The final color center distribution spans approximately 0.12µm with a tail extending approximately an additional 0.5µm. The distribution is much more Gaussian than in Figure 8.6, likely due to increased mobility in the vacancies. Note that the distribution of vacancies and interstitials is nearly uniform. 128 Chapter 9. Conclusion and Continuing Work 9.1 Conclusion The first part of this dissertation documents the development of a solidstate neutron detector based on semiconducting diamond using a rectifying Schottky barrier and reports the subse- quent manufacture and testing. The Positionally Opposed Schottky Semi-Metal (POSSM) neutron detector uses a sandwhich architecture consisting of two diamonds surrounding a 10µm thick gadolinium foil. The foil forms a rectifying contact with each diamond, becoming a shared cathode for the two diodes. Gadolinium was chosen for its large neutron cross sec- tion which enhances the detector’s efficiency. Preamplifier circuitry was designed following the typical photodiode transimpediance amplifier architecture. The diamond diodes were constructed using electronic grade HPHT diamonds purchased from New Diamond Technol- ogy. The semiconducting diamond layers were grown using the CVD method by personnel of the Engineering Research Center, Microwave and Plasma Processing Laboratory, located aboard Michigan State University. Sample cutting, preparation, cleaing, and imaging were performed by staff engineers at Fraunhofer USA - Center for Coatings and Diamond Tech- nology, co-located with the Engineering Research Center. The detector was assembled at Hillsdale College and tested at the Weapons Neutron Research Facility located at the Los Alamos Neutron Science Center. Testing demonstrated the overall success of the POSSM detector with peak detector efficiency achieved with neu- trons between 50 and 100MeV. Neutron detection at these enrgies surpasses the performance of existing diamond based detectors. However, during experimental setup, the preamplifier electronics were damaged rending one channel inoperable. Consequently, only one one chan- 129 nel of data was recorded and analyzed. The initial testing demonstrated that the preamplifier electronics were too slow, and a set of revised specifications was developed. A new design, dubbed Design Revision 1, was produced, though at the time of writing has not been con- structed. The revised design includes a more durable method of diamond diode construction, faster reacout electonics, over-voltage protection, and improved input/output. The second portion of this dissertation reports the development of a computer simu- lation of the formation of quantum color centers in diamond. The process is cast into a reaction-diffusion formalism. Stochastic reaction terms are derived resulting in a system of coupled nonlinear partial differential equations. A computer program was written to solve this system. The program, the “Color Center ANnealing And Reaction-Diffusion” (CCA- NARD) program, linearizes the system through the substitution of the Gateaux derivative and solves the linearized system using a finitie difference Crank-Nicolson scheme. To aide in the development of CCANARD, a comparative study of commonly available Monte Carlo simulation frameworks was undertaken. Monte Carlo simulation is used to simulate the process of ion implantation into diamond (or other) substrates. The simulated implants provide the initial conditions for the CCANARD program. The comparative study simulated H+ , 4 He+ , 14 N+ , and 82 Pb2+ ions implanted into the Group IV semiconductors at energies ranging from 100 keV to 1000 keV. The frameworks used were SRIM, GEANT4, and GEANT4 with the GEANT Crystal Object (GECO) classes allowing for the simulation of crystal structure in materials. After analysing the results and comparing the simulation results with literature, SRIM was used in the remainder of this work. The results reported in this dissertation demonstrate proof of concept and are promising for the future development of diamond based sensing devices. Nevertheless, this work rep- resents only a small step forward and continuing work will continue to refine the POSSM 130 technology and CCANARD program. 9.2 Continued Development of the POSSM Detector The next step in the continued development and proving of the POSSM prototype is to construct the entire device according to the Design Rev. 1 plans and specs. This includes rebuilding the diamond diodes and fabricating the preamplifier electronics as specified on 90x70 mm perfboard. Once the Design Rev. 1 prototype has been assembled it should be thoroughly tested. Given more time than available for the testing of the Design Rev. 0 prototype, a full battery of local tests should be conducted with charged-particle sources for response characterization of the device. The device will need to be tested in neutron beams, ideally at WNR or equivalent facilities where energies in excess of several hundred MeV can be achieved. High flux testing at lower energies, while not only more readily available through a variety of user facilities, will offer additional insight into the potential of the device to compete against extant technologies for use in beam line monitoring, reactor monitoring, and other similar environments, and should be undertaken. Diamond detectors have been demonstrated to perform well in these environments [2][38][40][41], and the POSSM device, though not explicitly designed for this purpose may serve admirably in the role adding capability such as self-vetoing to the application. It is likely that minor changes in the form of a Design Rev. 2 will be necessary following the completion of Design Rev. 1 testing. Once the requisite modifications have been made, Design Rev. 2 should be extended to include full SolidWorks or equivalent models of the POSSM diode package, and the preamplifier design package should be converted for PCB 131 fabrication. Additional considerations should be given to size – Surface Mount Technology will allow for a much smaller preamplifier footprint, and the size can be standardized to commercially available jiffy boxes. Interest has been expressed by researchers at the Army Research Laboratory in the POSSM device. In addition to its intended purpose, there is interest in the technology for charge collection and conversion, non-destructive testing of the integrity of mission-critical materiel, and for battlefield sensing applications. Whatever form these future collaborations take, and whatever revisions, modifications, and improvements will be made to the POSSM technology, this work has clearly demon- strated proof of concept and laid the foundation for the next generation of solid state neutron detection. 9.3 Continued Development of the CCANARD pro- gram Once experimental benchmarking has been completed, attention must be turned to the im- proving the accuracy of CCANARD. This could be as simple as adding a “calibration” factor to the reaction terms allowing CCANARD to be empirically calibrated to experimental data. This has the advantage of being easily implemented and accurate. This method is versatile and extensible allowing CCANARD to be recalibrated for a variety of possible applications. The drawback to this method is that it relies on extant experimental data which may or may not accurately represent other implant/host material combinations. Alternatively, a more theoretical approach could be taken in efforts to derive better reaction terms. This may result in a better description of the meso-scale color center formation process, leading 132 to increased accuracy across the periodic table. Additional future work on CCANARD will be driven by the needs of AFRL. Possible short-term work includes porting the code to Windows for wider availability on government systems. Such a port should be a relatively simple procedure requiring a change to system commands primarily in the user interface. However, the inclusion of the Windows Subsystem for Linux on later releases of Windows may negate the necessity of such a port. Again, this work will be conducted at the discretion of AFRL. Concurrency in the CCANARD program should be further exploited to reduce wall time. This could include extending arrayed data into Fortran coarrays [131][132], increasing GPU offload, and making GPU acceleration more general through the use of the OpenMP host- device model to make GPU acceleration both brand and compiler agnostic [139]. Maximum performance could be obtained through writing critical sections of GPU-accelerated code in a dedicated language. With proper memory management, the Krylov solver is expected be able to be run highly efficiently and massively parallel on GPU hardware, possibly reducing run time by two orders of magnitude on mid-tier consumer hardware. Future versions of CCANARD could include automatic selection of the iterative solver method based on anal- ysis of the eigenvalues of the linear system once created. An automatic selection tool and more advanced preconditioning schemes could further improve CCANARD’s computational efficiency. Longer-term development of CCANARD should focus on increasing the fidelity of the physics in the simulation, incorporating higher order effects on diffusion and color center conversion, such as pressure dependence (see Chapter 7), diffusion in the presence of electric and magnetic fields, planar crystallographic defects, and the simulation of more complex color centers. 133 For maximum end-user convenience, simulations should be able to be saved and loaded, a graphical user interface should be developed, and a batch processing capability (simple automation from an input file rather than requiring the user to set up each simulation) should be developed. This could ease the use of CCANARD on shared cluster or super- computer resources. Monte Carlo implant should be incorporated directly into CCANARD. Implementing an internal Monte-Carlo implant code could take the form of simply calling an existing code from the CCANARD user interface, incorporating the code of an open source MC framework directly into the CCANARD code base, or developing MC code from scratch. Tangentially useful to CCANARD would be a library of ion activation energies. This could be developed through a thorough literature search and a number of DFT simulations. While such a library could not be exhaustive including every ion-solid combination possible, it could include the majority of ions of interest into to common semi-conducting solids. This would, in many cases, reduce the burden of research on the end user. 134 BIBLIOGRAPHY [1] D. C. Palmer and S. E. Palmer, “CrystalMaker X,” 2019. [2] S. Almaviva, M. Marinelli, E. Milani, G. Prestopino, A. Tucciarone, C. Verona, G. Verona-Rinati, M. Angelone, and M. Pillon, “Improved performance in synthetic diamond neutron detectors: Application to boron neutron capture therapy,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 612, pp. 580–582, Jan. 2010. [3] M. Chadwick, M. Herman, P. Obloinsk, M. Dunn, Y. Danon, A. Kahler, D. Smith, B. Pritychenko, G. Arbanas, R. Arcilla, R. Brewer, D. Brown, R. Capote, A. Carlson, Y. Cho, H. Derrien, K. Guber, G. Hale, S. Hoblit, S. Holloway, T. Johnson, T. Kawano, B. Kiedrowski, H. Kim, S. Kunieda, N. Larson, L. Leal, J. Lestone, R. Little, E. Mc- Cutchan, R. MacFarlane, M. MacInnes, C. Mattoon, R. McKnight, S. Mughabghab, G. Nobre, G. Palmiotti, A. Palumbo, M. Pigni, V. Pronyaev, R. Sayer, A. Sonzogni, N. Summers, P. Talou, I. Thompson, A. Trkov, R. Vogt, S. van der Marck, A. Wallner, M. White, D. Wiarda, and P. Young, “ENDF/B-VII.1 Nuclear Data for Science and Technology: Cross Sections, Covariances, Fission Product Yields and Decay Data,” Nuclear Data Sheets, vol. 112, pp. 2887–2996, Dec. 2011. [4] J. Shive, The Properties, Physics, and Design of Semiconductor Devices. Princeton, NJ: D. Van Nostrand Company, Inc., 1959. [5] R. Kuate Defo, E. Kaxiras, and S. L. Richardson, “How carbon vacancies can affect the properties of group IV color centers in diamond: A study of thermodynamics and kinetics,” Journal of Applied Physics, vol. 126, p. 195103, Nov. 2019. [6] S. Pezzagna, B. Naydenov, F. Jelezko, J. Wrachtrup, and J. Meijer, “Creation efficiency of nitrogen-vacancy centres in diamond,” New J. Phys., vol. 12, p. 065017, June 2010. [7] T. Lhmann, N. Raatz, R. John, M. Lesik, J. Rdiger, M. Portail, D. Wildanger, F. Kleiler, K. Nordlund, A. Zaitsev, J.-F. Roch, A. Tallaire, J. Meijer, and S. Pez- zagna, “Screening and engineering of colour centres in diamond,” J. Phys. D: Appl. Phys., vol. 51, p. 483002, Dec. 2018. [8] V. M. Acosta, E. Bauch, M. P. Ledbetter, C. Santori, K.-M. C. Fu, P. E. Barclay, R. G. Beausoleil, H. Linget, J. F. Roch, F. Treussart, S. Chemerisov, W. Gawlik, and D. Budker, “Diamonds with a high density of nitrogen-vacancy centers for magnetom- etry applications,” Phys. Rev. B, vol. 80, p. 115202, Sept. 2009. [9] F. Waldermann, P. Olivero, J. Nunn, K. Surmacz, Z. Wang, D. Jaksch, R. Taylor, I. Walmsley, M. Draganski, P. Reichart, A. Greentree, D. Jamieson, and S. Prawer, “Creating diamond color centers for quantum optical applications,” Diamond and Re- lated Materials, vol. 16, pp. 1887–1895, Nov. 2007. 135 [10] “Diamond History and Lore,” 2023. [11] J. Robertson, “Electronic and atomic structure of diamond-like carbon,” Semiconduc- tor Science and Technology, vol. 18, p. 9, 2003. [12] K. Seeger, Semiconductor Physics, An Introduction. Heidelber, Fed. Rep. of Germany: Springer-Verlag, 4th ed., 1989. [13] L. Berry and B. Mason, Mineralogy: Concepts, Descritpions, Determinations. San Francisco, CA: W. H. Freeman and Company, 1959. [14] C. Wild and W. Eckhard, “cvd diamond booklet.pdf.” [15] G. Schmid, J. Koch, R. Lerche, and M. Moran, “A neutron sensor based on single crystal CVD diamond,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 527, pp. 554– 561, July 2004. [16] S. Bhagavantam and D. Narayana Rao, “Dielectric Constant of Diamond,” Nature, vol. 161, no. 729, 1948. [17] C. M. Breeding and J. E. Shigley, “The ”Type” Classification System of Diamonds and Its Importance in Gemology,” Gems & Gemology, vol. 45, pp. 96–111, July 2009. [18] S. Eaton-Magana, J. E. Shigley, and C. M. Breeding, “Observations on HPHT-Grown Synthetic Diamonds: A Review,” G&G, vol. 53, pp. 262–284, Nov. 2017. [19] I. Carmichael, F. Turner, and J. Verhoogen, Igneous Petrology. McGraw-Hill interna- tional series in the earth and planetary sciences, New York, NY: McGraw-Hill, 1974. [20] P. J. Wyllie, The Dynamic Earth: Textbook in Geosciences. New York, NY: John Wiley & Sons, 1971. [21] A. Butcher, “A Brief History of Lab-Grown Diamonds,” 2022. [22] U. F. D’Haenens-Johansson, K. Soe Moe, P. Johnson, S. Y. Wong, R. Lu, and W. Wang, “Near-Colorless HPHT Synthetic Diamonds from AOTC Group,” G&G, vol. 50, pp. 30–45, May 2014. [23] S. S. Nicley, THE BORON DOPING OF SINGLE CRYSTAL DIAMOND FOR HIGH POWER DIODE APPLICATIONS. PhD Thesis, Michigan State University, East Lansing, MI, 2015. [24] Y. Gu, J. Lu, T. Grotjohn, T. Schuelke, and J. Asmussen, “Microwave plasma reactor design for high pressure and high power density diamond synthesis,” Diamond and Related Materials, vol. 24, pp. 210–214, Apr. 2012. 136 [25] G. Dhanaraj, K. Byrappa, V. Prasad, and M. Dudley, Springer Handbook of Crystal Growth. Heidelberg: Spinger-Verlag, 2010. [26] H. Sternschulte, M. Schreck, B. Stritzker, A. Bergmaier, and G. Dollinger, “Lithium addition during CVD diamond growth: influence on the optical emission of the plasma and properties of the films,” Diamond and Related Materials, vol. 9, pp. 1046–1050, Apr. 2000. [27] M. Ovezmyradov, I. V. Magedov, L. V. Frolova, G. Chandler, J. Garcia, D. Bethke, E. A. Shaner, and N. G. Kalugin, “Chemical Vapor Deposition of Phosphorous- and Boron-Doped Graphene Using Phenyl-Containing Molecules,” [28] Y. Wu, J. Tong, L. Ruan, F. Luo, G. Liu, R. Zhang, X. Han, Y. Zhang, F. Tian, and X. Zhang, “N-type diamond semiconductor induced by co-doping selenium and boron,” Computational Materials Science, vol. 196, p. 110515, Aug. 2021. [29] C. J. Wort and R. S. Balmer, “Diamond as an electronic material,” Materials Today, vol. 11, pp. 22–28, Jan. 2008. [30] K. M. Gupta and N. Gupta, Advanced Semiconducting Materials and Devices. Engi- neering Materials, Switzerland: Springer International Publishing, 2016. [31] “Ge - Germanium: Electrical Properties,” 2021. [32] K. Zeljami, J. Gutierrez, J. P. Pascual, T. Fernandez, A. Tazon, and M. Boussouis, “CHARACTERIZATION AND MODELING OF SCHOTTKY DIODES UP TO 110 GHZ FOR USE IN BOTH FLIP-CHIP AND WIRE-BONDED ASSEMBLED ENVI- RONMENTS,” PIER, vol. 131, pp. 457–475, 2012. [33] Z. Horvth, V. Rakovics, B. Szentpli, S. Pspki, and K. dnsk, “InP Schottky junctions for zero bias detector diodes,” Vacuum, vol. 71, pp. 113–116, May 2003. [34] M. Angelone, D. Lattanzi, M. Pillon, M. Marinelli, E. Milani, A. Tucciarone, G. Verona-Rinati, S. Popovichev, R. Montereali, M. Vincenti, and A. Murari, “De- velopment of single crystal diamond neutron detectors and test at JET tokamak,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spec- trometers, Detectors and Associated Equipment, vol. 595, pp. 616–622, Oct. 2008. [35] M. Marinelli, E. Milani, G. Prestopino, M. Scoccia, A. Tucciarone, G. Verona-Rinati, M. Angelone, M. Pillon, and D. Lattanzi, “High performance Li6F-diamond thermal neutron detectors,” Appl. Phys. Lett., vol. 89, p. 4, 2006. [36] D. Wooldridge, A. Ahearn, and J. Burton, “Conductivity Pulses Induced in Diamond by Alpha-Particles,” Physical Review, vol. 71, no. 12, 1947. 137 [37] R. F. Kozlov, E. A. Konorova, M. I. Krapivin, V. A. Nadein, and V. G. Yudina, “USING DIAMOND DETECTORS AS IMMERSED a COUNTERS,” Atomic Energy, vol. 40, no. 6, pp. 574–575, 1976. [38] D. Jain, J. Nuwad, N. Manoj, and V. Sudarsan, “Diamond-Based Radiation Detectors for Applications in Highly Corrosive Solutions and High-Radiation Fields,” in Materials Under Extreme Conditions, pp. 683–715, Elsevier, 2017. [39] P. Bergonzo, A. Brambilla, D. Tromson, C. Mer, B. Guizard, F. Foulon, and V. Amosov, “CVD diamond for radiation detection devices,” Diamond and Related Materials, vol. 10, pp. 631–638, Mar. 2001. [40] T. Shimaoka, J. H. Kaneko, K. Ochiai, M. Tsubota, H. Shimmyo, A. Chayahara, H. Umezawa, H. Watanabe, S.-i. Shikata, M. Isobe, and M. Osakabe, “A diamond 14 MeV neutron energy spectrometer with high energy resolution,” Review of Scientific Instruments, vol. 87, p. 023503, Feb. 2016. [41] M. Majerle, “The response of single crystal diamond detectors to 17–34 MeV neutrons,” Nuclear Inst. and Methods in Physics Research, A, vol. 951, p. 6, 2020. [42] V. Zerkin, “Experimental Nuclear Reaction Data (EXFOR) Database.” [43] S. Agostinelli, J. Allison, K. Amako, J. Apostolakis, H. Araujo, P. Arce, M. Asai, D. Axen, S. Banerjee, G. Barrand, F. Behner, L. Bellagamba, J. Boudreau, L. Broglia, A. Brunengo, H. Burkhardt, S. Chauvie, J. Chuma, R. Chytracek, G. Cooperman, G. Cosmo, P. Degtyarenko, A. Dell’Acqua, G. Depaola, D. Dietrich, R. Enami, A. Fe- liciello, C. Ferguson, H. Fesefeldt, G. Folger, F. Foppiano, A. Forti, S. Garelli, S. Giani, R. Giannitrapani, D. Gibin, J. Gmez Cadenas, I. Gonzlez, G. Gracia Abril, G. Greeni- aus, W. Greiner, V. Grichine, A. Grossheim, S. Guatelli, P. Gumplinger, R. Hamatsu, K. Hashimoto, H. Hasui, A. Heikkinen, A. Howard, V. Ivanchenko, A. Johnson, F. Jones, J. Kallenbach, N. Kanaya, M. Kawabata, Y. Kawabata, M. Kawaguti, S. Kelner, P. Kent, A. Kimura, T. Kodama, R. Kokoulin, M. Kossov, H. Kurashige, E. Lamanna, T. Lampn, V. Lara, V. Lefebure, F. Lei, M. Liendl, W. Lockman, F. Longo, S. Magni, M. Maire, E. Medernach, K. Minamimoto, P. Mora de Freitas, Y. Morita, K. Murakami, M. Nagamatu, R. Nartallo, P. Nieminen, T. Nishimura, K. Ohtsubo, M. Okamura, S. O’Neale, Y. Oohata, K. Paech, J. Perl, A. Pfeiffer, M. Pia, F. Ranjard, A. Rybin, S. Sadilov, E. Di Salvo, G. Santin, T. Sasaki, N. Sav- vas, Y. Sawada, S. Scherer, S. Sei, V. Sirotenko, D. Smith, N. Starkov, H. Stoecker, J. Sulkimo, M. Takahata, S. Tanaka, E. Tcherniaev, E. Safai Tehrani, M. Tropeano, P. Truscott, H. Uno, L. Urban, P. Urban, M. Verderi, A. Walkden, W. Wander, H. We- ber, J. Wellisch, T. Wenaus, D. Williams, D. Wright, T. Yamada, H. Yoshida, and D. Zschiesche, “Geant4a simulation toolkit,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 506, pp. 250–303, July 2003. 138 [44] E. Bagli, M. Asai, A. Dotti, L. Pandola, and M. Verderi, “Allowing for crystalline structure effects in Geant4,” Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, vol. 402, pp. 304–307, July 2017. [45] B. Shu, “Q-Value Calculator (QCalc).” [46] E. V. Ryabeva, V. V. Kadilin, and G. L. Dedenko, “Neutron Detector Based on Polystyrene and Cadmium Layers,” J. Phys.: Conf. Ser., vol. 798, p. 012178, Jan. 2017. [47] G. Festa, F. Grazzi, A. Pietropaolo, A. Scherillo, and E. Schooneveld, “Thermal neu- tron radiative capture on cadmium as a counting technique at the INES beam line at ISIS: A preliminary investigation of detector cross-talk,” Applied Radiation and Isotopes, vol. 130, pp. 264–269, Dec. 2017. [48] “Mathematica.” [49] K. S. Krane, Introductory Nuclear Physics. New York, NY: John Wiley & Sons, 1988. [50] A. B. Smith, “Fast-neutrons incident on gadolinium,” Annals of Nuclear Energy, vol. 31, pp. 1813–1831, Nov. 2004. [51] J. C. Khong, D. Daisenberger, G. Burca, W. Kockelmann, A. S. Tremsin, and J. Mi, “Design and Characterisation of Metallic Glassy Alloys of High Neutron Shielding Capability,” Sci Rep, vol. 6, p. 36998, Nov. 2016. [52] “ASTAR.” [53] L. Ry, L. Dobrzaski, F. Dubecky, S. Jaboski, P. Parys, W. Sysz, and M. Rosiski, “De- velopment of x-ray and ion diagnostics of plasma obtained with a 10-TW femtosecond laser,” Phys. Scr., vol. 91, p. 074008, July 2016. [54] B. Zako, L. Hrubn, A. agtov, J. Osvald, P. Bohek, E. Kovov, Y. Halahovets, S. V. Ro- zov, and V. Sandukovskij, “Study of Schottky barrier detectors based on a high qual- ity 4H-SiC epitaxial layer with different thickness,” Applied Surface Science, vol. 536, p. 147801, Jan. 2021. [55] J. Isberg, J. Hammersberg, E. Johansson, T. Wikstrm, D. J. Twitchen, A. J. White- head, S. E. Coe, and G. A. Scarsbrook, “High Carrier Mobility in Single-Crystal Plasma-Deposited Diamond,” Science, vol. 297, pp. 1670–1672, Sept. 2002. [56] M. Schreck, P. ajev, M. Trger, M. Mayr, T. Grnwald, M. Fischer, and S. Gsell, “Charge carrier trapping by dislocations in single crystal diamond,” Journal of Applied Physics, vol. 127, p. 125102, Mar. 2020. 139 [57] W. F. Paxton, T. Wade, M. Howell, N. Tolk, W. P. Kang, and J. L. Davidson, “Thermionic emission characterization of boron-doped microcrystalline diamond films at elevated temperatures: Thermionic emission characterization of B-doped microcrys- talline diamond films,” Phys. Status Solidi A, vol. 209, pp. 1993–1995, Oct. 2012. [58] L. Diederich, O. Kttel, P. Aebi, and L. Schlapbach, “Electron affinity and work func- tion of differently oriented and doped diamond surfaces determined by photoelectron spectroscopy,” Surface Science, vol. 418, pp. 219–239, Nov. 1998. [59] S. M. Baker, G. R. Rossman, and J. D. Baldeschwieler, “Observation of surface charge screening and Fermi level pinning on a synthetic, borondoped diamond,” Journal of Applied Physics, vol. 74, pp. 4015–4019, Sept. 1993. [60] N. Rouger and A. Marchal, “Design of Diamond Power Devices: Application to Schot- tky Barrier Diodes,” Energies, vol. 12, p. 2387, June 2019. [61] R. Ramey and S. Katzberg, “Surface Work Function of Gadolinium,” Journal of Chem- ical Physics, vol. 53, no. 4, pp. 1347–1348, 1970. [62] N. Connor, “Gadolinium Electrical Resistivity and Electrical Conductivity,” 2020. [63] B. H. Gaury and P. M. Haney, “Sesame Documentation,” p. 69, 2019. [64] B. H. Gaury and P. M. Haney, “Sesame,” 2019. [65] D. Zhu, J. Xu, A. N. Noemaun, J. K. Kim, E. F. Schubert, M. H. Crawford, and D. D. Koleske, “The origin of the high diode-ideality factors in GaInN/GaN multiple quantum well light-emitting diodes,” Appl. Phys. Lett., vol. 94, p. 081113, Feb. 2009. [66] J. Baliga and A. Huang, “Y9.FS1.1: SiC Power Devices for SST Applications,” 2017. [67] N. A. Al-Ahmadi, “Metal oxide semiconductor-based Schottky diodes: a review of recent advances,” Mater. Res. Express, vol. 7, p. 032001, Mar. 2020. [68] D. Xue, K. Betzler, and H. Hesse, “Dielectric constants of binary rare-earth com- pounds,” J. Phys.: Condens. Matter, vol. 12, pp. 3113–3118, Apr. 2000. [69] “Photodiode Circuit Design Wizard.” [70] P. Horowitz and W. Hill, The Art of Electronics. New York, NY: Cambridge University Press, 3rd ed., 2015. [71] “Preamplifier Introduction.” [72] J. Ardizzoni, “A Practical Guide to High-Speed Printed-Circuit-Board Layout,” Analog Dialogue, vol. 39, no. 09, 2005. 140 [73] A. Hardy, “Diamond Detector Update,” 2022. [74] I. A. Prokhorov, V. G. Ralchenko, A. P. Bolshakov, A. V. Polskiy, A. V. Vlasov, I. A. Subbotin, K. M. Podurets, E. M. Pashaev, and E. A. Sozontov, “Analysis of synthetic diamond single crystals by X-ray topography and double-crystal diffractometry,” Crys- tallogr. Rep., vol. 58, pp. 1010–1016, Dec. 2013. [75] D. Howell, Quantifying stress and strain in diamond. Thesis (Doctoral), University College London, London, England, 2009. [76] F. Lawrence, C. Mallika, U. Kamachi Mudali, R. Natarajan, D. Ponraju, S. Se- shadri, and T. Sampath Kumar, “Radiation degradation in the mechanical properties of Polyetheretherketonealumina composites,” Journal of Nuclear Materials, vol. 420, pp. 338–341, Jan. 2012. [77] P. Wady, A. Wasilewski, L. Brock, R. Edge, A. Baidak, C. McBride, L. Leay, A. Grif- fiths, and C. Valls, “Effect of ionising radiation on the mechanical and structural properties of 3D printed plastics,” Additive Manufacturing, vol. 31, p. 100907, Jan. 2020. [78] “DIRECTIVE 2011/65/EU OF THE EUROPEAN PARLIAMENT AND OF THE COUNCIL of 8 June 2011 on the restriction of the use of certain hazardous substances in electrical and electronic equipment,” 2011. [79] “B8 Open-Design Diamond Detector,” 2023. [80] R. Brun and F. Rademakers, “ROOT An object oriented data analysis framework,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spec- trometers, Detectors and Associated Equipment, vol. 389, pp. 81–86, Apr. 1997. [81] “Los Alamos Neutron Science Center - Facilities.” [82] I. Harris, C. J. Ciccarino, J. Flick, D. R. Englund, and P. Narang, “Group III Quantum Defects in Diamond are Stable Spin-1 Color Centers,” Phys. Rev. B, vol. 102, p. 195206, Nov. 2020. arXiv:1907.12548 [cond-mat, physics:physics, physics:quant-ph]. [83] P. Dek, B. Aradi, M. Kaviani, T. Frauenheim, and A. Gali, “Formation of NV centers in diamond: A theoretical study based on calculated transitions and migration of nitrogen and vacancy related defects,” Phys. Rev. B, vol. 89, p. 075203, Feb. 2014. [84] D. Suter and F. Jelezko, “Single-spin magnetic resonance in the nitrogen-vacancy cen- ter of diamond,” Progress in Nuclear Magnetic Resonance Spectroscopy, vol. 98-99, pp. 50–62, Feb. 2017. 141 [85] G. Thiering and A. Gali, “Color centers in diamond for quantum applications,” in Semiconductors and Semimetals, vol. 103, pp. 1–36, Elsevier, 2020. [86] H. Bernien, B. Hensen, W. Pfaff, G. Koolstra, M. S. Blok, L. Robledo, T. H. Taminiau, M. Markham, D. J. Twitchen, L. Childress, and R. Hanson, “Heralded entanglement between solid-state qubits separated by three metres,” Nature, vol. 497, pp. 86–90, May 2013. [87] S. Ajisaka and Y. B. Band, “Decoherence of Nitrogen Vacancy Centers in Diamond,” Phys. Rev. B, vol. 94, p. 134107, Oct. 2016. arXiv:1511.09049 [cond-mat, physics:quant- ph]. [88] P. C. Humphreys, N. Kalb, J. P. J. Morits, R. N. Schouten, R. F. L. Vermeulen, D. J. Twitchen, M. Markham, and R. Hanson, “Deterministic delivery of remote entangle- ment on a quantum network,” Nature, vol. 558, pp. 268–273, June 2018. [89] R. U. A. Khan, P. M. Martineau, B. L. Cann, M. E. Newton, and D. J. Twitchen, “Charge transfer effects, thermo and photochromism in single crystal CVD synthetic diamond,” J. Phys.: Condens. Matter, vol. 21, p. 364214, Sept. 2009. [90] J. R. Maze, P. L. Stanwix, J. S. Hodges, S. Hong, J. M. Taylor, P. Cappellaro, L. Jiang, M. V. G. Dutt, E. Togan, A. S. Zibrov, A. Yacoby, R. L. Walsworth, and M. D. Lukin, “Nanoscale magnetic sensing with an individual electronic spin in diamond,” Nature, vol. 455, pp. 644–647, Oct. 2008. [91] F. Dolde, H. Fedder, M. W. Doherty, T. Nbauer, F. Rempp, G. Balasubramanian, T. Wolf, F. Reinhard, L. C. L. Hollenberg, F. Jelezko, and J. Wrachtrup, “Electric- field sensing using single diamond spins,” Nature Phys, vol. 7, pp. 459–463, June 2011. [92] P. Neumann, I. Jakobi, F. Dolde, C. Burk, R. Reuter, G. Waldherr, J. Honert, T. Wolf, A. Brunner, J. H. Shim, D. Suter, H. Sumiya, J. Isoya, and J. Wrachtrup, “High- Precision Nanoscale Temperature Sensing Using Single Defects in Diamond,” Nano Lett., vol. 13, pp. 2738–2742, June 2013. [93] P. Pringsheim and R. C. Voreck, “Farbzentren in Diamanten,” Z. Physik, vol. 133, pp. 2–8, Sept. 1952. [94] R. Jones, J. Goss, H. Pinto, and D. Palmer, “Diffusion of nitrogen in diamond and the formation of A-centres,” Diamond and Related Materials, vol. 53, pp. 35–39, Mar. 2015. [95] Y.-C. Chen, P. S. Salter, S. Knauer, L. Weng, A. C. Frangeskou, C. J. Stephen, S. N. Ishmael, P. R. Dolan, S. Johnson, B. L. Green, G. W. Morley, M. E. Newton, J. G. Rarity, M. J. Booth, and J. M. Smith, “Laser writing of coherent colour centres in diamond,” Nature Photon, vol. 11, pp. 77–80, Feb. 2017. 142 [96] E. Vagena, E. Androulakaki, M. Kokkoris, N. Patronis, and M. Stamati, “A compar- ative study of stopping power calculations implemented in Monte Carlo codes and compilations with experimental data,” Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, vol. 467, pp. 44–52, Mar. 2020. [97] E. Lagzdina, D. Lingis, A. Plukis, R. Plukien, M. Gasparinas, I. Matulaitien, V. Ko- valevskij, G. Niaura, and V. Remeikis, “Structural investigation of RBMK nuclear graphite modified by 12C+ ion implantation and thermal treatment,” Nuclear Instru- ments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, vol. 444, pp. 23–32, Apr. 2019. [98] J. F. Ziegler, J. P. Biersack, and M. D. Ziegler, SRIM The Stopping and Range of Ions in Matter. Chester, MD: SRIM Co., 2015. [99] S. Pezzagna and J. Meijer, “Single-Ion Implantation in Diamond with a High Lateral Resolution,” in Comprehensive Hard Materials, pp. 321–336, Elsevier, 2014. [100] E. Bagli, M. Asai, D. Brandt, A. Dotti, V. Guidi, and D. H. Wright, “A model for the interaction of high-energy particles in straight and bent crystals implemented in Geant4,” Eur. Phys. J. C, vol. 74, p. 2996, Aug. 2014. [101] “GEANT4 Physics Reference Manual,” 2020. [102] E. Bagli and V. Guidi, “DYNECHARM++: a toolkit to simulate coherent interac- tions of high-energy charged particles in complex structures,” Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, vol. 309, pp. 124–129, Aug. 2013. [103] E. Bagli, “G4 mc channeling/tree/master/data,” 2018. [104] T. Cokelaer, “fitter,” 2022. [105] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wil- son, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, . Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cim- rman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, SciPy 1.0 Contributors, A. Vijaykumar, A. P. Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C. N. Woods, C. Fulton, C. Masson, C. Hggstrm, C. Fitzgerald, D. A. Nicholson, D. R. Hagen, D. V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G. A. Price, G.-L. Ingold, G. E. Allen, G. R. Lee, H. Audren, I. Probst, J. P. Dietrich, J. Silterra, J. T. Webber, J. Slavi, J. Nothman, J. Buchner, J. Kulick, 143 J. L. Schnberger, J. V. de Miranda Cardoso, J. Reimer, J. Harrington, J. L. C. Ro- drguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. Thoma, M. Newville, M. Km- merer, M. Bolingbroke, M. Tartre, M. Pak, N. J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P. A. Brodtkorb, P. Lee, R. T. McGibbon, R. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T. J. Pingel, T. P. Robitaille, T. Spura, T. R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upad- hyay, Y. O. Halchenko, and Y. Vzquez-Baeza, “SciPy 1.0: fundamental algorithms for scientific computing in Python,” Nat Methods, vol. 17, pp. 261–272, Mar. 2020. [106] T. E. Derry, R. W. Fearick, and J. P. F. Sellschop, “Ion channeling in natural diamond. II. Critical angles,” Phys. Rev. B, vol. 26, pp. 17–25, July 1982. [107] V. Raineri, V. Privitera, G. Galvagno, F. Priolo, and E. Rimini, “Channeling implants in silicon crystals,” Materials Chemistry and Physics, vol. 38, pp. 105–130, July 1994. [108] S. T. Picraux and J. U. Andersen, “Measurements and Calculations of Critical Angles for Planar Channeling,” Phys. Rev., vol. 186, pp. 267–272, Oct. 1969. [109] F. V. Kaminsky, S. N. Shilobreeva, B. Y. Ber, and D. Y. Kazantsev, “Quantification of Hydrogen in Natural Diamond by Secondary Ion Mass Spectrometry (SIMS),” Dokl. Earth Sc., vol. 494, pp. 699–703, Sept. 2020. [110] J. Deasy, “ICRU report 49, stopping powers and ranges for protons and alph particles,” Medical Physics, vol. 21, no. 5, 1994. [111] R. Wilson, “Depth distributions and range parameters for elements implanted into single-crystal diamonds and chemically vapor-deposited polycrystal diamond films,” Surface and Coatings Technology, vol. 47, pp. 559–571, Aug. 1991. [112] R. Kalish, C. Uzan-Saguy, B. Philosoph, V. Richter, J. Lagrange, E. Gheeraert, A. Deneuville, and A. Collins, “Nitrogen doping of diamond by ion implantation,” Diamond and Related Materials, vol. 6, pp. 516–520, Mar. 1997. [113] C. W. Magee and C. P. Wu, “Hydrogen ion implantation profiles as determined by SIMS,” Nuclear Instruments and Methods, vol. 149, pp. 529–533, Feb. 1978. [114] I. Kyriakou, D. Sakata, H. N. Tran, Y. Perrot, W.-G. Shin, N. Lampe, S. Zein, M. C. Bordage, S. Guatelli, C. Villagrasa, D. Emfietzoglou, and S. Incerti, “Review of the Geant4-DNA Simulation Toolkit for Radiobiological Applications at the Cellular and DNA Level,” Cancers, vol. 14, p. 35, Dec. 2021. [115] G. Santin, V. Ivanchenko, H. Evans, P. Nieminen, and E. Daly, “GRAS: a general- purpose 3-D Modular Simulation tool for space environment effects analysis,” IEEE Trans. Nucl. Sci., vol. 52, pp. 2294–2299, Dec. 2005. 144 [116] J. F. Ziegler, M. Ziegler, and J. Biersack, “SRIM The stopping and range of ions in matter (2010),” Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, vol. 268, pp. 1818–1823, June 2010. [117] T. Tan, “Point Defects and Diffusion in Semiconductors,” in Handbook of Solid State Diffusion, Volume 1, pp. 239–319, Elsevier, 2017. [118] A. Fick, “Ueber Diffusion,” Ann. Phys. Chem., vol. 170, no. 1, pp. 59–86, 1855. [119] A. Shiryaev, D. Frost, and F. Langenhorst, “Impurity diffusion and microstructure in diamonds deformed at high pressures and temperatures,” Diamond and Related Materials, vol. 16, pp. 503–511, Mar. 2007. [120] M. S. Daw, W. Windl, N. N. Carlson, M. Laudon, and M. P. Masquelier, “Effect of stress on dopant and defect diffusion in Si: A general treatment,” Phys. Rev. B, vol. 64, p. 045205, June 2001. [121] C. B. Slawson, “TWINNING IN THE DIAMOND,” American Minerologist, vol. 35, no. 3-4, pp. 193–206, 1950. [122] J. Bernholc and A. Antonelli, “Mechanism of Self-Diffusion in Diamond,” PHYSICAL REVIEW LETTERS, vol. 61, no. 23, p. 4, 1988. [123] X. J. Hu, Y. B. Dai, R. B. Li, H. S. Shen, and X. C. He, “The diffusion of vacancies near a diamond (001) surface,” Solid State Communications, p. 4, 2002. [124] B. Zhang and X. Wu, “Calculation of self-diffusion coefficients in diamond,” Appl. Phys. Lett., vol. 100, p. 051901, Jan. 2012. [125] N. Gothard, “C interstitial activation energy in diamond?,” Nov. 2022. [126] J. Orwa, K. Ganesan, J. Newnham, C. Santori, P. Barclay, K. Fu, R. Beausoleil, I. Aharonovich, B. Fairchild, P. Olivero, A. Greentree, and S. Prawer, “An upper limit on the lateral vacancy diffusion length in diamond,” Diamond and Related Materials, vol. 24, pp. 6–10, Apr. 2012. [127] K. Blevins, A Comparison of Finite-Difference Methods for the Solution of the Tran- sient Heat Conduction Equation in Inhomogeneous Media. PhD thesis, Air Forse In- stitute of Technology, Wright-Patterson AFB, 1982. [128] J. Ramos, “A review of some numerical methods for reaction-diffusion equations,” Mathematics and Computers in Simulation, vol. 25, pp. 538–548, Dec. 1983. [129] R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equa- tions. Philadelphia, PA: SIAM, 2007. 145 [130] M. Adak, “Comparison of Explicit and Implicit Finite Difference Schemes on Diffusion Equation,” in Mathematical Modeling and Computational Tools (S. Bhattacharyya, J. Kumar, and K. Ghoshal, eds.), vol. 320, pp. 227–238, Singapore: Springer Singapore, 2020. Series Title: Springer Proceedings in Mathematics & Statistics. [131] M. Curcic, Modern Fortran: Building Efficient Parallel Applications. Shelter Island, NY: Manning Publications Co., 2020. [132] M. Metcalf, J. Reid, and M. Cohen, Modern Fortran Explained: Incorporating For- tran 2018. Numerical Mathematics and Scientific Computation, Oxford, UK: Oxford University Press, 2018. [133] “Information technology Programming languages Fortran Part 1: Base language,” Language Specification ISO/IEC 1539:2010, ISO/IEC, 2010. [134] Y. Saad, “SPARSEKIT: a basic tool kit for sparse matrix computations,” Technical Report 90.20, Research Institute for Advanced Computer Science NASA Ames Re- search Center, Columbia, MD, 1990. [135] J. W. Demmel, Applied Numerical Linear Algebra. Philadelphia, PA: SIAM, 1997. [136] S. Hajian, “spmat90,” 2017. [137] M. Rahmani, “ogpf,” 2018. [138] K. Moreland, “Diverging Color Maps for Scientic Visualization (Expanded),” in Pro- ceedings of the 5th International Symposium on Visual Computing, 2009. [139] T. G. Mattson, Y. He, and A. Koniges, The OpenMP Common Core: Making OpenMP Simple Again. Scientigic and Engineering Computation, 2019. [140] J. G. Feng, Y. P. He, and Q. M. Tao, “Evaluation of Compilers Capability of Automatic Vectorization Based on Source Code Analysis,” Scientific Programming, vol. 2021, pp. 1–15, Nov. 2021. [141] G. Ozen and G. Lopez, “Accelerating Fortran DO CONCURRENT with GPUs and the NVIDIA HPC SDK,” 2020. [142] A. Tratsevskiy, “Calculate Linux Scratch Edition,” 2022. [143] J. Martin, W. Grebner, W. Sigle, and R. Wannemacher, “Confocal microscopy of color center distributions in diamond,” Journal of Luminescence, vol. 83-84, pp. 493–497, Nov. 1999. 146 APPENDIX A. Design Drawings Rev. 0 Due to formatting requirements and publishing constraints, the mechanical drawings have been removed from this edition. They are available upon request from the author. 147 APPENDIX B. Design Drawings Rev. 1 Due to formatting requirements and publishing constraints, the mechanical drawings have been removed from this edition. They are available upon request from the author. 148 APPENDIX C. Circuit Schematics Due to formatting requirements and publishing constraints, circuit diagrams have been re- moved from this edition. They are available upon request from the author. 149 APPENDIX D. Implant and Vacancy Profiles Table D.1: Statistical fits of implant profiles Target Material Implant Simulation Energy Mean Ion Straggle Best-fit Fit Parameters (keV) Depth (µm) (µm2 ) Distribution Diamond alpha GEANT4 100 0.29 0.0 genlogistic (0.6, 0.31, 0.02) Diamond alpha G4/GECO 100 0.29 0.0 genlogistic (0.63, 0.31, 0.02) Diamond alpha SRIM 100 0.31 0.0 genlogistic (0.43, 0.33, 0.01) Diamond alpha GEANT4 150 0.39 0.0 genlogistic (0.57, 0.4, 0.02) Diamond alpha G4/GECO 150 0.38 0.0 genlogistic (0.55, 0.4, 0.02) Diamond alpha SRIM 150 0.41 0.0 genlogistic (0.44, 0.44, 0.01) Diamond alpha GEANT4 250 0.55 0.0 genlogistic (0.53, 0.57, 0.02) Diamond alpha G4/GECO 250 0.55 0.0 genlogistic (0.54, 0.57, 0.02) Diamond alpha SRIM 250 0.59 0.0 genlogistic (0.38, 0.62, 0.01) Diamond alpha GEANT4 500 0.9 0.0 genlogistic (0.5, 0.93, 0.02) Diamond alpha G4/GECO 500 0.9 0.0 genlogistic (0.53, 0.93, 0.02) Diamond alpha SRIM 500 0.97 0.0 genlogistic (0.4, 1.0, 0.01) Diamond alpha GEANT4 1000 1.61 0.0 genlogistic (0.52, 1.64, 0.02) Diamond alpha G4/GECO 1000 1.61 0.0 genlogistic (0.53, 1.64, 0.02) Diamond alpha SRIM 1000 1.72 0.0 genlogistic (0.4, 1.75, 0.02) Diamond H GEANT4 100 0.48 0.0 genlogistic (0.49, 0.5, 0.02) Diamond H G4/GECO 100 0.48 0.0 genlogistic (0.49, 0.5, 0.02) Diamond H SRIM 100 0.44 0.0 genlogistic (0.47, 0.46, 0.01) Diamond H GEANT4 150 0.68 0.0 genlogistic (0.48, 0.71, 0.02) Diamond H G4/GECO 150 0.68 0.0 genlogistic (0.5, 0.71, 0.02) Diamond H SRIM 150 0.65 0.0 genlogistic (0.46, 0.67, 0.01) Diamond H GEANT4 250 1.16 0.0 genlogistic (0.57, 1.19, 0.02) Diamond H G4/GECO 250 1.16 0.0 genlogistic (0.55, 1.19, 0.02) Diamond H SRIM 250 1.14 0.0 genlogistic (0.51, 1.16, 0.02) Diamond H GEANT4 500 2.86 0.01 genlogistic (0.69, 2.88, 0.04) Diamond H G4/GECO 500 2.86 0.01 genlogistic (0.7, 2.88, 0.04) Diamond H SRIM 500 2.85 0.01 genlogistic (0.58, 2.88, 0.03) Diamond H GEANT4 1000 7.99 0.07 genlogistic (0.64, 8.06, 0.09) Diamond H G4/GECO 1000 7.99 0.06 genlogistic (0.66, 8.06, 0.09) Diamond H SRIM 1000 8.04 0.03 genlogistic (0.77, 8.07, 0.08) 150 Table D.1 (cont’d) Diamond N GEANT4 100 0.11 0.0 gennorm (1.15, 0.11, 0.02) Diamond N G4/GECO 100 0.11 0.0 gennorm (1.17, 0.11, 0.02) Diamond N SRIM 100 0.11 0.0 genlogistic (0.39, 0.13, 0.01) Diamond N GEANT4 150 0.16 0.0 genlogistic (0.83, 0.16, 0.01) Diamond N G4/GECO 150 0.16 0.0 genlogistic (0.82, 0.16, 0.01) Diamond N SRIM 150 0.16 0.0 genlogistic (0.36, 0.18, 0.01) Diamond N GEANT4 250 0.24 0.0 genlogistic (0.62, 0.25, 0.01) Diamond N G4/GECO 250 0.24 0.0 genlogistic (0.68, 0.25, 0.02) Diamond N SRIM 250 0.25 0.0 genlogistic (0.34, 0.28, 0.01) Diamond N GEANT4 500 0.39 0.0 genlogistic (0.54, 0.41, 0.01) Diamond N G4/GECO 500 0.39 0.0 genlogistic (0.56, 0.41, 0.02) Diamond N SRIM 500 0.44 0.0 genlogistic (0.34, 0.47, 0.01) Diamond N GEANT4 1000 0.61 0.0 genlogistic (0.46, 0.64, 0.01) Diamond N G4/GECO 1000 0.62 0.0 genlogistic (0.5, 0.64, 0.01) Diamond N SRIM 1000 0.71 0.0 genlogistic (0.35, 0.74, 0.01) Diamond Pb GEANT4 100 0.07 0.0 lognorm (0.48, -0.02, 0.08) Diamond Pb G4/GECO 100 0.07 0.0 chi2 (4.98, -0.0, 0.01) Diamond Pb SRIM 100 0.03 0.0 gamma (82.0, -0.0, 0.0) Diamond Pb GEANT4 150 0.08 0.0 lognorm (0.63, -0.0, 0.07) Diamond Pb G4/GECO 150 0.07 0.0 lognorm (0.58, -0.0, 0.06) Diamond Pb SRIM 150 0.03 0.0 gamma (109.62, -0.01, 0.0) Diamond Pb GEANT4 250 0.09 0.0 chi2 (4.47, 0.02, 0.01) Diamond Pb G4/GECO 250 0.09 0.0 lognorm (0.58, 0.01, 0.07) Diamond Pb SRIM 250 0.05 0.0 gamma (164.82, -0.03, 0.0) Diamond Pb GEANT4 500 0.13 0.0 lognorm (0.49, 0.04, 0.08) Diamond Pb G4/GECO 500 0.13 0.0 lognorm (0.51, 0.05, 0.07) Diamond Pb SRIM 500 0.08 0.0 gamma (268.72, -0.09, 0.0) Diamond Pb GEANT4 1000 0.2 0.0 genlogistic (345.58, -0.03, 0.04) Diamond Pb G4/GECO 1000 0.2 0.0 genlogistic (351.37, -0.03, 0.03) Diamond Pb SRIM 1000 0.13 0.0 gamma (6322.36, -1.19, 0.0) Si alpha GEANT4 100 0.68 0.03 genlogistic (0.6, 0.76, 0.08) Si alpha G4/GECO 100 0.69 0.03 genlogistic (0.59, 0.77, 0.08) Si alpha SRIM 100 0.67 0.02 genlogistic (0.34, 0.77, 0.04) Si alpha GEANT4 150 0.88 0.04 genlogistic (0.54, 0.98, 0.08) 151 Table D.1 (cont’d) Si alpha G4/GECO 150 0.89 0.04 genlogistic (0.54, 0.98, 0.08) Si alpha SRIM 150 0.87 0.02 genlogistic (0.34, 0.98, 0.04) Si alpha GEANT4 250 1.23 0.04 genlogistic (0.53, 1.33, 0.08) Si alpha G4/GECO 250 1.23 0.04 genlogistic (0.53, 1.33, 0.08) Si alpha SRIM 250 1.21 0.02 genlogistic (0.3, 1.33, 0.04) Si alpha GEANT4 500 1.99 0.04 genlogistic (0.49, 2.1, 0.08) Si alpha G4/GECO 500 1.99 0.04 genlogistic (0.51, 2.1, 0.08) Si alpha SRIM 500 1.96 0.03 genlogistic (0.33, 2.08, 0.04) Si alpha GEANT4 1000 3.54 0.05 genlogistic (0.46, 3.67, 0.08) Si alpha G4/GECO 1000 3.55 0.05 genlogistic (0.47, 3.67, 0.08) Si alpha SRIM 1000 3.5 0.04 genlogistic (0.3, 3.64, 0.05) Si H GEANT4 100 0.79 0.01 genlogistic (0.4, 0.84, 0.03) Si H G4/GECO 100 0.79 0.01 genlogistic (0.36, 0.85, 0.02) Si H SRIM 100 0.86 0.01 genlogistic (0.37, 0.92, 0.03) Si H GEANT4 150 1.24 0.01 genlogistic (0.4, 1.3, 0.03) Si H G4/GECO 150 1.24 0.01 genlogistic (0.38, 1.3, 0.03) Si H SRIM 150 1.33 0.01 genlogistic (0.35, 1.39, 0.03) Si H GEANT4 250 2.33 0.02 genlogistic (0.45, 2.4, 0.04) Si H G4/GECO 250 2.33 0.02 genlogistic (0.47, 2.39, 0.04) Si H SRIM 250 2.43 0.02 genlogistic (0.42, 2.5, 0.04) Si H GEANT4 500 5.92 0.08 genlogistic (0.41, 6.06, 0.08) Si H G4/GECO 500 5.92 0.08 genlogistic (0.44, 6.06, 0.08) Si H SRIM 500 6.09 0.07 genlogistic (0.44, 6.21, 0.07) Si H GEANT4 1000 16.15 0.22 genlogistic (0.65, 16.3, 0.2) Si H G4/GECO 1000 16.14 0.29 genlogistic (0.59, 16.34, 0.19) Si H SRIM 1000 16.56 0.3 genlogistic (0.5, 16.79, 0.16) Si N GEANT4 100 0.3 0.02 gennorm (1.81, 0.3, 0.17) Si N G4/GECO 100 0.3 0.02 gennorm (1.83, 0.3, 0.17) Si N SRIM 100 0.24 0.0 genlogistic (0.41, 0.29, 0.02) Si N GEANT4 150 0.4 0.02 genlogistic (0.73, 0.44, 0.07) Si N G4/GECO 150 0.4 0.02 genlogistic (0.77, 0.43, 0.07) Si N SRIM 150 0.35 0.01 genlogistic (0.36, 0.42, 0.03) Si N GEANT4 250 0.58 0.02 genlogistic (0.59, 0.65, 0.07) Si N G4/GECO 250 0.58 0.02 genlogistic (0.6, 0.65, 0.07) 152 Table D.1 (cont’d) Si N SRIM 250 0.55 0.01 genlogistic (0.33, 0.64, 0.03) Si N GEANT4 500 0.92 0.03 genlogistic (0.51, 1.01, 0.07) Si N G4/GECO 500 0.92 0.03 genlogistic (0.49, 1.01, 0.07) Si N SRIM 500 0.93 0.02 genlogistic (0.3, 1.04, 0.04) Si N GEANT4 1000 1.38 0.04 genlogistic (0.45, 1.49, 0.07) Si N G4/GECO 1000 1.38 0.04 genlogistic (0.46, 1.49, 0.07) Si N SRIM 1000 1.45 0.02 genlogistic (0.28, 1.57, 0.04) Si Pb GEANT4 100 0.33 0.0 genlogistic (2.11, 0.28, 0.04) Si Pb G4/GECO 100 0.32 0.0 genlogistic (5.37, 0.23, 0.04) Si Pb SRIM 100 0.05 0.0 gamma (29.8, -0.01, 0.0) Si Pb GEANT4 150 0.39 0.02 gennorm (1.22, 0.39, 0.14) Si Pb G4/GECO 150 0.41 0.02 gennorm (1.22, 0.4, 0.14) Si Pb SRIM 150 0.06 0.0 gamma (33.8, -0.02, 0.0) Si Pb GEANT4 250 0.45 0.03 genlogistic (1.03, 0.44, 0.1) Si Pb G4/GECO 250 0.45 0.03 norm (0.45, 0.18) Si Pb SRIM 250 0.08 0.0 norm (0.08, 0.02) Si Pb GEANT4 500 0.56 0.04 norm (0.56, 0.2) Si Pb G4/GECO 500 0.55 0.04 gamma (57.97, -1.0, 0.03) Si Pb SRIM 500 0.14 0.0 norm (0.14, 0.03) Si Pb GEANT4 1000 0.73 0.05 gennorm (1.89, 0.73, 0.3) Si Pb G4/GECO 1000 0.72 0.05 norm (0.72, 0.22) Si Pb SRIM 1000 0.25 0.0 gennorm (2.15, 0.25, 0.07) Ge alpha GEANT4 100 0.49 0.03 genlogistic (0.57, 0.57, 0.08) Ge alpha G4/GECO 100 0.49 0.03 genlogistic (0.56, 0.58, 0.08) Ge alpha SRIM 100 0.49 0.02 genlogistic (0.33, 0.6, 0.04) Ge alpha GEANT4 150 0.67 0.04 genlogistic (0.5, 0.78, 0.08) Ge alpha G4/GECO 150 0.67 0.04 genlogistic (0.49, 0.79, 0.08) Ge alpha SRIM 150 0.68 0.03 genlogistic (0.32, 0.81, 0.05) Ge alpha GEANT4 250 1.0 0.05 genlogistic (0.42, 1.15, 0.08) Ge alpha G4/GECO 250 1.0 0.05 genlogistic (0.42, 1.15, 0.08) Ge alpha SRIM 250 1.0 0.04 genlogistic (0.28, 1.16, 0.05) Ge alpha GEANT4 500 1.72 0.06 genlogistic (0.39, 1.89, 0.08) Ge alpha G4/GECO 500 1.73 0.06 genlogistic (0.38, 1.9, 0.08) Ge alpha SRIM 500 1.72 0.05 genlogistic (0.27, 1.9, 0.05) 153 Table D.1 (cont’d) Ge alpha GEANT4 1000 3.09 0.09 genlogistic (0.33, 3.3, 0.08) Ge alpha G4/GECO 1000 3.09 0.08 genlogistic (0.33, 3.3, 0.08) Ge alpha SRIM 1000 3.06 0.07 genlogistic (0.26, 3.27, 0.06) Ge H GEANT4 100 0.62 0.01 genlogistic (0.33, 0.7, 0.03) Ge H G4/GECO 100 0.62 0.01 genlogistic (0.34, 0.7, 0.03) Ge H SRIM 100 0.7 0.01 genlogistic (0.32, 0.79, 0.03) Ge H GEANT4 150 0.99 0.02 genlogistic (0.32, 1.09, 0.04) Ge H G4/GECO 150 0.98 0.02 genlogistic (0.32, 1.09, 0.04) Ge H SRIM 150 1.08 0.02 genlogistic (0.29, 1.18, 0.03) Ge H GEANT4 250 1.83 0.03 genlogistic (0.32, 1.96, 0.05) Ge H G4/GECO 250 1.82 0.04 genlogistic (0.31, 1.96, 0.05) Ge H SRIM 250 1.92 0.04 genlogistic (0.3, 2.05, 0.05) Ge H GEANT4 500 4.39 0.12 genlogistic (0.3, 4.62, 0.08) Ge H G4/GECO 500 4.39 0.13 genlogistic (0.29, 4.62, 0.07) Ge H SRIM 500 4.51 0.11 genlogistic (0.29, 4.72, 0.07) Ge H GEANT4 1000 11.41 0.3 genlogistic (0.39, 11.73, 0.16) Ge H G4/GECO 1000 11.41 0.32 genlogistic (0.39, 11.74, 0.16) Ge H SRIM 1000 11.55 0.37 genlogistic (0.35, 11.89, 0.14) Ge N GEANT4 100 0.19 0.01 gennorm (2.78, 0.19, 0.16) Ge N G4/GECO 100 0.19 0.01 exponpow (1.6, -0.02, 0.31) Ge N SRIM 100 0.17 0.01 exponpow (2.33, -0.04, 0.28) Ge N GEANT4 150 0.26 0.01 exponpow (2.02, -0.04, 0.42) Ge N G4/GECO 150 0.26 0.01 exponpow (2.03, -0.04, 0.42) Ge N SRIM 150 0.24 0.01 exponpow (2.96, -0.08, 0.41) Ge N GEANT4 250 0.4 0.02 exponpow (2.87, -0.1, 0.64) Ge N G4/GECO 250 0.4 0.02 exponpow (2.81, -0.09, 0.64) Ge N SRIM 250 0.39 0.02 genlogistic (0.33, 0.49, 0.04) Ge N GEANT4 500 0.7 0.03 exponpow (4.42, -0.22, 1.09) Ge N G4/GECO 500 0.7 0.03 exponpow (4.54, -0.23, 1.09) Ge N SRIM 500 0.69 0.03 exponpow (5.84, -0.37, 1.22) Ge N GEANT4 1000 1.14 0.05 laplace (1.18, 0.16) Ge N G4/GECO 1000 1.14 0.04 laplace (1.18, 0.16) Ge N SRIM 1000 1.15 0.04 laplace (1.19, 0.15) Ge Pb GEANT4 100 0.16 0.01 genlogistic (445.7, -0.23, 0.06) 154 Table D.1 (cont’d) Ge Pb G4/GECO 100 0.17 0.01 lognorm (0.38, -0.06, 0.22) Ge Pb SRIM 100 0.03 0.0 rayleigh (0.01, 0.02) Ge Pb GEANT4 150 0.18 0.01 genlogistic (473.43, -0.34, 0.08) Ge Pb G4/GECO 150 0.19 0.01 lognorm (0.4, -0.07, 0.24) Ge Pb SRIM 150 0.04 0.0 gamma (11.07, -0.01, 0.0) Ge Pb GEANT4 250 0.21 0.01 rayleigh (-0.0, 0.17) Ge Pb G4/GECO 250 0.21 0.01 rayleigh (0.0, 0.16) Ge Pb SRIM 250 0.05 0.0 norm (0.05, 0.02) Ge Pb GEANT4 500 0.29 0.02 rayleigh (0.01, 0.22) Ge Pb G4/GECO 500 0.28 0.02 norm (0.28, 0.13) Ge Pb SRIM 500 0.09 0.0 norm (0.09, 0.03) Ge Pb GEANT4 1000 0.4 0.03 gamma (12.1, -0.16, 0.05) Ge Pb G4/GECO 1000 0.4 0.03 gennorm (1.89, 0.4, 0.22) Ge Pb SRIM 1000 0.16 0.0 gamma (30.87, -0.13, 0.01) Table D.2: Statistical fits of vacancy profiles Target Material Implant Simulation Energy Mean Vacancy Straggle Best-fit Fit Parameters (keV) Depth (µm) (µm2 ) Distribution Diamond alpha GEANT4 100 0.19 0.01 genlogistic (0.09, 0.31, 0.01) Diamond alpha G4/GECO 100 0.19 0.01 genlogistic (0.09, 0.3, 0.01) Diamond alpha SRIM 100 0.25 0.01 genlogistic (0.14, 0.33, 0.01) Diamond alpha GEANT4 150 0.26 0.02 genlogistic (0.07, 0.4, 0.01) Diamond alpha G4/GECO 150 0.26 0.02 genlogistic (0.07, 0.4, 0.01) Diamond alpha SRIM 150 0.33 0.01 genlogistic (0.11, 0.44, 0.01) Diamond alpha GEANT4 250 0.38 0.04 genlogistic (0.05, 0.57, 0.01) Diamond alpha G4/GECO 250 0.38 0.04 genlogistic (0.05, 0.57, 0.01) Diamond alpha SRIM 250 0.48 0.02 genlogistic (0.08, 0.62, 0.01) Diamond alpha GEANT4 500 0.66 0.1 gennorm (0.41, 0.86, 0.02) Diamond alpha G4/GECO 500 0.66 0.1 gennorm (0.41, 0.86, 0.02) Diamond alpha SRIM 500 0.83 0.04 genlogistic (0.06, 1.01, 0.01) Diamond alpha GEANT4 1000 1.23 0.3 gennorm (0.33, 1.57, 0.01) 155 Table D.2 (cont’d) Diamond alpha G4/GECO 1000 1.23 0.31 gennorm (0.32, 1.57, 0.01) Diamond alpha SRIM 1000 1.5 0.12 cauchy (1.67, 0.08) Diamond H GEANT4 100 0.22 0.05 lognorm (3.41, -0.0, 0.02) Diamond H G4/GECO 100 0.22 0.05 lognorm (3.42, -0.0, 0.02) Diamond H SRIM 100 0.36 0.01 genlogistic (0.1, 0.46, 0.01) Diamond H GEANT4 150 0.32 0.1 lognorm (3.62, -0.0, 0.08) Diamond H G4/GECO 150 0.32 0.1 lognorm (3.38, -0.0, 0.03) Diamond H SRIM 150 0.54 0.02 genlogistic (0.07, 0.67, 0.01) Diamond H GEANT4 250 0.57 0.29 lognorm (3.38, -0.0, 0.05) Diamond H G4/GECO 250 0.57 0.29 lognorm (3.39, -0.0, 0.05) Diamond H SRIM 250 0.97 0.06 genlogistic (0.05, 1.18, 0.01) Diamond H GEANT4 500 1.33 1.78 lognorm (3.28, -0.0, 0.13) Diamond H G4/GECO 500 1.33 1.77 lognorm (3.27, -0.0, 0.13) Diamond H SRIM 500 2.44 0.44 gennorm (0.36, 2.81, 0.01) Diamond H GEANT4 1000 2.49 3.2 exponpow (0.92, -0.0, 4.36) Diamond H G4/GECO 1000 2.49 3.19 exponpow (0.9, -0.0, 4.18) Diamond H SRIM 1000 6.86 3.69 cauchy (7.9, 0.28) Diamond N GEANT4 100 0.05 0.0 gennorm (3.78, 0.05, 0.05) Diamond N G4/GECO 100 0.05 0.0 gennorm (4.13, 0.05, 0.05) Diamond N SRIM 100 0.08 0.0 exponpow (2.45, -0.02, 0.13) Diamond N GEANT4 150 0.08 0.0 gennorm (5.21, 0.08, 0.08) Diamond N G4/GECO 150 0.08 0.0 gennorm (5.19, 0.08, 0.08) Diamond N SRIM 150 0.12 0.0 exponpow (3.16, -0.05, 0.21) Diamond N GEANT4 250 0.13 0.0 gennorm (6.94, 0.12, 0.12) Diamond N G4/GECO 250 0.13 0.0 gennorm (7.46, 0.12, 0.12) Diamond N SRIM 250 0.18 0.01 genlogistic (0.14, 0.27, 0.01) Diamond N GEANT4 500 0.23 0.01 genlogistic (0.1, 0.38, 0.01) Diamond N G4/GECO 500 0.23 0.01 genlogistic (0.1, 0.38, 0.01) Diamond N SRIM 500 0.33 0.01 cauchy (0.38, 0.06) Diamond N GEANT4 1000 0.39 0.03 genlogistic (0.06, 0.61, 0.01) Diamond N G4/GECO 1000 0.39 0.03 genlogistic (0.06, 0.61, 0.01) Diamond N SRIM 1000 0.55 0.03 cauchy (0.64, 0.08) Diamond Pb GEANT4 100 0.02 0.0 gamma (1.5, -0.0, 0.01) Diamond Pb G4/GECO 100 0.02 0.0 gamma (1.3, -0.0, 0.02) 156 Table D.2 (cont’d) Diamond Pb SRIM 100 0.01 0.0 norm (0.01, 0.01) Diamond Pb GEANT4 150 0.02 0.0 chi2 (2.87, -0.0, 0.01) Diamond Pb G4/GECO 150 0.02 0.0 lognorm (0.69, -0.0, 0.02) Diamond Pb SRIM 150 0.02 0.0 norm (0.02, 0.01) Diamond Pb GEANT4 250 0.02 0.0 gamma (2.15, -0.0, 0.01) Diamond Pb G4/GECO 250 0.03 0.0 gamma (2.02, -0.0, 0.01) Diamond Pb SRIM 250 0.03 0.0 norm (0.03, 0.01) Diamond Pb GEANT4 500 0.04 0.0 exponpow (1.16, 0.0, 0.07) Diamond Pb G4/GECO 500 0.04 0.0 exponpow (1.17, 0.0, 0.07) Diamond Pb SRIM 500 0.05 0.0 norm (0.05, 0.02) Diamond Pb GEANT4 1000 0.07 0.0 rayleigh (-0.01, 0.07) Diamond Pb G4/GECO 1000 0.07 0.0 exponpow (1.23, -0.0, 0.12) Diamond Pb SRIM 1000 0.08 0.0 norm (0.08, 0.04) Si alpha GEANT4 100 0.46 0.05 genlogistic (0.32, 0.66, 0.07) Si alpha G4/GECO 100 0.47 0.05 genlogistic (0.32, 0.66, 0.07) Si alpha SRIM 100 0.53 0.04 cauchy (0.58, 0.11) Si alpha GEANT4 150 0.61 0.08 genlogistic (0.22, 0.9, 0.07) Si alpha G4/GECO 150 0.62 0.08 genlogistic (0.21, 0.9, 0.06) Si alpha SRIM 150 0.69 0.06 cauchy (0.78, 0.13) Si alpha GEANT4 250 0.88 0.14 genlogistic (0.14, 1.27, 0.06) Si alpha G4/GECO 250 0.88 0.14 genlogistic (0.14, 1.27, 0.06) Si alpha SRIM 250 0.98 0.1 cauchy (1.1, 0.15) Si alpha GEANT4 500 1.49 0.33 genlogistic (0.09, 2.08, 0.05) Si alpha G4/GECO 500 1.49 0.33 genlogistic (0.09, 2.07, 0.05) Si alpha SRIM 500 1.64 0.21 cauchy (1.84, 0.17) Si alpha GEANT4 1000 2.78 0.97 genlogistic (0.05, 3.68, 0.05) Si alpha G4/GECO 1000 2.78 0.96 genlogistic (0.05, 3.68, 0.05) Si alpha SRIM 1000 3.0 0.57 cauchy (3.37, 0.21) Si H GEANT4 100 0.61 0.09 cauchy (0.77, 0.08) Si H G4/GECO 100 0.61 0.09 cauchy (0.77, 0.08) Si H SRIM 100 0.69 0.05 genlogistic (0.11, 0.91, 0.02) Si H GEANT4 150 0.97 0.2 cauchy (1.21, 0.1) Si H G4/GECO 150 0.97 0.2 cauchy (1.21, 0.1) Si H SRIM 150 1.08 0.1 genlogistic (0.07, 1.39, 0.02) 157 Table D.2 (cont’d) Si H GEANT4 250 1.84 0.68 cauchy (2.29, 0.14) Si H G4/GECO 250 1.84 0.69 cauchy (2.29, 0.13) Si H SRIM 250 1.99 0.35 genlogistic (0.05, 2.53, 0.03) Si H GEANT4 500 4.52 4.98 cauchy (5.86, 0.29) Si H G4/GECO 500 4.52 4.95 cauchy (5.86, 0.29) Si H SRIM 500 4.95 2.38 genlogistic (0.04, 6.31, 0.05) Si H GEANT4 1000 4.62 8.43 exponpow (1.12, -0.02, 7.35) Si H G4/GECO 1000 4.63 8.44 gennorm (3.64, 4.81, 4.87) Si H SRIM 1000 13.09 20.74 gennorm (0.37, 16.29, 0.13) Si N GEANT4 100 0.15 0.01 exponpow (1.16, -0.0, 0.23) Si N G4/GECO 100 0.15 0.01 exponpow (1.16, -0.0, 0.23) Si N SRIM 100 0.17 0.01 norm (0.17, 0.08) Si N GEANT4 150 0.2 0.01 exponpow (1.23, -0.0, 0.32) Si N G4/GECO 150 0.2 0.01 exponpow (1.23, -0.0, 0.32) Si N SRIM 150 0.25 0.01 norm (0.25, 0.12) Si N GEANT4 250 0.31 0.03 gennorm (5.24, 0.32, 0.31) Si N G4/GECO 250 0.3 0.03 gennorm (5.03, 0.31, 0.31) Si N SRIM 250 0.4 0.03 norm (0.4, 0.17) Si N GEANT4 500 0.49 0.07 gennorm (7.77, 0.49, 0.49) Si N G4/GECO 500 0.49 0.07 gennorm (7.77, 0.49, 0.49) Si N SRIM 500 0.69 0.07 norm (0.69, 0.27) Si N GEANT4 1000 0.76 0.17 norm (0.76, 0.41) Si N G4/GECO 1000 0.76 0.17 norm (0.76, 0.41) Si N SRIM 1000 1.1 0.15 laplace (1.22, 0.31) Si Pb GEANT4 100 0.08 0.0 gamma (1.31, -0.0, 0.06) Si Pb G4/GECO 100 0.08 0.0 exponpow (0.85, 0.0, 0.14) Si Pb SRIM 100 0.03 0.0 norm (0.03, 0.02) Si Pb GEANT4 150 0.12 0.01 gamma (1.29, 0.0, 0.09) Si Pb G4/GECO 150 0.12 0.01 gamma (1.35, -0.0, 0.09) Si Pb SRIM 150 0.04 0.0 norm (0.04, 0.02) Si Pb GEANT4 250 0.13 0.01 gamma (1.4, -0.0, 0.09) Si Pb G4/GECO 250 0.13 0.01 gamma (1.48, -0.0, 0.09) Si Pb SRIM 250 0.06 0.0 norm (0.06, 0.03) Si Pb GEANT4 500 0.17 0.01 rayleigh (-0.03, 0.16) 158 Table D.2 (cont’d) Si Pb G4/GECO 500 0.17 0.01 rayleigh (-0.03, 0.16) Si Pb SRIM 500 0.09 0.0 norm (0.09, 0.05) Si Pb GEANT4 1000 0.23 0.02 rayleigh (-0.03, 0.21) Si Pb G4/GECO 1000 0.23 0.02 rayleigh (-0.03, 0.21) Si Pb SRIM 1000 0.16 0.01 expon (0.0, 0.16) Ge alpha GEANT4 100 0.31 0.03 gennorm (5.39, 0.3, 0.31) Ge alpha G4/GECO 100 0.3 0.03 gennorm (5.49, 0.3, 0.31) Ge alpha SRIM 100 0.38 0.03 norm (0.38, 0.16) Ge alpha GEANT4 150 0.43 0.05 gennorm (7.05, 0.41, 0.42) Ge alpha G4/GECO 150 0.43 0.05 gennorm (7.2, 0.41, 0.42) Ge alpha SRIM 150 0.53 0.05 norm (0.53, 0.21) Ge alpha GEANT4 250 0.66 0.11 genlogistic (0.18, 1.03, 0.07) Ge alpha G4/GECO 250 0.66 0.11 genlogistic (0.17, 1.03, 0.07) Ge alpha SRIM 250 0.79 0.09 laplace (0.86, 0.23) Ge alpha GEANT4 500 1.19 0.31 genlogistic (0.09, 1.82, 0.06) Ge alpha G4/GECO 500 1.19 0.31 genlogistic (0.09, 1.82, 0.06) Ge alpha SRIM 500 1.39 0.21 laplace (1.53, 0.35) Ge alpha GEANT4 1000 2.23 0.95 genlogistic (0.05, 3.24, 0.05) Ge alpha G4/GECO 1000 2.24 0.94 genlogistic (0.05, 3.24, 0.05) Ge alpha SRIM 1000 2.52 0.58 laplace (2.81, 0.54) Ge H GEANT4 100 0.43 0.07 chi2 (0.92, -0.0, 0.63) Ge H G4/GECO 100 0.43 0.07 chi2 (0.94, -0.0, 0.63) Ge H SRIM 100 0.54 0.04 genlogistic (0.17, 0.75, 0.04) Ge H GEANT4 150 0.69 0.17 genlogistic (0.06, 1.13, 0.03) Ge H G4/GECO 150 0.69 0.17 powerlaw (0.55, -0.0, 1.33) Ge H SRIM 150 0.83 0.09 genlogistic (0.12, 1.15, 0.04) Ge H GEANT4 250 1.3 0.56 chi2 (1.03, -0.0, 1.39) Ge H G4/GECO 250 1.3 0.55 chi2 (1.03, -0.0, 1.39) Ge H SRIM 250 1.5 0.26 genlogistic (0.08, 2.04, 0.04) Ge H GEANT4 500 2.84 3.54 chi2 (1.08, -0.0, 2.52) Ge H G4/GECO 500 2.85 3.56 chi2 (1.09, -0.0, 2.5) Ge H SRIM 500 3.5 1.53 genlogistic (0.04, 4.75, 0.06) Ge H GEANT4 1000 3.26 4.1 gennorm (3.59, 3.39, 3.4) Ge H G4/GECO 1000 3.26 4.12 gennorm (3.58, 3.39, 3.4) 159 Table D.2 (cont’d) Ge H SRIM 1000 8.83 10.62 genlogistic (0.03, 12.05, 0.1) Ge N GEANT4 100 0.09 0.0 exponpow (1.08, -0.0, 0.15) Ge N G4/GECO 100 0.09 0.0 exponpow (1.09, 0.0, 0.15) Ge N SRIM 100 0.12 0.0 norm (0.12, 0.07) Ge N GEANT4 150 0.13 0.01 exponpow (1.16, -0.0, 0.21) Ge N G4/GECO 150 0.13 0.01 exponpow (1.16, -0.0, 0.21) Ge N SRIM 150 0.18 0.01 norm (0.18, 0.09) Ge N GEANT4 250 0.22 0.02 exponpow (1.24, -0.0, 0.34) Ge N G4/GECO 250 0.22 0.02 exponpow (1.25, -0.0, 0.34) Ge N SRIM 250 0.29 0.02 norm (0.29, 0.14) Ge N GEANT4 500 0.39 0.05 gennorm (6.56, 0.39, 0.39) Ge N G4/GECO 500 0.39 0.05 norm (0.39, 0.22) Ge N SRIM 500 0.52 0.05 norm (0.52, 0.22) Ge N GEANT4 1000 0.64 0.12 norm (0.64, 0.35) Ge N G4/GECO 1000 0.64 0.12 norm (0.64, 0.35) Ge N SRIM 1000 0.87 0.12 norm (0.87, 0.34) Ge Pb GEANT4 100 0.07 0.0 exponpow (0.76, 0.0, 0.13) Ge Pb G4/GECO 100 0.08 0.0 exponpow (0.79, 0.0, 0.15) Ge Pb SRIM 100 0.02 0.0 norm (0.02, 0.01) Ge Pb GEANT4 150 0.06 0.0 chi2 (1.73, 0.0, 0.04) Ge Pb G4/GECO 150 0.06 0.0 chi2 (1.68, 0.0, 0.04) Ge Pb SRIM 150 0.03 0.0 norm (0.03, 0.02) Ge Pb GEANT4 250 0.07 0.0 gamma (0.97, 0.0, 0.07) Ge Pb G4/GECO 250 0.07 0.01 expon (0.0, 0.07) Ge Pb SRIM 250 0.04 0.0 expon (0.0, 0.04) Ge Pb GEANT4 500 0.1 0.01 expon (0.0, 0.1) Ge Pb G4/GECO 500 0.09 0.01 gamma (1.25, -0.0, 0.08) Ge Pb SRIM 500 0.07 0.0 expon (0.0, 0.07) Ge Pb GEANT4 1000 0.14 0.01 rayleigh (-0.03, 0.14) Ge Pb G4/GECO 1000 0.15 0.01 rayleigh (-0.03, 0.15) Ge Pb SRIM 1000 0.12 0.0 expon (0.0, 0.12) 160 APPENDIX E. CCANARD Quick Start Guide The Color Center ANnealing and Reaction Diffusion (CCANARD) simulation is a computer program designed to solve a system of nonlinear coupled partial differential equations de- scribing the diffusion and interaction of point defects in a crystal lattice during the annealing process. This software was developed by Hank Thurston as a portion of his doctoral research and is available by request. No explicit or implicit guarantee or warranty is provided, to neither the functionality of the code nor accuracy of results. CCANARD was developed on Calculate Linux, a Gentoo based 64-bit Linux operating system. CCANARD should run without modification on any 64-bit Linux system with the GNU tool chain installed. Additional features of CCANARD may be available with other compiler suite and hardware combinations. E.1 Minimum System Requirements ˆ 64-bit Linux operating system ˆ GNU tool chain ˆ Python3.x installed on PATH ˆ Access to SRIM software ˆ 4GiB RAM 161 E.2 Installing CCANARD The CCANARD source code is hosted on GitLab. To clone the repository use the following command in a terminal emulator: $ git clone git@gitlab.com:HankT/ccanard.git This will clone the git repository which includes a source code directory named src, a shell script, build.sh, and a README. To build CCANARD: 1. cd into CCANARD directory: $ cd CCANARD 2. excute build.sh providing the desired compiler as an argument: $ ./build.sh There are four recognized arguments to the build.sh command: ˆ gfortran: Builds CCANARD using gfortran (gfortran must be on PATH) ˆ nvfortran: Builds CCANARD using nvfortran, the NVIDIA Fortran compiler (nvfortran must be on PATH) ˆ ifort: Builds CCANARD with ifort, the Intel Fortran compiler (ifort must be on PATH) ˆ clean: removes all object files, module files, binary files, and symlinks. This cleans the CCANARD environment for rebuilding. 3. execute CCANARD program: $ ./CCANARD (Or click the CCANARD symlink) E.3 Running a CCANARD Simulation The first step to running a CCANARD simulation is to generate a set of initial conditions. This provides the initial concentrations of point defects. In CCANARD version 00.01 the 162 input conditions are provided by SRIM2013 simulations. SRIM simulations must be run with the “Full Damage Cascade” option selected to generate the necessary data for CCA- NARD. CCANARD will use the ASCII files produced by SRIM, specifically EXYZ.txt and COLLISON.txt. Once a SRIM simulation has been run to generate initial conditions, the operator may launch the CCANARD program. The operator will be met with the main menu, shown in Figure E.1. Figure E.1: CCANARD main menu, visible when first entering the program. From the main menu the operator will be prompted to enter a number corresponding to the options presented in the terminal. Options 1-5 must be confirmed before a CCANARD simulation can begin. Option 0 will show the CCANARD “About” informational display. 163 E.3.1 Option 1) Simulation Control Figure E.2: Menu Option 1) Simulation Control. This menu allows the modification of administrative properties of a CCANARD simulation. E.3.1.1 Folders and Directories This section of the simulation control menu contains information about where simulation files are loaded and saved. The options are: ∗ Simulation name: The operator may choose to give the simulation a unique name for ease of later reference. By default the name will reference the date and time. The simulation name should not contain any spaces. ∗ Simulation directory: The operator may choose to specify the directory in which CCANARD will store output. By default it will be saved in a directory labeled with the date and time placed in the operator’s Home directory. The absolute path must be provided using Unix syntax and be enclosed in “ ”. 164 ∗ SRIM directory containing EXYZ.txt and COLLISON.txt: The operator must tell CCANARD the directory containing the SRIM ASCII output files EXYZ.txt and COL- LISON.txt. The absolute path must be provided using Unix syntax and be enclosed in “ ”. These files will be copies into the simulation directory. E.3.1.2 Runtime Parameters This section of the simulation control menu contains options which change the performance and communication of CCANARD. The options are: ∗ Console print interval: This is an integer value which specifies how frequently CCANARD will print the current simulation status to the console. Each interval is one time step. The default value is 100. ∗ Number of CPU cores: The operator may choose how many CPU cores will be used to parallelize CCANARD. The default value is 4. The integer entered here is a re- quest; CCANARD will try to provide the specified number of CPU cores, but cannot guarantee this number is provided. ∗ Use GPU acceleration: This option can be toggled to True or False, providing GPU acceleration to portions of the CCANARD program. The default value if False. This option can only be turned on if the specific compiler used and hardware available permits the use of GPU acceleration. This is currently only available with nvfortran and Nvidia GPUs. ∗ Verbosity level: This option can be toggled between values of 0 and 3, inclusive, specifying the amount of detail printed to the console and recorded in the log file during a CCANARD simulation. 0 turns all console and log output off. 1 provides 165 only the most important messages. 2 provides important messages and regular runtime updates. 3 outputs regular runtime updates and all messages, including debugging messages. The default value is 2. Once all settings have been set as desired, the operator will press [Enter] and be prompted to confirm the settings. After double checking the values, Y should be entered to confirm or N to deny. Once the settings have been confirmed, a check mark will appear next to the option on the main menu, indicating the settings have been confirmed. E.3.2 Option 2) Physical Domain Figure E.3: Menu Option 2) Physical Domain. This menu allows the modification of the physical and material properties modeled in a CCANARD simulation. E.3.2.1 Diffusion Properties ∗ Simulate color center formation: This option can be toggled to True or False. When True, CCANARD will model the pseudo-covalent bonding between vacancy and substitutional defects. The default value is True. 166 ∗ Simulate interstitial-vacancy recombination: This option can be toggled to True or False. When True, CCANARD will model the recombination of interstitial host material atoms into lattice vacancies. The default value is True. ∗ Reaction Range for color center formation: This value specifies the distance at which vacancy and substitution defects can interact. [reaction range divided by bond length] gives the number of coordination spheres in which the reaction can occur. The default value is 1.545 Angstroms, equal to a single coordination sphere in diamond. ∗ Reaction Range for vacancy-interstitial recombination: This value specifies the distance at which vacancy and host material interstitial atoms can interact. [reac- tion range divided by bond length] gives the number of coordination spheres in which the reaction can occur. The default value is 1.545 Angstroms, equal to a single coor- dination sphere in diamond. ∗ Maximal Diffusivity constant: This value, in units of (cm2 /s), is the theoreti- cal diffusion constant as temperature approaches infinity. The default value is D0 = 0.360E-05, the maximal diffusivity of nitrogen in diamond. ∗ Ion mobility activation energy: This is the energy, in units of eV, at which sub- stitutional implant ions become mobile. The default value is 5.600 eV, the mobility activation energy of nitrogen implants in diamond. ∗ Vacancy mobility activation energy: This is the energy, in units of eV, at which vacancy defects become mobile. The default value is 2.600 eV, the mobility activation energy of vacancies in diamond. ∗ Interstitial mobility activation energy: This is the energy, in units of eV, at which interstitial defects become mobile. The default value is 1.839 eV, the mobility 167 activation energy of carbon interstitial atoms in diamond. ∗ Color-center mobility activation energy: This is the energy, in units of eV, at which color centers become mobile. The default value is 10.000 eV, an arbitrarily high value to prevent the mobility of color centers by default. ∗ More Options...: Pressing 0 will display a second page of Physical Domain options. These additional options must also be confirmed before beginning the CCANARD simulation. Figure E.4: Additional physics options accessed through the More Options... selection of the Physical Domain menu. E.3.2.2 Material Properties ∗ Atomic Density: This is the number of atoms per cubic centimeter of host material. The default value is 0.1770E+24 (atoms/cc), the atomic density of diamond. ∗ Atomic bond length: This is the value, in Angstroms, of the bond length of the host material. The default value is 1.545 Angstroms, the bond length of diamond. 168 ∗ Length of sample in x (depth) direction: This is the length of the sample to be modeled in units of microns. The x-direction is defined as the direction of beam propagation in ion implantation. The default value is 5.00µm. ∗ Length of sample in y (lateral) direction: This is the length of the sample to be modeled in units of microns. The y-direction is defined as the direction parallel to the lab floor in a physical sample. The default value is 5.00µm. ∗ Length of sample in z (vertical) direction: This is the length of the sample to be modeled in units of microns. The z-direction is defined as the vertical direction in a physical sample. The default value is 5.00µm. E.3.2.3 Implant Procedure ∗ Implant beam spot size: This value is the area, in units of mm2 , of the ion implant beam that will be used in experiment or physical construction of the sample being modeled in CCANARD. The default value is 1mm2 ∗ Implant ion fluence: This value is the fluence, in units of ions/cm2 , of the ion implant beam that will be used in experiment or physical construction of the sample being modeled in CCANARD. The default value is 0.100E+16 ions/cm2 . ∗ Maximum ion concentration: this is the maximum concentration of implant ions in the sample. By default this value will be calculated using the provided Implant ion fluence. If the operator instead desires to provide a desired maximum ion concentration, the implant ion fluence will be calculated. Once all settings have been set as desired, the operator will press [Enter] and be prompted to confirm the settings. After double checking the values, Y should be entered to confirm 169 or N to deny. Once the settings have been confirmed, a check mark will appear next to the option on the main menu, indicating the settings have been confirmed. E.3.3 Option 3) Discretization and Solver Figure E.5: Menu Option 3) Discretization and Solver. This menu allows the operator to specify how CCANARD will discretize the domain and solve systems of coupled PDEs. E.3.3.1 System Discretization ∗ Number of nodes in x-direction: This option takes an integer value specifying the number of finite difference nodes in the x-direction to be used during system dis- cretization. The x-direction is defined as the direction of beam propagation in ion implantation. The default value is 21. ∗ Number of nodes in y-direction: This option takes an integer value specifying the number of finite difference nodes in the y-direction to be used during system discretiza- tion. The y-direction is defined as the direction parallel to the lab floor in a physical sample. The default value is 21. 170 ∗ Number of nodes in z-direction: This option takes an integer value specifying the number of finite difference nodes in the z-direction to be used during system discretiza- tion. This is the length of the sample to be modeled in units of microns. The z-direction is defined as the vertical direction in a physical sample. The default value is 21. E.3.3.2 Solver Parameters ∗ Solver method: This option specifies the Krylov subspace method which will be used to solve the linear system. The operator can toggle between CG (Conjugate Gradient) and PCCG (Preconditioned Conjugate Gradient). PCCG uses Jacobi preconditioning. The default value is CG which performs suitably well in all test cases, and outperforms PCCG in some. ∗ Solver convergence tolerance: this option specifies the tolerance which must be reached to consider the linear solved and terminate Krylov iteration. The default value is 0.1E-05. ∗ Maximum Krylov iterations: This option specifies the maximum number of Krylov iterations that will be attempted. If the solver has not reached convergence in this number of iterations, solving will terminate and return NaN. ∗ Time step ∆t: This option allows the operator to specify the time step, in seconds, to be used in the Crank-Nicolson solver. The default value is 5.00s. Once all settings have been set as desired, the operator will press [Enter] and be prompted to confirm the settings. After double checking the values, Y should be entered to confirm or N to deny. Once the settings have been confirmed, a check mark will appear next to the option on the main menu, indicating the settings have been confirmed. 171 E.3.4 Option 4) Annealing Options Figure E.6: Menu Option 4) Annealing options. This menu allows the operator to enter the annealing recipe which CCANARD will simulate. E.3.4.1 Plots This section allows the operator to specify the annealing recipe that will be used in the CCANARD simulation of diffusion and defect interaction. The annealing recipe is built with a series of individual steps, each specifying a start time, end time, starting temperature, and ending temperature. All times are given in seconds and all temperatures are given in Kelvin. A default recipe is provided consisting of a one-hour ramp from room temperature to 1073 Kelvin (800 deg C) followed by a one-hour hold at 1073K. The operator may choose to: 1) Add a subsequent step to the existing recipe 2) Remove the last step form the recipe 3) Clear the recipe 172 The first recipe step must begin at 0.0 seconds. Each subsequent recipe step must be continuous in both time and temperature. Discontinuous jumps are not permitted in the recipe. A recipe may hold up to 20 steps. Once the desired recipe has been built, the operator will press [Enter] and be prompted to confirm the settings. After double checking the values, Y should be entered to confirm or N to deny. Once the settings have been confirmed, a check mark will appear next to the option on the main menu, indicating the settings have been confirmed. E.3.5 Option 5) Output Options Figure E.7: Menu Option 5) Output Options. This menu allows the operator to specify and modify CCANARD’s visual output in the form of plots and animations. Up to 4 plots may be specified by the operator. Plots show the concentration of the specified defect at a specified plane in the sample. The operator can select to edit a plot by pressing [1]-[4]. Each plot may be turned on or off. The operator may specify a plane for each plot. The plane choices are the XY, YZ, 173 and ZX planes. The operator will also be prompted for the intercept value for the given plane. The value of x=0 for a YX plane corresponds the to the front surface of the sample. The values of y=0, z=0 correspond to planes bisecting the sample in the ZX, or XY plane, respectively. Finally the operator may choose the defect to be plotted. The default values are shown below in table E.1. Table E.1: Default plots generated by CCANARD. Plot Status Plane Plane intercept Species 1 on XY 0.000 ion 2 on XY 0.000 vac 3 on XY 0.000 int 4 on XY 0.000 ccs E.3.5.1 Output Method ∗ Create animation: This option can be toggled True or False. If True, an animation will be produced at the completion of the CCANARD simulation using ffmpeg. This requires ffmpeg to be on the PATH. This option also requires plots to be saved at the console print interval. The default value is False. ∗ Print plots to monitor: This option can be toggled True or False. If True, plots will be displayed on the monitor and be saved to disk for future review. If False, plots will be saved to disk, but will NOT be displayed to the screen. The default value is True. ∗ Save plots at print interval: This option can be toggled True or False. If True, plots will be made and recorded every time CCANARD alerts the operator with normal 174 runtime status updates. The frequency of these updates can be set in the Simulation Control menu. The default value is False. E.3.5.2 Color Map This option allows the operator to choose a color map for plots and animations. The options are shown below in Figure E.8. Figure E.8: These eight color maps are available in CCANARD. Once all settings have been set as desired, the operator will press [Enter] and be prompted 175 to confirm the settings. After double checking the values, Y should be entered to confirm or N to deny. Once the settings have been confirmed, a check mark will appear next to the option on the main menu, indicating the settings have been confirmed. E.3.6 Option 6) Begin Simulation This option will only be available once all the settings have been confirmed. Once the desired configuration has been set, pressing [6] will begin the simulations. Figure E.9 shows all the settings have been verified and the operator has the ability to begin the simualation. Figure E.9: All settings have been verified, indicated by check marks. The operator is able to begin the simulation. E.3.7 Option 0) About This option will display the “About” informational display. This will include the version of CCANARD, copyright information, and the compiler used to build the running binary. Press [Enter] to exit this window. 176 Figure E.10: The CCANARD “About” informational display reporting that it was compiled with nvfortran version 22. Note the ASCII ducks, the unofficial logo of CCANARD. (Canard is the French word for duck.) E.4 Output Files Within the simulation directory specified in the Simulation Control menu, CCANARD will create several files and folders as shown below in Figure E.11. 177 Figure E.11: The directory structure of a CCANARD simulation. The files are as follows: ∗ COLLISON.txt: A SRIM output file which CCANARD copies into the simulation di- rectory ∗ EXYZ.txt: A SRIM output file which CCANARD copies into the simulation directory ∗ interstital.txt: The (x,y,z) coordinates of the initial interstitial distribution ∗ ion.txt: The (x,y,z) coordinates of the initial ion distribution ∗ vacancy.txt: The (x,y,z) coordinates of the initial vacancy distribution ∗ : plots saved as .png image files ∗