\section{Introduction}\label{sec:intro}
Discovering hidden rules from sequential observations has been an essential topic in machine learning, with a large variety of applications such as physics simulation \citep{sanchez2020learning}, autonomous driving \citep{diehl2019graph}, ECG analysis \citep{golany2021ecg} and event analysis \citep{chen2020learning}, to name a few. A standard scheme is to consider sequential data at each timestamp to be holistic and homogeneous under some ideal assumptions (i.e., only the temporal behavior of one entity is involved in a sequence), under which data/observation is treated as a collection of slices at a different time from a unified system.
A series of sequential learning models fall into this category, including variants of recurrent neural networks (RNNs) \citep{cho2014properties,hochreiter1997long}, neural differential equations (DEs) \citep{chen2018neural,kidger2020neural,rusch2021unicornn,zhu2021neural} and spatial/temporal attention-based approaches \citep{vaswani2017attention,fan2019multi,song2017end}.
These variants fit well into the scenarios agreeing with the aforementioned assumptions and are proved effective in learning or modeling for relatively simple applications with clean data sources.

In the real world, a system may not only describe a single and holistic entity but also consist of several \emph{distinguishable} interacting but simple subsystems, where each subsystem corresponds to a physical entity. For example, we can think of the movement of a solar system as the mixture of distinguishable subsystems of the sun and surrounding planets, while interactions between these celestial bodies over time are governed by the laws of gravity. Back centuries ago, physicists and astronomers made enormous efforts to discover the rule of celestial movements from the records of every single body and eventually delivered the neat yet elegant differential equations (DEs) depicting principles of moving bodies and interactions therein. Likewise, nowadays, researchers also developed a series of machine learning models for sequential data with distinguishable partitions \citep{qin2017dual}. Two widely adopted strategies for learning the interactions between subsystems are graph neural networks \citep{iakovlev2020learning,ha2021unraveling,kipf2018neural,yildiz2022learning,xhonneux2020continuous} and attention mechanism \citep{vaswani2017attention,lu2019understanding,goyal2019recurrent}, while the interactions are typically encoded with ``messages'' between nodes and pair-wise ``attention scores'', respectively.

It is worth noting an even more difficult scenario:
\begin{itemize}
    \item \emph{The data/observation is so cluttered that cannot be readily distinguished into separate parts}.
\end{itemize}
This can be either due to the way of data collection (e.g., videos consisting of multiple objects) or because there are no explicit physical entities originally (e.g., weather time series).
To tackle this, a fair assumption can be introduced that complex observations can be decoupled into several relatively independent modules in the feature space, where each module corresponds to a \emph{latent entity}. Latent entities may not have exact physical meanings, but learning procedures can greatly benefit from such decoupling, as this assumption can be viewed as strong regularization to the system. This assumption has been successfully incorporated in several models for learning from \emph{regularly} sampled sequential data by emphasizing ``independence'' to some extent between channels or groups in the feature space \citep{li2018independently,yu2020rhyrnn,goyal2019recurrent,madan2021fast}.
Another successful counterpart in parallel benefiting from this assumption is transformer \citep{vaswani2017attention} which stacks multiple layers of self-attention and point-wise feedforward networks. In transformers, each attention head can be viewed as a relatively independent module, and interaction happens throughout the head re-weighting procedure following the attention scores. \citet{lu2019understanding} presented an interpretation from a dynamic point of view by regarding a basic layer in the transformer as one step of integration governed by differential equations derived from interacting particles. \citet{vuckovic2020mathematical} extended this interpretation with more solid mathematical support by viewing the forward pass of the transformer as applying successive Markov kernels in a particle-based dynamic system.

We note, however, despite the ubiquity of this setting, there is barely any previous investigation focusing on learning for \emph{irregularly sampled} and \emph{cluttered} sequential data. The aforementioned works either fail to handle the irregularity \citep{goyal2019recurrent,li2018independently} or neglect the independence/modularity assumption in the latent space \citep{chen2018neural,kidger2020neural}.
In this paper, inspired by recent advances of neural controlled dynamics \citep{kidger2020neural} and novel interpretation of attention mechanism \citep{vuckovic2020mathematical}, we take a step to propose an effective approach addressing this problem under the dynamic setting. To this end, our approach explicitly learned to decouple a complex system into several latent sub-systems and utilizes an additional meta-system capturing the evolution of interactions over time.
Specifically, taking into account the meta-system capturing interactions evolving in a constrained set (e.g., simplex), we further characterized such interactions using projected differential equations (ProjDEs) with neural-friendly projection operators. We argued our \textbf{contributions} as follows:
\begin{itemize}
    \item We provide a novel modeling strategy for sequential data from a system decoupling perspective;
    \item We propose a novel and natural interpretation of evolving interactions as a ProjDE-based meta-system, with insights into projection operators in the sense of Bregman divergence;
    \item Our approach is parameter-insensitive and more compatible with other modules and data, thus being flexible to be integrated into various tasks.
\end{itemize}
Extensive experiments were conducted on either regularly or irregularly sampled sequential data, including both synthetic and real-world settings. It was observed that our approach achieved prominent performance compared to the state-of-the-art on a wide spectrum of tasks.

\section{Related Work}

\paragraph{Sequential learning.} Traditionally, learning with sequential data can be performed using variants of recurrent neural networks (RNNs) \citep{hochreiter1997long,cho2014properties,li2018independently} under the Markov setting. While such RNNs are generally designed for regular sampling frequency, a more natural line of counterparts lies in the continuous time domain allowing irregularly sampled time series as input. As such, a variety of RNN-based methods are developed by introducing exponential decay on observations \citep{che2018recurrent,mei2017neural}, incorporating an underlying Gaussian process \citep{li2016scalable,futoma2017learning}, or integrating some latent evolution under ODEs \citep{rubanova2019latent,de2019gru}. A seminal work interpreting forward passing in neural networks as an integration of ODEs was proposed in \citet{chen2018neural}, followed by a series of relevant works \citep{liu2019neural,li2020scalable,dupont2019augmented}. As integration over ODEs allows for arbitrary step length, it is natural modeling of irregular time series and proved powerful in many machine learning tasks (e.g., bioinformatics \citep{golany2021ecg}, physics \citep{nardini2021learning} and computer vision \citep{park2021vid}). \cite{kidger2020neural} studied a more effective way of injecting observations into the system via a mathematical tool called Controlled differential Equation, achieving state-of-the-art performance on several benchmarks. Some variants of neural ODEs have also been extended to discrete structure \citep{chamberlain2021grand,xhonneux2020continuous} and non-Euclidean setting \citep{chamberlain2021beltrami}.

\paragraph{Learning with independence.} Independence or modular property serves as strong regularization or prior in some learning tasks under static setting \citep{wang2020learning,liu2020independence}. In the sequential case, some early attempts over RNNs emphasized implicit ``independence'' in the feature space between dimensions or channels \citep{li2018independently,yu2020rhyrnn}. As independence assumption commonly holds in vision tasks (with distinguishable objects), \citet{pang2020complex,li2020hoi} proposed video understanding schemes by decoupling the spatiotemporal patterns. For a more generic case where the observations are collected without any prior, \citet{goyal2019recurrent} devised a sequential learning scheme called recurrent independence mechanism (RIM),
and its generalization ability was extensively studied in \citet{madan2021fast}. \citet{lu2019understanding} investigated self-attention mechanism \citep{vaswani2017attention} and interpreted it as a nearly independent multi-particle system with interactions therein. \citet{vuckovic2020mathematical} further provided more solid mathematical analysis with the tool of Markov kernel. The study of such a mechanism in the dynamical setting was barely observed.

\textbf{Learning dynamics under constraints.} It is practically significant as a series of real-world systems evolve within some manifolds, such as fluid \cite{vinuesa2022enhancing}, coarse-grained dynamics \cite{kaltenbach2020incorporating}, and molecule modeling \cite{chmiela2020accurate}. While some previous research incorporates constraints from a physical perspective \cite{kaltenbach2020incorporating,linot2020deep}, an emerging line is empowered by machine learning to integrate or even discover the constraints \cite{kolter2019learning,lou2020neural,goldt2020modeling}. To ensure a system evolves in constraints, efficient projections or pseudo-projections are required, about which Bregman divergence provides rich insights \cite{martins2016softmax,krichene2015efficient,lim2016efficient}. Despite these results, to our best knowledge, there is barely any related investigation about neural-friendly projections.

\section{Methodology}

\subsection{Background}\label{sec:background}
In this section, we briefly review three aspects related to our approach. Our approach is built upon the basic sub-system derived from \emph{Neural Controlled Dynamics} \citep{kidger2020neural}, while the interactions are modeled at an additional meta-system analogous to \emph{Self-attention} \citep{lu2019understanding,vuckovic2020mathematical}, and further interpreted and generalized using the tool of \emph{Projected Differential Equations} \citep{dupuis1993dynamical}.

\paragraph{Neural Controlled Dynamics.} Continuous-time dynamics can be expressed using differential equations $\mathbf{z}'(t)=d\mathbf{z}/dt =f(\mathbf{z}(t),t)$, where $\mathbf{z}\in\mathbb{R}^d$ and $t$ are a $d$-dimension state and the time, respectively. Function $f:\mathbb{R}^d\times\mathbb{R}_+\rightarrow \mathbb{R}^d$ governs the evolution of the dynamics. Given the initial state $\mathbf{z}(t_0)$, the state at any time $t_1$ can be evaluated with:
\begin{equation}\label{eq: ode}
    \mathbf{z}(t_1)=\mathbf{z}(t_0)+\int_{t_0}^{t_1}f(\mathbf{z}(s),s)\mathrm{d}s
\end{equation}
In practice, we aim at learning the dynamics from a series of observations or controls $\{\mathbf{x}(t_k)\in\mathbb{R}^b | k=0,1,...\}$ by parameterizing the dynamics with $f_\theta(\cdot)$ where $\theta$ is the unknown parameter to be learned. Thus, a generic dynamics incorporating outer signals $\mathbf{x}$ can be written as:
\begin{equation}\label{eq: node}
    \mathbf{z}(t_1)=\mathbf{z}(t_0)+\int_{t_0}^{t_1}f_\theta(\mathbf{z}(s),\mathbf{x}(s),s)\mathrm{d}s
\end{equation}
Rather than directly injecting $\mathbf{x}$ as in Eq.~(\ref{eq: node}), Neural Controlled Differential Equation (CDE) proposed to deal with outer signals with a Riemann–Stieltjes integral \citep{kidger2020neural}:
\begin{equation}\label{eq: ncde}
    \mathbf{z}(t_1)=\mathbf{z}(t_0)+\int_{t_0}^{t_1} \mathbf{F}_\theta(\mathbf{z}(s))\mathbf{x}'(s)\mathrm{d}s
\end{equation}
where $\mathbf{F}_\theta:\mathbb{R}^d \rightarrow \mathbb{R}^{d\times b}$ is a learnable vector field and $\mathbf{x}'(s)=\mathrm{d}\mathbf{x}/\mathrm{d}s$ is the derivative of signal $\mathbf{x}$ w.r.t. time $s$, thus ``$\mathbf{F}_\theta(\mathbf{z}(s))\mathbf{x}'(s)$'' is a matrix-vector multiplication. During implementation, \citet{kidger2020neural} argued that a simple cubic spline interpolation on $\mathbf{x}$ allows dense calculation of $\mathbf{x}'(t)$ at any time $t$ and exhibits promising performance. In \cite{kidger2020neural}, it is also mathematically shown that incorporating observations/controls following Eq. (\ref{eq: ncde}) is with greater representation ability compared to Eq. (\ref{eq: node}), hence achieving state-of-the-art performance on several public tasks.

\paragraph{Self-attention.} It is argued in \citet{lu2019understanding,vuckovic2020mathematical} that a basic unit in Transformer \citep{vaswani2017attention} with self-link consisting of one self-attention layer and point-wise feedforward layer amounts to simulating a multi-particle dynamical system. Considering such a layer with $n$ attention-heads (corresponding to $n$ particles), given an attention head index $i\in\{1,2,...,n\}$, the update rule of the $i$th unit at depth $l$ reads:
\begin{subequations}\label{eq: att}
    \begin{alignat}{2}                  &\tilde{\mathbf{z}}_{l,i}=\mathbf{z}_{l,i}+\mathrm{MHAtt}_{W^l_{\mathrm{att}}}\left(\mathbf{z}_{l,i},\left[\mathbf{z}_{l,1},...,\mathbf{z}_{l,n}\right]\right) \label{eq: att_dif} \\
    &\mathbf{z}_{l+1,i}=\tilde{\mathbf{z}}_{l,i}+\mathrm{FFN}_{W^l_{\mathrm{ffn}}}\left(\tilde{\mathbf{z}}_{l,i}\right) \label{eq: att_conv}
    \end{alignat}
