## Abstract

An important challenge in particle-based modeling of electron–phonon interactions is the large difference in the statistical weight of the particles in the two simulated populations. Each change in the state of a simulated phonon during scattering is statistically representative of an interaction with multiple simulated electrons, which results in a large numerical burden accurately represent both populations. We developed two stochastic approaches to mitigate this numerical problem. The first approach is based on Poisson modeling of the scattering processes coupled with a thinning algorithm, which works effectively at steady-state, but it is prone to statistical errors in the energy during the transient regime. The second approach is based on point process (PP) modeling of the scattering, allowing stochastical book-keeping, which corrects the energy error. Here, we present a mathematical description of the problem and the two stochastic approaches along with the numerical results we obtained for the synchronous transient simulation of the electron and phonon populations.

## 1 Introduction

The increasing density and operation frequency of active and passive electron devices drives the need for more sophisticated heat management techniques. As a result, effective thermal management and controlled heat dissipation play a crucial role in device engineering because of their effects on both the device performance and reliability.

The fundamental equation governing semiclassical transport of electrons is Boltzmann's transport equation (BTE) [1], while the heat in a crystal lattice can be modeled in terms of phonons [2] described by a formulation of the BTE derived by Peierls [3]. These two equations describe the interaction of free charge carriers with the crystal in terms of absorption and emission of phonons [4].

There are a number examples in current electronics and photonics in which both the electronic and phonon populations are strongly coupled and driven far from equilibrium, necessitating a full solution of the BTE. The current problem in power dissipation in integrated circuit technology has already slowed Moore's law in terms of the scaling of clock speeds. For small dimension transistors and/or high power applications, the driving fields in the device are quite high, approaching MV/cm scale for wide bandgap technologies such as GaN, where electronic transport is quasi-ballistic and highly nonstationary. Such nonequilibrium electrons dissipate their energy on the drain side creating a phonon hot spot [5], where a localized far from equilibrium population of optical phonons is established, and diffusive transport of these nonequilibrium optical phonons as well as their harmonic decay to acoustic phonons are required to dissipate heat away from the device and to the environment. A second example is that of quantum cascade lasers based on intersubband emission of light in multiquantum well semiconductor structures, in which optical phonon emission plays a crucial role in the performance of such technologies. It is increasing recognized that the build-up of a nonequilibrium population of phonons in the active lasing region of III-V material-based quantum cascade lasers lead to detrimental effects in the lasing performance [6,7].

In general, a closed-form solution of the BTE cannot be obtained without severe approximations being made on the nature of the collision integral (e.g., the relaxation time approximation [8]), or on the distribution function itself (e.g., the near-equilibrium approximation). For this reason, both deterministic and stochastic integration schemes [5,918] have been developed to numerically determine the distribution function. Deterministic approaches are often based on lower order terms of the moment expansion of the distribution function [19], with approximations required to obtain closure. Furthermore, numerical representation of the energy dispersion cannot be included in such models, which have a relatively low computational cost. However, they usually require reliable closure relations in terms of phenomenological transport parameters [20,21], with their scope being limited [22] by the accuracy of such parameters. Alternatively, the particle-based ensemble Monte Carlo (EMC) method [23] provides a space-time stochastic solution of the full BTE, which is exact in statistical terms [9].

Previous coupled EMC approaches treated electrothermal interactions by employing the net energy lost via electron scattering to calculate a heat generation rate [24]. This rate was then used to obtain a temperature map by solving the heat transport equation either analytically [24,25] or stochastically [26,27]. This work expands on the particle-based EMC approach for phonon transport presented in Ref. [18] by directly modeling phonon transport within a particle-based framework. Both charge and heat are modeled simulating a statistically significant ensemble of particles, with each simulated particle representing an arbitrary number of physical particles [28]. In particular, each simulated charge carrier represents Wc physical electrons while each simulated phonon represents Wp physical phonons. Therefore, each simulated electron scattering represents the scattering of Wc physical electrons, and each simulated phonon scattering represents the scattering of Wp physical phonons. In this context, the direct modeling of the emission and absorption processes corresponds to the addition and removal of a simulated phonon from the simulation domain. Therefore, in the ideal case Wp = Wc, where one simulated phonon is either added or removed with each scattering event. Unfortunately, this approach is not feasible since, in a typical semiconductor, the density of phonons is orders of magnitude higher than the density of free charges.

For example, at room temperature, the phonon density is typically comparable to the density of atoms, about $1028 (m3)$, while in a highly doped semiconductor, the electron density is at most $1025 (m3)$. That means that a particle-based simulation with Wp = Wc requires at least a thousand simulated phonons for each simulated electron. Under these conditions, a Monte Carlo simulation requiring tens of thousands of simulated electrons would require 10 × 106 of simulated phonons. Moreover, the number of simulated electrons is typically in the tens (or hundreds) of thousands, since their number must be high enough to accurately represent the charge dynamics of the device; thus, this number cannot be arbitrary reduced. It is then obvious that imposing Wp = Wc is usually not a practical approach.

The next logical choice of allowing Wp > Wc complicates the numerical simulation of interactions between them. Since one physical electron interacts with one physical phonon, the scattering of one simulated electron represents the interaction of Wc physical electrons with Wc physical phonons. However, each simulated phonon represents Wp physical phonons; therefore, each scattering of a simulated electron must be treated as a scattering with a fraction of a simulated phonon.

In Sec. 2, we present a stochastic solution to this problem based on the traditional (for Monte Carlo simulation) premise of dealing with processes that obey Poisson statistics. We first start with a simple assumption of stationarity to describe the phonon emission (absorption) process with homogeneous Poisson statistics. This assumption allows us to monitor the creation (destruction) of simulated phonons to statistically represent the correct rate of emission (absorption) of physical phonons by employing a technique based on the so-called thinning algorithm [29], discussed in Sec. 3.

While this first approach results in a statistical error in the conservation of the system energy, during steady-state simulation, this error averages to zero. However, the equivalent approach (based on nonhomogeneous Poisson) for the transient, nonstationary regime results in a net error in the conservation of the system energy. For this reason, we developed another approach based on the Hawkes point process (PP), that is presented in 2 Sec. 4 of the paper.

This second approach uses the point process formalism to address the statistical error resulting from the simulated electron–phonon (Markovian) interaction during a transient by introducing a non-Markovian treatment of the error. The goal of this approach is to closely replicate the Poisson behavior in the stationary regime, while employing a stochastic form of book keeping capable of preserving (in statistical terms) the distribution of the generated and absorbed phonons, with negligible computational and memory requirement.

Section 5 presents numerical tests for the implementation of the presented algorithms for the case of a semiconductor system consisting of electrons and phonons. In particular, we show how the point process approach results in much smaller error in the system energy conservation.

## 2 The Monte Carlo Procedure

The particle-based Monte Carlo approach solves the semiclassical Boltzmann transport equation stochastically by modeling complicated transport phenomena as elementary interactions. In this approach, the many-body system is modeled by representing a statistically significant sample of its components as individual particles carrying physical quantities in phase space and time. Elementary interactions are treated as collisions, where the duration of scattering events is assumed much shorter than the particle mean free flight times. If the probability of the scattering event follows a Poisson distribution, the collision integral can be solved using the stochastic Monte Carlo procedure, that requires only the average time between scattering events. In this approach, the properties of the many-body system are obtained from the average properties of the simulated distribution.

The particle-based Monte Carlo method presented in this work solves the electron–phonon transport problem, defined by the coupled Boltzmann and Peierls–Boltzmann transport equations:
$∂fc∂t+∇k(j)fc·∂k(j)∂t+∇rfc·∂r∂t=(∂fc∂t(fp))Coll$
(1)
$∂fp∂t+∇rfp·∂r∂t=(∂fp∂t(fc))Coll$
(2)

where fc and fp are the distributions of electrons and phonons, respectively, $r$ is the position-space vectors while $k(j)$ represents a state in the jth band with momentum $k$. The terms on the right side of the equations represent the change of states due to collisions.

The reciprocal dependencies of each distribution function on the other have been explicitly written to highlight the coupled nature of these transport phenomena. Finally, for the sake of notation clarity, from here on the band index will be dropped, and the simple momentum notation will be used to represent both a state and the corresponding momentum vector.

### 2.1 Scattering Algorithm.

The cellular Monte Carlo method (CMC) extended via the acceptance/rejection technique is used in this work to reduce the number of calculations required by each scattering and maintain the accuracy of the EMC technique. In this approach, the scattering procedure [18] can be separated into four steps. First, a transition is selected from a precomputed look-up table. Then, new occupation numbers are obtained by sampling the run-time distribution of the simulated particles. A new scattering probability is calculated by using the new occupation numbers. Finally, the acceptance/rejection technique is used to sample the new transition probability from the precomputed one.

