An Electro-Optical Simulation Methodology for the Analysis of Single-Event Radiation Effects in Photonic Devices

A multiphysics simulation methodology is presented for the analysis of transient radiation effects in photonic integrated devices and circuits. Electrical and optical simulators are coupled by computing steady-state optical solutions at various time intervals during transient charge transport simulations. This approach allows for the simulation of transient effects, enabling new studies of radiation effects in photonic integrated circuits (PICs), currently limited to analysis of single optical devices. This methodology will also allow for the development of compact models for lumped parameter simulation of radiation effects at the circuit and system levels, which will aid the design and analysis of PICs for operation in space or other harsh radiation environments.

x device. A power absorption profile is computed, and this is passed to stage B) where an impulse of charges is generated, diffused, and recombined to capture the electronic effects of the ion strike. A set of discrete-time 3D charge profiles are then generated, and these are passed to stage C) where each profile is used to compute the shifts in optical phase and power loss that represent the overall device response to the ion strike. While this figure shows Lumerical  14 The object tree in Lumerical FDTD with the tools required for the ion strike simulation. The objects named "waveguide", "buried oxide", and "surface oxide" are the shapes used to define the dimensions and materials of the waveguide in this example. "FDTD" is the FDTD solver object that defines the simulation region and mesh and the object called "mesh" is the mesh override object that is defining a finer mesh directly surrounding the ion strike location. Unique to this simulation is the Gaussian light source (called "Strike") and the power absorption monitor 20 The gaussian source geometry settings. It is set to be a "z-normal" source so that it injects light in the z direction and is set to be a diameter of about 50 nm with "x span" and "y span". The z coordinate is set to 220 nm to place it on top of the waveguide .. 32 Lumerical DEVICE electrical contact settings. To define a shape in the simulation region as an electrical contact, (a) was set to a solid geometry and (b) was then set to the specific solid name that the user wishes to define as the contact (here it is set to the "contact" as seen in Figure 24). (c1) was set to initial time and voltage of the contact, while (c2) was set to the simulation end time and voltage. Both voltages were set to 0 to ensure that the electrical contact does not affect the passive diffusion simulation in this example .. 35 The settings of the charge monitor. Both "record electrons" and "record holes" must be selected. "Save data" allows the charge monitor to automatically save the computed charge profile to a file named in the "filename" field. "Integrate total charge" computes the total number of electrons and holes within the charge monitor which is convenient for estimating LET . 40 The FDTD simulation region dimensions for the parameter extraction simulation. The length of this simulation region (y-direction) is set to 51 μm so that it is slightly longer than the imported charge profile. The other dimensions are identical to the simulation region in Figure 15 making it a 1. other radiation harsh environments. More recently, research has shown that photonic integrated circuits (PICs) may also be affected by ionizing radiation [1,2,3,4,5]. This discovery has opened up a whole new field of radiation effects research in photonic devices that is still in its infant stages.
Investigating the effects of radiation on PICs is further motivated by their growing list of potential applications in space systems, including communications, sensing, autonomous navigation, and imaging [6].
Radiation effects are broken down into two primary subclasses of study: Total ionizing dose (TID) and single-event effects (SEEs). TID is caused by continued exposure of circuits to ionizing radiation and may cause the physical parameters of devices to degrade, such as transistor threshold voltages and transconductance, due to trapped charge in the transistor oxides [7]. SEEs on the other hand are the result of the interactions of single ionizing particles with semiconductor junctions, which may cause current and voltage fluctuations (called a single-event transient or SET) capable of corrupting data stored in the circuit. The severity of a SET depends on the charge generated and collected within a semiconductor junction and the response of the circuit to that transient current [8].
Until now, radiation effects studies on PICs have been primarily limited to studies of TID [1,3,4,5], however, recent studies have demonstrated that SEEs may also be a significant area of study for PICs [2]. In particular, there is currently a lack of ability to simulate and model the time state of research on radiation effects in PICs is then described to illustrate the reasons for developing a simulation methodology for SEEs in PICs. Finally, the theory behind electronic charge transport and optical FDTD simulations is described.

Photonic Integrated Circuits
Conventional electronics have facilitated the transmission and processing of electrical signals for more than a century. More recently, however, photonic technology has allowed for the transmission and processing of optical (light) signals that yield much higher data bandwidths, reaching into the tens of Gbps per channel with the possibility of scaling into the range of Tbps if parallel channels are implemented [9]. In addition to the commercial applications, photonics is also expected to see increased use in space systems to augment the performance of conventional electronic systems [6]. Conveniently, the same silicon (Si) technology that allows the fabrication of CMOS transistors, also facilitates the fabrication of Si photonic integrated circuits. This may allow for full photonic/electronic systems to be implemented on the same chip, assuming that the challenge of fabricated lasers on a chip can be overcome.
The primary component of any photonic system is the waveguide. Waveguides act as the means of transporting information in an optical system. Any optical waveguide requires at least two materials to function: A core material transparent to the carrier wavelength, and a bounding material that will reflect the carrier to continuously guide it within the waveguide core. With current Si photonic circuit technology, it is possible to fabricate a basic waveguide as a slab of Si  Although SiO2 is also transparent to these wavelengths, the theory of total internal reflection allows the light to still be reflected back into the Si at the boundary between the Si and SiO2 [10] (Figure 1). This occurs due to the difference in refractive index between the materials where the Si has refractive index 3.48 at 1550 nm [11] and the SiO2 has refractive index of 1.44 at 1550 nm [12]. Note that the refractive index of the bounding material (SiO2) is lower than that of the core material (Si). This is necessary for total internal reflection to guide the light through the core material.
As the wave reflects between the boundaries of the waveguide, standing wave patterns are created along the directions perpendicular to the direction of propagation. These standing wave patterns (illustrated in Figure 2), referred to as modes of propagation, are derived as solutions to the wave equation resulting from Maxwell's equations. The shape and type of mode is dependent upon the physical dimensions of the waveguide and the wavelength of the carrier signal that is being propagated down the waveguide. The mode profile may be used as a means of determining how well guided the mode is, that is how much of the signal energy falls within the waveguide itself, and how much gets radiated and lost in the SiO2. The more the mode profile falls outside of the waveguide, the more optical energy is lost per unit length of waveguide.
One way to categorize modes is by their order. A 1 st order mode is excited in the waveguide if the mode profile has a single peak at the center of the waveguide ( Figure 2). This 1 st order mode is also called the fundamental mode. A 2 nd order mode may be excited, which consists of two peaks of opposite polarity. A 3 rd order mode may be excited that has three peaks in the field. This trend of modes may in theory continue indefinitely; however, higher order modes tend to be less well guided and thus subject to higher power losses. Higher order modes may often be excited as artifacts of the excitation of the fundamental mode in an imperfect waveguide [13] or they may be intentionally excited via a device called a mode transducer that converts one mode to another [14].
These modal orders are analogous to the vibrational modes of vibrating strings that are tied at each end. A waveguide should be designed such that each of its physical dimensions are equal to a multiple of half the carrier wavelength as this allows for the modes to have minimal losses.