\end{subequations}
where $\mathrm{MHAtt}_{W^l_{\mathrm{att}}}$ and $\mathrm{FFN}_{W^l_{\mathrm{ffn}}}$ are multi-head attention layer and feedforward layer with parameters $W^l_{\mathrm{att}}$ and $W^l_{\mathrm{ffn}}$, respectively. Eq.~(\ref{eq: att}) can then be interpreted as an interacting multi-particle system:
\begin{equation}\label{eq: att_overall}
    \frac{\mathrm{d}\mathbf{z}_i(t)}{\mathrm{d}t}=F(\mathbf{z}_i(t),[\mathbf{z}_1(t),...,\mathbf{z}_n(t)],t)+G(\mathbf{z}_i(t))
\end{equation}
where function $F$ corresponding to Eq.~(\ref{eq: att_dif}) represents the diffusion term and $G$ corresponding to Eq.~(\ref{eq: att_conv}) stands for the convection term. Notably, the attention score obtained via $\mathrm{softmax}$ in Eq.~(\ref{eq: att_dif}) is regarded as a Markov kernel.
Readers are referred to \citet{lu2019understanding,vuckovic2020mathematical} for more details.

\paragraph{Projected DEs.} It is a tool depicting the behavior of dynamics where solutions are constrained within a (convex) set. Concretely, given a closed polyhedral $\mathcal{K}\subset \mathbb{R}^n$ and a mapping $H:\mathcal{K}\rightarrow \mathbb{R}^n$, we can introduce an operator $\Pi_{\mathcal{K}}: \mathbb{R}^n\times\mathcal{K} \rightarrow \mathbb{R}^n$ which is defined by means of directional derivatives as:
\begin{equation}\label{eq: derivative}
    \Pi_{\mathcal{K}}(\mathbf{a},H(\mathbf{a}))=\lim_{\alpha\rightarrow 0_+}\frac{P_{\mathcal{K}}(\mathbf{a}+\alpha H(\mathbf{a}))-\mathbf{a}}{\alpha}
\end{equation}
where $P_{\mathcal{K}}(\cdot)$ is a projection onto $\mathcal{K}$ in terms of Euclidean distance:
\begin{equation}\label{eq:proj-l2}
    \lVert P_{\mathcal{K}}(\mathbf{a})-\mathbf{a}\rVert _2=\inf_{\mathbf{y}\in\mathcal{K}}\lVert\mathbf{y}-\mathbf{a}\rVert_2
\end{equation}
Intuitively, Eq.~(\ref{eq: derivative}) pictures the dynamics of $\mathbf{a}$ driven by function $H$, but constrained within $\mathcal{K}$. Whenever $\mathbf{a}$ reaches beyond $\mathcal{K}$, it would be projected back using Eq.~(\ref{eq:proj-l2}). By extending Eq.~(\ref{eq: derivative}), \cite{dupuis1993dynamical,zhang1995stability} considered the projected differential equations as follows:
\begin{equation}
    \frac{\mathrm{d}\mathbf{a}(t)}{\mathrm{d}t}=\Pi_{\mathcal{K}}(\mathbf{a},H(\mathbf{a}))
\end{equation}
which allows for discontinuous dynamics on $\mathbf{a}$.

\subsection{Learning to decouple}\label{sec:method}
Our method is built upon the assumption that cluttered sequential observations are composed of several relatively independent sub-systems and, therefore, explicitly learns to decouple them as well as to capture the mutual interactions with a meta-system in parallel.
Let the cluttered observations/controlls be $\mathbf{c}(t)\in\mathbb{R}^k$ at time $t$ for $t=1,...,T$, where $T$ is the time horizon. We employ $k$ distinct mappings with learnable parameters (e.g., MLP) to obtain respective controls to each sub-system: $\mathbf{x}_i(t)=p_i(\mathbf{c}(t))\in\mathbb{R}^m$ for $i=1,...,n$. A generic dynamics of the proposed method can be written as:
\label{sec:method}
\begin{subequations}\label{eq:dec}
    \begin{alignat}{2}
        \frac{\mathrm{d}\mathbf{z}_i(t)}{\mathrm{d}t}&=f_{i}\left(\mathbf{z}_i(t),\left[\mathbf{z}_1(t),...,\mathbf{z}_n(t)\right],\mathbf{x}_i(t),\mathbf{a}(t)\right) \label{eq:dec_sub}\\
        \frac{\mathrm{d}\textbf{a}(t)}{\mathrm{d}t}&=\Pi_\mathcal{S}\left(\mathbf{a}(t),g(\mathbf{a}(t),\left[\mathbf{z}_1
        (t),...,\mathbf{z}_n(t)\right])\right) \label{eq:dec_att}
    \end{alignat}
\end{subequations}
where Eq.~(\ref{eq:dec_sub}) and Eq.~(\ref{eq:dec_att}) refer to the $i$th sub-system describing the evolution of a single latent entity and meta-system depicting the interactions, respectively. $\mathbf{z}_i(t)\in\mathbb{R}^q$ is the hidden state for the $i$th subsystem, and $\mathbf{a}$ is a tensor governs the dynamics of the interactions. Here $\Pi_\mathcal{S}(\cdot)$ is a projection operator, which projects the evolving trajectory into set $\mathcal{S}$. We introduce such an operator as it is assumed that interactions among latent entities should be constrained following some latent manifold structure. $f_i(\cdot)$ and $g(\cdot)$ are both learnable functions and also the essential roles for capturing the underlying complex dynamics.

\begin{remark}
It is seen the projection operator $\Pi_\mathcal{S}(\cdot)$ and the set $\mathcal{S}$ play important roles in Eq.~(\ref{eq:dec_att}). For $\Pi_\mathcal{S}(\cdot)$, while previous works of ProjDEs only consider L2-induced projection, we propose novel interpretation and extension under Bregman divergence. For $\mathcal{S}$, we consider a probabilistic simplex following the setting in \citet{lu2019understanding,vuckovic2020mathematical}, though it can be any polyhedral.
\end{remark}

According to Eq.~(\ref{eq:dec}), we fully decouple a complex system into several components. Although we found some decoupling counterparts in the context of RNNs \citep{li2018independently,yu2020rhyrnn} and attention-like mechanism \citep{lu2019understanding,goyal2019recurrent}, their decoupling could not be applied to our problem.
We elaborate on the details of implementing Eq.~(\ref{eq:dec}) in the following.

\paragraph{Learning sub-systems.} Sub-systems corresponding to the latent entities seek to model relatively independent dynamics separately. Specifically, we employ the way of integrating $\mathbf{x}_i$s into Eq.~(\ref{eq:dec_sub}) in a controlled dynamical fashion as in the state-of-the-art method \citep{kidger2020neural}:
\begin{equation}\label{eq: sub_generic}
    \mathrm{d}\mathbf{z}_i(t)= \mathbf{F}_i \left( \mathbf{z}_i(t),\mathbf{a}(t),\left[\mathbf{z}_1(t),...,\mathbf{z}_n(t) \right] \right) \mathrm{d}\mathbf{x}_i(t)
\end{equation}
where $\mathbf{F}_i(\cdot)\in\mathbb{R}^{q\times m}$ is a learnable vector field. Concretely, if we let $\mathbf{z}(t)=\left[\mathbf{z}_i(t),...,\mathbf{z}_n(t) \right]$ be the tensor collecting all sub-systems, the $i$th sub-system in a self-attention fashion reads:
\begin{equation}\label{eq: sub}
    \mathrm{d}\mathbf{z}_i(t)= \mathbf{F}( \left[ \mathbf{A}(t) \cdot \mathbf{z}(t) \right] _i) \mathrm{d}\mathbf{x}_i(t)
\end{equation}
where $[\cdot]_i$ takes the $i$th slice from a tensor. Note timestamp $t$ can be arbitrary, resulting in irregularly sampled sequential data. To address this, we follow the strategy in \citet{kidger2020neural} by performing cubic spline interpolation on $\mathbf{x}_i$ over observed timestamp $t$, resulting in $\mathbf{x}_i(t)$ at dense time $t$. Note that for all sub-systems, different from Eq.~(\ref{eq: sub_generic}) we utilize an identical function/network $\mathbf{F}(\cdot)$ as in Eq.~(\ref{eq: sub}), but with different control sequence $\mathbf{x}_i(t)=p_i(\mathbf{c}(t))$. Since in our implementation, $p_i(\cdot)$ is a lightweight network such as MLP, this can significantly reduce the parameter size.

\paragraph{Learning interactions.} In our approach, interactions between latent entities are modeled separately as another meta-system. This is quite different from some related methods \citep{lu2019understanding, vuckovic2020mathematical} where sub-systems and interactions are treated as one holistic step of forward integration. For the meta-system describing the interactions in Eq.~(\ref{eq:dec_att}), two essential components are involved: domain $\mathcal{S}$ and the projection operator $\Pi$. In the context of ProjDEs, a system is constrained as $\mathbf{a}(t)\in\mathcal{S}$ for any $t$. In terms of interactions, a common choice of $\mathcal{S}$ is the stochastic simplex which can be interpreted as a transition kernel \citep{vuckovic2020mathematical}. We allow follow this setting by defining $\mathcal{S}$ be a row-wise stochastic $(n-1)$-simplices:
\begin{equation}\label{eq: simplex}
    \mathcal{S}\triangleq\{\mathbf{A}\in\mathbb{R}^{n\times n}|\mathbf{A1=1},\mathbf{A}_{ij}\geq 0\}
\end{equation}
where $\mathbf{1}$ is a vector with all $1$ entries. $\mathbf{A}=\mathrm{mat}(\mathbf{a})$ is a $n\times n$ matrix. In the sequel, we will use the notation $\mathbf{A}$ throughout.
Thus the meta-system capturing the interactions can be implemented as follows:
\begin{equation}\label{eq: att_implement}
    \frac{d\mathbf{A}(t)}{dt}=\Pi_\mathcal{S}\left(\mathbf{A}(t),g(\mathbf{A}(t),\left[\mathbf{z}_1(t),...,\mathbf{z}_n(t)\right])\right)
\end{equation}
For the projection operator, we consider two versions shown in Eq.~(\ref{eq:projs}). In Eq.~(\ref{eq:pj_soft}), we give a row-wise projection onto the $(n-1)$-simplex with entropic regularization \citep{amos2019differentiable}, which has a well-known closed-form solution $\mathrm{softmax}(\cdot)$ appearing in attention mechanism. In Eq.~(\ref{eq:pj_sparse}), we adopt a standard L2-induced projection identical to Eq.~(\ref{eq:proj-l2}), which leads to sparse solutions \cite{wainwright2008graphical}. Intuitively, the projection of a point onto a simplex in terms of L2 distance tends to lie on a facet or a vertex of a simplex, thus being sparse.
\begin{subequations}\label{eq:projs}
    \begin{equation}\label{eq:pj_soft}
    P_\mathcal{S}^{\text{soft}}(\mathbf{A}_{j,:})=\argmin_{\mathbf{B}\in\mathcal{S}}\mathbf{A}_{j,:}^\top \mathbf{B}_{:,j}-\mathbb{H}^{\text{entr}
    }(\mathbf{B}_{:,j})
    \end{equation}
    \begin{equation}\label{eq:pj_sparse}
    \begin{split}
        P_\mathcal{S}^{\text{sparse}}(\mathbf{A}_{j,:})&=\argmin_{\mathbf{B}\in\mathcal{S}}\mathbf{A}_{j,:}^\top \mathbf{B}_{:,j}-\mathbb{H}^{\text{gini}}(\mathbf{B}_{:,j}) \\
        &=\argmin_{\mathbf{B}\in\mathcal{S}}|\mathbf{A}_{j,:}-\mathbf{B}_{:,j}|^2
    \end{split}
    \end{equation}
\end{subequations}
where $\mathbb{H}^{\text{entr}}(\cdot)$ and $\mathbb{H}^{\text{gini}}(\mathbf{y})= \frac{1}{2} \sum_i\mathbf{y}_i (\mathbf{y}_i - 1)$ are the standard entropy and the gini-entropy, respectively. $\mathbf{A}_{j,:}$ and $\mathbf{B}_{:,j}$ are the $i$th row and column of $\mathbf{A}$ and $\mathbf{B}$, respectively.
While the solution to Eq.~(\ref{eq:pj_soft}) is $\mathrm{softmax}(\mathbf{A})$, Eq.~(\ref{eq:pj_sparse}) also has closed-form solution shown in Appendix \ref{app:soft_sparse}.

