title: Learning Compact Physics-Aware Delayed Photocurrent Models Using Dynamic Mode Decomposition
authors: Joshua Hanson Pavel Bochev Biliana Paskaleva
text: 
  
  [
      Shijie Gu
      September 19th, 2021.
  =========================
  
  
  
  
  
  Radiation-induced photocurrent in semiconductor devices can be simulated using complex physics-based models, which are accurate, but computationally expensive. This presents a challenge for implementing device characteristics in high-level circuit simulations where it is computationally infeasible to evaluate detailed models for multiple individual circuit elements. In this work we demonstrate a procedure for learning compact delayed photocurrent models that are efficient enough to implement in large-scale circuit simulations, but remain faithful to the underlying physics. Our approach utilizes Dynamic Mode Decomposition (DMD), a system identification technique for learning reduced order discrete-time dynamical systems from time series data based on singular value decomposition. 
  
  
  To obtain physics-aware device models, we simulate the excess carrier density induced by radiation pulses by solving numerically the Ambipolar Diffusion Equation, then use the simulated internal state as training data for the DMD algorithm. Our results show that the significantly reduced order delayed photocurrent models obtained via this method accurately approximate the dynamics of the internal excess carrier density – which can be used to calculate the induced current at the device boundaries – while remaining compact enough to incorporate into larger circuit simulations.
  
  
  
  Dynamic Mode Decomposition, Ambipolar Diffusion Equation, delayed photocurrent, machine learning, data-driven compact models.
  
  
  
  
  
  § INTRODUCTION
  
  
  Ionizing radiation can affect operation of electronics in multiple application contexts including space, terrestrial and manmade nuclear environments. Specifically, fluence of ionizing radiation generates electron-hole pairs within semiconductor devices, resulting in excess currents flowing through the devices. Such excess current is not present during device operation in a normal environment and its impact on electronic systems can be catastrophic, ranging from instantaneous interruptions in service, loss of stored memory, and even burnout of the entire system. For example spacecraft depend on electronic components that must perform reliably over missions measured in years and decades and space radiation is a primary source of degradation, reliability issues, and potentially failure for these electronic components. Physics-based modeling and simulation of radiation effects on electronic systems in various radiation environments can facilitate understanding of the mechanisms governing the radiation response of the electronic materials, parts, and systems, and can be used to devise ways to mitigate radiation effects, and create new materials and devices that are resilient to radiation exposure. 
  
  Thus, computational analysis of radiation effects on electronic systems has utility ranging from guiding the initial designs of systems, setting up the design of experiments, and final qualification. At a device level the excess carrier behavior can be accurately modeled by the Drift-Diffusion equations (DDE)<cit.> given by
  
      ∇· (ϵ𝐄)    = q (p - n + C) 
      ∂ n/∂ t   = ∇· (n μ_n 𝐄 + D_n ∇ n) - R + g 
      ∂ p/∂ t   = ∇· (p μ_p 𝐄 + D_p ∇ p) - R + g
  
  where n and p are the concentrations of the electrons and holes, respectively, μ_n, μ_p, D_p, and D_n are the carrier mobilities and diffusivities, C is the doping concentration, 𝐄 is the electric field, q is the electron charge, R is the recombination rate, and g is the generation rate.
  
  Numerical solution of (<ref>)–(<ref>) using, e.g., a multi-dimensional finite element discretization of the device, forms the basis of the so-called Technology Computer-Aided Design (TCAD) simulators. However, (<ref>)–(<ref>) is a coupled system of nonlinear Partial Differential Equations (PDEs) whose accurate numerical solution can be very time consuming. As a result, although TCAD device simulators based on DDE could in principle be coupled to circuit simulators, their high computational cost makes the analysis of all but very small circuits computationally intractable. 
  
  At the other end of the spectrum are the so-called compact device models, which are computationally efficient but rely on empirical approximations and/or simplified analytic solutions to the semiconductor transport equations. Often such models must be recalibrated for different operating regimes and inclusion of new physics may force redevelopment of the model from scratch. 
  
  Data-driven techniques present an opportunity to automate the development of compact models by learning them directly from laboratory measurements and/or suitable synthetic data. For such models to be practically useful though, they must be able to correctly predict the device response when the device is integrated into a circuit and exposed to a wide range of stimuli. The caveat is that in a laboratory settings one can only apply a limited type of signals and directly measure the device response, leading to “sparse” training sets. In contrast, traditional Machine Learning (ML) applications, such as natural language processing and image classification 
  <cit.> operate in “big data” environments. As a result, a compact device model learned solely from such “sparse” data may fail to generalize to all relevant analysis conditions because it will learn salient patterns in the dataset rather than the causal physics underpinning the device operation. 
  
  The latter, i.e., the fact that device behavior is governed by strict physics laws can be used to counter the lack of big data by incorporating the physics knowledge into the model development. Recent work on scientific Machine Learning <cit.> suggests this strategy is effective and can lead to generalizable models. 
  
  
  The main goal of this work is to establish the viability of physics-aware machine learning for the development of compact data-driven photocurrent models that are efficient enough to allow large circuit simulations, while remaining faithful to the basic physics principles embodied by models such as (<ref>)–(<ref>). 
  
  To that end we adopt a setting that has been used extensively in the past five decades at Sandia National Laboratories for the development of compact analytic photocurrent models. Although this setting uses a simplified version of the DDE equations, it provides a stepping stone towards future developments of data-driven models based on the fully coupled system (<ref>)–(<ref>). In so doing we are able to leverage a wealth of experiences and physics knowledge accumulated through the use of the existing compact models, as well as provide a reference point for evaluation of our data-driven models. 
  
  We have organized the paper as follows. Section <ref> reviews the current state-of-the art in compact photocurrent models and establishes the physics basis for our data-driven model. Section <ref> briefly summarizes the finite element discretization of the physics model used to generate synthetic training data. The core of the paper is Section <ref> where we use Dynamic Mode Decomposition ideas <cit.> to develop our model. Section <ref> contains computational results highlighting the performance of the model as well as comparison with published results in <cit.>. We summarize our findings and discuss future research in Section <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  § PHYSICS-BASED COMPACT ANALYTIC PHOTOCURRENT MODELS
  
  
  
  
  
  Although the DDE system (<ref>)-(<ref>) accurately describes the behavior of the excess carriers, it does not lend itself to exact analytical solution. To enable analytic approximation of the governing equations, most photocurrent models in use today follow the same basic approach as in the classic paper <cit.> and split each PN-junction within a device into a depletion region and quasi-neutral P- and N-regions; see Fig. <ref>.
  
  Carriers in the depletion region are quickly converted to photocurrent and yield the so-called prompt photocurrent I_prompt. Carriers in the P and N regions have a delayed response and produce the delayed photocurrents I_p and I_n, respectively. As a result, the total photocurrent is given by
  
      I_total = I_prompt + I_p + I_n  .
  
  To calculate these currents one makes additional simplifying assumptions. The first is the electrical neutrality approximation, which stipulates that the excess electron and hole densities are equal throughout the entire device. The second is the congruence assumption, which states that the electron and hole fluxes into or out of any region must be equal; see, e.g., <cit.>. Under these assumptions the DDE model (<ref>)-(<ref>) in the P and N regions can be replaced by the Ambipolar Diffusion Equation (ADE)
  
  
      ∂ u/∂ t = D_a ∇^2 u - μ_a 𝐄·∇ u - 1/τ_a u + g
  
  
  where u, D_a, μ_a and τ_a are the excess carrier density (electrons or holes), the ambipolar diffusion coefficient, the ambipolar mobility and the carrier lifetime, respectively. In general, these parameters may depend on the excess carrier density and, as a result, (<ref>) is still a nonlinear PDE. However, for moderate radiation dose rates these coefficients can be approximated by constant values and the ADE becomes a linear parabolic PDE. The final assumption is that the depletion region width W is not affected by the excess carriers and is a constant. 
  
  Under these assumptions, the prompt photocurrent is modeled as 
  
      I_prompt = qgAW ,
   
  where A is the effective area of the PN-junction, while I_p and I_n are modeled by the ADE (<ref>). 
  
  At this point traditional compact photocurrent model development proceeds with deriving analytic approximations for the solutions of (<ref>) in one dimension and using them to obtain expressions for the delayed photocurrents. Early work <cit.> considered unbounded P and N regions and negligible electric fields. The resulting Wirth-Rogers model tends to overestimate the photocurrent as it neglects the effects of an ohmic contact at a finite distance from the depletion region <cit.>. Subsequent work <cit.> relaxed these conditions, assumed that 𝐄 is constant, and used approximate Laplace transforms to obtain analytic expressions for I_p and I_n. However, approximation of the Laplace transform results in a model that yields unphysical current estimates when 𝐄 exceeds roughly 10; see <cit.>. The latter work used Fourier analysis techniques to develop a more accurate  photocurrent model that avoids these drawbacks. The model in <cit.> was further improved in <cit.> by using a transformation of (<ref>) into an inhomogeneous heat equation and solving the latter exactly by Fourier techniques. 
  Another popular photocurrent model is the Fjeldly Model <cit.>. Unlike the Axness-Kerr model <cit.>, which solves the time-dependent equation (<ref>) exactly, the Fjeldly Model uses a steady-state solution of ADE combined with an RC delay circuit to achieve time dependence. 
  
  In this work we adopt the above setting and focus on the development of physics-aware compact data-driven models for the delayed photocurrents in the P- and N-regions that can be used as "plug-and-play" substitutes for conventional models. The core idea is to replace the simplified analytic solutions of the ADE, comprising the basis of most standard compact models, by a discrete-time dynamical system approximating the flow map of the ADE, i.e. the “solution operator” that returns the internal device state for any given external input. 
  
  To that end we apply a Dynamic Mode Decomposition (DMD) <cit.> approach to samples of the internal state (carrier density) of the device obtained by solving (<ref>) numerically by a Finite Element Method (FEM). In this paper we assume that all device parameters such as diffusion coefficients, carrier lifetimes and doping concentrations are known. In this case ADE already contains all the necessary physics information and development of physics-aware compact data driven models can be done entirely from synthetic data[In this setting the resulting compact model can be interpreted as a discrete reduced order model of the flow map generated by the ADE.].
  
  
  A more general setting occurs when one or more of the ADE parameters are either unknown or have large uncertainties. In this case the development of the compact model must also include a parameter identification step to refine the ADE, which requires laboratory data. Since the purpose of this paper is to demonstrate the viability of DMD as an effective tool for the generation of compact device models, detailed discussion of this more general setting is beyond the scope of this paper and will be addressed in a forthcoming work.
  
  
  
  
  
  
  Regardless of the particular setting though the physics-awareness of our models stems from the fact that they represent approximations of the flow map engendered by the physics model, and built from simulated internal states of this model. 
  
  These states contain physics information that cannot be obtained by laboratory instruments, which can typically only measure the currents at the device terminals. The latter may not be "rich enough" for a traditional ML approach, as well as for DMD, to obtain a reliable model of the underlying causal physics.
  
  
  
  
  
  
  § NUMERICAL SOLUTION OF THE ADE
   
  
  We consider the ADE on the space-time domain Ω := 𝒳×𝒯⊂ℝ^2, where 𝒳 := (0,ℓ) and 𝒯 := (0,t_final). Without loss of generality, one may assume that 𝒳 is the N-region of the device; see Fig. <ref> and so we set ℓ=ℓ_n. 
  
  To obtain a well-posed problem, we augment equation (<ref>) with homogeneous initial and boundary conditions, i.e., 
  
      u(x,0) = 0  ∀ x∈𝒳    
      u(0,t) = u(ℓ,t) = 0  ∀ t ∈𝒯 .
  
  The homogenous initial condition corresponds to the fact that at t=0 there are no excess carriers present in the device. The boundary condition choice corresponds to assuming infinite carrier recombination velocities and ohmic contacts at x=0 and the boundary of the depletion region x=ℓ; see <cit.>. More general boundary conditions not requiring these assumptions can also be considered but are not necessary for the purpose of this work. 
  
  We will simulate the internal state of the device by using the method of lines to solve the ADE numerically. For the spatial discretization we consider a standard Galerkin finite element method and then solve the resulting system of Ordinary Differential Equations (ODEs) using an implicit numerical integration scheme. For completeness we briefly review the discretization process below. 
  
  Let 𝒳^h denote a uniform[Utilizing a variable mesh spacing with increased node density near the boundary may provide some potential advantages such as more accurate gradient estimation at the edge points, however in the context of this work we restrict our attention to uniform meshes.] partition of 𝒳 into n+1 elements κ_i with vertices {x_i}_i=0^n+1, i.e., κ_i = [x_i,x_i+1], i=0,…,n, x_0 = 0 and x_n+1 = ℓ. The mesh parameter is given by h=ℓ/n+1. 
  
  
  
  
  
  As usual, L^2(𝒳) denotes the space of all square integrable functions on 𝒳 with norm and inner product denoted by ·_0 and (·,·)_0, respectively, and H^1_0(𝒳) is the Sobolev space of order one whose elements are constrained to vanish at the boundary points. The weak variational form of the ADE is then given by seek u∈ H^1_0(𝒳) such that
  
  
      (u_t,v)_0 + Q(u,v) = (g,v)_0  ∀ v ∈ H^1_0(𝒳) .
  
   
   The bilinear form Q(·,·):H^1_0(𝒳)× H^1_0(𝒳) ↦ℝ is defined as 
   
   
      Q(u,v) = D_a(u_x,v_x)_0 + μ_a 𝐄(u_x,v)_0 + 1/τ_a(u,v)_0  .
  
  
  To discretize (<ref>) in space we consider a nodal (Lagrangian) conforming finite element subspace V^h_0⊂ H^1_0(𝒳); see, e.g., <cit.>. Let {v_i}_i=1,n be the standard nodal basis having the property that v_i(x_j) = δ_ij. We then seek an approximate solution of (<ref>) as
  
  
      u_h(x,t) = ∑_i=1^n u_i(t) v_i(x)  ,
  
  
  where 𝐮(t) = (u_1(t),…,u_n(t)) ∈ℝ^n is a vector of unknown solution coefficients. 
  Inserting (<ref>) into the weak form (<ref>) and restricting the test space to V^h_0 then yields the system of ODEs
  
  
      M 𝐮̇(t) + K𝐮(t) = 𝐠
  
  
  where M,K∈ℝ^n× n are the (consistent) finite element mass and stiffness matrices with elements
  
      M_ij = (v_i,v_j)_0
          K_ij = Q(v_i,v_j) ,
  
  respectively, and 𝐠(t)∈ℝ^n is a discrete source term with g_i(t) = (g,v_i)_0. 
  
  The ODE system (<ref>) can be solved by any standard time-integration scheme. However, in general (<ref>) is stiff and an implicit scheme is preferred. In this paper, we use an implicit multi-step variable-order routine based on a backward differentitation formula (BDF) for approximating the state derivative; for more details see <cit.>. This method is included by default in the  submodule within the SciPy v1.5.1 package for Python 3.
  
  
  
  
  
  
  
  
  
  
  § A DYNAMIC MODE DECOMPOSITION COMPACT PHOTOCURRENT MODEL
   
  
  
  Assume that an approximate solution of the ODE (<ref>) is available at uniformly spaced time steps t_k:=kΔ t, k=0,…,m, for some sampling interval Δ t > 0. Then, the approximate numerical solution of the ADE can be represented as a linear discrete-time dynamical system acting on samples of the state 𝐮 and the input 𝐠, i.e.,
  
      𝐮_k+1 = A 𝐮_k + B 𝐠_k
  
  where
  
      𝐮_k   := [u_1(k Δ t) ⋯ u_n(k Δ t)]^𝖳
      𝐠_k   := g(k Δ t)
  
  and A, B ∈ℝ^n × n are linear maps. Expressing the system in this form facilitates the use of many familiar system identification techniques. In particular, Dynamic Mode Decomposition (DMD) is a data-driven method for learning the maps A and B from a time series of state and input measurements {𝐮_k, 𝐠_k}_k=0^m, with the goal of identifying a small number of dominant dynamic modes. Below we summarize the DMD algorithm adapted for control inputs as described in <cit.>. By organizing the samples into the following matrices
  
      X  = [     |           |;   𝐮_0     ⋯ 𝐮_m-1;     |           | ] ;   
          X'  = [   |       |; 𝐮_1   ⋯ 𝐮_m;   |       | ] ;    
          G = [     |           |;   𝐠_0     ⋯ 𝐠_m-1;     |           | ],
  
  we can express the linear relationships within the data as
  
      X' = AX + BG = [ A B ][ X; G ] =: [ A B ] S.
  
  Therefore the maps A and B can be approximated by
  
      [ A B ]≈[ A̅ B̅ ] := X' S^†
  
  where † indicates the Moore-Penrose pseudoinverse. An efficient and accurate algorithm for estimating the pseudoinverse of a rectangular matrix is realized via truncated singular value decomposition. The matrix of samples S can be factored as
  
      S = U Σ V^𝖳 = [      Ũ Ũ_trun ][      Σ̃       0;       0 Σ̃_trun ][      Ṽ^𝖳; Ṽ_trun^𝖳 ]≈ŨΣ̃Ṽ^𝖳
  
  where U ∈ℝ^n × n, Ũ∈ℝ^n × p, Σ∈ℝ^n × m, Σ̃∈ℝ^p × p, V^𝖳∈ℝ^m × m, and Ṽ^𝖳∈ℝ^p × m. Here _trun denotes the m-p truncated singular values, so that the p greatest singular values are kept. In practice, one sets an error tolerance and truncates all singular values below this threshold. Now the pseudoinverse of S is naturally approximated by
  
      S^† = ṼΣ̃^-1Ũ^𝖳
  
  Substituting equation (<ref>) into (<ref>) gives
  
      A    ≈A̅≈ X' ṼΣ̃^-1Ũ_̃1̃^𝖳∈ℝ^n × n
      
          B    ≈B̅≈ X' ṼΣ̃^-1Ũ_̃2̃^𝖳∈ℝ^n × n,
  
  where Ũ = [ Ũ_1^𝖳 Ũ_2^𝖳 ]^𝖳 with Ũ_1 ∈ℝ^n × p, Ũ_2 ∈ℝ^n × p.
  
  We can achieve a more compact model by incorporating the projection 𝐮̃ = P𝐮 of the state onto the canonical dynamic mode coordinates. In the same manner as equation (<ref>), we factor the matrix X' ≈ÛΣ̂V̂^𝖳, where the truncations are chosen to preserve the r greatest singular values. Then the projection onto dynamic mode coordinates is given simply by P = Û^𝖳∈ℝ^r × n, hence the desired reduced order model obtained via DMD is
  
      𝐮̃_k+1   = Ã𝐮̃_k + B̃𝐠_k
      𝐮̃_0   = Û^𝖳𝐮_0
      𝐮_k   ≈Û𝐮̃_k
  
  where
  
      Ã    := Û^𝖳 A̅ Û        = Û^𝖳 X' Ṽ Σ̃^-1 Ũ_̃1̃^𝖳 Û        ∈ℝ^r ×r 
  
      B̃    := Û^𝖳 B̅        = Û^𝖳 X' Ṽ Σ̃^-1 Ũ_̃2̃^𝖳        ∈ℝ^r ×n.
  
  
  The parameter p represents the number of dynamic modes to fit to the data, which controls the model precision. The parameter r represents the number of modes to project onto, that is, the order of the final reduced-order model, which controls the model compactness. The case where r > p usually results in diminished performance; r = p retains exactly the same number of modes fit to the data in the compactified model; r < p results in a more compact model, but ignores the p - r least significant modes fit to the data, which can result in slightly reduced accuracy. In this work we will use r = p.
  
  
  
  
   §.§ Training of the compact model
  
  
  Construction of the DMD model (<ref>) requires training samples representing time series of state and input measurements {𝐮_k, 𝐠_k}_k=0^m. To generate such samples we recall the assumption in Section <ref> that all device parameters are known. In particular, here we consider a generic PN-junction device characterized by the parameters in Table <ref>. The values in this table are adapted from <cit.> to enable a direct comparison with a published compact photocurrent model.
  
  
  
  
  
  We then select a suitable set of generation density functions {g^k_train}_k=1^M and use the computational scheme in Section <ref> to solve (<ref>) numerically with homogeneous initial and Dirichlet boundary conditions. Selection of the inputs g^k_train depends on the type of the anticipated testing input(s) g_test for the model and will be revisited in Section <ref>. 
  
  
  Obtaining a numerical solution for the ADE requires proper scaling and non-dimensionalization of the governing equations. For convenience we scale the computational domain 𝒳 so that the length of the N-region, i.e., ℓ, becomes a unit of length and t_final becomes a unit of time. After rescaling the domain and the equations the computational domain 𝒳 becomes the unit square, and the non-dimensional ADE coefficients are given by
  
  
      [ D_a = 0.063; μ_a = 0.063 ]
      ;  [ L_p    = 0.5;   τ_a = 3.92;              ];    [ 𝐄 = -23.84   0;       g = 2.15 ] .
  
  
  
  We highlight that the numerical solution of the ADE constitutes the "physics-based" element of our procedure. Specifically, the physics information is incorporated by using trusted a priori dynamics models – which are calibrated or driven by experimental measurements – to generate training data from simulating the unobservable internal state of a device using robust numerical techniques. 
  
  
  In this work we use the ADE as the physics basis for the data-driven model because it has been used to develop almost all compact photocurrent models in use today, i.e., it is an example of model that is trusted based on decades-long practical experiences. 
  
  However, we emphasize that the DMD algorithm is also suitable for more complex physics-based models, such as the full drift-diffusion equations or detailed molecular dynamics simulations, with the capability of producing dramatically reduced order approximations that are feasible to implement in high-level circuit simulators but remain faithful to the underlying physics.
  
  
  
  
  
  
  § SIMULATION RESULTS
  
  
  
  Assuming that a DMD model (<ref>) has been trained according to the procedure in Section <ref>, we test its performance as follows. 
  
  Let g(x,t) be a target generation density for which we seek the response of our device. We sample g(x,t) in space using the vertices defining the finite element mesh 𝒳^h, and in time using a desired time step Δ t for a total of m time steps.
  
  This sampling produces the inputs 𝐠_k to the DMD model. We then set 𝐮_0=0 and use (<ref>) to recover the internal state of the device, i.e., the excess carrier concentration, at the mesh nodes {x_i}_i=1^n for every t_k = kΔ t:
  
      𝐮̃_k+1 = Ã𝐮̃_k + B̃𝐠_k ,
         k=0,…,m-1 .
  
  Each vector 𝐮̃_k induces a C^0 finite element function 
  
      u_h^DMD(x,t_k) = ∑_i=1^nũ_i(t_k) v_i(x)
  
  which is the predicted internal carrier density. Using u_h^DMD and taking into account the homogeneous Dirichlet boundary condition u_h^DMD(x_0,t_k)=u_h^DMD(x_n+1,t_k) = 0, we define the approximate boundary photocurrent at t=t_k as 
  
      J_h^DMD(x,t_k) = D_a ∂_x u_h^DMD(x,t_k) = D_a ∑_i=1^nũ_i(t_k) ∂_x v_i(x),
  
  
  for x=x_0 and x=x_n+1. Note that owing to the local support of the basis functions v_i(x), J_h^DMD(x_0,t_k) and J_h^DMD(x_0,t_k) only include contributions from the basis functions supported on elements κ_0 and κ_n, respectively. 
  
  We then compare the predicted DMD flux with simulated experimental measurements of the current out of the device terminals. These measurements are obtained using the same numerical procedure as in Section <ref>, i.e., by computing a finite element solution u_h^FEM(x,t_k) and setting
  
  
      J_h^FEM(x_0,t_k)    = D_a ∂_x u_h^FEM(x_0,t_k) 
      
      	J_h^FEM(x_n+1,t_k)    = D_a ∂_x u_h^FEM(x_n+1,t_k).
  
  
  
  
  
   §.§ Manufactured solution test
  
  
  
  We first evaluate the ability of the DMD model to reproduce an artificially manufactured solution:
  
      u^MNF(x,t) := t e^-2tsin(3 π x/ℓ).
  
  Observe that this solution satisfies homogenous initial and boundary conditions. Since (<ref>) is expressed in closed-form, we can directly compute the boundary photocurrent of the manufactured solution:
  
      J^MNF(x,t) = D_a ∂_x u_h^MNF(x,t) = D_a t e^-2t(3 π/ℓ) cos(3 π x/ℓ).
  
  
  Thus, for this test we will compare J^DMD with the known manufactured solution current J^MNF, rather than with the simulated current J^FEM in (<ref>)–(<ref>).
  
  Substitution of u^MNF(x,t) into the governing equation (<ref>) yields a generation density
  
      g^MNF(x,t)    = ∂ u^MNF/∂ t(x,t) - D_a ∂^2 u^MNF/∂ x^2(x,t) + μ_a 𝐄(x) ∂ u^MNF/∂ x(x,t) + 1/τ_a u^MNF(x,t) 
         = (1 - 2t) e^-2tsin(3 π x/ℓ) + D_a t e^-2t(3 π/ℓ)^2 sin(3 π x/ℓ) 
            + μ_a E(x) t e^-2t(3 π/ℓ) cos(3 π x/ℓ) + 1/τ_a t e^-2tsin(3 π x/ℓ)
  
  such that when the ADE (<ref>) is driven by g^MNF(x,t), its solution will exactly match the desired manufactured solution u^MNF(x,t). The inputs 𝐠_k to the DMD model are obtained by sampling g^MNF(x,t) according to the method described earlier.
  
  Since the input g^MNF(x,t) corresponding to the manufactured solution is spatially irregular, a spatially uniform training input will generally result in poor performance. To address this, we design a sequence of localized pulses which will excite different regions of the device. This sequence will then be used as the training input for the DMD model. We choose to use a Gaussian profile that has been windowed by a cosine function as the spatial envelope for the input pulses, where the window function is applied to restrict the support of the envelope to a compact interval. This profile is consistent with experimentally viable radiation doses; other reasonable choices include Lorenz or Voigt profiles, which reflect different radiation broadening mechanisms. The windowed Gaussian profile with center x_i and support [x_i-w/2,x_i + w/2] is given by
  
      ρ_i(x) := cos(πx - x_i/w) exp(-16(x - x_i/w)^2)    if x_i - w/2≤ x ≤ x_i + w/2
      
      	0    if x < x_i - w/2 or x > x_i + w/2
  
  where x_i = iℓ/N_pulses - 1 for i = 0,…,N _pulses-1, and w = ℓ/N_pulses - 2. The profiles ρ_i(x) are illustrated in Figure <ref>. Combining the spatial envelopes defined above with a square temporal envelope gives
  
      g_i(x,t) =
      gρ_i(x)     0 ≤ t ≤0.5
      
      0          ,
  
  where the value of g is defined in (<ref>).
  
  Now we can express the training input as
  
      g_train = ∑_i=0^N_pulses - 1 g_i(x,t - t_i)
  
  where t_i = i(5.0) for i = 0,…,N_pulses-1. For the manufactured solution (<ref>), we choose N_pulses = 10.
  
  To generate the training samples, we solve (<ref>) with the source (<ref>) and the parameters (<ref>) on a mesh 𝒳^h comprising 512 uniform elements, and sample the solution in time at Δ t = 0.005 increments.
  
  
  
  Using the the manufactured photocurrent (<ref>) and carrier density (<ref>) and approximate photocurrent (<ref>) and approximate carrier density (<ref>) derived from the DMD model, we define the following error quantities:
  
      E_J(t)    := | J^DMD_h(t) - J^MNF(t) | 
      
          E_u(t)    := | u^DMD_h(t) - u^MNF(t) |
  
  which will characterize the ability of the DMD model to reproduce manually selected dynamic modes. For simplicity, in this work the electric field is not considered an input to the system, so the DMD model must be trained separately for each unique electric field strength.
  
  We report results for the DMD model using several different values for the parameter p, which defines the number of the greatest singular values kept in the truncated SVD decomposition (<ref>).
  
  Figure <ref> compares the boundary photocurrent from the manufactured solution and from the FEM and DMD models subject to the manufactured input (<ref>). Figures <ref> and <ref> show the carrier density from the manufactured solution and from the FEM and DMD models, with respect to time t or position x, respectively. Figures <ref>, <ref>, and <ref> show the same, but in the absence of an electric field.
  
  
  
  
  
  
  
  
  
  
  
   §.§ Verification test
  
  
  
  
  
  In this test we compare our compact DMD model with the Axness-Kerr <cit.> compact analytic model. We consider the case of a lightly doped diode as described in <cit.> for which the length of the N-region equals 2 diffusion lengths (case ξ_p = 2 in <cit.>.) 
  
  The device is irradiated by a 1.0 step pulse
  
  
      g_test(x,t) =
      g    0 ≤ t ≤1.0
      
      0
  
  
  with the value of g rescaled as in (<ref>). This example from <cit.> corresponds to the parameters in <ref> with the exception of a few corrections to account for typographical errors in that paper. 
  
  In contrast to the manufactured solution test where the desired input g^MNF(x,t) to the model is spatially irregular, now the target generation density, defined in (<ref>), is spatially constant. As a result, the training input does not have to excite different regions of the device and can be chosen to be spatially constant as well. Thus, for this example we choose the training input g_train(x,t) to be a constant in space and a discontinuous piecewise linear in time function such that
  
  
  
  
  
  
  
  
  
  
  
  
  
      g_train(x,t) = 
      
      0                4t    < 0     2 ≤ 4t    < 5    4t ≥7
      g(4t)       0 ≤ 4t    < 1
      g    1 ≤ 4t    < 2   5 ≤ 4t    < 6
      g(7-4t)     6 ≤ 4t    < 7
      ,
  
  
  See Figure <ref> for an illustration of g_train and g_test. 
  
  We then solve (<ref>) with the source (<ref>) and the parameters (<ref>) on a mesh 𝒳^h comprising 1024 uniform elements, and sample the solution in time at Δ t = 0.0025 increments. The compact DMD photocurrent model for the device is now defined according to (<ref>).
  
  
  Since the input is applied uniformly across the entire device (and thus the entire state space), we could reduce the input dimension of the DMD model (<ref>) from N to 1, however the higher dimensional input permits spatially irregular excitations such as localized radiation pulses or non-transversal plane waves as required for the manufactured solution test in Section <ref>.
  
  
  To determine an appropriate dimension for the reduced-order DMD model, it is informative to inspect the relative magnitudes of the singular values from the decompositions of the sample matrices.
  
  For the training input (<ref>), Figure <ref> illustrates the magnitude roll-off in the singular values for the state and state-input sample matrix decompositions.
  
  Observe that even for tight error thresholds, only a few modes are necessary to construct an accurate approximation of the state transition and input matrices.
  
  We wish to emphasize that including too many modes in the reduced-order model often leads to a realization which is unstable over a long time horizon.
  
  This phenomenon occurs due to unstable modes corresponding to small singular values. 
  
  Including these low-magnitude modes will typically lead to a better fit for the state dynamics resulting from the training input over the original training time horizon, but may lead to model divergence over a longer time horizon, signifying an example of overfitting.
  
  
  
  Based on the photocurrent (<ref>)-(<ref>) and carrier density (<ref>) derived from the finite-elements model and the approximate photocurrent (<ref>) and approximate carrier density (<ref>) derived from the DMD model, we redefine the following error quantities:
  
      E_J(t)    := | J^DMD_h(t) - J_h^FEM(t) | 
      
          E_u(t)    := | u^DMD_h(t) - u_h^FEM(t) |
  
  which will facilitate the performance evaluation of the DMD model for typical, spatially uniform input functions. As before, the DMD model must be trained separately for each electric field strength.
  
  Again, we show the results for the DMD model using multiple values for the parameter p. However, observe that for the case where the generation density is spatially uniform, fewer modes are necessary to achieve good performance than when the input magnitude varies along the length of the device, as in the manufactured solution test.
  
  Figures <ref> and <ref> compare the boundary photocurrent produced by the FEM and DMD models for the training input (<ref>) and test input (<ref>), respectively. Figures <ref> and <ref> show the same, but in the absence of an electric field. Figures <ref> and <ref> illustrate the simulated excess carrier density (i.e., the internal state) from the FEM and DMD models for both the training and test inputs. 
  
  Several conclusions can be drawn from these results. First, the photocurrent plots at the left boundary of the N-region shown in Figures <ref> and <ref> are in an excellent agreement with the results reported in <cit.> for ξ_p=2.
  
  Second, the error plots on the bottom rows of the figures quantify the differences between the reference FEM solution and its DMD approximation as a function of the number p of selected dynamical modes. These results reveal that, as expected, the error decreases with increase of the number of dynamic modes; however, even with just 6 modes selected the DMD photocurrent model yields excellent accuracy.
  
  Overall, these results suggest that a data-driven approach is indeed a viable and effective alternative to traditional analytic model development that can be used to quickly develop accurate and computationally efficient photocurrent models directly from data.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § CONCLUSIONS
  
  
  We have developed a compact data-driven delayed photocurrent model given by a low-dimensional discrete-time dynamical system, which approximates the flow map of the Ambipolar Diffusion Equation. 
  
  To obtain the approximate flow map we use the Ambipolar Diffusion Equation to reconstruct numerically the internal state of the device, which is not directly observable through a laboratory measurement, and then apply Dynamic Mode Decomposition to the simulated internal state samples.
  
  In doing so physics knowledge is incorporated into the model development, which allows us to obtain models from sparse data sets that accurately approximate the dynamics of the excess carrier density.
  
  This in turn allows us to accurately estimate the induced current at the device boundaries, which is the quantity required for circuit simulations.
  
  
  Our results confirm that such physics-aware data-driven models are a viable alternative to traditional analytic compact models that use simplified analytic solutions of the governing equations and often must undergo recalibration and/or redevelopment to include new physics effects. 
  
  Our future work will consider extension of the approach to include an additional parameter identification step, and to the fully coupled DDE system (<ref>)-(<ref>). The latter will allow us to model the total photocurrent in the device and eliminate the need to split it into three separate regions.
  
  
  Applying nonlinear observable functions to the state of the DDE would allow us to model the nonlinear problem using the same DMD algorithm, requiring little to no additional computational cost for running and training the DMD model (besides applying the nonlinearities to the measurement data, which is inexpensive). We also plan to incorporate and test our models in circuit simulators to demonstrate their utility for circuit design and analysis tasks.
  
  
  
  
  
  
  
  
  § ACKNOWLEDGMENTS
  
  
  This work was supported by the Sandia National Laboratories (SNL) Laboratory-directed Research and Development (LDRD) program, and the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under Award Number DE-SC-0000230927 and under the Collaboratory on Mathematics and Physics-Informed Learning Machines for Multiscale and Multiphysics Problems (PhILMs) project. 
  
  The work of J. Hanson was performed as part of Sandia Summer Student Program Internship at the Computer Science Research Institute (CSRI).
  
  
  Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energys National Nuclear Security Administration contract number DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. 
  
  
  This work benefited from numerous discussions and interactions with our colleagues Eric Keiter, Suzey Gao and Larry Musson who shared their expertise with compact model development and TCAD simulations. We are grateful for their help and advice throughout the preparation of this paper.
  
  
  
  
  
  § INTRODUCTION
   
  I like chocolate pudding.  For that matter everyone else does too. This relationship is captured in (<ref>),
  
      f(x) = x^100-d,
  
  where f is happiness in utils, x is the quantity of chocolate pudding consumed in ounces and d is the age of the
  subject in years.
  
  It has been argued that that chocolate pudding consumption by age is shown in Figure <ref>.  However,
  we propose an alternative pudding consumption model, shown in Figure <ref>.
  
  
  
  
  
  We also want to show these figures side by side, to better illustrate our comparison.  We can see this in Figures
  <ref> and <ref>.
  
  
  
  
  
  § ACTUAL CONTENT
  
  
  As we saw in Section <ref>, we all like chocolate pudding. This is where I wish
   worked. It would fill the page with meaningless technobabble so I could illustrate this
  package. Instead, I'll talk about how to use quotations in latex. "Never use these quotations." “Always use these,
  instead.”
  
  
  
  § CONCLUSIONS
  
  Herein, we repeat the abstract in past tense.
  
  Unlike many other baked goods, chocolate pudding is subject to a myriad of interesting (and unique) effects on both the
  meso and nano scales.  Understanding these phenomena is critical, not only to America's restaurant industry, but to
  children everywhere.  We have examined these effects and have proposed new potential models which accurately capture
  the material structure of chocolate pudding.
  
  
  
  siam
  
  
  
  
  
  
  
  