Figure 2
Transmission modes of optical waveguides. Modes may take on various shapes that depend on the waveguide dimensions, the optical properties of the waveguide material, the carrier wavelength. Left illustrates a transverse magnetic (TM) mode and right illustrates a transverse electric (TE) mode. A TM mode has its magnetic field transverse (or perpendicular) to the direction of propagation and its electric field aligned in the direction of propagation. A TE mode is the opposite having its electric field transverse and its magnetic field in the direction of propagation Modes are also categorized by their transverse field [15]. A transverse electric (TE) mode will have its electric field polarized such that it is transverse (perpendicular) to the direction of propagation, and its magnetic field will be polarized in the direction of propagation. A transverse magnetic (TM) mode is the opposite of TE, with its magnetic field polarized transverse to the direction of propagation and its electric field polarized in the direction of propagation. A pure TE or TM mode is only possible in theory. In a physical waveguide, a mode manifests itself as some superposition of both TE and TM, though it is possible to get close to pure TE or TM.
PICs include a variety of components other than waveguides. Passive components include waveguides, ring resonators, interferometers, and couplers, for example. There are also active components such as modulators, optical amplifiers, lasers, and photodiodes. Though there are numerous types of optical components, most of them are fabricated as different structures of waveguides, such as ring resonators, modulators, directional couplers, and many more. Therefore, to make this research generalizable as possible to all optical components, this thesis focuses solely on the response of waveguides to ionizing radiation.

Radiation Environments
The design of electronic and photonic devices that operate within the earth's troposphere generally requires little concern for radiation effects. Some of the limited exceptions to this would be systems that operate within or near radiation test facilities and high energy particle accelerators [1]. However, as systems leave the earth's atmosphere, they enter a space environment filled will radiation from a variety of sources ( Figure 3) [11].
One reason for the lower earth atmosphere having limited exposure to radiation is that the earth has a magnetic field that traps ionizing radiation that comes from solar flares, effectively shielding the earth from being struck with ionizing particles. The regions in which the radiation becomes trapped are known as the Van Allen Belts. These radiation belts pose a problem to electronic systems that pass through them to travel beyond the earth as well as for satellites that orbit the earth and must travel through the belts continually.
Further away from the earth is an environment susceptible to radiation that comes directly from the sun. During a solar flare, it is likely that a system will be struck by many protons and heavy ions. Another radiation concern in this environment comes in the form of galactic cosmic rays (GCR) which are high energy particles that are traveling at significant fractions of the speed of light. GCRs originate from outside of the solar system from events elsewhere in the galaxy such as supernovas. [16] Figure 3 The near-earth radiation environment consists of radiation from various sources, including the belts of trapped particles surrounding the earth, solar activity and events, and galactic cosmic rays (GCRs) that occur outside of the solar system [17]

Radiation Effects
When ionizing particles strike semiconductor material, they transfer energy to the system in the form of generating free charges in the material. There are two primary subclasses of effects that these free charges may have on a system. One type occurs when a circuit is under long-term exposure to ionizing radiation and extraneous positive charges get trapped in the device oxides.
This effect is due to the low mobility of holes in the oxides. In CMOS technology, these trapped excess charges have the effect of degrading device parameters such as transistor threshold voltages and gate capacitances. The severity of the effect is proportional to the amount of radiation exposure in both time and intensity and is classified by the total dose which is measured in units of kRads (1 kRad = 10 J/kg). As such, these types of effects are called total ionizing dose (TID) effects.
The second subclass of radiation effect is the instantaneous transient effect that occurs when a single ionizing particle strikes a component (known as single-event effects or SEEs). In CMOS technology, the free charge that is generated by the ion strike can be collected at nearby PN junctions through the electric field in the depletion region ( Figure 4) [12]. This results in a current pulse at each of the affected circuit nodes. The severity of this effect depends on the charge that is generated by the ion strike, the total charge collected at the circuit node, and the response of the circuit itself.
The amount of charge that is generated in the device depends on the total energy that is transferred from the ion to the semiconductor material. In the study of SEEs, this transfer is quantified by the linear energy transfer (LET) of a particle (measured in units of MeV cm 2 /mg).
Due to the generation of larger amounts of free carriers, higher LET values equate to a larger effect on the overall circuit.
The current pulse that is generated in the device junctions in response to the ion strike causes a perturbation in the voltage at the affected nodes that is known as a single-event transient (SET). The intensity of a SET is highly dependent on the capacitance at the affected node, lower node capacitances being more susceptible to SETs due to them requiring less total charge to perturb the voltage. If a SET is severe enough to cause a bit to flip in a circuit that results in information corruption, this is known as a single-event upset (SEU).

Figure 4
An ion striking a Si PN junction causes the generation of free carriers that are swept away by the electric field in the depletion region (left) and causes a current transient at the circuit nodes connected to the PN junction (right) [18]

Single-Photon Absorption and Two-Photon Absorption
Laser pulses are often injected into electronic devices to emulate the effect of an ion striking the device [19,20]. This is because photons absorb into silicon to generate EHPs similarly to how the transfer of energy from ionizing particles may generate EHPs during a single event.
There are two distinct processes by which lasers are used to induce charges in silicon: single-photon absorption (SPA) and two-photon absorption (TPA). In SPA, a single photon interacts with the material to generate a single EHP. In TPA, two photons interact with the material simultaneously to generate a single EHP.