Comparing Eq.~(\ref{eq:pj_soft}) to the standard Euclidean projection in Eq.~(\ref{eq:pj_sparse}), we note the entropic regularization $\mathbb{H}(\cdot)$ in Eq.~(\ref{eq:pj_soft}) allows for a smoother trajectory by projecting any $\mathbf{A}$ into the interior of $(n-1)$-simplex.
We visualize the two versions of projections in Eq.~(\ref{eq:projs}) onto $1$-simplex from some random points in Fig.~\ref{fig:softVSsparse}. One can readily see that Eq.~(\ref{eq:pj_sparse}) is an exact projection such that points far from the simplex are projected onto the boundary. However, $\mathrm{softmax}$ is smoother by projecting all points onto a relative interior of 1-simplex without sudden change. In the context of Bregman divergence, different distances can facilitate efficient convergence under different ``L-relative smoothness'' \cite{dragomir2021fast}, which can potentially accelerate the learning of dynamics. We leave this to our future work.

We further discuss some neural-friendly features of Eq.~(\ref{eq:pj_soft}) and (\ref{eq:pj_sparse}) facilitating the neural computation:

\textbf{(1)} First, the neural computational graph can be simplified using projection Eq.~(\ref{eq:pj_soft}). Though Eq.~(\ref{eq: att_implement}) using projection Eq.~(\ref{eq:pj_soft}) defines a projected dynamical system directly on $\mathbf{A}$, we switch to update the system using $\mathbf{L}$ as follows, which is considered to further ease the forward integration. This is achieved by instead modeling the dynamics of the feature before fed into $\mathrm{softmax}(\cdot)$:
\begin{subequations}\label{eq: transform}
\begin{small}
    \begin{alignat}{2}
        \mathbf{A}(t) &= \text{Softmax}(\mathbf{L}(t)) \\
        \mathbf{L}(t) &= \mathbf{L}(0) + \int_0^t \frac{\mathrm{d}}{\mathrm{d}s} \frac{\mathbf{Q}(\mathbf{z}(s)) \cdot \mathbf{K}^\top(\mathbf{z}(s))}{\sqrt{d_k}} \mathrm{d} s, \\
        \mathbf{L}(t + \Delta t) &= \mathbf{L}(t) + \Delta t \cdot \frac{\mathrm{d}}{\mathrm{d}s} \frac{\mathbf{Q}(\mathbf{z}(s)) \cdot \mathbf{K}^\top(\mathbf{z}(s))}{\sqrt{d_k}} \bigg|_{s = t}
    \end{alignat}
\end{small}
\end{subequations}
where $\mathbf{Q}(\cdot)$ and $\mathbf{K}(\cdot)$ correspond to the query and key in the attention mechanism, respectively. $\mathbf{L}(0) = \mathbf{Q}(\mathbf{z}(0)) \cdot \mathbf{K}^\top(\mathbf{z}(0)) / \sqrt{d_k}$. We show that updating the dynamic of $\mathbf{L}$ following Eq.~(\ref{eq: transform}) is equivalent to directly updating $\mathbf{A}$ in Appendix \ref{app:equi}.

\textbf{(2)} Second, both the solution to projection Eq.~(\ref{eq:pj_sparse}) and its gradient w.r.t. $\mathbf{A}$ are in closed form. See Proposition \ref{prop:sol_sparse} and Proposition \ref{prop:gradient_sparse} in Appendix \ref{app:soft_sparse} for more details. This, in turn, eases the computational flow in the neural architecture with high efficiency and stability.

Though only two versions of projections are discussed under Bregman divergence, we believe they are sufficiently distinguishable for analyzing the behavior of ProjDEs. For generic neural-friendly projections, we leave them to our future work.

\begin{figure}[tb]
  \centering
    \includegraphics[width=.30\textwidth]{figs/soft_vs_sparse.pdf}
  \caption{Comparsion of softmax and L2 projection onto a simplex. We see that the softmax projection trends to project onto the ``center'' of the simplex while the L2 projection trends to project onto the corner.}
  \label{fig:softVSsparse}
\end{figure}

\textbf{Integration.} We employ the standard Euler's discretization for performing the forward integration by updating $\mathbf{z}$ and $\mathbf{A}$ simultaneously with a sufficiently small time step. We term our approach a \textbf{d}ecoupling-based \textbf{n}eural \textbf{s}ystem (\textbf{DNS}) using projection Eq.~(\ref{eq:pj_soft}) and \textbf{DNS\textsubscript{G}} using projection Eq.~(\ref{eq:pj_sparse}), respectively.

\begin{figure*}[tb]
\captionsetup[subfigure]{labelformat=empty}
    \centering
    \foreach \t in {5,...,12}{
        \begin{subfigure}{.11\textwidth}
            \centering
            \includegraphics[width=.98\textwidth, trim=100 100 100 100,clip]{figs/softmax_att/trajectory_\t.pdf}\\
            \includegraphics[width=.98\textwidth, trim=100 100 100 100,clip]{figs/softmax_att/attention_\t.pdf}
            \caption{$t=$\t}
        \end{subfigure}
    }
  \vspace{-8pt}
  \caption{A figure showing the corresponding three-body trajectory (on the top), as well as the evolution over time on interactions (at the bottom) between three \textbf{latent sub-systems} in a Three-Body environment. Timestamp from 5 to 12.}
  \label{fig: att}
\end{figure*}
\section{Experiments}
We evaluate the performance of DNS on multiple synthetic and real-world datasets. More details about the dataset and implementation details can be found in Appendix \ref{app: experiment details} and \ref{app: dataset}. Throughout all the tables consisting of the results, ``-'' indicates ``not applicable'' since RIM cannot handle irregular cases.

\begin{remark}
    In all the experiments, the input feature is treated holistically without any distinguishable parts. For example, in the Three Body dataset, the input is a 9-dimensional vector, with every 3 dimensions (coordinates) from a single object. However, this prior is not fed into any models in comparison. Thus, we do not compare to models integrated with strong prior such as \cite{kipf2018neural}.
\end{remark}

\textbf{Baselines.} We compare DNS with several selected models capturing interactions or modeling irregular time series, including \textbf{CT-GRU} \citep{mozer2017discrete} using state-decay decay mechanisms, \textbf{RIM} \citep{goyal2019recurrent} updating almost independent modules discretely, and \textbf{NeuralCDE} \citep{kidger2020neural} which reports state-of-the-art performance on several benchmarks.
\begin{table}[t]
\begin{center}
\begin{small}
\begin{sc}
\caption{\textbf{Trajectory prediction}. MSE loss of the three body dataset ($\times 10^{-2}$).}
\label{tab: three_body}
\begin{tabular}{lccr}
\toprule
    Model  & Regular & Irregular \\
\midrule
CT-GRU    & 1.8272   & 2.4811 \\
NeuralCDE & 3.3297   & 5.0077 \\
RIM       & 2.4510   & -      \\
\midrule
DNS       & \textbf{1.7573}  & \textbf{2.2164} \\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\textbf{Adapting DNS to the noisy case.}
To allow DNS fitting to noisy and uncertain circumstances, we create a variant by slightly modifying it.
This variant is obtained by replacing cubic spline interpolation over $\mathbf{x}_i(t)$ with natural smoothing spline \citep{green1993nonparametric}, in consideration of incorporating smoother controls and alleviating data noise. This version is termed as \textbf{DNS\textsubscript{S}}.

\subsection{Three body}
The three-body problem is characterized by a chaotic dynamical system for most randomly initial conditions. A small perturbation may cause drastic changes in the movement. Taking into account the problem's complexity, it is particularly suitable for testing our approach.
In this experiment, we consider a trajectory predicting problem given the noisy historical motion of three masses, where gravity causes interactions between them.
Therefore, models need to (implicitly) learn both Newton's laws of motion for modeling sub-system dynamics and Newton's law of universal gravitation to decouple the latent interaction. This dataset consists of 50k training samples and 5k test samples. For each sample,
8 historical locations for the regular setting and 6 historical locations (randomly sampled from 8) for the irregular setting in the 3-dimensional space of three bodies are given to predict 3 subsequent locations. To equip with the cluttered setting, the correspondence between dimensions and bodies will not be fed into the learning models, hence a 9-dimensional observation at each time stamp. Models' performance is summarized in Table~\ref{tab: three_body}.
\begin{figure}[t]
    \centering
    \includegraphics[width=.4\textwidth, trim = {100 40 100 0, clip}]{figs/three-body-focus.pdf}
    \caption{A figure showing the focus of 3 sub-systems on 9-dimensional input of Three Body. The strength of focus is reflected by the thickness of the lines.}
    \label{fig:body_focus}
\end{figure}
We can conclude that DNS outperformed all the selected counterparts in both regular and irregular settings. Notably, although our method is built on NeuralCDE, with the decoupling, the performance can be significantly improved. See Table~\ref{tab: appendix three body} in Appendix \ref{app:three body} for more detailed results.
\begin{table*}[t]
\begin{center}
\begin{small}
\begin{sc}
\caption{\textbf{Link prediction}. Accuracy on Spring ($\%$). \textsc{Clean}, \textsc{Noisy}, and \textsc{Short} correspond to settings with clean, noisy, and short portion data, respectively. Detailed results for \textsc{Clean} and \textsc{Noisy} are separately summarized in Tab.~\ref{tab:spring app} and Tab.~\ref{tab: appendix noisy spring} in the appendix.}
\label{tab: spring}
\begin{tabular}{lcc|cc|cc}
\toprule
    \multirow{2}{*}{Model} & \multicolumn{2}{c}{Clean} & \multicolumn{2}{c}{Noisy} & \multicolumn{2}{c}{Short}\\
    \cline{2-7}
    & Regular & Irregular & Train\&Test & Test & 50\% & 25\% \\
    \midrule
    CT-GRU & 92.89$\pm$0.52   &  88.47$\pm$0.34 & 92.71$\pm$0.55 & 92.80$\pm$0.53 & -& -\\
    NeuralCDE & 92.47$\pm$0.06   &  89.74$\pm$0.18 & 90.76$\pm$0.08 & 89.61$\pm$0.09 & -& -\\
    RIM & 89.73$\pm$0.07   &  -  &89.65$\pm$0.14 & 89.64$\pm$0.10  & -& - \\
    \midrule
    DNS\textsubscript{G} & 94.31$\pm$0.48 & \textbf{94.25$\pm$0.29} & - & - & 92.58 & 92.31 \\
    DNS\textsubscript{S} & - & - & \textbf{93.67}$\pm$\textbf{0.57} & \textbf{92.99}$\pm$\textbf{1.30} & -& -\\
    DNS & \textbf{94.44$\pm$0.69}  & 93.60$\pm$1.21 & 93.42$\pm$1.05 & 89.56$\pm$0.42 & 91.11 & 92.13\\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table*}

\textbf{Visualization and analysis.} We visualize dynamics $\mathbf{A}$ of DNS along the movements of three body system. See Fig. \ref{fig: att} for results. We set the time stamps starting from 5 to 12 to make visualization more informative. It is seen in the beginning ($t=5,6$ or even earlier), $\mathbf{A}$ remains stable as the three bodies are apart from each other without intensive interactions. At $t=7$, $\mathbf{A}$ demonstrates obvious change when two bodies start to the coil. Another body joins in this party at $t=8$, yielding another moderate change of $\mathbf{A}$. When flying apart, one body seems more independent, while another two keep entangled together. These are well reflected via the meta-system $\mathbf{A}$. To further see how the holistic 9-dimensional input is decoupled into sub-systems $\mathbf{z}_i$, we visualize the sub-system focus in Fig.~\ref{fig:body_focus} (also see Appendix \ref{app:cam}). Interestingly, latent entities (sub-systems) do not correspond to physical entities (three bodies). Instead, the first sub-system puts more focus on the whole input, but the remaining two sub-systems concentrate on the x-axis and y/z-axis, respectively. Though counterintuitive, this unexpected decoupling exhibits good performance. We will investigate how to decouple out physical entities from cluttered observations in our future work.
\begin{table}[tb]
\centering
\begin{small}
\begin{sc}
\caption{\textbf{Link prediction}. Ablation study. ($\%$).
}
\label{tab:abl}
\begin{tabular}{lc}
\toprule
Control & Accuracy ($\%$) \\
\midrule
No encoding       & 91.57                \\
MLP(2$\times$input)     & 91.51                \\
MLP(16$\times$input)      & 91.17                \\
DNS (8$\times$MLP(2$\times$input))      & \textbf{95.38}                \\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{table}

\begin{figure*}[tb]
\captionsetup[subfigure]{labelformat=empty}
\centering
    \foreach \t in {1,2,3,4,8,13,19,28}{
        \begin{subfigure}{.11\textwidth}
        \centering
        \includegraphics[width=.98\textwidth, trim={25 25 25 0}]{figs/spring_att/traj_\t.pdf}\\ \vspace{-5pt}
         \includegraphics[width=.98\textwidth, trim={25 25 25 0}]{figs/spring_att/soft_\t.pdf}\\ \vspace{-5pt}
        \includegraphics[width=.98\textwidth, trim={25 25 25 0}]{figs/spring_att/sparse_\t.pdf}
        \vspace{-15pt}
        \caption{$t=$\t}
        \end{subfigure}
    }
    \vspace{-10pt}