To compute the look-up table, we assume that all states $k$ in a small enough region of momentum space, $Ωk0$, centered around the state $k0$ have a transition probability similar to the one of $k0$. The total transition rate from $k0$ to a small region of space $Ωk0′$ centered around the final state $k0′$ is computed using Fermi's golden rule
$Γ(k0,Ωk0′)=2πℏ||2Dos[Ωk0′](E(k0))$
(3)

where M is the matrix element, E(k0) is the energy of the state $k0$, and $Dos(Ωk0′)$ is the density of states in the cell $Ωk0′$. Explicit selection rules are not necessary since they arise from the interaction of the dispersion/band structure with momentum and the energy conservation. Momentum conservation is explicitly enforced in the matrix element, while the energy conservation is implicitly accounted for in the Dos operator.

With these considerations, the momentum space can be subdivided into a finite number of cells. In this space, allowed transitions from states of each cell to all other cells can be computed and stored in the look-up table. The scattering procedure can then simply pick one particular transition from the precomputed table that is loaded at the start of a simulation. This approach exchanges fast access memory (RAM) usage for the computations used in each EMC scattering to: search for all the possible final states, calculate all the scattering rates, and compute the scattering partial sums.

The simulations presented in this paper employ full-band electron band structures computed via the empirical pseudo-potential method, and a phonon dispersion computed via a shell model. Details regarding the matrix elements can be found in Refs. [18] and [30].

The second step requires computing the occupation number, that is the distribution function of a particular state. This task is achieved by setting the occupation number of a state equal to its average value around that state. In other words, the occupation number in the position r of the state $k$ is calculated as the ratio of the number of simulated particles in a phase-space volume centered around $(r,k)$ and the number of states in the same volume. The sampling volume itself is defined by a component in position space, $Vr$, centered around r, and a component in momentum space, $Ωk$, centered around $k$.

The final step employs the acceptance/rejection method, which is a basic statistical technique used to generate a distribution, difficult to obtain directly, by rejection sampling from a different, more convenient, distribution. In this approach, the acceptance/rejection method provides the flexibility and the accuracy of the EMC method while keeping the computational demands comparable to that of the CMC method.

### 2.2 Numerical Considerations.

Compared to a more standard CMC method [30], there are two additional aspects impacting the computational demands of this approach: the acceptance/rejection technique, and the run-time calculation of the occupation number. In this work, the latter is the largest contributor to the computational complexity. Indeed, each occupation number requires searching over a portion of the simulated particles resulting in a computational complexity of $O(Np)$, where Np is the number of simulated phonons. Therefore, getting a new occupation number with each scattering event gives a total computational complexity for the scattering algorithm of the order of $O(N2)$, where N is the total number of simulated particles. While this method is computationally intensive, and further optimizations may be used, it allows altering each scattering rate according to the local run-time phase-space distribution.

In contrast, the simple acceptance/rejection method requires a constant number of calculations for each scattering event. When using less demanding methods to obtain the occupation number, e.g., an analytic distribution, this technique results in an increase of the wall-clock simulation time of a few percent.

Finally, in position-space/drift simulations, the use of a position-space component of the sampling volume, $Vr$, requires a tradeoff between $Vr$ itself and N. Smaller $Vr$ requires an increased number of particles to reduce the statistical noise of the occupation number. Larger $Vr$ limits the spatial resolution and may result in nonphysical interactions between particles too distant from each other. The values of $Vr$ and N will depend on the phenomena of interest and material properties. For example, phenomena strongly localized in phase-space will require smaller $Vr$ and a larger number of particles, while phenomena requiring low spatial resolution may use larger $Vr$ and smaller N.

## 3 Phonon Absorption and Emission as a Poison Distribution

### 3.1 Poisson Model.

When the statistics of the scattering modeling for the electron–phonon emission process is stationary, we can make the following assumptions:

1. The number of emissions occurring during any time interval is an integer.

2. The emissions occur independently.

3. The scattering rate is not time-dependent.

4. The number of scattering events during any time interval is proportional to the length of the interval.

5. The probability of two emission events happening in the exact same time is vanishing small.

Under these conditions, the phonon emission can be simply modeled as a homogeneous Poisson distribution [31,32]
$Pe(dηe(q,t)=k)=(dtλe(q))ke−dtλe(q)k!$
(4)

where P represents a probability, the superscript e is used to indicate an emission process, $dηe(q,t)$ is the total number of (physical) charges emitting a phonon into the state $q$ during the time interval $[t,t+dt ]$, k is an integer number, and the parameter $λe(q)$ is called the rate or intensity of the process, and, since we assumed a homogeneous Poisson distribution, is constant in time.

It follows from the previous equation that [31,32] we can write the expectation of $dηe$ as follows:
$E[dηe(q,t)]=dt ηcΓe(q)=dt λe(q)$
(5)

where the symbol E represents the expectation operator [31] such that the expected value, $E(x)$, of the random variable x is the mean value of x, ηc represents the total number of (physical) charged particles (either holes or electrons), and $Γe$ is the charge-phonon emission rate.

In the context of particle-based Monte Carlo simulations, the number of simulated charges undergoing scattering through phonon emission in a time interval dt centered around the time $t, dNce$, satisfy the following relation:
$WcE[dNce(q,t)]=Wc dt λce(q)=E[dηe(q,t)]$
(6)

where N is used to represent the number of simulated particles and the superscript c is used to indicate charged particles. Consequently, $dNce$ represents the number of simulated charged electrons, each one representing Wc physical particles undergoing phonon emission scattering. Finally, $λce(q)$ is the intensity of the Poisson process modeling the simulated scattering.

A similar approach is employed for the number of emitted simulated phonons, $dNpe$
$WpE[dNpe(q,t)]=Wp dt λpe(q)=E[dηe(q,t)]$
(7)

where the superscript p is used to indicate a simulated phonon representing Wp physical particles, and $λpe(q)$ is the intensity of the Poisson process modeling the emission (creation) of simulated phonons.

The simple approach of generating a new simulated phonon with each scattering of a simulated charge corresponds to a Poisson process with intensity $λce(q)$ that generates an average of $dt λce(q)$ new simulated phonons in the interval dt. However, it is clear from Eqs. (6) and (7) that the following is true:
$dtλpe(q)=dtλce(q)WcWp$
(8)

Consequently, $dt λce(q)>dt λpe(q)$ since $Wc. Therefore, this approach would generate a number of simulated phonons larger than the expected value calculated in Eq. (7). A similar treatment in the case of phonon absorption leads to an equivalent conclusion; consequently, this simple approach cannot be used.

An intuitive solution to this problem is to create a new simulated phonon as result of only a fraction of the scattering events, while preserving the Poisson nature of the phonon generation process. Indeed, this goal can be achieved by employing a technique called thinning [29]. A theoretical overview and the implementation details of this approach are presented in the next Secs. 3.2 and 3.3.

### 3.2 Thinning Algorithm.

The thinning technique generates events obeying a particular Poisson distribution, which will be referred to as the target distribution, by down-sampling events according to a different Poisson distribution, which will be referred to as the sampling distribution. In other words, the algorithm reproduces the target distribution by randomly selecting a subset of events from a suitable sampling distribution. Obviously, this technique requires the events of the sampling distribution to be more frequent than the events of the target distribution; the sampling Poisson distribution is therefore said to be stochastically dominating [33] the target distribution. The requirement of stochastic domination results in a lower bound for the Poisson intensity of the sampling distribution $λs(t)$
$λs(t)>λt(t) ∀t$
(9)

where $λt(t)$ is the intensity of the target distribution. The events are then retained with probability $λt(t)/λs(t)$ or rejected with probability $1−λt(t)/λs(t)$. Finally, $λs(t)$ is typically chosen to be a constant in order to simplify the algorithm that generates the sampling distribution. In the specific case of phonon emission, $λpe(q)$ represents the intensity of the target distribution and the intensity of the dominating distribution is chosen to be $λce(q)$.