Figure 5
The optical absorption coefficients for single-photon absorption (obtained from [21]) and twophoton absorption (obtained from [22]) in Si The equation describing the combined effect of SPA and TPA is shown below [19].
In Equation 2.1, N is the concentration of generated EHPs in cm -3 , ( ) is the wavelength dependent SPA coefficient in cm -1 , ( ) is the wavelength dependent TPA coefficient in cm/GW, and is the irradiance. The values of ( ) (obtained from Jones et. al. [21]) and ( ) (obtained from Bristow et. al. [22]) are plotted versus wavelength in Figure 5. The coefficients of SPA and TPA intersect at a wavelength of about 1150 nm which is about the wavelength of a photon whose energy equals the Si bandgap of 1.1 eV. To the left of this point, for small wavelengths, the SPA mechanism dominates since the energy of each photon is greater than the band gap of Si. For wavelengths greater than 1150 nm, the TPA mechanism dominates because two photons are required to efficiently generate EHPs.
According to Buchner et. al. [19], both SPA and TPA are useful for emulating ion strikes in silicon devices, however, SPA doesn't allow for the injection of carriers at variable depth, which can be a problem when there are other fabrication layers blocking the laser from reaching the device. This issue may be solved by using a laser with a TPA wavelength (above 1150 nm) so that it can propagate through the other fabrication layers and focusing it so that the irradiance becomes very strong at the point where the carrier generation is desired. Carriers may be generated at this focal point because of the square relationship for the irradiance of the TPA term in Equation 2.1.

Research of Single-Event Effects in Photonics
Radiation effects have been studied in the context of CMOS technology for decades, however, PICs have received limited treatment of radiation effect studies. Primarily, studies have been focused on TID effects in PICs. Some work has been accomplished in observing the behavior of photonic devices due to TID radiation exposure as in Bhandaru et. al. [3] Other research has been done in identifying the mechanisms that cause degradation of the performance of a Mach-Zehnder modulator [4] and to model the effects caused by the TID [1]. A limited set of work has studied SEEs in PICs In the field of optics, beams of light have been used to modulate the transmissive properties of Si waveguides in a similar manner to ionizing radiation [23,24]. Since light is itself a form of radiation, pulsed laser light is often used to emulate an ion strike in electronics [20] and may also be used in a similar manner for photonic devices. A recent paper by Goley et. al. [2] has demonstrated that SEEs may be of some significance to the operation of PICs in a space environment. This furthers the incentive to develop a comprehensive understanding of the mechanisms that may be involved in SEEs in PICs.
In the Goley paper [2], simulations were performed in the Lumerical Tool Suite for a waveguide that was given a profile of ion induced charge with an associated LET. For each charge profile, one optical simulation was performed by using the carrier densities of the profile as the input to a charge-to-index conversion model called the Soref and Bennett model [25] represented in Figure 6 [2]. The charge-to-index conversion model allows the simulation software to shift the refractive index of the Si in the waveguide in response to the charges induced by the single-event.
Now that the perturbation of the waveguide's optical properties can be simulated, an optical simulation is used to compute the optical power losses and phase shifts seen at the output of the waveguide for each of the given charge profiles. A plot of the transmission loss versus LET with a trend curve is shown in Figure 7 [2].
The Soref and Bennett charge to index conversion model is represented in [25] by the where Δ and Δ represent the change in refractive index and optical absorption respectively from their non-excited values at a specified optical wavelength. The values Δ and Δ are the change in the Si's hole and electron concentrations in cm -3 , and the values _ , _ , _ , _ , _ , _ , _ , and _ is a set of coefficients that depend on the optical wavelength.
The values of these coefficients for a wavelength of 1550 nm are given in Table 1.

Figure 6
A graph of Si refractive index versus EHP density for the Soref and Bennet charge-to-index conversion model. As EHP density increases in Si, decreases in refractive index are observed. As the EHP density rises to about 1.83e21 cm -3 the refractive index goes below 1 and starts to become reflective and absorbing like a metal [2] Figure 7 Trend of waveguide transmission losses versus LET with trendline from [2]. These values were computed in Lumerical FDTD for static charge profiles in a slab waveguide of cross-sectional dimensions of 220 nm wide by 450 nm high The mechanism of single-event radiation effects in a Si waveguide. Charges induced at the ion strike location cause shifts in the Si's refractive index. Significant free charges cause the Si to become reflective and absorbing In Equation 2.5, is the electrostatic potential, 0 is vacuum permittivity, is the charge of an electron, is hole density, is electron density, is donor density, and is acceptor density.
Poisson's equation relates the electrostatic potential (and thus the electric field) to the local charge density. This equation allows the charge transport solver to obtain the electric field which is required to solve the drift-diffusion equations.
The drift-diffusion equations are given by where and are the electron and hole current density vectors respectively, and are the material electron and hole mobilities respectively, and are the electron and hole diffusivities respectively, and is the electric field vector. The diffusivities are given by Einstein's relation

The Physics of Optical Simulations
The most common method of computation for 3D optical simulations is called the finitedifference time-domain (FDTD) method. This method computes Maxwell's equations of electric and magnetic fields by discretizing the problem in both space and time. Maxwell's equations for non-magnetic materials are given in [26] as where 0 is the permittivity of free space, is relative permittivity, 0 is the permeability of free space, is local charge density, is current density, is the electric field, and is the magnetic field. All optical devices and effects in may be sufficiently simulated with these equations. In an FDTD solver, the derivatives in these equations are replaced with finite differences to facilitate computation.
Due to the relation between a material's refractive index and its permittivity and permeability the shifts in refractive index from the Soref and Bennett model (Equations 2.2 and 2. 3) may be captured as shifts in relative permittivity.

Chapter Conclusions
The fundamental component of PICs is the waveguide. Basic waveguide transmission theory shows that optical waves may be confined to a waveguide through the law of total internal reflection. The confinement of the optical waves in the waveguide gives rise to various modes of oscillation that depend heavily on the waveguide dimensions and carrier wavelength. Some of these modes may be used to carry information, however, the fundamental mode is the most commonly used as it is the most convenient and efficient mode to excite.
There are several different radiation sources that are encountered by space systems. One source being the Van Allen belts that are belts of particles trapped in the magnetic field around the earth. This poses a threat to satellites orbiting the earth as well as systems leaving the earth that must pass through the belts. Solar activity is another radiation source which includes solar flares where the sun may eject large volumes of heavy ions into the space environment of electronic systems. The final primary radiation source is GCRs that are highly energetic rays that come from outside the solar system.
SEEs are the effects in electronic systems caused by the interaction of a single high energy particle with a device in the system. If a heavy ion strikes a semiconductor junction, energy is The methodology starts with step A in Figure 9 with a simulation that emulates the effect of the ion strike on the device in the optical domain by injecting a Gaussian laser light source into the device at the desired strike location and angle. With this simulation, a 3D profile of power absorption per unit volume is created that corresponds to the generation rate of electrons and holes at each point in the device. For each strike location and angle the user wishes to simulate, it will be necessary to perform one of these simulations to generate a unique power absorption profile corresponding to that location and angle.