\caption{Visualization of the evolution of the meta-systems of DNS and DNS\textsubscript{G} on Spring dataset. On each time stamp $t$, from top to bottom, we show the trajectory of the 5 balls, the meta-system state of DNS, and the meta-system state of DNS\textsubscript{G}, respectively.}
\label{fig:vis_spring}
\end{figure*}

\subsection{Spring}
We experiment with the capability of DNS in decoupling the independence in complex dynamics controlled by simple physics rules. We use a simulated system in which particles are connected by (invisible) springs \citep{kuramoto1975self,kipf2018neural}. Each pair of particles has an equal probability of having an interaction or not. Our task is to use observed trajectory to predict whether there are springs between any pair of two particles, which is analogous to the task of link prediction under a dynamical setting. This can be inferred from whether two trajectories change coherently. The spring dataset consists of 50k training examples and 10k test examples. Each sample has a length of 49. We test a variety of combinations of the number of sub-systems and dimensions of the hidden state. Experimental results are in Table~\ref{tab: spring}. To test the models' noise resistance, we add Gaussian noise to the spring dataset and obtain the noisy spring dataset. We set two scenarios, ``Train\&Test'' and ``Test'', corresponding to injecting noise at both training and test phases and only at testing phases, respectively. Experimental results are in Table~\ref{tab: spring}.

\textbf{Clean Spring.} From \textsc{Clean} part of Table \ref{tab: spring}, we see variants of DNS stably outperform all the selected counterparts by a large margin. Especially, under the irregularly sampled data, DNS and DNS\textsubscript{G} have a remarkable performance gap with all other methods and maintain reliability as in the regular setting. We believe this is significant since learning from irregularly sampled data is typically much more difficult than learning from normal data.

\textbf{Noisy Spring.} According to \textsc{Noisy} part of Table \ref{tab: spring}, DNS\textsubscript{S} is quite reliable in noisy cases. It seems a smoothing procedure on the controls can be helpful under massive uncertainty. Also, we see that adding noise tends to damage the performance of all methods. This also raises one of our future research directions to investigate how to handle different controls. Without applying a smooth cubic spline, DNS can still have a good performance, which indicates that by decoupling, the model focuses on learning latent interaction patterns, and patterns are less susceptible to noise.

\textbf{Visualization and analysis.}  We also visualize state $\mathbf{A}$ of meta-systems over time in Fig.~\ref{fig:vis_spring} for Spring. From top to bottom, the first, second and third rows correspond to the trajectory of particles, meta-system state of DNS, and meta-system state of DNS\textsubscript{G}. One interesting thing we note is that the interactions in DNS\textsubscript{G} almost concentrate on the starting portion of all the time stamps. At $t=8$ and after, there is no interaction at all. Though not obvious, this also happens to DNS in the sense that $\mathbf{A}$ tends to be diagonal. We suppose this is because DNS and DNS\textsubscript{G} only need a portion of data from the start to determine the existence of a link rather than looking into all the redundant time stamps.

\textbf{Short Spring.} We thus verify this by training and testing both variants with 50\% and 25\% of data cropped from the starting time stamp and summarize results in \textsc{Short} part of Table \ref{tab: spring}. It is seen that incomplete data in this task only slightly impact the performance. And this can be surprisingly reflected in the evolution of meta-systems. This also aligns with the intuition that \emph{Link prediction} needs fewer data than \emph{Trajectory prediction} as in Three Body.
\begin{table}[tb]
\begin{center}
\begin{small}
\begin{sc}
\caption{\textbf{Video classification}. Accuracy of the human actions dataset ($\%$). \textsc{Norm} and \textsc{Unnorm} refer to normalized and unnormalized inputs, respectively. Detailed results with superscript $^\dagger$ and $^\ddagger$ are in Tab.~\ref{tab:action app} and Tab.~\ref{tab:v_unnorm}, respectively.}
\label{tab: video}
\begin{tabular}{lccc}
\toprule
    \multirow{2}{*}{Model} &  \multicolumn{1}{c}{Norm} & \multicolumn{2}{c}{Unnorm}\\
    \cline{2-4}
     & Irreg & Reg & Irreg \\
\midrule
CT-GRU    & 67.30$\pm$6.19$^\dagger$   & 60.33$^\ddagger$ & 66.67$^\ddagger$\\
NeuralCDE & 89.73$\pm$3.38$^\dagger$   & 70.33$^\ddagger$ & 59.17$^\ddagger$\\
RIM       & -      & 55.50$^\ddagger$ & - \\
\midrule
DNS       & \textbf{91.35}$\pm$\textbf{3.48}$^\dagger$  & \textbf{97.00}$^\ddagger$ & \textbf{95.33}$^\ddagger$\\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\textbf{Ablation study.} Since our method merely incorporates an extra meta-system and a control encoder for modeling the interaction compared to standard NeuralCDE, we conduct experiments under different settings to see how different encoders and hidden state dimensions can contribute to improving NeuralCDE. To ensure fairness, we cast a 2-layer MLP with different output sizes (2 and 16 times of input size) as in DNS to obtain varying sizes of controls. Results are summarized in Table \ref{tab:abl} (detailed in Tab.~\ref{tab: appendix spring abl}). We see that with an extra control encoder, there is no obvious performance difference among these settings. However, once the interaction meta-system is imposed, DNS can achieve quite significant performance gain. This, in turn, shows the necessity of the proposed meta-system for explicitly modeling the evolving interactions.

\subsection{Human actions}
The recognition of human actions dataset contains three types of human actions, which are hand clapping, hand waving, and jogging \citep{schuldt2004recognizing}. For this dataset, we consider the limbs of the character as subsystems. When the character does one kind of action, subsystems interact in a specific pattern. We test the performance of all the selected models with the learnable backbone Resnet18 \citep{he2016deep}. We also test the compatibility of all methods with different dynamical ranges: \textsc{Norm} and \textsc{Unnorm} indicate pixel value in $[0,1]$ and $[0,255]$, respectively. Experimental results are summarized in Table \ref{tab: video}. DNS consistently outperforms all other methods, and exhibits strong compatibility to drastically changed ranges under \textsc{Unnorm} setting. Thus it is potentially more flexible to be integrated into various tasks with a large dynamical range (e.g., earthquake).

To view how the decoupling works for video recognition tasks, we visualize the strength of the learned parameters by mapping the 128-D feature into 6 latent sub-systems in Figure \ref{fig: cnn top feature} with re-ordered indices for better view.
It can be seen that there are some obvious latent structures in the grouping of the parameters 128-D control to the system. Each sub-system mainly focuses on a small portion of the control, based on which we can infer that each sub-system models different components in inputted images.

\begin{figure}[tb]
  \centering
    \includegraphics[width=.45\textwidth, trim=0 20 0 20, clip]{figs/top40feature.pdf}
    \vspace{-10pt}
  \caption{A figure showing the importance of each feature vector entry for subsystems}
  \label{fig: cnn top feature}
\end{figure}

\section{Conclusion}
In this paper, we propose a method for modeling cluttered and irregularly sampled sequential data. Our method is built upon the assumption that complex observation may be derived from relatively simple and independent latent sub-systems, wherein the interactions also evolve over time. We devise a strategy to explicitly decouple such latent sub-systems and a meta-system governing the interaction. Inspired by recent findings of projected differential equations and the tool of Bregman divergence, we present a novel interpretation of our model and pose some potential future directions. Experiments on various tasks demonstrate the prominent performance of our method over previous state-of-the-art methods.


\section{Introduction}

Electron transport in semiconductors is of fundamental interest and of high relevance for microelectronic devices \cite{sze2007, Lundstrom_2000, Hartnagel_2001}. The upper limit for the mobility of a semiconductor is governed by scattering of electrons by phonons. Early studies of charge transport properties employed a semi-empirical description of the band structure and electron-phonon scattering. The introduction of the Monte Carlo method allowed for the numerical simulation of transport with fewer approximations \cite{reggiani_1983}. Later, full-band MC tools capable of simulating realistic device geometries were developed, but the treatment of the e-ph scattering rates in general required some fitting parameters.  \cite{ Hess:1991,Fischetti:1991, fischetti_1991_2}.  The development of the ab-initio description of the electron-phonon interactions based on density functional theory (DFT), density functional perturbation theory (DFPT) and Wannier interpolation has enabled the parameter-free computation of low-field charge transport properties such as mobility  \cite{Bernardi_2016, Giustino_2017}. These methods have now been applied to a range of semiconductors, including Si \cite{restrepo2009, Ponce_2018}, GaN \cite{ponce2019gan}, GaAs \cite{Zhou_2016}, two-dimensional materials \cite{borysenko2010, kaasbjerg2012, 2dmobi2020}, and others. Recent methodology developments, including two-phonon scattering \cite{lee2020}, quadrupole interactions \cite{brunin2020,quadprb2020}, and $GW$ corrections \cite{Ponce_2018,gwpt_2019}, have facilitated a rigorous comparison of the accepted level of theory with experiment.

The accuracy of the ab-initio calculations has been mainly based on comparison to low-field mobility values. Recent works have extended these calculations to beyond low-field transport and noise properties. \cite{hatanpaa2022two,warmelectrons,cheng2022high,catherall2023high}
In GaAs, it was found that although the qualitative shape of the drift velocity versus electric field curve was predicted correctly compared to experiment, the magnitude of the drift velocity was overpredicted by about 50\%.\cite{cheng2022high} The inclusion of the next-leading-order term of scattering involving two phonons (2ph)  yielded a low-field mobility and drift velocity with substantially improved agreement. However, the characteristic nonmonotonic trend of PSD with electric field was not predicted even with the 2ph theory. Owing to the high cost of the ab-initio calculations, the treatment of 2ph processes in that work required several approximations, such as the neglect of off-shell 2ph scattering processes.  Whether these neglected processes or other numerical considerations can account for the PSD discrepancy remains unknown.

Here, we introduce a semi-analytical model for both 1ph and 2ph e-ph scattering via the Fröhlich mechanism, allowing the full 2ph scattering term to be treated over the wide range of energies needed for high-field transport while introducing error on the order of only a few percent. We find that the transport and noise properties are qualitatively unchanged compared to the ab-initio calculations. The calculated scattering rates agree with those obtained from photoluminescence experiments to  within 25\% for the $\Gamma$ valley rates, and a factor of two for the $\Gamma$-L intervalley rates. Despite this degree of agreement, the qualitative discrepancy observed previously for the PSD remains. We consider the remaining approximations in the semi-analytical model and find that they are unlikely to account for the PSD discrepancy. Therefore, we conclude that the characteristic peak in the PSD with electric field arises from the formation of space charge domains rather than intervalley scattering as has been assumed in the literature. This finding has implications for the use of transport measurements to study intervalley scattering.

\section{Theory}

\subsection{Overview of formalism for charge transport and noise properties}

We first review the ab-initio treatment of electron transport and electronic noise using the BTE \cite{hamaguchi2010basic}.
For a spatially homogeneous system with electric field $\vec{\mathcal{E}}$, the electron distribution function $f_{m\kwv}$
is governed by

\begin{equation}\label{eq:BTE}
    \frac{\partial f_{m\boldsymbol{k}}}{\partial t} + \frac{e\vec{\mathcal{E}}}{\hbar}\cdot \nabla_{\boldsymbol{k}}f_{m\kwv} = \mathcal{I}[f_{m\kwv}],
\end{equation}

where $f_{m\boldsymbol{k}}$ is the electron occupation in the state with band index $m$ and wave vector $\boldsymbol{k}$, $e$ is the fundamental charge, and $\mathcal{I}[f_{m\kwv}]$ is the collision term, which describes the scattering of electrons with phonons \cite{Bernardi_2016}.
In this work, index $m$ is dropped in all the following derivations for simplicity as only one band is relevant for electron transport in GaAs in the range of electric fields considered.

For non-degenerate electrons, the collision term can be linearized as \cite{cheng2022high}:

\begin{equation} \label{eq:linear_BTE}
\mathcal{I}[f_{\boldsymbol{k}^\prime}] \approx  - \sum_{\boldsymbol{k}}  \Theta_{\boldsymbol{k}^\prime, \boldsymbol{k}} f_{\boldsymbol{k}},
\end{equation}
where $\Theta_{\boldsymbol{k}^\prime, \boldsymbol{k}}$ is e-ph collision matrix. The diagonal elements $\Theta_{\boldsymbol{k}, \boldsymbol{k}}$ are equal to the total scattering rates as $\Theta_{\boldsymbol{k}, \boldsymbol{k}} = \Gamma_{\boldsymbol{k}} = - \sum_{\boldsymbol{k}^\prime\neq \boldsymbol{k}} \Theta_{\boldsymbol{k}^\prime, \boldsymbol{k}}$. With this linearization and a finite difference representation of the  derivative operator $\nabla_{\boldsymbol{k}}$ \cite{warmelectrons}, Eq.~\ref{eq:BTE} becomes a linear partial-differential equation which can be solved by numerical linear algebra. The equation for the steady distribution function  $f_{\boldsymbol{k}}^s$ is given by Eq.~6 of Ref. \onlinecite{warmelectrons}. Steady-state mean transport properties such as drift velocity can be calculated  with the appropriate Brillouin zone sum using this distribution.