We implemented the thinning procedure into the algorithm shown in Fig. 1. After a simulated charge scatters, a rejection probability is computed as
$Prej=λpe(q)λce(q)=WcWp$
(10)
Fig. 1
Fig. 1
Close modal
Then, a random number, $ϱu$, is generated from the uniform distribution $U(0;1)$. Finally, if $ϱu, a new simulated phonon is created in the same position of the simulated charged particle; otherwise, nothing occurs. The validity of this approach is proven by comparing the value calculated in Eq. (7) with the expected number of simulated phonons, $NThinpe$, created by the thinning procedure up to time t
$E[dNThinpe(q,t)]=E[dNce(q,t)|ϱu
(11)

Where $|$ is the conditional expectation operator such that $E[A|B]$ is the expected value of the event A computed with respect to the conditional probability distribution function $p[A|B]$ [31]. Equation (11) shows that the average number of simulated phonons created by the thinning algorithm corresponds to the expected value as calculated in Eq. (7).

A similar treatment of the mechanism of phonon absorption leads to an algorithm employing the same rejection probability as in Eq. (10). In this case, the simulated phonon to be removed is chosen according to the procedure described in Ref. [18].

### 3.3 Transient Regime.

The limits of the technique presented in Sec. 3.2 become evident when analyzing the expected rate at which energy is lost and gained by the two simulated populations at steady-state and in the transient regime. Specifically, the expected energy loss by the simulated charge carriers in a differential time interval $dt , dEce$, can be computed by adding together the contribution of each scattering occurring during dt over all the possible states in the first Brillouin zone (BZ1) [34]
$E[dEce(t)]=E[∑q⊂BZ1ℏωqWcdNce(t,q)]$
(12)
where $ωq$ is the frequency of the state $q$, Wc is the statistical weight of the simulated charges, and $dNce(t,q)$ is as defined in Eq. (6) The linearity of the expectation operator and the properties of the Poisson distribution allow us to rewrite the previous equation as
$E[dEce(t)]=∑q⊂BZ1ℏωqWcE[dNce(t,q)]=∑q⊂BZ1ℏωqWc dt λce(q,t)$
(13)

where the product $dt λce(t,q)$ is the expectation, or the average value, of $dNce(t,q)$.

In a similar way, the expected energy gained by the population of simulated phonons during $dt , dEpe(t)$, is calculated by adding the contribution of each newly created particle
$E[dEpe(t)]=∑q⊂BZ1ℏωqWpE[dNpe(t)]=∑q⊂BZ1ℏωqWp dt λpe(t,q)$
(14)

where $dNpe$ is defined in Eq. (7), and, as usual, Wp is the number of physical phonons represented by each simulated phonon. Finally, the product $dt λpe(t,q)$ is the expectation of $dNpe(t,q)$.

Since the thinning algorithm requires that $λpe(t,q)=λce(t,q) Wc/Wp$, it is clear that the average energy lost by the population of simulated charge carriers equals the energy gained by the population of simulated phonons.

A similar approach can be applied to analyze the thinning algorithm applied to the absorption process. In this case, the equations show that the average energy gained by the population of simulated charge carriers via absorption, $dEca$, equals the energy lost because of the removal of simulated phonons, $dEpa$.

The corresponding central variances for Eqs. (12) and (14) are, respectively,
$var[dEce(t)]=∑q⊂BZ1( dt ℏωqWc)2λce(q)var[dEpe(t)]=∑q⊂BZ1( dt ℏωqWp)2λpe(q)=∑q⊂BZ1( dt ℏωq)2WpWcλce(q)$
(15)

Since Wp > Wc (typically, $Wp≫Wc$), the variance of the energy gained by the simulated phonon population is larger than the variance of the energy lost by the simulated charge population: $var[dEpe(t)]>var[dEce(t)]$. In the case of phonon absorption, the same argument gives $var[dEpa(t)]>var[dEca(t)]$.

In the steady-state regime, instantaneous differences between the energy gained by the simulated phonon population, $dEpe(t)$, and the energy lost by the simulated charged carriers, $dEce(t)$, have no lasting effect since, from Eq. (11), both quantities have the same average value. The total energy of the charge-phonon system is therefore conserved on the average, and the effect of the difference in variance results in larger noise in the energy of the simulated phonon population as compared to the simulated charge population. As a consequence, the relative difference between the energy of the two populations existing before reaching the steady-state regime will remain unaltered. More explicitly, the conditions present at the end of a transient are preserved in the steady-state regime, as they should be. Following the same argument, any random fluctuation of the difference in energy $dEpe(t)−dEce(t)$ present at the end of a transient will be preserved in the steady-state regime. In other words, the statistical fluctuation of the transient regime carries over in the steady-state regime. This creates an irremovable discrepancy between the energy gained by the simulated crystal and the energy lost by the simulated charge particles. The larger the fluctuation, the larger the error in the conservation of energy. In other words, the thinning algorithm is not fast enough to ensure the system energy conservation in the presence of fast transients.

It is, in principle, possible to employ a simple book-keeping scheme that accounts for the error in the conservation of the energy exchanged between the two populations. This approach would compensate for any error by adding or removing the appropriate number of simulated phonons, when the energy discrepancy justifies it. However, a simple book-keeping scheme requires somewhat arbitrary decisions on the state of the simulated phonon to be removed and in which state a new simulated phonon should be created. In order to avoid such arbitrary decisions, a technique based on the probabilistic model called the point process [35] is used in place of the Poisson distribution approach presented in Sec. 3.1. The proposed framework allows modeling the creation and removal of simulated phonons in a way that statistically enforces the conservation of the system energy. The mathematical basis for the modeling of absorption and emission processes and the algorithm employed in this work will be discussed in Sec. 4.

## 4 Phonon Absorption and Emission as Point Processes

### 4.1 Point Process.

In general terms, a PP [36] can be described as a sequence of points randomly located in some space. However, in the context of this work, the only type of point process considered is defined on the one-dimensional time space, represented by the letter t. In this simple case, a PP can be thought of as being a collection of points (in our case the scattering events) randomly placed on the time line, and the specific location patterns and these points constitute the process itself. In other words, the term process will denote here a collection of events occurring in time, and the term stochastic process will identify a process that develops in time according to some probabilistic rule. A point process is commonly described in terms of its counting measure $M(t)$ [35,37,38], which represents the cumulative count of events in a process. This means that for two points in times t1 and t0, where $t1>t0$, the value $M(t1)−M(t0)$ is a positive integer representing the number of events occurring in the time interval $[t1,t0]$. In particular, $M(t)$ is called a random counting measure when the counted events are parts of a stochastic process.

More formally, a function $M(t)$ defined for time $t>0$ is said to be a counting measure if [35,37]

1. $M(0)=0$

2. $M(t)$ is a positive integer $∀t>0$

3. $M(t)$ is nondecreasing, right-continuous step function with jumps of finite size.

The counting function can be expressed as a summation of step functions [31]
$M(t)=∑i=1∞Si(t)$
(16)
where the step function $Si(t)$ is defined as
$Si(t)={1 t≥ti0t
(17)

where $ti$ is the time of the ith event.

The PP of interest in this work can be defined in terms of the probability for new events to occur
$lim dt →0+P[M(t+dt )−M(t)=j]{o( dt )j>1λ(t) dt+o(dt)j=11−λ(t) dt+o(dt)j=0$
(18)
where $P[A]$ represents the probability of the event A, $M$ is the random counting measure of the PP, λ(t) is referred to as intensity, $o( dt )$ is the Bachmann–Landau (or little o) notation [35], which simply indicates something much smaller than dt, this term explicitly represents the vanishingly small probability of two events happening at the same instant in time. A point process with this last property is said to be simple [35,36], and results in a straightforward relationship between the intensity of the PP, the probability, and the expected value [35,38] of new events occurring in a time interval dt. In particular, for small values of $dt, o(dt)=0$ and the previous quantities are simply related by
$dt λ(t)=P[M(t+dt )−M(t)=1]=E[M(t+dt )−M(t)]$
(19)

Figure 2 shows a graphical representation of a PP described by Eq. (18). The dots on the time line mark the corresponding events, while solid line represents the PP intensity as a function of time. A shorter (average) distance between events corresponds to a higher probability of a new event happening. Since Eq. (19) states that the probability and intensity are proportional to each other, it follows that a higher density of events results in a higher intensity while a more sparse distribution of events results in a lower intensity. In other words, the average distance between events, centered around a specific time t, has an inverse relationship with the corresponding intensity, λ(t). This particular case shows the alternation between high and low intensity and the corresponding clustering of events centered around the peak of intensity.

Fig. 2
Fig. 2
Close modal

The point process framework provides extensive flexibility in the definition of the intensity λ(t) of Eq. (18). When λ(t) is a deterministic function, the number of events in an interval dt obeys a Poisson distribution with intensity λ(t). In other words, the quantity $M(t+dt)−M(t)$ is a nonhomogeneous Poisson random variable with intensity λ(t) [31,35]. In this general case, the PP is referred to as a nonhomogeneous Poisson point process. Obviously, if λ(t) is a constant, the quantity $M(t+dt )−M(t)$ is a homogeneous Poisson random variable, and the PP is referred to as a homogeneous Poisson point process. As one would expect, this definition of Poisson PP in terms of a counting measure allows the derivation of all the properties of the Poisson distribution [31,32]. If λ(t) is not a deterministic function, but another independent Poisson random variable, the PP is referred to as a Cox Doubly stochastic Poisson process [32,39,40], or simply a Cox process, which is considered a generalization of the Poisson process where the intensity itself is a random variable. The term Cox process is also commonly used in the literature to describe any PP with a Markovian [35] intensity λ(t). In the case of a non-Markovian intensity, the PP is commonly referred to as a Hawkes process [32,35,39]. In the same way that the Cox process can be considered a generalization of the Poisson process, the Hawkes process can be thought of as a non-Markovian extension of the Poisson PP. In a sense, the Hawkes PP is a Poisson PP which accounts for the effect of past events. For this reason, an approach based on the Hawkes PP seems a natural extension to the model based on the Poisson distribution presented in Sec. 3.2. It may seem unusual to use such a non-Markovian approach in the context of a stochastic method to solve the collision integral of the BTE that assumes Markovian interactions. The reasons behind this approach can be understood by comparing the PP-based approach and a basic book-keeping scheme. Specifically, both cases keep track of the error, then in the book-keeping scheme, particles are added or removed outside of the scattering process. However, the PP approach transfers this second part of the book-keeping scheme directly into the scattering algorithm by altering the probability of generating or removing particles according to the value of the tracked error. As such, this PP approach can be thought of as a statistical form of book-keeping. More explicitly, since the statistical weight of the simulated phonon is much higher than the statistical weight of the simulated electron, as such particles cannot be simply created or removed with each scattering event, and a choice must be made on when to alter the number of simulated phonons. In this work, we decided to make such choice by a non-Markovian PP that compensates for the error in the energy conservation. In other words, it is not the Markovian electron–phonon interaction that is modeled via a non-Markovian PP, but how the simulated phonons are created and removed as a result of scattering.

### 4.2 Hawkes Model.

The Hawkes [39] point process is typically used to model systems where each new event impacts the whole stochastic process. More explicitly, the chance of a new event has some dependence on the sequence of events that preceded it. This type of process is typically employed to model systems where the observed events tend to cluster in time. For example, seismic events present short-term clustering; an earthquake main-shock is both preceded and followed by intermediate-sized events [41]. Hawkes processes have been employed to model criminal activity such as property crime, gang violence, and terrorism. Each initial event increases the likelihood of more events, as the successful offender tends to replicate the crime in a nearby location in the following weeks [42]; or, in the case of gang violence, each violent event tends to ignite retaliation [43].

In all these examples, each event makes the occurrence of the subsequent events more probable for some period of time; the corresponding types of PPs are referred to as self-exciting Point Processes. Formally, this behavior can be expressed in terms of the counting measure as
$cov(M(t1)−M(t0),M(t2)−M(t1))>0$
(20)

for any $t0, where cov denotes the covariance of the two quantities. When the covariance is negative, the PP is called self-correcting or self-regulating; in that case each event makes the occurrence of the following events less probable for some period of time. It is clear that by simultaneously modeling two such opposite processes, that the energy error can be both positive and negative. Consequently, the employed PP can be neither self-exciting nor self-correcting, but it should account for both positive and negative covariance.

The Hawkes process is defined by relating the counting measure to the collection of previously occurred events, which will be represented by the symbol $H(t)$ and referred to as the history1. In particular, the Hawkes process is defined by expressing the conditional probability of observing a particular counting function given the history
$lim dt →0+P(M(t+dt )−M(t)=m|H(t))={o( dt )m>1λ(t,H(t)) dt+o( dt )m=11−λ(t,H(t)) dt+o( dt )m=0$
(21)
where $|$ represents the conditional probability operator2 [31], and $λ(t,H(t))$ is referred to as the intensity of the process. For example, in the case studied by Hawkes [39], the intensity takes the form
$λ(t,H(t))=λ0+∫0taeb(t−s)dM(s)$
(22)
where $λ0$ is the background intensity of the process, a and b are positive real constants, and the term $dM(s)$ is
$dM(t)=M′(t) dt=∑i=1∞δ(t−ti) dt$
(23)
where ti is the time of the ith event, $M$ is replaced by the definition in Eq. (16), and the delta function is the derivative of the step function. The previous equation allows us to rewrite Eq. (22) as
$λ(t,H(t))=λ0+∑ti
(24)

It is now clear that each event instantaneously increases the arrival intensity by A, and that this influence decays exponentially over time.

In more general terms, the intensity of a nonlinear Hawkes process can be written as [35,44]
$λ(t,H)=ϑ(ζ(t)+∫−∞th(t,s)dM(s))$
(25)
where the function $ϑ$ is such that $ϑ:R→ R+$, and the function $ζ$ accounts for the evolution of the stochastic process in the absence of events; in other words, it provides a background intensity to the PP. The term h is referred to as the exciting function, and can be rewritten by using the definition of $dM(s)$ in Eq. (23) as
$∫0th(t,s)dM(s)=∑i=1∞∫0th(t,s)δ(s−ti)ds$
(26)
where the order of the summation and integration has been reversed. The final integration can be split in two cases
$∫0th(t,s)δ(s−ti)ds=h(t,ti) ∀ti≤t∫0th(t,s)δ(s−ti)ds=0∀ti≥t$
(27)
These last two equations allow us to rewrite Eq. (25) as
$λ(t,H)=ϑ(ζ(t)+∑ti
(28)

This expression clearly shows the relationship between the exciting function, past events, and the intensity. More explicitly, each past event influences the future intensity $λ(t,H)$ thought its weight h.

### 4.3 Multivariate Hawkes Point Process.

In this work, the creation and removal of simulated phonons are modeled in terms of a multivariate [44] Hawkes PP in order to minimize the effects of the statistical fluctuations of the energy conservation. The definition of the PP must ensure that the variation in energy of the phonon population closely mirrors the variation in the population of simulated charges, in spite of the difference between the statistical weights of the particles in the two populations. In other words, in order to correctly model transients, the model must alter the rate of generation or removal of simulated phonons to achieve conservation of the system energy.

In general, the effect of all other variates on the ith variate of a multivariate Hawkes PP is represented in the exciting function
$λi(t,H)=ϑ(ζ(t)+∑jNv∫hij(t,s)dMj(s))$
(29)

where Nv is the number of variates, $hij$ defines how the jth variate affects the ith variate, and $Mj$ is the random counting measure of the jth variate.

The electron–phonon interactions are modeled by employing four types of variates to represent all the possible energy exchanges happening between the simulated populations. In particular, each phonon state with momentum $q$ is connected to four variates

• Electron emission of a phonon in $q$, represented by the superscript ce.

• Electron absorption of a phonon in $q$, represented by the superscript ca.

• Phonon in $q$ created by phonon emission, represented by the superscript pe.

• Phonon in $q$ removed via phonon absorption, represented by the superscript pa.

These four types of variates allow the expression of the energy exchange between the two populations in terms of an energy balance equation
$dE (t)=dEpe(t)+dEpa(t)−dEce(t)−dEca(t)$
(30)
where the notation $dE (t)$ represents $E(t+dt)−E(t)$. In particular, expressions for $dEce$ and $dEpe$ have been explicitly defined in Eqs. (13) and (14) in the case of phonon emission. The quantities $dEca(t)$ and $dEpa(t)$ are similarly defined for the case of phonon absorption
$dEca(t)=∑q⊂BZ1ℏωqWcdNca(q,t)dEpa(t)=∑q⊂BZ1ℏωqWpdNpa(q,t)$
(31)

The time integral of $dE (t)$ is a measure of the error in the conservation of energy of the simulated interacting populations.

Moreover, the expected value, or expectation, of $dE (t)$ when generation and removal of simulated phonons are modeled as a Hawkes PP can be written as
$E[ dE (t)|H]=∑qℏωqWp(E[dNpa(q,t)|H]−E[dNpe(q,t)|H])−ℏωqWc(E[dNca(q,t)|H]−E[dNce(q,t)|H])$
(32)
where the sum has been separated because of the linearity of the expectation operator. The definition of the PP intensity and Eq. (19) allow rewriting the previous equation in terms of the intensities
$E[ dE (t)|H]=dt ∑qℏωqWp(λpa(q,t)−λpe(q,t))−ℏωqWc(λca(q,t)−λce(q,t))$
(33)

where $λca(q,t)$ and $λce(q,t)$ are the intensities corresponding to the scattering of simulated charged particles as defined in Sec. 3.3. The quantities $λpa(q,t)$ and $λpe(q,t)$ are the intensities corresponding to the removal and generation of simulated phonons, respectively. A simple algebraic substitution shows that the thinning algorithm presented in Sec. 3.2 produces an expected value of zero for $dE (t)$ at steady-state.

When the creation and the removal of simulated phonons are modeled as Hawkes processes, the intensities $λpe(q,t)$ and $λpa(q,t)$ are designed to minimize the error in the conservation of energy between the simulated populations, which is represented by the integral $∫ dE (t) dt$. In particular, the desired properties for the intensity $λpa(q,t)$ are as follows:

• $λpa(q,t)$ has a background intensity of $λca(q,t)*Wc/Wp$.

• $λpa(q,t)$ has an inverse dependence on $∫ dE (t)$.

• $λpa(q,t)$ has a maximum intensity of $λca(q,t)$.

The first property ensures that, in absence of an error, the process behaves like the Poisson model described in Sec. 3.2; the second property imposes a negative feedback on the value of $∫ dE (t)$, while the final property allows the simplification of the algorithm by ensuring that each simulated charged particle interacts with no more than one simulated phonon.

In a similar fashion, the desired properties for the intensity λpe are

• $λpa(q,t)$ has a background intensity of $λca(q,t)*Wc/Wp$.

• $λpa(q,t)$ has a forward dependence on $∫ dE (t)$.

• $λpa(q,t)$ has a maximum intensity of $λca(q,t)$.

There are a variety of definitions of λpe and λpa that can satisfy these properties, and each definition has its own merits and drawbacks. The most basic approach would employ a step-like intensity
$λpe(q,t)={λce(q,t)WcWp∫ dE (t)≤0λce(q,t)∫ dE (t)>0$
(34)
The implementation of this approach would result in a numerically large burden on the algorithm by continuously imposing an unnecessary amount of removal and creation of simulated phonons. Moreover, the intensity reduces to the Poisson model presented in Sec. 3.2 only when $∫ dE (t)=0$, which is not a likely event. This means that at steady-state, the removal and creation of simulated phonons will largely differ from the expected Poisson distribution. This problem could be mitigated by using some constant $E max∗$ as a threshold
$λpe(q,t)={λce(q,t)WcWp∫ dE (t)
(35)

This approach does not burden the algorithm and it reproduces the desired Poisson-like distribution behavior for $∫ dE (t). However, there is no effect on the intensity until $∫ dE (t)≥E max∗$, and this can result in large delay of the algorithms response to the error in $∫ dE (t)$. This effect results in a statistical difference between the momentum distribution of the created phonons and the phonon states which created the error by scattering with the simulated charge carriers. In other words, the momentum distribution of the generated phonons does not correspond to the variation in the momentum of the scattered charged particles.

This issue can also be mitigated by employing some type of transition function between the maximum and minimum value of the intensity. The simplest approach is to employ a linear relationship in the definition of intensity
$λpe(q,t)=λce(q,t)WcWp+a∫ dE (t)$
(36)

where the proportionality constant a determines the sensitivity of the intensity. Small values of a closely reproduce the Poisson-like behavior presented in Sec. 3.2 but allow for large values of the error $∫ dE (t) dt$. In contrast, while large values of a limit the magnitude of the error, the behavior of the intensity would largely differ from the Poisson case presented in Sec. 3.2.

This work employs a log-linear model to balance the efficiency of the algorithm, the need for an instantaneous response, and the need to limit the error $∫ dE (t) dt$. This approach provides an immediate and modulated response to the error represented by $∫ dE (t) dt$. In particular, the effect is minimal and the rate of removal and creation of simulated phonons does not diverge much from the Poisson model presented in Sec. 3.2. As the error grows, the effect on the intensity becomes more pronounced, ensuring the desired small delay between the increase of the error value and the algorithm response.

### 4.4 Log-Linear Model of the Intensity.

In general, a log-linear approach models functions whose logarithm equals a linear combination of parameters. In this case, it is the logarithm of the intensity that is proportional to a linear function. This is simply achieved by requiring $ϑ$ in Eq. (28) to be an exponential function. In particular, $λpe(q,t)$ and $λpa(q,t)$ are respectively given by
$λpe(q,t,H)=exp(Lpe(q,t)+∑q∫−∞thpe(q,t,s)dM(s))λpa(q,t,H)=exp(Lpa(q,t)+∑q∫−∞thpa(q,t,s)dM(s))$
(37)
In the absence of the exciting function h, the intensity should reduce to the Poisson model presented in Sec. 3.2. This allows the derivation of expressions for the terms $Lpe(q,t)$ and $Lpe(q,t)$ by imposing
$WcWpλce=λpe(q,t,H)=exp(Lpe(q,t))WcWpλca=λpa(q,t,H)=exp(Lpa(q,t))$
(38)
which leads us to rewrite Eq. (37) as
$λpe(q,t,H)=WcWpλce exp(∑q∫−∞thpe(q,t,s)dM(s))λpa(q,t,H)=WcWpλca exp(∑q∫−∞thpa(q,t,s)dM(s))$
(39)
The exciting functions must be proportional to the measure of error presented in Eq. (12), and the two intensities should have opposite behavior
$hpe(q,t,s)=−hpa(q,t,s)=aq(dEpa(q,t)−dEpe(q,t)−dEca(q,t)+dEce(q,t))$
(40)
This equation allows us to explicitly write the arguments of the exponentials in Eq. (39)
$∑q∫−∞thpe(q,t,s)dM(s)=+a∫−∞t dE (t)∑q∫−∞thpa(q,t,s)dM(s)=−a∫−∞t dE (t)$
(41)

where the factor $aq=a$ is assumed constant in q, the summation and integration order have been reversed, and $dE (t)$ is the same as in Eq. (30).

Finally, the intensities for the creation and removal of simulated phonons are, respectively,
$λpe(q,t,H)=min(WcWpλcee+a∫ dE (t),λce)λpa(q,t,H)=min(WcWpλcae−a∫ dE (t),λca)$
(42)

where the minimum function accounts for the limit imposed on the maximum intensity. The factor a has the dimensions of the inverse of an energy; intuitively, this factor determines the sensitivity of the intensity to changes in $∫ dE (t)$.

In the specific case of
$λpa=WcWpλcae−a∫ dE (t)λpe=WcWpλceea∫ dE (t)$
(43)
the expectation in Eq. (33) can be simplified as follows:
$E[ dE (t)|H]=dt ∑qℏωqWc(λca(q,t)(e−a∫ dE (t)−1)−λce(q,t)(1−ea∫ dE (t)))=dt (e−a∫ dE (t)−1) E ca(t)−dt (ea∫ dE (t)−1) E ce(t)$
(44)
These equations allow us to write the expected values for the integral of $dE (t)$ in terms of a recurrence relation [45]
$E[∫−∞ti+1 dE (t)|H]=Ei+1=Ei+(e−aEi−1) dt E ca−(eaEi−1) dt E ce$
(45)

where ti is the time of the ith event. A simple substitution proves that $Ei+1=0$ in the special case of Ei = 0. At the start of the simulation we can assume no error in the energy exchange between the simulated populations, i.e., we impose $E0=0$. Therefore, Ei = 0 for each n, and the expected value for the integral of $dE (t)$ must be zero for each t. Since $∫ dE (t)$ is a measure of the error in the conservation of energy in scattering between the simulated populations, the energy gained by the simulated phonon population must equal the energy lost by the population of simulated charged particles.

Finding the analytical expression for the variance of $∫ dE (t)$ is challenging and beyond the scope of this work. However, an upper bound to the central variances can be computed by employing the Popoviciu's inequality [46]. In particular, the variance of a random variable x is bound by
$var[x]≤(max(x)−min(x))2/4$
(46)

where $max(x)$ and $min(x)$ represent the maximum and minimum value of x, respectively.

In order to compute the extrema of $∫ dE (t)$, we consider separately the cases of small and large values of a. For large values of a, the intensities defined in Eq. (42) oscillate between zero and the maximum value. In this case, the integral $∫ dE (t)$ is bounded by the discrete energy of the simulated phonons
$−Wpℏωmax<∫ dE (t)<+Wpℏωmax$
(47)

where $ωmax$ is the maximum phonon energy.

For small values of a, the extremum of $∫ dE (t)$ can be obtained from the definitions of intensity in Eq. (42). In particular, for positive values of $∫ dE (t)$, the upper bound is computed by imposing $λpe≤λce$
$∫ dE (t)<+1alnWpWc$
(48)
Similarly, for negative values of $∫ dE (t)$, the lower bound is computed by setting $λpa≥λca$
$∫ dE (t)>−1alnWpWc$
(49)
The inequalities in Eqs. (47)(49) provide the final bound to $∫ dE (t)$
$∫ dE (t)<+1alnWpWc+Wpℏωmax∫ dE (t)>−1alnWpWc−Wpℏωmax$
(50)
Finally, the upper bound of the central variances is computed from Eq. (46) as
$var[∫ dE (t)]<(1aln(WpWc)+Wpℏωmax)2$
(51)

This equation confirms the intuitive notion that the variance depends on the ratio of Wp and Wc. Moreover, the equation provides a practical approach for the selection of the parameter a by choosing the magnitude of the variance.

### 4.5 Hawkes Model Implementation.

The literature proposes various sampling algorithms [4750] to generate samples of point processes. This work employs a technique based on Ogata's [49] thinning approach.

Ogata frames the thinning algorithm in terms of a target and a sampling PP when both PPs are simple, as defined in Sec. 4.1. In this context, the rejection probability is expressed as a ratio of the PP intensity [35,49]. If the intensities of both PPs are deterministic, as in a Poisson PP, the thinning procedure follows the basic thinning algorithm [29]. However, when at least one of the intensities is a random variable, as in a Hawkes processes, the algorithm must account directly [35] for the inequality in Eq. (9). That is, at any time the intensity of the sampling PP must be larger than the intensity of the target PP. Ogata achieves this goal by employing a sequential variant of the thinning algorithm. In this sequential approach, the inequality is ensured by recomputing the intensity of the sampling PP after each event alongside to the value of the intensity of the target PP.

In the case of the log-linear model presented in this work, the intensities of the Hawkes PP as defined in Eq. (42) are bounded random variables. This condition allows to largely simplify the sampling algorithm by selecting an appropriate sampling PP that ensures the inequality in Eq. (9). The definition in Eq. (42) provides that the intensity for the generation of simulated phonons is always smaller than the intensity for the emission scattering of simulated charge carriers. Similarly, the intensity for the removal of simulated phonons is always smaller than the intensity for the absorption scattering of simulated charge carriers. These two relationships allow the implementation of the thinning procedure by simply following the procedure detailed in Sec. 3.2, but employing the following rejection probabilities:
$Preja=λpa(q,t,H)λca(q,t)=WcWp e−a∫ dE (t)Preje=λpe(q,t,H)λce(q,t)=WcWp ea∫ dE (t)$
(52)

where $Preja$ and $Preje$ represent the probability of removing and creating a simulated phonon, respectively, a is the same as in Eq. (42), and $∫ dE (t)$ is defined in Eq. (30).

## 5 Simulation Results

In this section, we provide a validation of the approach described in Sec. 4. These simulations employ the full-band Monte Carlo method presented in Refs. [18] and [30] with the addition of the techniques presented in this work for the creation and removal of simulated phonon as result of electron–phonon scattering. Moreover, both types of particles are simulated in a synchronous heterogeneous ensemble. This approach employs numerically calculated electron band structure (or the full dispersion for phonons) to calculate all the possible interactions and their respective scattering rates. All this data is precomputed and stored in a compressed (but rather large) look-up table that is loaded at the start of the simulation. This approach allows one to simply choose a transition from the look-up table without the need to search for all the possible final states with each scattering event. Since more than one phonon can occupy a state, the absorption (emission) rate for a phonon state depends on the state occupation number. This effect is accounted for by employing stochastic acceptance/rejection technique to adapt the scattering rate to the run-time simulation distributions. This technique requires to first store the expected scattering maximum rate in the look-up table; then the scattering algorithm extracts the state occupation number from the local phonon distribution at the time of scattering. Finally, the occupation number is used in a stochastic procedure to decide whether to proceed or not with the scattering event.

### 5.1 Conservation of Energy.

This first set of simulations tested the influence of the scattering algorithm and the weight of the simulated particles on energy conservation in the system. It can be inferred from Eqs. (15) and (51) that the error in the energy conservation strongly depends on the values of Wc and Wp. Moreover, these values impact the simulation performances since smaller values require a higher number of particles. It is clear that there must be a tradeoff between the computational burden of the simulation and the error in the conservation of energy $dE$. As explained in the introduction, the value of $Wc$ is determined by the characteristic of the electron simulation; consequently, it is convenient to think about the computational performances and noise level of electrothermal simulation in terms of Wc/Wp. When $Wc/Wp=1$ scattering happens one to one and dE always equal zero. Since typical Monte Carlo simulations have tens of thousands of simulated electron, imposing $Wc/Wp=1$ would requires a number of simulated phonons of the order of 10 × 106. Even with modern high performance computing (HPC) systems, simulating such a large number of particles extremely impractical. Decreasing the value of $Wc/Wp<1$ reduces the number of simulated particles while increasing the value of dE, and hence some sort of trade off must be made between the value dE and the computational demands of the simulation.

The simulated system is an isolated Si cube with an electron-hole pair population initialized with a Maxwellian distribution shifted in energy, and a simulated phonon population that is initially in thermal equilibrium at 100 K, i.e., Bose–Einstein distribution. At time $t=0$, the system is allowed to evolve according to the simulation dynamics as described at the start of this section and in Ref. [18]. During the relaxation process, the instantaneous particle distribution is recorded at intervals of 0.1 fs. Figures 35 show the time evolution of the system in terms of the total energy lost by the charges (dashed line) and the energy gained by the simulated phonons (solid line) in three different cases. Since scattering events are supposed to conserve energy, the two lines should superimpose perfectly in each case. The ratios Wc/Wp are obtained from the simulation output after calculating both Wc and Wp as the ratio between the calculated number of physical particles and the imposed number of simulated particles.

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal

Figure 3 shows what happens when the charge-phonon transient is simply treated in terms of the Poisson probability distribution, as described in Sec. 3.2. This simulation employs a rather large ratio, Wc/Wp, but the difference between the energy gained by the simulated phonons and the energy lost by the simulated charges is a sizable fraction of the total energy exchange.

In Figs. 4 and 5, the electron–phonon interaction is treated as a Hawkes process with different Wc/Wp ratios, both smaller than in the case presented in Fig. 3. In particular, the simulation plotted in Fig. 4 has a ratio Wc/Wp that is a fifth of the case presented in Fig. 3 and is half compared to the case presented in Fig. 5. The importance of this parameter is shown in Eq. (51), which suggests that the average distance between the two curves should scale with the inverse square of Wc/Wp. This means that simulations with a larger Wc/Wp ratio should perform better. Indeed, compared to the case of Fig. 4, the plot in Fig. 5 shows a much closer agreement between the energy gained by the simulated phonons and the energy lost by the simulated charges.

The case in Fig. 3 presents a ratio Wc/Wp at least twice as large than in the cases modeled as Hawkes processes presented in Figs. 4 and 5 and an overall larger error. Indeed, this is the expected behavior predicted in Sec. 4.3. These plots show the importance of using an appropriate scattering algorithm and a sufficiently large ratio Wc/Wp.

### 5.2 Electron Decay.

The full algorithm for the coupled charge-phonon dynamics has been tested by simulating the thermalization process of electron and hole in an isolated cube of Si that is initialized out of equilibrium. In particular, the simulated electrons-hole pairs are initialized with a Maxwellian distribution that is shifted in energy, and the simulated phonons are initialized at 100 K thermal equilibrium before the system is left to relax. The focus of this test is to observe the energy exchange between charges and phonons when modeling absorption and emission in terms of a Hawkes process, as discussed in Sec. 4.4.

Figures 6 and 7 show the evolution of the energy distribution of holes and electrons at time 0.1 fs and 30.1 fs, respectively. The ratio Wc/Wp used is $5.42×10−4$, which shows good energy conservation performance in Fig. 5. These plots represent the evolution of the charge population in the early transient period of the simulation. In Fig. 6, few carriers have scattered and the majority of the charges are still at near the position of the initial energy distribution. It is interesting to note that secondary peaks are already visible for both holes and electrons.

Fig. 6
Fig. 6
Close modal
Fig. 7
Fig. 7
Close modal

The distributions in Fig. 7 show a clear pattern in the carrier decay. In particular, the picture shows evenly spaced density peaks for both electrons and holes. In order to better observe this behavior, Fig. 8 shows an enlarged representation of the electron distribution where the distance between the minor vertical grid lines is set to 0.06 eV, which is roughly the average energy of an optical phonon. From this picture, it is clear that the density peaks are evenly spaced, with the distance between consecutive peaks corresponding to the optical phonon energy. Indeed, this behavior is expected since the emission of optical phonons is the dominant scattering mechanism, and hence the energy relaxation process, for high energy carriers out of equilibrium.

Fig. 8
Fig. 8
Close modal

### 5.3 Electron Thermalization.

In the final set of tests, all the simulated electrons in the isolated Si cube are initialized with a delta function distribution with kinetic energy of 20 meV, while the simulated phonons are initialized at 300 K. This setup tests the code's capability to correctly thermalize the simulated population of particles by both absorption and emission for different ratios of Wc/Wp. Figure 9 compares an equilibrium distribution of simulated electron (circles) used as a reference against distributions obtained after 0.6 fs (dashed line) and 300 fs (solid line) with a ratio $Wc/Wp=5.42×10−3$. At the 0.6 fs mark, few electrons have scattered and the majority of the carriers are still close to the initial energy distribution (dashed line). From the picture, it is clear that the distribution after 300 fs corresponds to the reference distribution, while the initial population was out of equilibrium. In Sec. 5.2, optical phonons were the predominant scattering mechanism for high energy electrons; however, in the present simulation, the acoustic phonons play the critical role in defining the final shape of the electron distribution function.

Fig. 9
Fig. 9
Close modal

Figure 10 compares the time evolution of the electron distribution for the simulation in Fig. 9 (dashed line) and for two other ratios of Wc/Wp corresponding to $2.71×10−3$ (solid line) and $2.71×10−4$ (dotted line). In this picture, the distributions are expressed in terms of their distance from the reference distribution by employing the Jensen–Shannon divergence [51]. We can see in this figure small differences in the transient regime: all three cases reach steady-state at the same time, around the 150 fs mark. After this point, the divergence value is dominated by the statistical noise.

Fig. 10
Fig. 10
Close modal

While all three cases reach steady-state at the same time, one can note that the solid and dashed lines are clearly less noisy during the 150 fs transient. This is due to by the fact that the three simulations have the same Wc but different Wp. Since all simulations have the same starting temperature, larger Wp corresponds to fewer simulated phonons. This causes a difference in the statistical noise when computing the momentum-dependent occupation number.

## 6 Conclusion

This work presented two stochastic approaches to tackle the challenges in modeling electrothermal phenomena within an all particle-based Monte Carlo method. The first approach models the electron–phonon interactions n terms of Poisson processes. This method excels in the steady-state regime but it results in energy conservation issues in the transient regime. This statistical error in energy conservation is addressed by modeling the electron–phonon interactions via the point process formalism. This formalism effectively implements a form of statistical book-keeping, which was found to mitigate the energy conservation issues. Compared to a full book-keeping scheme, this approach presents negligible computational and memory requirement, but it is capable of preserving (in statistical terms) the distribution of the generated and absorbed phonons. We implemented these algorithms within the framework of the particle-based full-band cellular Monte Carlo, and provided some numerical results obtained from the simulation of synchronous and interacting populations of electron and phonon.

## Funding Data

• Air Force Office of Scientific Research (AFOSR) (Grant No. FA9550-16-1-0406; Funder ID: 10.13039/100000181).

• Air Force Research Laboratory (AFRL) (Grant No. FA8650-14- 1-7418; Funder ID: 10.13039/100006602).

## Nomenclature

• a =

A constant number

•
• BZ1 =

first Brillouin zone

•
• $E$ =

energy

•
• E =

expectation

•
• f =

distribution function

•
• h =

the collection of previously occurred events of a process

•
• $H$ =

the collection of previously occurred events of a system

•
• kB =

Boltzmann constant

•
• $M$ =

counting measure

•
• N =

number of simulated particles

•
• $ϱu$ =

random number obtained from the uniform distribution

•
• p =

probability distribution function

•
• P =

Poisson statistic

•
• P =

probability

•
• $R$ =

set of real numbers

•
• S =

step function

•
• t =

time

•
• Var =

variance

•
• W =

number of physical particles represented by a simulated particle

•
• x =

random variable

•
• $Γ$ =

scattering rate

•
• ζ =

real function

•
• η =

number of physical particles

•
• λ =

intensity of a point process

•
• σ =

thermal conductivity

•
• $ϑ$ =

function from R to R+

•
• ω =

frequency

•
• $()a$ =

value for absorption

•
• $()c$ =

value for charges (electrons)

•
• $()e$ =

value for emission

•
• $()p$ =

value for phonons

## Footnotes

1

In measure theory, $H$ is called a filtration or history [35], that is an increasing sequence of σ-algebras on the measurable space.

2

$P(A|B)$ is the probability of the event A happening given the event B has already occurred.

## References

1.
Cercignani
,
C.
,
1972
, “
On the Boltzmann Equation for Rigid Spheres
,”
Transp. Theory Stat. Phys.
,
2
(
3
), pp.
211
225
.10.1080/00411457208232538
2.
Srivastava
,
G. P.
,
1990
,
The Physics of Phonons
,
Taylor & Francis Group
,
New York
.
3.
Peierls
,
R.
,
1955
,
Quantum Theory of Solids. International Series of Monographs on Physics
,
Clarendon Press
,
Oxford, UK
.
4.
Aksamija
,
Z.
,
Hahm
,
H.-S.
, and
Ravaioli
,
U.
,
2008
, “
Emission and Absorption of Phonons in Silicon
,”
Phys. Status Solidi C
,
5
(
1
), pp.
90
93
.10.1002/pssc.200776563
5.
Sinha
,
S.
,
Pop
,
E.
,
Dutton
,
R.
, and
Goodson
,
K.
,
2006
, “
Non-Equilibrium Phonon Distributions in Sub-100 nm Silicon Transistors
,”
ASME J. Heat Transfer
,
128
(
7
), pp.
638
647
.10.1115/1.2194041
6.
Iotti
,
R. C.
,
Rossi
,
F.
,
Vitiello
,
M. S.
,
Scamarcio
,
G.
,
Mahler
,
L.
, and
Tredicucci
,
A.
,
2010
, “
Impact of Nonequilibrium Phonons on the Electron Dynamics in Terahertz Quantum Cascade Lasers
,”
Appl. Phys. Lett.
,
97
(
3
), p.
033110
.10.1063/1.3464977
7.
Shi
,
Y.
, and
Knezevic
,
I.
,
2014
, “
Nonequilibrium Phonon Effects in Midinfrared Quantum Cascade Lasers
,”
J. Appl. Phys.
,
116
(
12
), p.
123105
.10.1063/1.4896400
8.
Lundstrom
,
M.
,
1990
,
Fundamentals of Carrier Transport
, Vol.
10
,
,
.
9.
Jacoboni
,
C.
, and
Reggiani
,
L.
,
1983
, “
The Monte Carlo Method for the Solution of Charge Transport in Semiconductors With Applications to Covalent Materials
,”
Rev. Mod. Phys.
,
55
(
3
), pp.
645
705
.10.1103/RevModPhys.55.645
10.
Joshi
,
A.
, and
Majumdar
,
A.
,
1993
, “
Transient Ballistic and Diffusive Phonon Heat Transport in Thin Films
,”
J. Appl. Phys.
,
74
(
1
), pp.
31
39
.10.1063/1.354111
11.
Chai
,
J. C.
,
Lee
,
H. S.
, and
Patankar
,
S. V.
,
1994
, “
Finite Volume Method for Radiation Heat Transfer
,”
J. Thermophys. Heat Transfer
,
8
(
3
), pp.
419
425
.10.2514/3.559
12.
Mazumder
,
S.
, and
Majumdar
,
A.
,
2001
, “
Monte Carlo Study of Phonon Transport in Solid Thin Films Including Dispersion and Polarization
,”
ASME J. Heat Transfer
,
123
(
4
), pp.
749
759
.10.1115/1.1377018
13.
Chen
,
G.
,
2002
, “
Ballistic-Diffusive Equations for Transient Heat Conduction From Nano to Macroscales
,”
ASME J. Heat Transfer
,
124
(
2
), pp.
320
328
.10.1115/1.1447938
14.
Lacroix
,
D.
,
Joulain
,
K.
, and
Lemonnier
,
D.
,
2005
, “
Monte Carlo Transient Phonon Transport in Silicon and Germanium at Nanoscales
,”
Phys. Rev. B
,
72
(
6
), p.
064305
.10.1103/PhysRevB.72.064305
15.
Kukita
,
K.
, and
Kamakura
,
Y.
,
2013
, “
Monte Carlo Simulation of Phonon Transport in Silicon Including a Realistic Dispersion Relation
,”
J. Appl. Phys.
,
114
(
15
), p.
154312
.10.1063/1.4826367
16.
Mei
,
S.
,
Maurer
,
L.
,
Aksamija
,
Z.
, and
Knezevic
,
I.
,
2014
, “
Full-Dispersion Monte Carlo Simulation of Phonon Transport in Micron-Sized Graphene Nanoribbons
,”
J. Appl. Phys.
,
116
(
16
), p.
164307
.10.1063/1.4899235
17.
Landon
,
C. D.
, and
,
N. G.
,
2014
, “
Deviational Simulation of Phonon Transport in Graphene Ribbons With Ab Initio Scattering
,”
J. Appl. Phys.
,
116
(
16
), p.
163502
.10.1063/1.4898090
18.
Sabatti
,
F. F.
,
Goodnick
,
S. M.
, and
Saraniti
,
M.
,
2016
, “
Simulation of Phonon Transport in Semiconductors Using a Population-Dependent Many-Body Cellular Monte Carlo Approach
,”
ASME J. Heat Transfer
,
139
(
3
), p.
032002
.10.1115/1.4035042
19.
Vecchi
,
M. C.
, and
Rudan
,
M.
,
1998
, “
Modeling Electron and Hole Transport With full-Band Structure Effects by Means of the Spherical-Harmonics Expansion of the BTE
,”
IEEE Trans. Electron Devices
,
45
(
1
), pp.
230
238
.10.1109/16.658836
20.
Murthy
,
J. Y.
,
Narumanchi
,
S. V.
,
José
,
P.-G.
,
Wang
,
T.
,
Ni
,
C.
, and
Mathur
,
S. R.
,
2005
, “
Review of Multiscale Simulation in Submicron Heat Transfer
,”
Int. J. Multiscale Comp. Eng.
,
3
(
1
), p.
5
.10.1615/IntJMultCompEng.v3.i1.20
21.
Zhao
,
Y.
,
Watling
,
J.
,
Kaya
,
S.
,
Asenov
,
A.
, and
Barker
,
J.
,
2000
, “
Drift Diffusion and Hydrodynamic Simulations of Si/SiGe P-Mosfets
,”
Mater. Sci. Eng.: B
,
72
(
2–3
), pp.
180
183
.10.1016/S0921-5107(99)00508-5
22.
Rupp
,
K.
,
Jungemann
,
C.
,
Hong
,
S.-M.
,
Bina
,
M.
,
Grasser
,
T.
, and
Jüngel
,
A.
,
2016
, “
A Review of Recent Advances in the Spherical Harmonics Expansion Method for Semiconductor Device Simulation
,”
J. Comput. Electron.
,
15
(
3
), pp.
939
958
.10.1007/s10825-016-0828-z
23.
Jacoboni
,
C.
, and
Lugli
,
P.
,
1989
,
The Monte Carlo Method for Semiconductor Device Equations
,
Springer–Verlag
,
Wien, New York
.
24.
Rowlette
,
J. A.
, and
Goodson
,
K. E.
,
2008
, “
Fully Coupled Nonequilibrium Electron-Phonon Transport in Nanometer-Scale Silicon Fets
,”
IEEE Trans. Electron Devices
,
55
(
1
), pp.
220
232
.10.1109/TED.2007.911043
25.
Raleva
,
K.
,
Vasileska
,
D.
,
Goodnick
,
S. M.
, and
Nedjalkov
,
M.
,
2008
, “
Modeling Thermal Effects in Nanodevices
,”
IEEE Trans. Electron Devices
,
55
(
6
), pp.
1306
1316
.10.1109/TED.2008.921263
26.
Kamakura
,
Y.
,
Mori
,
N.
,
Taniguchi
,
K.
,
Zushi
,
T.
, and
Watanabe
,
T.
,
2010
, “
Coupled Monte Carlo Simulation of Transient Electron-Phonon Transport in Nanoscale Devices
,”
International Conference on Simulation of Semiconductor Processes and Devices
, Bologna, Italy, Sept. 6–8, pp.
89
92
27.
Mohamed
,
M.
,
Aksamija
,
Z.
,
Vitale
,
W.
,
Hassan
,
F.
,
Park
,
K.-H.
, and
Ravaioli
,
U.
,
2014
, “
A Conjoined Electron and Thermal Transport Study of Thermal Degradation Induced During Normal Operation of Multigate Transistors
,”
IEEE Trans. Electron Devices
,
61
(
4
), pp.
976
983
.10.1109/TED.2014.2306422
28.
Fischetti
,
M. V.
, and
Laux
,
S. E.
,
1988
, “
Monte Carlo Analysis of Electron Transport in Small Semiconductor Devices Including Band-Structure and Space-Charge Effects
,”
Phys. Rev. B
,
38
(
14
), p.
9721
.10.1103/PhysRevB.38.9721
29.
Lewis
,
P. A. W.
, and
Shedler
,
G. S.
,
1979
, “
Simulation of Nonhomogeneous Poisson Processes by Thinning
,”
Nav. Res. Logistics Q.
,
26
(
3
), pp.
403
413
.10.1002/nav.3800260304
30.
Saraniti
,
M.
, and
Goodnick
,
S.
,
2000
, “
Hybrid Full-Band Cellular Automaton/Monte Carlo Approach for Fast Simulation of Charge Transport in Semiconductors
,”
IEEE Trans. Electron Devices
,
47
(
10
), pp.
1909
1915
.10.1109/16.870571
31.
Ross
,
S. M.
,
1995
,
Stochastic Processes
,
Wiley
,
New York
.
32.
Streit
,
R. L.
,
2010
,
The Poisson Point Process
,
Springer
,
Boston, MA
.
33.
,
J.
, and
Russell
,
W. R.
,
1969
, “
Rules for Ordering Uncertain Prospects
,”
Am. Econ. Rev.
,
59
(
1
), pp.
25
34
.http://www.jstor.org/stable/1811090
34.
Ashcroft
,
N. W.
, and
Mermin
,
N. D.
,
1981
,
Solid State Physics
,
Holt–Saunders International Editions
,
Tokyo, Japan
.
35.
Daley
,
D. J.
, and
Vere-Jones
,
D.
,
2007
,
An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods
,
Springer
,
New York
.
36.
Cox
,
D.
, and
Isham
,
V.
,
1980
,
Point Processes. Monographs on Statistics & Applied Probability
,
Chapman and Hall/CRC
,
London
.
37.
Aleksandr
,
I.
, and
Khinchin
,
A.
,
1949
,
Mathematical Foundations of Statistical Mechanics
,
Courier Corporation
,
New York
.
38.
Daley
,
D. J.
, and
Vere-Jones
,
D.
,
2007
,
An Introduction to the Theory of Point Processes: Volume II: General Theory and Structure
,
Springer-Verlag
,
New York
.
39.
Hawkes
,
A. G.
,
1971
, “
Spectra of Some Self-Exciting and Mutually Exciting Point Processes
,”
Biometrika
,
58
(
1
), pp.
83
90
.10.1093/biomet/58.1.83
40.
Cox
,
D. R.
,
1955
, “
Some Statistical Methods Connected With Series of Events
,”
J. R. Stat. Soc. Ser. B (Methodological)
,
17
(
2
), pp.
129
164
.10.1111/j.2517-6161.1955.tb00188.x
41.
Ogata
,
Y.
,
1988
, “
Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes
,”
J. Am. Stat. Assoc.
,
83
(
401
), pp.
9
27
.10.1080/01621459.1988.10478560
42.
Mohler
,
G. O.
,
Short
,
M. B.
,
Brantingham
,
P. J.
,
Schoenberg
,
F. P.
, and
Tita
,
G. E.
,
2011
, “
Self-Exciting Point Process Modeling of Crime
,”
J. Am. Stat. Assoc.
,
106
(
493
), pp.
100
108
.10.1198/jasa.2011.ap09546
43.
Tita
,
G.
, and
Ridgeway
,
G.
,
2007
, “
The Impact of Gang Formation on Local Patterns of Crime
,”
J. Res. Crime Delinquency
,
44
(
2
), pp.
208
237
.10.1177/0022427806298356
44.
Brémaud
,
P.
, and
Massoulié
,
L.
,
1996
, “
Stability of Nonlinear Hawkes Processes
,”
Ann. Probab.
,
24
(3), pp.
1563
1588
.https://www.jstor.org/stable/2244985
45.
Bilotta
,
S.
,
Pergola
,
E.
,
Pinzani
,
R.
, and
Rinaldi
,
S.
,
2015
,
Recurrence Relations, Succession Rules, and the Positivity Problem
,
Springer International Publishing
,
Cham, Switzerland
, pp.
499
510
.
46.
Sharma
,
R.
,
Gupta
,
M.
, and
Kapoor
,
G.
,
2010
, “
Some Better Bounds on the Variance With Applications
,”
J. Math. Inequalities
, (
3
), pp.
355
363
.10.7153/jmi-04-32
47.
Jensen
,
J. L.
, and
Möller
,
J.
,
1991
, “
Pseudolikelihood for Exponential Family Models of Spatial Point Processes
,”
Ann. Appl. Probab.
,
1
(
3
), pp.
445
461
.10.1214/aoap/1177005877
48.
Propp
,
J. G.
, and
Wilson
,
D. B.
,
1996
, “
Exact Sampling With Coupled Markov Chains and Applications to Statistical Mechanics
,”
Random Struct. Algorithms
,
9
(
1–2
), pp.
223
252
.10.1002/(SICI)1098-2418(199608/09)9:1/2<223::AID-RSA14>3.0.CO;2-O
49.
Ogata
,
Y.
,
1981
, “
On Lewis' Simulation Method for Point Processes
,”
IEEE Trans. Inf. Theory
,
27
(
1
), pp.
23
31
.10.1109/TIT.1981.1056305
50.
Billera
,
L. J.
, and
Diaconis
,
P.
,
2001
, “
A Geometric Interpretation of the Metropolis-Hastings Algorithm
,”
Stat. Sci.
,
16
(
4
), pp.
335
339
.10.1214/ss/1015346318
51.
Lin
,
J.
,
1991
, “
Divergence Measures Based on the Shannon Entropy
,”
IEEE Trans. Inf. Theory
,
37
(
1
), pp.
145
151
.10.1109/18.61115