Figure 9
Flowchart of the proposed simulation methodology that makes it possible to simulate the full time-domain effect of an ion strike in any arbitrary optical device. In A) an optical domain simulation is performed to emulate the effect of an ion striking the device. A power absorption profile is computed, and this is passed to stage B) where an impulse of charges is generated, diffused, and recombined to capture the electronic effects of the ion strike. A set of discrete-time 3D charge profiles are then generated, and these are passed to stage C) where each profile is used to compute the shifts in optical phase and power loss that represent the overall device response to the ion strike. While Figure 9 uses the example of the Lumerical tool suite (FDTD and DEVICE) that this simulation work was performed in, this simulation methodology should be compatible with any simulation tools of a similar feature set. One example of a software substitution that may be compatible with this methodology would be substituting Lumerical DEVICE with Synopsys Sentaurus [28], which may be used in conjunction with Lumerical FDTD. Unlike Lumerical DEVICE, Sentaurus natively includes radiation effects tools that could prove to be beneficial to research performed with this methodology.

Methods
The waveguide evaluated in this paper is a Si slab waveguide shown in Figure 10. A Gaussian beam with a wavelength of 500 nm and diameter of 50 nm is used to emulate the ion striking the waveguide. In Si, a wavelength of 500 nm causes SPA to occur since its photon energy is above the bandgap in Si (1.12 eV) [19]. SPA is often used experimentally to interrogate the detailed spatial dependence on the SEE vulnerability, not possible with broadbeam heavy ion irradiation. In SPA testing, the photon energy is greater than the semiconductor bandgap, thus a photon with a wavelength of less than 1.1 μm may be used to generate a single electron-hole pair (EHP) [14]. In this thesis, the simulated ion strikes are assumed to be striking the waveguide normally from the upper surface (as shown in Figure 12) and alternate strike angles are not considered here.

Figure 10
A cross-sectional slice of the simulated Si slab waveguide. The waveguide is 700 nm wide by 220 nm high. This Si slab is surrounded entirely by SiO2 to create the reflective barrier required for guiding the light waves. The SiO2 is split into separate a separate surface oxide and buried oxide Figure 11 The normalized electric field intensity of the fundamental TE mode of the simulated waveguide shown in Figure 10 Figure 12 The 3D model of the Si waveguide in Lumerical FDTD. The white arrow represents the location and direction of the simulated ion strike using SPA For each simulation (both optical and electrical), mesh constraints were added that generate a finer mesh in the region directly surrounding the ion strike (as shown in Figure 13). This is necessary to attain accurate simulation results and to prevent the divergence of the simulations caused by fast changing gradients in the charge density at the location of the ion strike. In general, having a mesh edge length that is longer than radius of the ion strikes charge generation profile (in this case about 25 nm) is not sufficient to capture the details of the ion strike in the simulation.
For both the optical and charge transport simulations in this thesis, the max mesh length has been set to 5 nm in the region surrounding the ion strike location so that a sufficiently fine mesh may be generated.

Figure 13
The generated mesh surrounding the ion strike location in Lumerical DEVICE. The mesh gets finer in the regions closer to the center of the strike location to prevent fast changing spatial gradients in the charge density from causes the simulation to diverge

Stage A -Emulation of Heavy Ion Strike in Optical Domain
The goal of the first simulation step (step A in Figure 9) in this methodology is to generate a power absorption profile that represents the charge deposition rate of an ion striking the optical device. The mechanism by which this happens is nearly identical in result to the mechanisms of the photoelectric effect, and thus lasers may be injected into devices to emulate and ion strike [19,20]. These tools are leveraged in this methodology to emulate the effect of an ion striking a photonic device. The output of this simulation is a 3D power absorption profile that represents the EHP generation rate at each point in the simulated device.

Device Setup
The device dimensions and materials were constructed in the solver using the software's built in shapes. In Lumerical, when a new shape or simulation tool is added to the simulation environment, it is put into the "Objects Tree". The objects tree that is necessary for this first simulation is shown in Figure 14.

Figure 14
The object tree in Lumerical FDTD with the tools required for the ion strike simulation. The objects named "waveguide", "buried oxide", and "surface oxide" are the shapes used to define the dimensions and materials of the waveguide in this example. "FDTD" is the FDTD solver object that defines the simulation region and mesh and the object called "mesh" is the mesh override object that is defining a finer mesh directly surrounding the ion strike location. Unique to this simulation is the Gaussian light source (called "Strike") and the power absorption monitor (called "Pabs") In Figure 14, the device dimensions and materials are defined by the object's "waveguide", "buried oxide", and "surface oxide". The "waveguide" object represents the waveguide's Si slab (that is the material that the light propagates through represented in blue in Figure 10). The "buried oxide" and "surface oxide" objects represent the SiO2 material encapsulating the waveguide slab (represented in red in Figure 10) that establishes the reflective barrier for the waveguide.
Separating the device's oxide into the "buried oxide" and "surface oxide" is not strictly necessary; however it may make it straightforward to define different materials for the surface and buried oxides if it is desired. The geometries of each these objects are set by right-clicking them and choosing "properties" from the drop-down menu. For this example, the waveguides crosssectional dimensions are set to the dimensions shown in Figure 10.

Simulation Region Setup
An FDTD solver simulation region must also be defined (represented in Figure 14 by "FDTD"). This region is generally made to be slightly smaller than the region of SiO2 represented by the surface oxide and buried oxide to prevent any part of the simulation region from not having any materials defined for it. For this simulation, the waveguide does not have to be set to be particularly long along the length of the waveguide (the y dimension) since this is not a simulation measuring typical light propagation through the waveguide. The Si's high absorption of the 500 nm wavelength light prevents the light from travelling very far along the y direction. Therefore, to save on computation time, the simulation region is defined to only be 5 μm long. The full dimensions have been defined as 1.5 μm x 5 μm x 0.44 μm for this example as shown in Figure   15.
The default meshing scheme is set up within the properties of the FDTD solver object as shown in Figure 16. The default mesh in the x and y directions is set to be defined by a "maximum mesh step" of 20 nm. The mesh in the z direction has been specially defined in this case to be defined by 10 "mesh cells per wavelength". The reason for defining it this way is due to the fact that the ion strike is occurring along the z axis in this example and defining by "mesh cells per wavelength" automatically changes the meshing resolution in the z direction if a different wavelength was chosen to represent the emulated ion strike. The "minimum mesh step" setting has been set to 5 nm which is an order of magnitude smaller than the diameter of the light source emulating the ion strike (50 nm). The remainder of the mesh settings in Figure 16 are set to their default values.