The current power spectral density (PSD) is used to characterize fluctuations in occupation about the mean distribution. The PSD is defined as the Fourier transform of the autocorrelation of the current density fluctuations (Eq. 12 of Ref.~\onlinecite{warmelectrons}).
Following Ref.~\onlinecite{warmelectrons}, the current PSD at frequency $\omega$ can be calculated as
\begin{equation}
\label{noisesum}
    S_{j_{\alpha}j_{\beta}}(\omega) =
    2 \bigg(\frac{2 e}{\mathcal{V}_0}\bigg)^2 \Re \left[ \sum_{\mathbf{k}} v_{\boldsymbol{k}\alpha} g_{\boldsymbol{k}\beta} \right],
\end{equation}
where $v_{\boldsymbol{k}\alpha}$ is the group velocity of the electron at state $\boldsymbol{k}$ along axis $\alpha$, $\mathcal{V}_0$ is the supercell volume, and $g_{\boldsymbol{k}\beta}$ is the effective distribution function \cite{rustagi2014} which satisfies the following equation:

\begin{equation} \label{eq:linear_g}
    \sum_{\boldsymbol{k}} (\Theta_{\boldsymbol{k}^\prime, \boldsymbol{k}} + \frac{e \vec{\mathcal{E}}}{\hbar} \cdot \vec{D}_{\boldsymbol{k}^\prime, \boldsymbol{k}} + i\omega \delta_{\boldsymbol{k}^\prime, \boldsymbol{k}}) g_{\boldsymbol{k} \alpha} = f_{\boldsymbol{k}^\prime}^s (v_{\boldsymbol{k}^\prime,\alpha} - V_{\alpha}),
\end{equation}
where $D$ is the finite difference representation of $\nabla_{\boldsymbol{k}}$ specified in Eq.~\ref{eq:BTE}, $f_{\boldsymbol{k}}^s$ is the steady-state occupation for the state at wave vector $\boldsymbol{k}$, and $V_\alpha$ is the drift velocity along axis $\alpha$.

The e-ph collision matrix is obtained from perturbation theory in orders of the e-ph interaction strength. For 1ph scattering of non-degenerate electrons,  the non-diagonal scattering matrix elements are given by \cite{cheng2022high}:
\begin{equation} \label{eq:Theta_1ph}
    \Theta_{\boldsymbol{k}^\prime = \boldsymbol{k}+\boldsymbol{q}, \boldsymbol{k}}^{\mathrm{(1ph)}} = -\frac{2 \pi}{\hbar} \frac{1}{N} \sum_{\nu}  \left|g_{\nu}(\boldsymbol{k},\boldsymbol{q})\right|^2
     \bigg(\delta(\epsilon_{\boldsymbol{k}}-\hbar \omega_{\nu \boldsymbol{q}}-\epsilon_{\boldsymbol{k}+\boldsymbol{q}})N_{\boldsymbol{q}} + \delta(\epsilon_{\boldsymbol{k}} + \hbar \omega_{\nu\boldsymbol{q}} -\epsilon_{\boldsymbol{k}+\boldsymbol{q}}) (N_{\boldsymbol{q}}+1) \bigg)
\end{equation}
for $\boldsymbol{k} \neq \boldsymbol{k}^\prime$, where $g_{\nu}(\boldsymbol{k},\boldsymbol{q})$ is the e-ph scattering matrix element, $\epsilon_{\boldsymbol{k}}$ is the energy of the electronic state $\boldsymbol{k}$, $\omega_{\nu \boldsymbol{q}}$ is the frequency of phonon with mode $\nu$ and wave vector $\boldsymbol{q}$, and $N$ is the total number of $\boldsymbol{k}$ points in the Brillouin zone, $N_{\boldsymbol{q}}$ is the phonon occupation according to the Bose-Einstein statistics. The two delta functions are energy conservation conditions for the emission and absorption subprocesses, respectively.

\subsection{Two-phonon scattering}

The collision integral for 2ph scattering was derived in Ref.~\onlinecite{lee2020}, and the linearized form is given in Ref.~\onlinecite{cheng2022high}. 2ph scattering can be divided into the two-emission (2e), one-emission-one-absorption (1e1a) and two-absorption (2a) subprocesses. Here, we rewrite the original formalism in Ref.~\onlinecite{lee2020} to facilitate the derivations in the next section.
Assuming non-degenerate electron statistics in Eq.~4 in Ref.~\onlinecite{lee2020}, we rewrite Eq.~12 in Ref.~\onlinecite{cheng2022high} by splitting $W^{(i)}$ to two parts:
\begin{equation}
\begin{aligned} \label{eq:old_W}
    W^{(i)} &= |\tilde{W}_{\boldsymbol{k}, \boldsymbol{q}, \boldsymbol{p}, \alpha_2} + \tilde{W}_{\boldsymbol{k}, \boldsymbol{p}, \boldsymbol{q}, \alpha_1}|^2
    \\
    &= \left[|\tilde{W}_{\boldsymbol{k}, \boldsymbol{q}, \boldsymbol{p}, \alpha_2}|^2 + \text{Re}(\tilde{W}_{\boldsymbol{k}, \boldsymbol{q}, \boldsymbol{p}, \alpha_2} \tilde{W}_{\boldsymbol{k}, \boldsymbol{p}, \boldsymbol{q}, \alpha_1}^*)\right] + \left[ |\tilde{W}_{\boldsymbol{k}, \boldsymbol{p}, \boldsymbol{q}, \alpha_1}|^2 + \text{Re}(\tilde{W}_{\boldsymbol{k}, \boldsymbol{p}, \boldsymbol{q}, \alpha_1}\tilde{W}_{\boldsymbol{k}, \boldsymbol{q}, \boldsymbol{p}, \alpha_2}^* ) \right],
\end{aligned}
\end{equation}
where
\begin{equation} \label{eq:W_tilde}
    \tilde{W}_{\boldsymbol{k}, \boldsymbol{q}, \boldsymbol{p}, \alpha} = \frac{ g_{\nu}(\boldsymbol{k},\boldsymbol{q})g_{\mu}(\boldsymbol{k}+\boldsymbol{q},\boldsymbol{p}) }
    { \epsilon_{\boldsymbol{k}'}-\epsilon_{\boldsymbol{k}+\boldsymbol{q}}+\alpha \hbar \omega_{\nu\boldsymbol{p}} + i\eta - i\hbar \Gamma_{\boldsymbol{k}+\boldsymbol{q}} / 2 },
\end{equation}
where $\eta$ is a positive infinitesimal and $\Gamma_{\boldsymbol{k} + \boldsymbol{q}} = \Gamma^{(\mathrm{1ph})}_{\boldsymbol{k} + \boldsymbol{q}} + \Gamma^{(\mathrm{2ph})}_{\boldsymbol{k} + \boldsymbol{q}}$ is the total scattering rate of the intermediate state $\boldsymbol{k} + \boldsymbol{q}$.
For the 2e and 2a subprocesses, the two terms give the same contribution after the summation.
For the 1e1a subprocess, the two terms in Eq.~\ref{eq:old_W} physically represent the emission-then-absorption (a-e) and absorption-then-emission (e-a) subprocesses, respectively.
Finally, we exchange the summation order of $\boldsymbol{q}$ and $\boldsymbol{p}$ in Eq.~9 in Ref.~\onlinecite{cheng2022high} for the second term of  Eq.~\ref{eq:old_W} and arrange the equations to obtain:

\begin{equation} \label{eq:Theta_2ph}
    \Theta^{(\mathrm{2ph})}_{\boldsymbol{k}^\prime, \boldsymbol{k}} = - \frac{2\pi}{\hbar}\frac{1}{N^2} \sum_{\alpha_1 = \pm 1} \sum_{\alpha_2 = \pm 1} \sum_{\boldsymbol{q} + \boldsymbol{p} = \boldsymbol{k}' -\boldsymbol{k}}\sum_{\nu\mu} \tilde{\Theta}^{(\alpha_1, \alpha_2)}_{\boldsymbol{k},\boldsymbol{q}\nu,\boldsymbol{p}\mu}
\end{equation}

for $\boldsymbol{k} \neq \boldsymbol{k}^\prime$, where $\alpha_1$ and $\alpha_2$ indicate whether the first and second phonon is emitted ($\alpha_{1,2}=1$) or absorbed ($\alpha_{1,2}=-1$), so that the four combinations of $\alpha_{1,2}=\pm 1$ describe the four subprocesses.

The term $\tilde{\Theta}^{(\alpha_1, \alpha_2)}_{\boldsymbol{k},\boldsymbol{q}\nu,\boldsymbol{p}\mu}$ in Eq.~\ref{eq:Theta_2ph} is defined as

\begin{equation} \label{eq:Theta_2ph_subprocess}
    \tilde{\Theta}^{(\alpha_1, \alpha_2)}_{\boldsymbol{k},\boldsymbol{q}\nu,\boldsymbol{p}\mu} = (N_{\boldsymbol{q}} + \delta_{\alpha_1, 1}) (N_{\boldsymbol{p}} + \delta_{\alpha_2, 1})  \left[|\tilde{W}_{\boldsymbol{k}, \boldsymbol{q}, \boldsymbol{p}, \alpha_2}|^2 + \text{Re}(\tilde{W}_{\boldsymbol{k}, \boldsymbol{q}, \boldsymbol{p}, \alpha_2} \tilde{W}_{\boldsymbol{k}, \boldsymbol{p}, \boldsymbol{q}, \alpha_1}^*)\right] \delta
    (\epsilon_{\boldsymbol{k}} - \epsilon_{\boldsymbol{k}'} - \alpha_{1}\omega_{\nu\boldsymbol{q}} - \alpha_{2}\omega_{\mu\boldsymbol{p}}).
\end{equation}

As the 2ph scattering rates depend on the intermediate state rates, the 2ph scattering must be calculated iteratively.
In each iteration, the intermediate state rates of the previous iteration is used to update the 2ph scattering matrix $\Theta^{(\mathrm{2ph})}_{\boldsymbol{k}^\prime, \boldsymbol{k}}$ and scattering rates $\Gamma^{(\mathrm{2ph})}_{\boldsymbol{k}}$.

Compared with 1ph scattering, 2ph scattering is much more computationally expensive, particularly for high-field transport which requires a larger energy window than for low-field transport. In Ref.~\onlinecite{cheng2022high}, several approximations were made to make the computation feasible, including limiting the number of self-consistent iterations of the 2ph rates to three, restricting the off-shell extent $|\epsilon_{\boldsymbol{k}'}-\epsilon_{\boldsymbol{k}+\boldsymbol{q}}+\alpha \hbar \omega_{\nu\boldsymbol{p}}|$ to 25 meV, and neglecting the interference term $\text{Re}(\tilde{W}_{\boldsymbol{k}, \boldsymbol{q}, \boldsymbol{p}, \alpha_2} \tilde{W}_{\boldsymbol{k}, \boldsymbol{p}, \boldsymbol{q}, \alpha_1}^*)$ term in Eq.~\ref{eq:Theta_2ph_subprocess}.  Additionally, the maximum grid density that could be used was $200 \times 200 \times 200$. The effect of these approximations on the observable transport and noise properties was not assessed. In particular, the on-shell approximation neglects off-shell processes and thus underestimates the scattering rates. These approximations were mentioned as possible reasons for the PSD discrepancy in Ref. \cite{cheng2022high}.

\subsection{Semi-analytical model for 1ph and 2ph $\mathbf{\Gamma}$-$\mathbf{\Gamma}$ scattering} \label{sec:semi-analytical_model}

In this section, we introduce a semi-analytical model to treat 1ph and 2ph $\Gamma-\Gamma$ intravalley scattering by the Fröhlich interaction that retains the accuracy of the ab-initio formalism to within a few percent while allowing the approximations described above to be lifted.
This model is based on the fact that over the range of wavevectors considered in this study, the $\Gamma$ valley in GaAs is nearly spherically symmetric, and $\Gamma\text{-}\Gamma$ scattering can be accurately described by using only the Fröhlich interaction \cite{verdi2015frohlich}.
The model is valid only for $\Gamma$ intravalley scattering because $\Gamma$-L intervalley scattering lacks an analytic description of similar accuracy.

The semi-analytical model uses the following approximations. First, the band structure is described using the Kane model \cite{kane1957band} for a spherically symmetric, non-parabolic band. This description is accurately satisfied for the $\Gamma$ valley, with the Kane model bands deviating from the ab-initio band structure by at most 7\% over the range of wave vectors considered ($\sim 0.1 G$, where  $G$ is the reciprocal lattice constant). Second, prior works have shown that $\Gamma\text{-}\Gamma$ e-ph scattering in GaAs is dominated by longitudinal optical (LO) phonons via the Fröhlich interaction \cite{frohlich1954electrons}. We therefore neglect  scattering processes involving other phonon branches and scattering mechanisms. The computed matrix elements $g_\text{LO}(\boldsymbol{k}, \boldsymbol{q})$  for Fröhlich scattering are found to exhibit negligible anisotropy so that $g_\text{LO}(\boldsymbol{k}, \boldsymbol{q}) = g_\text{LO}(k)$, enabling an analytic form of $g_\text{LO}(k)$ to be fit to the ab initio values as described  in Sec.\ref{sec:details}. In the range of wave vectors considered, this approximation is satisfied to within 3\% \cite{verdi2015frohlich}. Third, we take the LO phonon frequency to be a constant $\omega_\text{LO} = 35$ meV. In the range of phonon wave vectors $k \in (0, 0.2G)$ considered here, this assumption is satisfied to within less than 0.3\%.

We now discuss the treatment of 1ph and 2ph e-ph scattering based on these simplifications. The summation in Eq.~\ref{eq:Theta_2ph} may be rewritten as an integral in the Brillouin zone over the intermediate wave vector $k_m$ by letting $\boldsymbol{q} \rightarrow k_m - \boldsymbol{k}$ and $\boldsymbol{p} \rightarrow \boldsymbol{k}^\prime - k_m$.
Additionally, we exploit spherical symmetry to rewrite all the quantities in spherical coordinates as: $\Theta_{\boldsymbol{k}^\prime, \boldsymbol{k}} = \Theta(k, k^\prime, \theta_{\kwv, \kwv^\prime})$ and $\Gamma_{\boldsymbol{k}} = \Gamma_{k}$. After some simplifications, we obtain the 1ph and 2ph collision matrices as:

\begin{equation} \label{eq:Theta_1ph_integral}
    \Theta^{\mathrm{(1ph)}}(k, k^\prime, \theta_{\kwv, \kwv^\prime}) = \frac{2 \pi}{\hbar} \frac{1}{\Omega_\text{BZ}}  |g_{\text{LO}}(|\boldsymbol{k}^\prime-\boldsymbol{k}|)|^2
    \sum_{\alpha=\pm1} A_{\alpha}
     \delta(\epsilon_{k}-\alpha \hbar \omega_\text{LO}-\epsilon_{k^\prime})
\end{equation}
and
\begin{equation} \label{eq:Theta_2ph_integral}
    \Theta^{(\mathrm{2ph})}(k, k^\prime, \theta_{\kwv, \kwv^\prime}) = \frac{2\pi}{\hbar}\frac{1}{\Omega_{\text{BZ}}^2} \sum_{\alpha_1 = \pm 1} \sum_{\alpha_2 = \pm 1}
    A_{\alpha_1}A_{\alpha_2} \delta  (\epsilon_k - \epsilon_{k^\prime} - (\alpha_{1} + \alpha_{2}) \hbar \omega_\text{LO}) I^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime}),
\end{equation}
where $\Omega_{\text{BZ}}$ is the Brillouin zone volume. $A_{\alpha}$ is the phonon occupation factor defined as

\begin{equation} \label{eq:phonon_occupation}
    A_{\alpha}=N_{\text{LO}} + \delta_{\alpha, +1},
\end{equation}
 where $N_{\text{LO}}=(\exp(\hbar \omega_{\text{LO}}/k_B T)-1)^{-1}$ is the LO phonon occupation, and $I^{(\alpha_1, \alpha_2)} = I_\text{1}^{(\alpha_1, \alpha_2)} + I_\text{2}^{(\alpha_1, \alpha_2)}$ is decomposed to the non-interference part $I_\text{1}^{(\alpha_1, \alpha_2)}$ and the interference part $I_\text{2}^{(\alpha_1, \alpha_2)}$:

\begin{equation} \label{eq:int_I}
\begin{aligned}
    I_\text{1}^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime}) &=
    \int |\tilde{W}_{\boldsymbol{k}, k_m - \boldsymbol{k}, \boldsymbol{k}^\prime-k_m, \alpha_2}|^2   \mathrm{d}^3 k_m
    \\
    &= \int \left|\frac{ g_\text{LO}(|k_m-\boldsymbol{k}|)g_\text{LO}(|\boldsymbol{k}^\prime-k_m|) }{ \epsilon_{k^\prime}-\epsilon_{k_m}+\alpha_2 \omega_\text{LO} + i\eta - \hbar \Gamma_{k_m} / 2 }\right|^2 \mathrm{d}^3 k_m,
\end{aligned}
\end{equation}
and
\begin{equation} \label{eq:int_I_itf}
\begin{aligned}
    I_\text{2}^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime}) &=
    \int_{\boldsymbol{q} + \boldsymbol{p} = \boldsymbol{k}^\prime - \boldsymbol{k}} \text{Re}(\tilde{W}_{\boldsymbol{k}, \boldsymbol{q}, \boldsymbol{p}, \alpha_2} \tilde{W}_{\boldsymbol{k}, \boldsymbol{p}, \boldsymbol{q}, \alpha_2}^{*})   \mathrm{d}^3 \boldsymbol{p} \mathrm{d}^3 \boldsymbol{q}
\end{aligned}
\end{equation}

Equation \ref{eq:int_I} can be further simplified by writing the integration in spherical coordinates and separating the radius and angular part:

\begin{equation} \label{eq:int_I_spherical}
    I_\text{1}^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime}) = \int \frac{ \tilde{I}^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime}, k_m) }{ |\epsilon_{k^\prime}-\epsilon_{k_m}+\alpha_2 \omega_\text{LO} + i\eta - \hbar \Gamma_{k_m} / 2|^2 } k_m^2 \mathrm{d} k_m ,
\end{equation}
where $\tilde{I}_\text{1}^{(\alpha_1, \alpha_2)}$ is the angular part defined as:
\begin{equation} \label{eq:I_tilde}
     \tilde{I}_\text{1}^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime}, k_m) = \int  |g_\text{LO}(|k_m-\boldsymbol{k}|)g_\text{LO}(|\boldsymbol{k}^\prime-k_m|)|^2 \sin \theta_{k_m} \mathrm{d} \theta_{k_m} \mathrm{d} \phi_{k_m},
\end{equation}
where $\theta_{k_m}$, $\phi_{k_m}$ are the polar angle and azimuthal angle defining the intermediate wave vector $k_m$, respectively.
Since $\tilde{I}_\text{1}^{(\alpha_1, \alpha_2)}$ is independent of the band structure and the self-energy, recomputation of this term in each 2ph iteration is not required.
In practice, to significantly reduce the computational cost, $\tilde{I}_\text{1}^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime}, k_m)$ are precomputed on a grid of $k$, $k^\prime$, $\theta_{\kwv, \kwv^\prime}$ and $k_m$ before the 2ph iteration. Here $k$, $k^\prime$, $\theta_{\kwv, \kwv^\prime}$ are not independent of  each other due to the energy conservation condition.
Once $\tilde{I}_\text{1}^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime}, k_m)$ is computed on a predefined grid, $\Theta^{(\mathrm{2ph})}(k, k^\prime, \theta_{\kwv, \kwv^\prime})$ can be calculated according to Eqs.~\ref{eq:Theta_2ph_integral}, \ref{eq:int_I_spherical} and \ref{eq:int_I_itf}.
We note that such separation of the radius and spherical part is not valid for $\tilde{I}_\text{2}^{(\alpha_1, \alpha_2)}$, so an expensive iterative update is required.
However, since $\tilde{I}_\text{2}^{(\alpha_1, \alpha_2)}$ is generally much smaller compared with $\tilde{I}_\text{2}^{(\alpha_1, \alpha_2)}$, we update $\tilde{I}_\text{2}^{(\alpha_1, \alpha_2)}$ every 10 iterations to decrease the cost of the self-consistent calculations.

To complete the 2ph iteration, the last quantities to be computed are the total 1ph and 2ph scattering rates:
\begin{equation} \label{eq:rates_int}
\Gamma_{k}^{(\text{type})} = \int \Theta^{(\text{type})}(k, k^\prime, \theta_{\kwv, \kwv^\prime}) \mathrm{d}^3 \boldsymbol{k}^\prime = \int \Theta^{(\text{type})}(k, k^\prime, \theta_{\kwv, \kwv^\prime}) 2\pi k^{\prime 2} \mathrm{d} k^\prime \sin \theta_{\kwv, \kwv^\prime} \mathrm{d} \theta_{\kwv, \kwv^\prime}
\end{equation}

where $\text{type}=\mathrm{1ph},\ \mathrm{2ph}$ indicates the type of scattering. We also perform the radius integration  over $k^\prime$ analytically to integrate the delta functions in Eq.~\ref{eq:Theta_1ph_integral} and ~\ref{eq:Theta_2ph_integral}.
The angular integrations in Eq.~\ref{eq:rates_int} are performed numerically; details are provided in  Sec.~\ref{sec:details}.

The  computational flow of the semi-analytical model is as follows.
First, we generate a grid of $k$, $\theta_{\kwv, \kwv^\prime}$, $k_m$ and  calculate the corresponding $k^\prime$ from the energy conservation conditions for each subprocess.
Second, we calculate $\Theta^{\mathrm{(1ph)}}(k, k^\prime, \theta_{\kwv, \kwv^\prime})$ by Eq.~\ref{eq:Theta_1ph_integral} and $\Gamma_{k}^{(\mathrm{1ph})}$ by Eq.~\ref{eq:rates_int}.
Then, we calculate $\tilde{I}_\text{1}^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime}, k_m)$ by Eq.~\ref{eq:I_tilde}.
Finally, we perform the self-consistent 2ph iterations through Eq.~\ref{eq:int_I_spherical}, \ref{eq:int_I_itf}, \ref{eq:Theta_2ph_integral} and \ref{eq:rates_int} until convergence, where Eq.~\ref{eq:int_I_itf} is calculated every 10 iterations.

From the perspective of computational cost, the semi-analytical model reduces the number of integration variables in the 2ph scattering rate calculation from 9 (integration over $\boldsymbol{k}$, $\boldsymbol{k}^\prime$, $k_m$) to 5 ($k$, $\theta_{\kwv, \kwv^\prime}$, $k_m$) due to the spherical symmetry, and avoids the recomputation of $I_\text{1}^{(\alpha_1, \alpha_2)}$ in the 2ph self-consistent iterations due to the separation of radius and angular integration in Eq.~\ref{eq:int_I_spherical}.  This reduction allows for the use of a  denser grid for the intermediate state integration and thereby  reduces the discretization error. Therefore, the total scattering rates can be calculated with negligible discretization error compared with the ab-initio calculation.
However, since the semi-analytical model is only for $\Gamma$-$\Gamma$ scattering, a discretized scattering matrix ($\Theta_{\boldsymbol{k}^\prime, \boldsymbol{k}}$) is still needed to compute the drift velocity and current PSD which are affected by $\Gamma$-L intervalley scattering.
Therefore, the discretization error in the final state integration cannot be avoided for the present calculations.
Nevertheless, the semi-analytical model still decreases the discretization error of the intermediate state integration and treats the full 2ph scattering term. The differences between the ab-initio calculation and the semi-analytical model for $\Gamma$ intravalley scattering are summarized in Table~\ref{tab:table}.

\begin{table}[htbp]
    \begin{ruledtabular}
    \begin{tabular}{ccc}
         & Ab-initio calculation  & Semi-analytical model \\ \hline
      \multirow{2}{*}{Final state integration} & \multirow{2}{*}{200$\times$200$\times$200} & 200$\times$200$\times$200 for observables \\
                  &                   & Effectively exact for scattering rates \\
      Intermediate state  integration &  200$\times$200$\times$200 & Effectively exact \\
      Processes & On-shell only & All processes included \\
      Two-phonon iterations & 3 & Iterate until convergence
      \\
      Interference term & Not included & Included
      \\
      Computational time (CPU hours) & 50000 & 40  \\
    \end{tabular}
    \end{ruledtabular}
    \caption{Comparison between the ab-initio calculation and the semi-analytical model for $\Gamma-\Gamma$ intravalley scattering.  The semi-analytical model improves upon the ab-initio model in all respects except the final state integration grid density for observables, for which the same grid is used.}
    \label{tab:table}
\end{table}

\section{Computational details} \label{sec:details}