Figure 15
Defined dimensions of the FDTD simulation region for the ion strike simulation. This region encapsulates the waveguide of Figure 10. The simulation region has been defined as 1.5 x 5 x 0.44 μm 3 Figure 16 The mesh settings of the ion strike simulation. The maximum mesh step is set to 20 nm and the minimum mesh step is set to 5 nm. The mesh in the z direction has been set to be defined by 10 "mesh cells per wavelength" instead of by "maximum mesh step" since the ion strike is propagating in the z direction for this example Figure 17 The boundary condition settings in Lumerical FDTD. The boundary condition at each boundary ("x min bc", "x max bx", "y min bc", etc.) is set to perfectly matched layer (PML). The PML settings are set to their defaults At the edges of the simulation region are the simulation boundaries. So that the reaction of light as it hits the simulation boundary may be defined, boundary conditions must be set as in Figure 17. For these simulations, the boundary conditions are all set to perfectly matched layer (PML) as this allows the light to react to the simulation boundary as though the material at the boundary is extended to a theoretical infinity. This allows light to pass through the boundaries where Si is present and allows the simulator to act as though it is simulating a theoretically infinitely long waveguide.

Figure 18
Mesh override settings for the region directly surrounding the ion strike. This defines the maximum mesh step size to be 5 nm along all axes. This guarantees that the generated power absorption profile will have a sufficiently fine resolution A mesh override region is also defined directly around the location of the ion strike in the device. This mesh override region has been set to a maximum mesh size of 5 nm along all axes to guarantee that the mesh is sufficiently fine to capture any details in the region of the ion strike.
The dimensions of this mesh override region span the height of the waveguide along the z direction (0 nm < z < 220 nm) and a 250 nm by 250 nm square centered around the region of the ion strike location as shown in Figure 19.

Figure 19
The location and dimensions of the mesh override region

Gaussian Light Source
A Gaussian light source object has been placed directly on top of the waveguide (a z coordinate of 220 nm as shown in Figure 20) and set to "z-normal" so that it injects light downward along the z axis. This Gaussian source is set to have a "x span" and "y span" of 50 nm to represent an ion strike with a diameter of about 50 nm. In Figure 21 the wavelength of this source has been defined as 500 nm by selecting "override global source settings" and "set frequency/wavelength", and then defining "wavelength start" and "wavelength stop" to both be 500 nm. This ensures that the source wavelength will be as close to a pure 500 nm wave pulse as is possible for the simulator.

Figure 20
The gaussian source geometry settings. It is set to be a "z-normal" source so that it injects light in the z direction and is set to be a diameter of about 50 nm with "x span" and "y span". The z coordinate is set to 220 nm to place it on top of the waveguide Figure 21 The gaussian light source wavelength settings. "Wavelength start" and "wavelength stop" have been set as to the same wavelength of 500 nm to achieve as close as possible to a pure 500 nm light signal

Power Absorption Monitor
A power absorption monitor object is found in Lumerical FDTD's object library as shown in Figure 22. This object is added to the simulation and is given the same physical dimensions as the simulation region in Figure 15. This object contains a prewritten script that computes the total absorbed power per unit volume profile that is necessary for this simulation methodology.

Figure 22
The power absorption monitor can be found in Lumerical FDTD's Object Library under "Optical Power". The non-advanced option is all that is necessary for these simulations

Simulation Execution
Once all the necessary objects and tools have been added to the simulation region, the simulation may be run. When the simulation is complete, the power absorption profile will be stored in the power absorption monitor and may be saved to a MATLAB (.mat) file. A crosssection of the power absorption profile at the center of the strike (y = 0) computed for the waveguide described in Figure 10 with a laser wavelength of 500 nm representing the ion strike is shown in Figure 23. The power absorption computed here is directly proportional to the generation rate of EHPs in the Si waveguide, and as such, the absorption profile may act as a generation source of EHPs in stage B of this methodology.

Figure 23
A cross-section of the waveguide showing the optical power absorbed frome the light source at the point where it is at full intensity. This power absorption profile gets converted to an electronhole pair generation rate profile for use in the Lumerical DEVICE simulation to act as an emulation of an ion strike in the electrical domain simulation

Stage B -Charge Transport Simulation
Following the simulation of the heavy ion strike via analysis of optical absorption, the power absorption profile is imported into Lumerical DEVICE to act as an EHP generation source in the charge transport simulation, step B of Figure 9, and excite a profile of EHPs in the device.
DEVICE then simulates the diffusion and recombination behavior of these free charges in the device to capture the electrical response of the ion strike. The output of this stage is a set of progressively diffusing discrete-time charge profiles that characterize the charge state of the waveguide over a set of discrete-time points. Once this time-dependent charge profile is computed, it may be imported into stage C of this methodology to characterize to optical response of the photonic device to an ion strike.

Figure 24
The Objects Tree in Lumerical DEVICE for the charge transport simulation The device dimensions and material in this simulation are designed identical to the device in the previous simulation (shown in Figure 10) with the exception, in this example, being the waveguide length along the y direction, which is set to about 50 um for this simulation. This allows the generated charges to have some room along the length of the waveguide to diffuse.
Otherwise, the diffusing charges would hit the hard simulation boundaries at the end of the waveguide and have nowhere else to diffuse. If the waveguide is not made long enough, the charges will eventually hit the two ends of the waveguide and prevent the charges from diffusing any furture due to the lack of boundary conditions in the charge transport simulation. In some experimentation, electrical contact boundary conditions have been added to the ends of the waveguide in an attempt to make the waveguide appear to be theoretically infinitely long to the simulator, similar to the PML boundary condition in the FDTD simulator. This however did not have the desired result and instead interfered with the simulation by giving an outlet for the charge to diffuse and recombine faster than expected in a passive waveguide.
As in the ion strike simulation in FDTD, the material of the waveguide slab is set to Si and the material of the material encompassing the slab is SiO2 to keep the device consistent with the previous stage.

Charge Solver Setup
After the device dimensions and materials are specified, a simulation region must be setup.
In Figure 24, the solver called "CHARGE" specifies the dimensions of the simulation region as well as the type of simulation to be carried out and other relevant settings.

Figure 25
The general charge solver settings. This set to a 3D transient simulation at standard room temperature Figure 25 shows the generale charge solver settings for this charge transport simulation.
This simulation was set to a "solver geometry" three-dimensions. The "solver mode" setting was set to transient to allow this simulation to compute the temporal diffusion of charge in the device.
The temperature settings in this thesis are left to their defaults, however, these may be useful in future research where the SEEs in PICs are determined in different temperature environments.
The charge solver's mesh settings are shown in Figure 26. The default mesh edge length has been specified to the solver to be between 0.1 μm and 10 μm. This allows the solver to generate a mesh that contains the mesh edges between those lengths. These values make for a relatively coarse mesh, however a mesh constraint region has been placed directly around the strike location that sets the maximum mesh edge length in a square of 700 nm around the strike location to 5 nm ( Figure 27). Having this mesh constraint region around the ion strike location is of particular importance in this charge transport simulation, as the charge transport algorithms are particularly sensitive to diverging due to fast changing (both spatially and temporally) charge gradients. Such charge gradients are seen in these simulations when generating the impulse of charges that represents the ion strike.

Figure 26
The charge solver mesh settings. The minimum mesh edge length is 0.1 μm and the maximum mesh edge length is set to 10 μm. A long max edge length allows the simulator to generate a mesh that is very coarse in regions far away from the ion strike location to save computational resources. "Max refine steps" has been set to 200000 so that the mesh generator can thoroughly refine the mesh Figure 27 The dimensions and settings for the mesh constraint region directly surrounding the ion strike location. This region prevents simulation divergence by making the mesh finer where the charge gradients are expected to be the steepest Shown in Figure 28(b) are the settings that determine the timing and shape of the shutter used to control the charge generation source. The shutter mode chosen here is "pulse on" as it allows the user to turn the generation source on for a short time to emulate the instantaneous charge generation of an ion strike. The setting of "shutter ton", the difference between "shutter toff" and "shutter ton, and "shutter tslew" must be greater than or equal to the "min time step" setting in Figure 28(a). Failure to follow this rule will throw an error because the time settings of the global source shutter cannot physically be less than the smallest possible simulation time step. "Shutter tslew" determines the time that it takes for the source shutter to go from zero intensity to maximum intensity. "Shutter ton" determines the time at which the shutter begins to allow charges to generate and "shutter toff' is the time the shutter begins to ramp the charge generation down.
"Shutter slew function" determines the shape of the ramp as the shutter turns on and off with choices of exponential and linear. The setting in Figure 28(c) is a global source intensity scaler that multiplies the entire charge generation profile by its value and is useful for LET tuning. These settings as they are shown in Figure 28(b) generate a triangular pulse starting at 1 ps in the simulation, hits peak generation at 2 ps, and ramps back down to no charge generation at 3 ps.