\subsection{Ab-initio calculations}
The ab-initio calculation parameters are identical to those in our previous work \cite{cheng2022high}. In brief,  electronic structure and e-ph matrix elements are computed using Density Functional Theory (DFT) and Density Functional Perturbation Theory (DFPT) with \textsc{Quantum Espresso} \cite{giannozzi_qe_2009, giannozzi_qe_2017} with an  $8 \times 8 \times 8$ coarse grid, 72 Ryd plane wave energy cutoff, a relaxed lattice parameter of 5.556 $\textrm{\AA}$, and a non-degenerate carrier concentration of $10^{15}$ cm$^{-3}$.
Following our previous work \cite{cheng2022high}, we apply a band structure correction for both the $\Gamma$ valley and the L valley.
For the $\Gamma$ valley, we use a spherically symmetric Kane model band structure \cite{conwell_1968} with an experimental effective mass of 0.067$m_e$ and a non-parabolicity of 0.64 Ryd$^{-1}$ \cite{Lundstrom_2000}. For the L valley, we shift the energy by 50 meV to achieve the experimental $\Gamma$-L valley separation of 300 meV.
Wannier interpolation in \textsc{Perturbo}~\cite{perturbo_2021} is then applied to interpolate the e-ph matrix elements to a finer grid of $200 \times 200 \times 200$.
After the e-ph matrix elements are obtained, the 1ph and 2ph scattering matrices are computed according to Eqs.~\ref{eq:Theta_1ph} and \ref{eq:Theta_2ph}.
The delta functions in Eqs.~\ref{eq:Theta_1ph} and \ref{eq:Theta_2ph_subprocess} are approximated by a Gaussian function  with a standard deviation of 5 meV.

During the computation of scattering matrices, a phonon frequency cutoff of 2 meV is applied to neglect phonons with low frequencies.
An energy cutoff of 360 meV above the conduction band minimum is used to reduce the number of k points in the Brillouin zone integration.
The 2ph calculation applies the on-shell approximation by restricting the off-shell extent to 25 meV.
Following Ref.~\cite{warmelectrons}, the linear systems of equations representing the Boltzmann equation are solved by the generalized minimal residual method (GMRES).

\subsection{Semi-analytical model}
The band structure used in the semi-analytical model is the same as that in the ab-initio calculation.
The  LO phonon energy is taken to be $\omega_\text{LO} = 35$ meV. The function $g_{\text{LO}}(k)$ is obtained by a weighted averaged of $g_\text{LO}(\boldsymbol{k}_i, \boldsymbol{p}_i)$:
\begin{equation}
g_{\text{LO}}(k)= \frac{\sum_i g_\text{LO}(\boldsymbol{k}_i, \boldsymbol{p}_i)\exp(-\frac{(k-p_i)^2}{2\sigma^2}{})}
{\sum_i \exp(-\frac{(k-p_i)^2}{2\sigma^2}{})}
\end{equation}
where the standard deviation $\sigma=5\times 10^{-4}$ Ryd, the summation is over all the on-shell processes, and the $g_\text{LO}(\boldsymbol{k}_i, \boldsymbol{p}_i)$ are calculated by the Wannier interpolation.
The spherical coordinates integration in Eq.~\ref{eq:I_tilde} is defined such that the the $\theta_{\boldsymbol{k}_m}=0$ direction is orthogonal with both $\boldsymbol{k}$ and $\boldsymbol{k}^\prime$.
200 grid points are used for each $\theta_{\boldsymbol{k}_m}$ and $\phi_{\boldsymbol{k}_m}$ integration.
The radius integration in Eq.~\ref{eq:int_I_spherical} is transformed into the integration of $\epsilon_{k_m}$ and performed using an adaptive integration range with 120 grid points.
The integration range is $(\epsilon_{k^\prime} + \alpha_2 \hbar \omega_0 - 6 \hbar \Gamma_{k_m}, \epsilon_{k^\prime} + \alpha_2 \hbar \omega_0 + 6 \hbar \Gamma_{k_m})$, corresponding to a width of $12 \hbar \Gamma_{k_m}$ around the center of the Lorentzian function in the denominator of Eq.~\ref{eq:int_I_spherical}.
The relative residual error from this choice of integration limits is estimated as $1/12^2 \approx 0.7\%$.
In the final state integration of both 1ph and 2ph, the angular integrations in Eq.~\ref{eq:rates_int} are performed with  200 grid points. All the above numerical integrations are performed on uniform grids using the midpoint rule.
The 2ph calculation is performed with 20 iterations, with the interference term Eq.~\ref{eq:int_I_itf} updated every 10 iterations.
The relative difference between the 10th and 20th iterations is less than 1\%, indicating convergence of the iterative process. The discretization of $\Theta(k, k^\prime, \theta_{\kwv, \kwv^\prime})$ to $\Theta_{\boldsymbol{k}^\prime,\boldsymbol{k}}$ is performed by the regular grid interpolation provided in scipy \cite{2020SciPy-NMeth}.

\section{Results}

\subsection{ $\mathbf{\Gamma}$ valley scattering rates at 300 K} \label{sec:rates}

\begin{figure}
    \centering
    {
    \phantomsubcaption\label{fig:rate_300K-a}
    \phantomsubcaption\label{fig:rate_300K-b}
    \includegraphics[width=1\textwidth]{figures/rates_300K_all.png}
    }
    \caption{ (a) 1ph and 2ph $\Gamma$-valley  scattering rates at 300 K obtained from ab-initio calculations (symbols and shaded regions) and semi-analytical model (lines).  Note that intervalley scattering is excluded from these rates. Due to the variations of the ab-initio rates, we apply a Gaussian smearing and plot the shaded region to indicate the region within a standard deviation. (b) Decomposition of  semi-analytical 2ph rates (black solid line) to four contributing subprocesses (blue, orange, green, and red solid lines). In both (a) and (b), vertical dashed lines indicating energies of $\hbar \omega_\text{LO}$ and $2\hbar \omega_\text{LO}$ are plotted.}
    \label{fig:rate_300K}
\end{figure}

We first present the 1ph and 2ph $\Gamma$ valley scattering rates versus energy  obtained by ab-initio calculation and the semi-empirical model for GaAs at 300K in Fig.~\ref{fig:rate_300K-a}. For both 1ph and 2ph rates, the  ab-initio calculations and the semi-analytical model are in quantitative agreement. Specifically, we observe a rapid increase of the 1ph and 2ph rates at $\hbar \omega_\text{LO} \approx 35$ meV followed by a nearly constant trend. The degree of agreement between  the semi-analytical model and the mean values of the ab-initio calculation is notable considering the semi-analytical model includes off-shell processes and the interference term, both of which are neglected in the ab-initio calculations. The agreement can be attributed in part to the cancellation of errors between the limitation on the iteration number and the on-shell approximation in the ab-initio calculation. The third iteration of the 2ph rates yields values that are overestimated from the converged value by about 9\%, while the on-shell approximation and the non-interference approximation are found to underestimate the rates by around 3\% and 5\%, respectively. These approximations offset each other to yield good agreement between the two approaches. Overall, the contribution of off-shell processes is found to make only a minor contribution to the $\Gamma$ intravalley 2ph scattering rates even up to energies of 300 meV.

The major difference between the ab-initio and the semi-analytical rates is the variation of the individual rates in the ab-initio calculation in a given energy range, which is due to the relatively low grid density used in the ab-initio calculation.
As explained in Sec.~\ref{sec:semi-analytical_model}, the semi-analytical model uses a significantly finer grid, leading to negligible variations in individual scattering rates in the same energy range.
Although anisotropy could in principle lead to similar variations of the ab-initio rates, this contribution is negligible owing to the high spherical symmetry of the band structure and e-ph matrix elements  (about 3\% as mentioned in Sec.~\ref{sec:semi-analytical_model}).

The high grid density in the semi-analytical model enables features in the scattering rates to be observed that cannot be discerned in the ab-initio calculations, including the previously mentioned rapid increase of 1ph and 2ph rates at $\hbar \omega_\text{LO}$, and also a small but evident kink at $2\hbar \omega_\text{LO}$ (about 70 meV).

We now analyze each of these observations. For the 1ph rates, the increase at $\hbar \omega_\text{LO}$ is because LO phonon emission from an electron may only occur above an energy of  $\hbar \omega_\text{LO}$.
For the 2ph rates, the situation is more complicated due to the existence of four subprocesses (2e, e-a, a-e, 2a) in 2ph scattering.
To better understand the features in the 2ph scattering rates, the scattering rates of the four subprocesses are plotted separately in Fig.~\ref{fig:rate_300K-b}.
We observe that the increase of the total 2ph rates at $\hbar \omega_\text{LO}$ can be attributed to the e-a subprocess due to a similar reason with the emission subprocess in 1ph process, namely that the e-a subprocess requires the electron to have energy exceeding $\hbar \omega_\text{LO}$.
The kink at $2\hbar \omega_\text{LO}$ comes from the cancellation between the increase of the 2e rates and the decrease of e-a rates.
The increase of the 2e rates is due to the emission of 2 LO phonons at energies higher than $2\hbar \omega_\text{LO}$.
The decrease of the e-a rates is due to the increase of the intermediate state rates in the denominator of Eq.~\ref{eq:W_tilde} (or Eq.~\ref{eq:int_I_spherical}).
Specifically, for an e-a subprocess with $\epsilon_{\boldsymbol{k}}=2\hbar \omega_\text{LO}$, the corresponding intermediate state has the band energy $\epsilon_{\boldsymbol{k}+\boldsymbol{q}}=\hbar \omega_\text{LO}$, where an increase of both $\Gamma^{(\mathrm{1ph})}$ and $\Gamma^{(\mathrm{2ph})}$ occurs as explained above.
In fact, the effect of the intermediate state self-energy leads to change at any integer multiple of $\hbar \omega_\text{LO}$, but the magnitudes are decaying with increasing energy such that they cannot be observed at high energy.

Another observation in Fig.~\ref{fig:rate_300K-b} is the difference between the e-a rates and a-e rates asymptotically decreases to zero with increasing energy.
This can be understood by analysis of Eq.~~\ref{eq:Theta_2ph_integral}.
Specifically, in the process $\boldsymbol{k}_i \rightarrow \boldsymbol{k}_m \rightarrow \boldsymbol{k}_f$, the relative differences between norms of the state vectors $k_i$, $k_m$, $k_f$ become small at  high energies (above a few multiples the LO phonon energy), so that the factor $I^{(\alpha_1, \alpha_2)}(k, k^\prime, \theta_{\kwv, \kwv^\prime})$ becomes insensitive to the subprocess type.
Thus, their magnitudes are fully determined by the phonon occupation factor $A_{\alpha_1}A_{\alpha_2}$ defined in Eq.~\ref{eq:phonon_occupation} which satisfy $\Gamma^{(2a)} / A_{-1}^2 = \Gamma^{(e-a)} / A_{-1}A_{+1} = \Gamma^{(a-e)} / A_{-1}A_{+1} = \Gamma^{(2e)} / A_{+1}^2$.
From this relationship, we find $\Gamma^{(e-a)} = \Gamma^{(a-e)}$.
Furthermore, a common ratio $\Gamma^{(2a)} / \Gamma^{(e-a)} = \Gamma^{(e-a)} / \Gamma^{(2e)} = A_{-1}/A_{+1}$ can also be obtained for the subprocess rates at the high energy region.
This relationship is observed in Fig.~\ref{fig:rate_300K-b}.

\subsection{Drift velocity and current PSD at 300 K}

We now examine the transport and noise properties from each model.  For the semi-analytical model results, the $\Gamma\text{-}\Gamma$ block of the scattering matrix is calculated by the semi-analytical model in Eq.~\ref{eq:Theta_1ph_integral} and \ref{eq:Theta_2ph_integral}, while the $\Gamma$-L and L-L blocks are those of the ab-initio calculation. Figures \ref{fig:drift_psd-a} and \ref{fig:drift_psd-b} display the drift velocity and normalized current PSD, respectively, from the ab-initio calculations and the semi-analytical model. The experimental measurements are also plotted for comparison.
In Fig.~\ref{fig:drift_psd-a}, the ab-initio calculation and the semi-analytical model give a similar prediction of the drift velocity versus electric field up to 5 kVcm$^{-1}$.
Both give the low-field mobility of around 7000 cm$^2$V$^{-1}$s$^{-1}$, which agrees with the experimental value at about 8000 cm$^2$V$^{-1}$s$^{-1}$ to within around 15\%. The similarity between the ab-initio and semi-analytical results is expected due to the agreement of their scattering rates as discussed in Sec.~\ref{sec:rates}.

In Fig.~\ref{fig:drift_psd-b}, the PSD obtained from different experimental measurements reveal a non-monotonic pattern characterized by an initial decrease, followed by a marked rise around the commencement of negative differential mobility, and a subsequent decrease. The origin of this trend was explained in Ref.~\cite{cheng2022high}.
However, both the ab-initio calculation and the semi-analytical model predict the PSD to be nearly independent of electric field and thus fail to predict the characteristic PSD peak at about 3 kVcm$^{-1}$.