Figure 29
Adding an optical generation source to the simulation region A generation source may be added to the simulation area by the "Import Optical Generation" tool under "Source" as shown in Figure 29. This is the simulation object allows a power absorption profile from the previous simulation ( Figure

Figure 30
Si recombination settings of the charge solver. The various type of recombination in Si may be enabled and disabled here. The default settings of each type of recombination are already set appropriately for Si Accurate charge transport simulation in semiconductors requires that recombination be taken into account. The material settings for semiconductors in Lumerical DEVICE have a tab for controlling each of the recombination mechanisms within the semiconductor (Figure 30). The red box in Figure 30 shows where each recombination mechanism may be enabled or disabled. The default settings within each mechanism should already be set to values appropriate to model recombination within Si, it is simply up to the user to determine whether each mechanism is significant enough to enable for each application. If a given mechanism is determined to have a negligible effect on the simulation, then it may be disabled here to save computational resources.
In this work, the enabled mechanisms have been set to Shockley-Read-Hall (SRH), radiative, and Auger recombination.

Figure 31
Dimensions and placement of the aluminum contact. It's width in the x direction is 600 nm and it is placed directly under the strike location Lumerical DEVICE requires that at least one electrical contact of a metal material be specified within the simulation region. Since this is an example of a passive optical waveguide, that has no inherent electrical contacts, a dummy contact must be created in a location where it will not have an effect on the natural diffusion of charges. An aluminum contact has been placed in the center of the simulation region (x = y = 0) with the dimensions as shown in Figure 31. Notice that the contact has been placed below the waveguide with 800 nm of SiO2 between them. This ensures that the aluminum contact is a sufficient distance away from the waveguide as to not affect the diffusion of the charges throughout the waveguide.

Figure 32
Lumerical DEVICE electrical contact settings. To define a shape in the simulation region as an electrical contact, (a) was set to a solid geometry and (b) was then set to the specific solid name that the user wishes to define as the contact (here it is set to the "contact" as seen in Figure 24). (c1) was set to initial time and voltage of the contact, while (c2) was set to the simulation end time and voltage. Both voltages were set to 0 to ensure that the electrical contact does not affect the passive diffusion simulation in this example

Figure 33
The settings of the charge monitor. "Record electrons", "record holes", "save data", and "integrate total charge" were all checked. A filename was also specified so that the charge profiles recorded by the monitor would be saved The charge profile information from this simulation may be exported from Lumerical DEVICE by manually exporting it from the CHARGE simulation object once the simulation is complete. However, an easier way to export the charge profiles by using a charge monitor was discovered. In the charge monitor settings, Figure 32, there is a setting to automatically save the contents of the charge monitor to a MATLAB (.mat) file once the simulation is complete. The settings "record electrons" and "record holes" were selected. "Save data" was also selected so that the data was automatically saved from the charge solver into the ".mat" file entered into the "filename" box. "Integrate total charge" was also selected to save a graph of total charge in the device overtime for the purpose of estimating the ion strike's LET. The charge monitor has been set to the exact dimensions of the waveguide so that the motion of the charges may be captured in all parts of the waveguide.

Simulation Execution
Once the simulation has been executed, the charge profile results will be stored in the CHARGE solver and the charge monitor. If the charge monitor was also set to save the results automatically, then the charge profiles will be available in a MATLAB file within the same directory that the simulation file resides. An example of a set of charge profiles, generated in one of these simulations, is shown in Figure 34. The color scale of the charge profiles in Figure 34 is renormalized at each time step to better visualize the diffusion of the charge throughout the device.
These charge profiles are used in the following FDTD simulation to compute perturbations in the refractive index of the waveguide and ultimately compute the overall effect on the transmission and phase shift in the waveguide.

Figure 34
An example set of charge profiles computed by Lumerical DEVICE. As time progresses, the charge spreads out through the waveguide and the peak charge density decreases. These discrete-time charge profile steps are used in the following simulation to compute the transmissive parameters of the waveguide for each time step. Note that the color scale is renormalizing for each time step in this image. This helps to better visualize the diffusive nature of the charge

Estimating and Tuning LET
The LET of the ion strike was estimated via analysis of the charge profiles, and tuned by adjusting the generation source intensity setting shown in Figure 28c. One method of estimating the LET is to enable the "integrate total charge" setting in Figure 35. This setting stores an array of results that can be displayed as total electrons and total holes in the waveguide over time. The

Figure 35
The settings of the charge monitor. Both "record electrons" and "record holes" must be selected. "Save data" allows the charge monitor to automatically save the computed charge profile to a file named in the "filename" field. "Integrate total charge" computes the total number of electrons and holes within the charge monitor which is convenient for estimating LET One problem encountered with the "integrate total charge" function in Lumerical DEVICE's charge monitors, is that it is only accessible to charge profiles that were just computed.
If the Lumerical DEVICE software is closed or a new simulation is started, then a given charge profile can no longer be put through DEVICE's charge integrator to compute the total charge. If estimating the LET was forgotten before closing the software, then it could not be done for the charge profile that was just computed. An external charge integration script was developed in MATLAB to fix this problem and is shown in Appendix A.1. The DEVICE generated output charge profile may be imported into this script to compute the total charge separately from the Lumerical software. This script has the added benefit that it not only computes the total charge but includes a line of code to automatically estimate the LET from the given charge profile.

Mesh Analysis
A mesh analysis was performed to ascertain the effect of the mesh size on the results of the DEVICE simulation. This was done by running a set of example simulations that kept all parameters constant except for maximum mesh length of the mesh constraint region in Figure 27.
This parameter sweep reran the simulation for maximum mesh lengths of 2.5 nm, 5 nm, 10 nm, 20 nm, and 40 nm. The LET of each resulting charge profile was estimated with the MATLAB script in Appendix A.1 and plotted in Figure 36.
The graph in Figure 36 appears to converge to an LET of about 97.7 MeV-cm 2 /mg as the maximum mesh length gets smaller. This convergence indicates that the simulation is getting more accurate by decreasing the mesh size surrounding the ion strike location. As mesh size increases, the LET becomes increasingly over-estimated. With a 5 nm maximum mesh length, the estimated LET is very close to converging and the mesh size is not so small cost a lot in processing power.
The maximum mesh length of 5 nm has therefore been deemed sufficient to generate an accurate mesh for these simulations.

Figure 36
Estimated LET versus maximum mesh length for an example charge transport simulation. The estimated LET converges to an LET of about 97.7 MeV-cm 2 /mg as maximum mesh length is reduced, indicating that the simulation becomes more accurate with reduced mesh size. LET seems to be over-estimated at max mesh lengths larger than 5 nm

Stage C -Extraction of Losses, Phase Shift, and Optical S-Parameters
The The Objects Tree in Lumerical FDTD for the s-parameter extraction simulations. The waveguide and mesh is setup identically to the device in Figure 14. There is now an "np density" object for including the carrier profiles in FDTD. The FDTD simulation object now includes input and output ports that facilitate the extraction of device's s-parameters Figure 38 The properties of the "np density" object. This object allows the user to import a set of charge profiles generated by Lumerical DEVICE through the "Import data…" button The objects tree from this simulation stage, and the simulation objects necessary for extracting the s-parameters are shown in Figure 37. The first object of interest that is unique to this stage is the "np density" object which provides a means of importing the charge profiles generated in the previous simulation stage into FDTD. The np density object supports timevarying charge profiles by including the "t_index" setting seen in Figure   The Si material settings for modeling the shifts in optical properties. A new Si material has been created "Si -IP" and has been set to a material type of "index perturbation" and set to inherit the optical properties of the "Si (Si) -Palik" material. To enable Si's refractive index dependency on free carriers, "include np density" is checked and the chosen model is the "Soref and Bennett model" with 1550 nm coefficients The device dimensions here are setup identical to the dimensions of the device in the first stage to keep the simulations consistent. An exception to this is the waveguide length (y-direction) which is set to ~51 μm as it is set to in the second stage. Although the dimensions remain mostly the same, the definition of the Si material must be changed to take the refractive index of Si's dependence on the free carriers in the np density object into account. A new material must be defined in FDTD's material database (Figure 39) that inherits the base characteristics of the Si material but computes the refractive index perturbations based upon free carriers. The material type "index perturbation" has this functionality and includes a model for the Soref and Bennett charge-to-index conversion model.

Simulation Region Setup
The simulation region dimensions and settings ( Figure   The FDTD simulation region dimensions for the parameter extraction simulation. The length of this simulation region (y-direction) is set to 51 μm so that it is slightly longer than the imported charge profile. The other dimensions are identical to the simulation region in Figure 15 making it a 1.5 μm x 51 μm x 0.44 μm simulation region

S-Parameter Extraction Setup
The device's s-parameters cannot be extracted without some form of monitors measuring the device input and output. In Lumerical FDTD, these monitors come in the form of ports. This simulation uses two ports (one input and one output), the dimensions of which are shown in Figure   41. These ports have been made larger than the waveguide cross-section so that they sufficiently span the waveguide's modal profile. The input port is placed at one end of the np density profile (y = -25 μm) and the output port is place at the opposing end (y = 25 μm) to allow the light injected by the input port to come into contact with all induced charges before reaching the output.

Figure 41
S-Parameter extraction input and output port dimensions. These ports are placed at opposing ends of the waveguide and are made slightly smaller than the overall simulation region but still larger than the cross section of the waveguide so that the ports can generate the whole waveguide mode Figure 42 shows the modal properties of the input port that determine the mode injected by the port as well as the direction, amplitude, and phase of the signal. So that the light is injected down the length of the waveguide, the "injection axis" is set to "y-axis". Having the "direction" set to "forward" causes the light to be injected in the positive y direction. The default amplitude of 1 V/m and phase of 0 are fine as this allows the port to output a normalized wave. The port is set to operate under the waveguide's fundamental mode (TE) over one frequency point that represents a 1550 nm wavelength (the waveguide's operating wavelength). The output ports modal properties are identical except that its direction is set to "backward". Setting all ports such that they are facing "into" the device is important for keeping the sign of the s-parameters consistent.

Figure 42
The modal properties of the input port. This port is set to inject forward along the y-axis so that it injects in the positive y direction. The amplitude and phase are set to 1 and 0 respectively to simply act as a normalized electromagnetic wave. The source is operating in the waveguide's fundamental mode over 1 frequency point (1550 nm)

Charge Profile Time Sweep Setup
All of the tools necessary for extracting a single set of s-parameters from the device have been set up so far. However, the goal of this methodology is to generate a series of s-parameters that represent the changing state of a device over time in response to a single event strike. It is necessary to set up an FDTD parameter sweep to accomplish this. The parameter to be swept is

Simulation Execution and Result Processing
Once the parameter sweep has been appropriately setup, the sweep may be executed. This will result in several simulations running consecutively (one for each np density time index specified in the parameter sweep). This will result in two sets of results to be stored in the parameter sweep object once it is complete (one for each port). The s-parameters obtained from the output port in this example represent the forward gain (the gain from the input port to the output port) of the waveguide. The s-parameters obtained from the input port represent the reflection coefficient, which is the amount of light that gets reflected back to the input.
To convert the s-parameters to the transient transmission loss and phase shift, the sparameters in the output port are first normalized with respect to the s-parameters at time t = 0.
Since no charges have been induced into the device at that time, this normalizes the transmission and phase shift with respect to an unperturbed device. The transmission losses are then obtained by = 1 − √ ( ( )) 2 + ( ( )) 2 (3.2) and phase shift by where ( ) is the time-dependent complex-valued set of s-parameters.

Chapter Conclusions
The simulation methodology developed in this work has been described in detail by walking through an example simulation of a slab waveguide with cross-sectional dimensions of 220 nm by 700 nm (Figure 10) in the Lumerical Tool Suite [27]. This methodology contains three stages. In the first stage, the effect of a heavy ion striking the device is emulated in Lumerical pulse width has not been captured in these transients due to technical issues in attempts at repeating these computations for longer (~1 mm) lengths of waveguide. Further work will be required to address these technical issues. However, the results shown here in Figure 44 do show a characterization of the waveform peaks.
The computed transient curves of transmission loss and phase shift in Figure 44 represent the photonic equivalent of the current transients induced in electronic circuits by ion strikes like the one in Figure 4. Figure 44 shows the effects of particles of different LET on the waveguide's transmissive properties over time. These transient curves provide a basic characterization of the overall effect of the ion strike on a signal passing through the device. It is possible that these transients could be used as a means of developing compact models for the analysis of the effect of ion strikes on entire photonic circuits instead of just individual devices.

Figure 44
The transient transmission characteristics of the waveguide in Figure 10 for three LET values computed via the simulation methodology presented in this work. As LET increases, the intensity of the transients becomes more severe, reaching as high as 1% transmission loss and -0.015 radians for an LET of 40 MeV-cm 2 /mg Figure 45 shows the peak values of the transients in Figure 44 plotted versus LET. This appears to show a linear relationship between the peak induced transmission loss and the LET, and the phase shift and LET. Transients induced by the ion strike become more intense as LET increases with the peak transmission loss reaching almost 1% and the peak phase shift reaching about -0.015 radians for an LET of 40 MeV-cm 2 /mg. From these results it seems that, at least for a simple straight waveguide structure and for the limited set of LETs studied here, the transmission loss is the dominantly affected property that could potentially lead to information loss in a photonics system. This is consistent with the findings in Goley et. al. [2].

Figure 45
The peak values of transmission loss and phase shift versus LET. This graph appears to show a somewhat linear relationship between LET, peak transmission loss, and peak phase shift in the waveguide One use of this methodology may be to aide in the discovery of techniques for hardening photonic devices to the effects of ion strikes. Although not analyzed in this work, one potential method of hardening waveguides to radiation may be to design a rib waveguide doped as a diode junction with grounded electrical contacts on each side. Such a device would be equivalent to a photonic modulator like the one in Figure 46 [4] with the two driving electrical contacts grounded.
The electrical contacts would provide an electrical outlet for the charges generated by the ion strike and the diode junction doping provides an electric field that aides in sweeping the generated charges out of the waveguide. Such a device may also be subject to the TID radiation effects of optical modulators described in [4] where the ability for the junction to sweep charges out of the waveguide degrades as holes build up in the oxides around the modulator. Studying the effect of adding the diode junction and electrical contacts on the time constant of transients like the ones in Figure 44 may be the focus of future research, however, this work has focused primarily on fully passive waveguides to understand the baseline effect of ion strikes on these devices.

Figure 46
An example of a silicon photonic modulator. A similar structure may be useful for draining single event induced charges out of a waveguide and therefore reducing the time constant of the effect induced by these charges. [4] CHAPTER V CONCLUSIONS This work has presented a simulation methodology for analyzing the behavior of SETs in PICs and computing shifts in the transmissive properties of a photonic device as they shift over time in response to an ion strike. This methodology shows promise as a utility in characterizing the response of devices to various forms of ionizing radiation as well as discovering methods for improving the radiation hardness of photonic devices. Further work to make this methodology more automated may allow it to become more user-friendly and reduce the learning curve that is required to run these simulations.
This work has used this methodology to generate a set of transient curves for various LETs in a waveguide, showing a proportionality between transmission loss and LET as well as between phase shift and LET. Further work to validate the accuracy of this methodology, including comparison to the physical experimentation of devices under laser or radiation source, will allow this methodology to be tuned if any discrepancies are discovered.

Charge Profile Integration and LET Estimation MATLAB Script
This MATLAB script integrates the total charge in Lumerical DEVICE charge profile and plots the total number of electrons over time. This script assumes that the number of electrons and holes remain approximately the same over time and therefore assumes the maximum amount of electrons is equal to the total number of EHPs. The maximum number of electrons is used to estimate the LET of the ion strike.