\begin{figure}
    \centering
    {
    \phantomsubcaption\label{fig:drift_psd-a}
    \phantomsubcaption\label{fig:drift_psd-b}
    \includegraphics[width=1\textwidth]{figures/observables.png}
    }
    \caption{Drift velocity and normalized PSD versus electric field for the (1+2)ph results obtained by ab-initio calculation (dashed blue line) and semi-analytical model (solid red line) at a temperature of 300 K. (a) Drift velocity versus electric field: The ab-initio calculation and semi-analytical model qualitatively agree with the measurements of Ruch et al. \cite{ruchkino1967} (filled circles) and Ashida et al. \cite{ashida1974} (open circles)  (b) Normalized PSD versus electric field. Values calculated using two approaches show consistent discrepancies compared with the results obtained from noise temperature and differential mobility measurements (filled circles, Ref.~\cite{gasquet1985} and open circles,  Ref.~\cite{Bareikis_1980}), and from time of flight experiments (triangles, Ref.~\cite{ruchkino1968}).}
    \label{fig:drift_psd}
\end{figure}

\subsection{Comparison of cryogenic $\mathbf{\Gamma}$ and $\mathbf{\Gamma}$-L scattering rates to experiment} \label{sec:helium}

\begin{figure}
    \centering
    {
    \phantomsubcaption\label{fig:helium-a}
    \phantomsubcaption\label{fig:helium-b}
    \includegraphics[width=1\textwidth]{figures/rates_helium_all.png}
    }
    \caption{Experimental (black symbols, Ref. \onlinecite{fasol1990}), ab-initio 1ph (orange symbols) and (1+2)ph (blue symbols), semi-analytical 1ph (green lines) and (1+2)ph results (red lines) of (a) $\Gamma$ valley scattering rates and (b) $\Gamma$-L intervalley scattering rates at helium temperatures. The (1+2)ph calculations agree better with experiment compared to the 1ph calculations in all cases. The (1+2)ph scattering rates agree with experiment to within 25\% and a factor of two for $\Gamma$ rates and $\Gamma$-L intervalley rates, respectively.}
    \label{fig:helium}
\end{figure}

The lifetimes of photoexcited hot electrons in GaAs have been experimentally measured at 10 K from an analysis of the linewidths of peaks from continuous-wave luminescence spectroscopy. \cite{fasol1990} In this section, we compare the ab-initio, semi-analytical and experimental scattering rates at cryogenic temperatures. In Fig.~\ref{fig:helium-a} and \ref{fig:helium-b}, we show the $\Gamma$ and $\Gamma$-L scattering rates obtained by ab-initio calculation, the semi-analytical model, and experiment. The calculations were performed at cryogenic temperatures to enable comparison with experiment. The experimental scattering rates and error bars are converted from the corresponding lifetimes and error bars directly reported in Ref.~\cite{fasol1990}.
Since the semi-analytical model is only valid for $\Gamma\text{-}\Gamma$ scattering, only the experimental and ab-initio results are shown in Fig.~\ref{fig:helium-b}.

Figure \ref{fig:helium-a} shows that the experimental and theoretical values for the $\Gamma$-valley scattering all yield a nearly constant value  between 100 meV and 325 meV. The experimental rates are about 8 ps$^{-1}$ in this energy range, while the semi-analytical 1ph and (1+2)ph calculations give about 3.5 ps$^{-1}$ and 6 ps$^{-1}$, respectively. This result affirms that the 2ph scattering makes a  non-negligible contribution to electron scattering in GaAs. Similarly with Sec.~\ref{sec:rates}, the ab-initio calculations give  mean values of the scattering rates that agree quantitatively with the semi-analytical model but with substantial scatter about the mean. Such observation again suggests that the approximations in the ab-initio calculations do not result in qualitative deviations.

Figure \ref{fig:helium-b} shows the $\Gamma$-L intervalley scattering rates of the ab-initio calculations and the photoluminescence experiments.
According to Ref.~\onlinecite{fasol1990}, the experimental $\Gamma$-L rates are obtained by
\begin{equation} \label{eq:rates_GL}
    \Gamma_{\Gamma-L}=\Gamma_{\text{tot}}-\Gamma_{\Gamma-\Gamma},
\end{equation}
where $\Gamma_{\Gamma-\Gamma}$ is taken as a constant estimated by fitting the data in Fig.~\ref{fig:helium-a}.
To make a consistent comparison with the experimental estimation, the ab-initio $\Gamma$-L intervalley scattering rates in Fig.~\ref{fig:helium-b} are also calculated by Eq.~\ref{eq:rates_GL} instead of being directly calculated from the scattering matrix.
It is found that the experimental intervalley rates is about 6 ps$^{-1}$ in the energy range from 340 meV and 400 meV, while the ab-initio 1ph and (1+2)ph give around 2 ps$^{-1}$ and 3 ps$^{-1}$, respectively.
Although the additional 2ph calculation decreases the deviation from experiment results, an underestimation of a factor of two is still observed. This discrepancy could be attributed to the need for an off-shell extent larger than 25 meV in the 2ph calculation, owing to the larger intermediate state scattering rates (appearing in the denominator of Eq.~\ref{eq:W_tilde} and Eq.~\ref{eq:int_I_spherical}) at energies above 300 meV. However, at present a larger off-shell extent is computationally infeasible.

\section{Discussion}

The semi-analytical model treats the full scattering term for $\Gamma$ intravalley 2ph scattering but does not qualitatively alter transport and noise properties. In particular, the marked discrepancy with the experimental PSD remains. We now examine alternate possible origins for the discrepancy.

\subsection{Underestimated $\mathbf{\Gamma}$-L intervalley scattering rates}

A comparison of our computed cryogenic rates with those measured from photoluminescence experiments indicates that the  $\Gamma$ rates agree to within 25\%, but the $\Gamma$-L intervalley rates are underestimated by around a factor of two.  Despite this underestimate, prior work suggests that this effect is unlikely to reconcile the PSD discrepancy. Specifically, Monte Carlo simulations of electron transport in GaAs with a three-valley $\Gamma$-L-X model have found that increased intervalley scattering suppresses the PSD feature (see Fig. 7 of Ref.~\cite{povzela1980electron}). Therefore, although the possibility cannot be definitively excluded at present, including intervalley scattering processes beyond those treated already is not expected to yield improved agreement with the PSD.

\subsection{Contribution from simultaneous electron-two-phonon interaction} \label{sec:e-2ph}

According to Ref.~\onlinecite{cheng2022high}, another possibility is that the contribution of electron-two-phonon (e-2ph) interaction \cite{kocevar_1980} is not considered.
Here, we make a qualitative estimation of the magnitude of this effect based on the Fröhlich mechanism for electron scattering.
According to Ref.~\onlinecite{kocevar_1980}, the e-ph Hamiltonian up to the second order can be written as
\begin{equation} \label{eq:eph-Hamiltonian}
    H_{\text{e-ph}} = \sum_{\boldsymbol{R}\kappa} \boldsymbol{u}_{\boldsymbol{R}\kappa} \cdot \nabla V(\boldsymbol{r}-\boldsymbol{R}_\kappa) + \frac{1}{2} \sum_{\kappa} \boldsymbol{u}_{\boldsymbol{R}\kappa} \cdot \nabla \nabla V(\boldsymbol{r}-\boldsymbol{R}_\kappa) \cdot \boldsymbol{u}_{\boldsymbol{R}\kappa},
\end{equation}
where $\kappa$ is the index of atom in a unit cell, $\boldsymbol{R}$ is the unit cell position, $\boldsymbol{u}_{\boldsymbol{R}\kappa}$ is the corresponding phonon-induced displacement, and $V$ is the electron potential. In the long wavelength limit, the electric potential for the Fr{\"o}hlich interaction can be obtained by assigning a point dipole to each atom \cite{verdi2015frohlich}. For acoustic phonons with the same displacements for atoms in the same unit cell, the net dipole moment will be zero and no scattering will occur.
Similarly, the electric potential for e-2ph interaction can be obtained by assigning a point quadrupole to each atom.
Following the same logic, the net quadrupole moment will be zero  if the quadrupoles are induced by two acoustic phonons or two optical phonons, which means that the simultaneous e-2ph interaction based on the Fr{\"o}hlich interaction can only be induced by an acoustic and optical phonon.

We estimate the order of magnitude of such simultaneous e-2ph interaction involving an acoustic and optical phonon. A full derivation can be found in the Appendix.
The final estimated scattering rates in the $\Gamma$ valley is
\begin{equation} \label{eq:e-2ph_rates}
\begin{aligned}
    \Gamma^{(\mathrm{e-2ph})}(k) &\sim
    \frac{8\pi^2}{\Omega_\text{BZ} \hbar}  \frac{k^3}{\epsilon_{k}} \left(\frac{eZ }{\Omega\epsilon_\infty} \sqrt{\frac{\hbar}{2M\omega_\text{A}}} \sqrt{\frac{\hbar}{2M\omega_\text{O}}}\right)^2 (N_\text{A} + 1) (N_\text{O} + 1),
\end{aligned}
\end{equation}
where $\omega_\text{O}$ and $\omega_\text{A}$ are frequencies of optical and acoustic phonons at the edge of Brillouin zone, $N_\text{O}$ and $N_\text{A}$ are the corresponding phonon occupations, $\Omega_\text{BZ}$ is the Brillouin zone volume, $M$ is the average atom mass in a unit cell, $\epsilon_\infty$ is the high-frequency permittivity, and $Z$ is the Born effective charge of a single atom.
For a typical $k$ such that $\epsilon_{k}=200$ meV, the e-2ph scattering rates at 300 K can be estimated as $10^{-2.5}$ ps\textsuperscript{-1}, which is about 3.5 orders of magnitude smaller than the ab-initio or semi-analytical 2ph scattering rates obtained in this work. Thus, we conclude that the effect of the simultaneous e-2ph interaction based on the Fr\"{o}hlich interaction can be neglected.

\subsection{Space charge domains and experimental non-idealities}

Finally, we consider an alternate mechanism for the PSD peak which does not rely solely on intervalley scattering. The earliest studies of negative differential resistance in GaAs arose from the observation of current instabilities at electric fields approaching a threshold value of around  3 kVcm$^{-1}$. \cite{gunn1964} These instabilities were  attributed to the formation of space charge domains associated with the negative differential resistance. The typical Boltzmann formalism used to describe charge transport from first principles would not include the contribution of such effects because it neglects real-space gradients and space charge effects which are essential to the instability.

Space charge instabilities manifest as current fluctuations, and therefore the nucleation of space-charge domains could explain the PSD peak around the threshold field. However, an inconsistency with this explanation is that the increase in PSD begins at a field below the threshold value for NDR as in Fig~\ref{fig:drift_psd}. This inconsistency can be accounted for by considering the possibility that the local electric field exceeds the threshold even though the average field does not. Such a possibility was investigated theoretically by McCumber and Chynoweth \cite{mccumber_1966}, who found that the dipole layer generation process was sensitive to inhomogeneities such as doping fluctuations that would arise from purely random Poisson statistics. The dipole layer was found to propagate  even if the average uniform field was less than the nominal threshold field based on the velocity-field characteristics.

These considerations support the explanation of the PSD peak in terms of  instabilities associated with the local electric field exceeding the threshold field for negative differential resistance space charge domain formation. This finding has implications for the use of transport studies to determine intervalley scattering strength. In particular, Monte Carlo methods have been used for decades for this purpose in GaAs by interpreting transport and noise measurements, and in those simulations, noise was assumed to arise from solely intervalley scattering. Our result indicates that this approach to studying intervalley scattering is not valid  because the physical origin of noise differs from that assumed in the model. Instead, photoluminescence methods which directly provide an electronic lifetime as in Ref.~\cite{fasol1990} should be employed.

\section{Conclusions}
We have introduced a semi-analytical model of 1ph and 2ph $\Gamma$ intravalley scattering for electrons in GaAs which allows for prior approximations in the treatment of the 2ph term for $\Gamma$ scattering to be lifted while incurring errors of  a few percent. We find that the calculated transport and noise properties are qualitatively unchanged from the ab-initio values. The computed drift velocity agrees with experiment to within 15\%, while agreement with measured cryogenic scattering rates are within 25\% for the $\Gamma$ valley scattering rates and a factor of two for  $\Gamma$-L intervalley scattering. However, the qualitative discrepancy for the  PSD  is not improved with the semi-analytical model. Considering these observations and prior work, we conclude that the PSD peak  arises from space charge domain formation rather than partition noise associated with intervalley scattering, as has been assumed for many decades. This result implies that care must be taken when  interpreting transport and electrical noise measurements in terms of intervalley scattering. Our findings highlight the insights into charge transport that can be obtained from a first-principles treatment of high-field charge transport and noise properties.

\section{Acknowledgements}

This work was supported by AFOSR under Grant Number FA9550-19-1-0321.
The authors thank P.~Cheng, B.~Hatanp\"a\"a, D.~Catherall, S-N.~Sun and T.~Esho for helpful discussions.