System identification through learning approaches is emerging as a promising strategy for understanding and simulating dynamical systems, which nevertheless faces considerable difficulty when confronted with power systems modeled by differential-algebraic equations (DAEs). This paper introduces a neural network (NN) framework for effectively learning and simulating solution trajectories of DAEs. The proposed framework leverages the synergy between Implicit Runge-Kutta (IRK) time-stepping schemes tailored for DAEs and NNs (including a differential NN (DNN)). The framework enforces an NN to cooperate with the algebraic equation of DAEs as hard constraints and is suitable for the identification of the ordinary differential equation (ODE)-modeled dynamic equation of DAEs using an existing penalty-based algorithm. Finally, the paper demonstrates the efficacy and precision of the proposed NN through the identification and simulation of solution trajectories for the considered DAE-modeled power system.

Differential-algebraic equations (DAEs), neural network (NN), power systems, differential NN (DNN).

§ INTRODUCTION

System identification for power networks is essential in the domain of power systems for many reasons: First, accurate models are necessary for effective design, analysis, and control, offering insights into system behaviors. Identifying model parameters allows engineers to further study power dynamics performance and stability, for example. Second, identification can produce predictive models for planning or optimizing the behaviors of power systems, which is important in advanced power models with renewables and emerging technologies. The identified power models obtained through identification techniques facilitate control strategies for enhancing the efficiency and resilience of the controlled systems.

Differential-algebraic equations (DAEs) are essential in power systems due to their ability to model dynamic behaviors and interconnections. They offer a solid mathematical framework capturing differential equations for dynamic responses and algebraic constraints between causal components. For instance, DAEs enable a well-formed representation of interconnected power system components, facilitating the modeling of generators and transformers. This advantage is beneficial as power grids evolve with renewable energy and advanced integration.

§.§ Existing works

The literature pays considerably less attention to system identification challenges in DAEs compared to, for instance, nonlinear ODE models (refer to <cit.> and <cit.> for more details). Existing methods for DAEs are primarily designed for scenarios of addressing numerical aspects within the optimization problem and involving environments with or without disturbances, under the condition that all inputs are known. In such cases, these methods typically have straightforward formulations, with a focus on addressing numerical aspects within the optimization problem (see, for instance, <cit.> and <cit.>, along with relevant references). Also, there exist some works studying the identification of DAEs with process disturbances, e.g., in <cit.>. However, the modeling of DAEs frequently encounters challenges in accurately representing nonlinearities through traditional methods. Neural networks (NNs), famous for their usefulness in capturing complex nonlinear interconnections, exhibit high performance, indicating they are ideal for accurately representing the intricate dynamics in DAE-modeled power systems.

NNs have extensively been investigated for learning dynamical systems that are described by ordinary differential equations (ODEs), considerably in nonlinear ODEs (NODEs), providing us with
an effective alternative to conventional numerical methods with high costs <cit.>.
Although NNs have succeeded in learning solution trajectories or simulating dynamic behaviors for ODEs <cit.>, developing an NN-based framework for learning and simulating solution trajectories in the context of nonlinear DAEs (NDAEs) remains as an unresolved challenge. The difficulty arises from the fact that DAEs exhibit infinite stiffness property <cit.>, which can result in gradient pathologies <cit.> and optimization problems with poor conditioning. These issues often lead to the failure of many trainings. To that end, observing the general forms of NDAEs (especially the dynamic system), in this paper, we introduce a differential NN (DNN) <cit.>, which is a continuous variant of an artificial NN (ANN). DNNs incorporate a feedback element, imparting a memory effect, enabling them to operate with the historical data of a process—a feature not feasible for static algorithms. The significant attribute of a DNN lies in its ability to identify unknown systems by incorporating established ODEs.

§.§ Main Contributions

The main contributions of the current paper can be summarized as follows:

* We present an NN scheme for identifying the parameters of an NDAE-modeled power system in which a DNN is included to approximate the dynamic equation of the system. Our approach can also be well-adapted to the general NDAEs with an index of $1$, which is the inherent physical property of the considered power model.

* To guarantee that the proposed NN architecture is effective, under some standard assumptions (imposed on nonlinearities), we further analyze the error dynamic between the true system and the proposed DNN and give the conditions for bounding the identification error at the infinite-time horizon. This is demonstrated by a numerical simulation. We also illustrate that the imposed conditions for identification error analysis are milder than the existing ones in the literature.

§.§ Paper organization

The rest of this paper is organized as follows: Section <ref> presents the considered power system and an overview of the system identification problem. In Section <ref>, we first give the details of the proposed NN scheme and the numerical methods for training and testing and then present the boundedness conditions for the identification errors. To illustrate the efficacy of the identification method, we present a case study of the considered DAE-modeled power system in Section <ref>. The paper summary and future research directions in Section <ref> conclude the paper. The used notation is provided next:

§.§ Notation

The set of real numbers (respectively non-negative real numbers) is
denoted by $\mathbb{R}$ (respectively $\mathbb{R}_{\geq 0}$). The symbol $\rVert \cdot \rVert$ denotes the Euclidean norm on
$\mathbb{R}^{n}$ (and the induced matrix norm $\rVert A \rVert$ for a matrix
$A\in\mathbb{R}^{m\times n}$). We use
$I_{n}$ to denote the $n\times n$
identity matrix. The minimum eigenvalue of a symmetric matrix $P$ is denoted by $\lambda_{\min}(P)$. To represent the positive definiteness and negative semi-definiteness, we utilize $\succ 0$ and $\preceq 0$, respectively. For $t_{1}$, $t_{2}\in \mathbb{R}$, with $t_{1}<t_{2}$, we
denote by $C^{1}_{n}([t_{1},t_{2}])$ the Banach space of continuously differentiable
functions $\psi:[t_{1},t_{2}] \to\mathbb{R}^{n}$ with the norm
$\Vert\psi\Vert_{[t_{1},t_{2}]}=\sup_{ r \in [t_1,t_2]} \|\psi(r)\| + \sup_{ r \in [t_1,t_2]}\|\frac{d \psi(r)}{d r}\| < +\infty$. For brevity, we denote a time-dependent signal $s(t)$ by $s$ wherever needed (e.g., $x_d$ instead of $x_d(t)$).

§ PROBLEM STATEMENT

In this paper, we consider the system identification problem for the standard $4^{\text{th}}$-order power system model <cit.>, which is represented via the following set of equations:

\begin{align}
\begin{split} \label{main_DAE_system_1}
\textit{Generator ODEs}: \;\;\;\;\; E_d \dot{x}_d &= A_d x_d + C_d f(x_d,x_a) +  B u + hw_0, \end{split}\\
\begin{split}  \label{main_DAE_system_2}
\textit{Power flow equations}: \;\;\;\;\;\; 0 & = A_a x_a + C_a g(x_d,x_a),
\end{split}
\end{align}

where $x_d \in\mathbb{R}^{n_d}$ contains the dynamic states, $x_a \in\mathbb{R}^{n_a}$ contains the algebraic variables, $A_d$, $C_d$, $B$, $A_a$, and $C_a$ are constant system matrices with appropriate dimensions, $u \in \mathbb{R}^m$ contains the control inputs while the functions $f: \mathbb{R}^{n_d} \times \mathbb{R}^{n_a} \to \mathbb{R}^{n_{f}}$ and $g: \mathbb{R}^{n_d} \times \mathbb{R}^{n_a} \to \mathbb{R}^{n_{g}}$ (it is continuously differentiable) are the nonlinearities in dynamic and algebraic equations, respectively; $E_d = I_{n_d}$ holds and $w_0 \in \mathbb{R}$ is the synchronous speed; the vector $h$ satisfies $h \in \mathbb{R}^{n_d}$. In the above power system model, Eq. (<ref>) lumps the dynamic equations (set of ODEs) for all the synchronous generators while (<ref>) includes the power flow/balance equations (the algebraic constraints) of the power grid. Due to space limitations, the detailed description/explanation of the power system model is not included in this work and can be found in the power system literature, e.g., see the references <cit.>.

That being said, w.l.o.g., we let $f(0,0)=0$ and $g(0,0)=0$ and assume that the DAE (<ref>) is with index $1$ (this is a physical property of the power system (<ref>)), which means that the Jacobian $\frac{\partial \tilde{g}}{\partial x_a}$ is invertible for $\tilde{g}(x_d,x_a) :=  A_a x_a + C_a g(x_d,x_a)$. This implies that, by Implicit Function Theorem <cit.>, for the algebraic Eq. (<ref>), there exists a unique solution: $x_a = \ell(x_d)$, where $\ell$ is a function. If we substitute $\ell(x_d)$ into $x_a$ in (<ref>), we then get $E_d \dot{x}_d = \dot{x}_d = A_d x_d + C_d f(x_d,\ell(x_d)) +  B u + hw_0$, which is a standard NODE, and the corresponding system identification techniques for NODEs can be applied in this case.

Throughout the paper, we deal with the system identification problem for (<ref>), specifically, the black-box identification. To collect data sets for identifications, we also consider the output of the NDAE system (<ref>) as follows:

\begin{align} \label{eq:output_data}
y &= \begin{bmatrix}
x_d^T & x_a^T
\end{bmatrix}^T.
\end{align}

For system identification, the model (<ref>) is further assumed to be smoothly parametrized by a
constant unknown parameter vector $\theta \in \mathbb{R}^{n_\theta}$, which is to be estimated.  Notice that in the implementation of identification, the input $u(t)$ is known and given by the human operator, while the output $y(t)$ has to be sampled at each time instant denoted by $t_k$. Then, we have the following complete data collection:

\begin{align*}
N_{\eta,T} &:= \{ (y(t_{k}), u(s)) \mid k=1,\dots,\eta, \; s \in [0,T] \},
\end{align*}

for some $T>0$. The main goal of this work is to estimate the parameter $\theta$ by using the data set $N_{\eta,T}$. To that end, we first propose an NN architecture for identifying the dynamic and algebraic models and then analyze the corresponding identification error between the true output and the predicted output generated by the proposed NN.

§ PROPOSED METHOD

This section presents the proposed method for identifying the dynamic and algebraic equations of the NDAE (<ref>). For the ODE-modeled dynamic system, a DNN scheme is employed, while for the algebraic one, a standard feedforward neural network (FNN) is utilized.

§.§ Implicit Runge-Kutta (IRK) method

In the system identification of (<ref>), we assume that the data of $\{x_d(t_k),x_a(t_k)\}_{k=1}^\eta$ is accessible to the human operator in discrete time. Also, numerical analysis methods can be used to approximate solutions of (<ref>). In the current paper, we utilize the Implicit Runge-Kutta (IRK) method with $\nu$ stages expressed as follows:

\begin{align}
\begin{split} \label{Runge-Kutta_1}
\alpha_j&= x_d^n + \Delta \sum_{i=1}^\nu b_{j,i} \tilde{f}\left(\alpha_i, \beta_i, u\right), \quad  0  = \tilde{g}(\alpha_j, \; \beta_j), \end{split}  \\
\begin{split}  \label{Runge-Kutta_3}
x_d^{n+1} &= x_d^n + \Delta \sum_{j=1}^\nu c_{j} \tilde{f} \left(\alpha_j, \beta_j,  u\right), ~\quad 0  = \tilde{g}(x_d^{n+1},\; x_a^{n+1}),
\end{split}
\end{align}

where $j = 1,\dots,\nu$; $\Delta>0$ is sufficiently small, $x_d^{n} :=x_d(t_{n})$,

\begin{align*}
\alpha_j &:= x_d\left(t_n+ \Delta \sum_{k=1}^\nu b_{j,k}\right),  \quad  \beta_j := x_a\left(t_n+ \Delta \sum_{k=1}^\nu b_{j,k}\right),  \\ \tilde{f}(x_d,x_a,u) &:= A_d x_d + C_d f(x_d,x_a) +  B u + hw_0, \quad
\tilde{g}(x_d,x_a) := A_a x_a + C_a g(x_d,x_a).
\end{align*}

Note that $\{b_{j,i},c_j\}$ are known parameters and since $E_d =  I_{n_d}$ holds, $E_d$ is omitted.

For the IRK scheme (<ref>), after parametrization (with the denotation $\theta$; for example, $\alpha_j$ is changed to $\alpha_j^\theta$), we can denote the dynamic states and algebraic variables by

\begin{align*}
[\alpha_1^\theta,\; \dots, \; \alpha_\nu^\theta, \; x_d^{n+1,\theta}],\quad
[\beta_1^\theta,\; \dots, \; \beta_\nu^\theta, \; x_a^{n+1,\theta}].
\end{align*}

These can be used in the training or testing process.

§.§ Architecture of neural networks

A schematic design of the proposed system identification technique.

In this section, we consider an architecture of NNs for identifying the DAE system (<ref>), in which the DNN is for identifying the dynamic Eq. (<ref>) and another FNN is for approximating the algebraic constraint (<ref>). This scheme is depicted in Fig. <ref>, in which the expression of the DNN is: $ \dot{x}_{nn} = A_{nn} x_{nn} + B_{nn} \hat{\rho}(x_{nn}) + h \omega_0 + C_{nn} \gamma(u_{nn})$, where the function $\hat{\rho}$ represents the upper NN; the matrices $A_{nn}, B_{nn}, C_{nn}$, the operator $\gamma$, the input $u_{nn}$ are defined as in (<ref>), and $x_{nn}$ is the state of the DNN (it approximates $x_d$). Also, in the identification processes, the matrices or functions $A_{nn}, B_{nn}, C_{nn}, \gamma, u_{nn}, h, \omega_0$ are assumed to be known.

§.§ Estimation method

In the identification process, the two constraints have to be satisfied (i.e., the NNs are restricted to satisfy the differential and algebraic equations described by (<ref>)). At the beginning, we select the initial conditions of NNs in the whole considered domain (the initial conditions are randomly distributed). Also, the output training or testing data (i.e., $x_d^n, x_a^n$) are measured by the human operator so that we have two training sets of size $\eta>0$ for the two NNs (one for the dynamic system (<ref>), and another one for the algebraic Eq. (<ref>)), respectively, that is

\begin{align*}
D_{\eta}^d := \{ x_d^{n,1},  \dots, x_d^{n,\eta} \},\quad D_{\eta}^a := \{ x_a^{n,1},  \dots, x_a^{n,\eta} \}.
\end{align*}

Then, to analyze the difference between the NNs and the IRK scheme dynamics (<ref>), we consider the following loss function:

\begin{equation*}
\mathcal{L}({\theta,\eta}) :=  l(\mathcal{L}_d({\theta,\eta}), \mathcal{L}_a({\theta,\eta})),
\end{equation*}

where $l:\mathbb{R}_{\geq 0} \times \mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}$ is an operator with $l(0,0)=0$, and the functions $l(s,\cdot), l(\cdot,s)$ increase as the variable $s \in \mathbb{R}_{\geq 0}$ increases, where $\cdot $ can be any fixed non-negative real number.

In this paper, $\mathcal{L}_d({\theta,\eta})$ and $\mathcal{L}_a({\theta,\eta})$ are specifically defined as

\begin{align*}
\mathcal{L}_d({\theta,\eta}) &:=  \frac{1}{\eta (\nu+1) } \sum_{x_d^n \in D^d_{\eta}} \; \sum_{j=1}^{\nu+1} \left \| x_d^n - \hat{x}_d^{n,j}(\theta) \right\|^2, \\
\mathcal{L}_a({\theta,\eta}) &:=  \frac{1}{ \nu+1}  \tilde{\ell} \left(\sum_{j=1}^{\nu}  \| \tilde{g}(\alpha_j^\theta,\beta_j^\theta) \| + \left( \| \tilde{g}(x_d^{n+1,\theta},x_a^{n+1,\theta}) \| \right) \right),
\end{align*}

\begin{align*}
& \hat{x}_d^{n,j}(\theta) := \alpha_j^\theta - h\sum_{i=1}^\nu b_{j,i} \tilde{f}\left(\alpha_i^\theta, \beta_i^\theta, u \right), \quad
\hat{x}_d^{n,\nu+1}(\theta) := x_d^{n+1,\theta} - h \sum_{j=1}^\nu c_{j} \tilde{f}(\alpha_j^\theta, \beta_j^\theta, u),
\end{align*}

and the function $\tilde{\ell}(s)$ increases as the variable $s \in \mathbb{R}_{\geq 0}$ increases.

Next, we use the two NNs (as depicted in Fig. <ref>) to identify the dynamic and algebraic equations. We train the parameters of the considered NNs by minimizing the loss function using a gradient-based optimizer <cit.>, which can be expressed as

\begin{align*}
\hat{\theta} &= \mathop{\arg  \min_{\theta} } \; \mathcal{L} ({\theta,\eta}).
\end{align*}

Note that the selection of the explicit form of the operator $\mathcal{L}$ can be tricky since there are two main constraints originating from the DAE system (<ref>), i.e., the solution of (<ref>) should satisfy the dynamic system (<ref>) and stay within the manifold defined as follows:

\begin{align*}
M &:= \{ (x_d,x_a) \mid \tilde{g}(x_d,x_a) = 0 \},
\end{align*}

which is deduced from (<ref>). For that purpose, we select a recurrent penalty method in the sequel to get a well-conditioned loss function. Also, for simplicity, we select the loss function as

\begin{align*}
\mathcal{L}({\theta,\eta}) &=  l(\mathcal{L}_d({\theta,\eta}), \mathcal{L}_a({\theta,\eta})) := w_d \mathcal{L}_d({\theta,\eta}) + w_a \mathcal{L}_a({\theta,\eta}).
\end{align*}

It is noteworthy that we employ an existing algorithm in the literature (Algorithm $2.1$ in <cit.> or Algorithm $1$ in <cit.>) for determining the values of $w_d, w_a >0$. Due to space limitations, we are unable to provide the details of such an algorithm in the current paper.

§.§ Identification error analysis of DNN

In this section, we consider the following structure of general DNN (the DNN shown in Fig. <ref> is a special case) for identifying (<ref>):

\begin{align} \label{eq:DNN}
\dot{x}_{nn} &= A_{nn} x_{nn} + B_{nn} \rho(x_{nn}) + C_{nn} \gamma(u_{nn}),
\end{align}

where $A_{nn} \in \mathbb{R}^{n \times n} $ is the state matrix of the DNN; $B_{nn} \in \mathbb{R}^{n \times n_B}, C_{nn} \in \mathbb{R}^{n \times n_C} $ are also matrices; the function $\rho: \mathbb{R}^{n} \to \mathbb{R}^{n_B}$ is nonlinear, the operator $\gamma: \mathbb{R}^{n_u} \to \mathbb{R}^{n_C}$, and the input $u_{nn} := u$ (thus, $n_u = m$), here $u$ contains the known inputs given in (<ref>). The subscript $``nn"$ represents the NN. Then, we analyze the identification error after the NN
parameters have been trained using the previous
learning method.

In the identification, one can assume the system state can be measured, and the dimensions of the input and states (i.e., $x_d, x_a$) of the true system are equal to those of the DNN. Under these mild assumptions, we then study the dynamic behavior of the error

\begin{align*}
e &:= x_d - x_{nn}.
\end{align*}

The resulting error dynamics can be expressed as follows:

\begin{align} \label{eq:error_system}
\dot{e} &= Ae + \tilde{\phi}(e,x_d,x_a,u,t),
\end{align}

where $A$ is a Hurwitz matrix, the function $\tilde{\phi}(e,x_d,x_a,u,t) := A_d x_d + C_d f(x_d,x_a) +  B u + hw_0 -  A_{nn} x_{nn} - B_{nn} \rho(x_{nn}) - C_{nn} \gamma(u_{nn}) - Ae$. In the expression of $\tilde{\phi}(e,x_d,x_a,u,t)$, we have the variable $x_a$, which is estimated by the previously-proposed FNN (see Fig. <ref> for clarification). Furthermore, there exists a function $\ell$ such that $x_a = \ell(x_d)$ holds under the condition that the function $\tilde{g} $ is continuously differentiable and its Jacobian is invertible (or the DAE system (<ref>) is with index $1$). Then, the variable $x_a$ in $\tilde{\phi}(e,x_d,x_a,u,t)$ can be substituted by $\ell(x_d)$, thus obtaining

\begin{align}\label{eq:error_system_replace}
\phi(e,x_d,u,t) := ~& A_d x_d + C_d f(x_d,\ell(x_d)) +  B u + hw_0 \nonumber \\
& -A_{nn} x_d - (A-A_{nn})e - B_{nn} \rho(x_d-e) - C_{nn} \gamma(u).
\end{align}

For the system (<ref>), several assumptions are imposed in this work to guarantee that the error $e$ is bounded. We first introduce the following assumption on the nonlinearity $\phi$:

There exist matrices $L \succ 0$ and $K \succ 0$ such that the following inequality:

\begin{align*}
& \phi(e,x_d,u,t)^\top  L \phi(e,x_d,u,t) \leq c_0(x_d,u,t)  + c_1(x_d,u,t) e^\top K e,
\end{align*}

holds where the real numbers $c_0, c_1 >0$ are bounded by

\begin{align*}
& \sup \; c_0(x_d,u,t) = c^0 < +\infty, \quad \sup \; c_1(x_d,u,t) = c^1 < +\infty,  \\
& \forall x_d \in C^1_{n_d}(\mathbb{R}_{\geq 0}), \;
u \in C^1_m(\mathbb{R}_{\geq 0}), \; t \geq 0.
\end{align*}

Note that Assumption <ref> is non-restrictive, since the right-hand side terms of (<ref>) can respectively be upper bounded by

\begin{align*}
& (A_d -  A_{nn}) x_d + C_d f(x_d,\ell(x_d)) +  B u + hw_0- B_{nn} \rho(x_d) - C_{nn} \gamma(u)  \leq \; c_0(x_d,u,t), \\
&- (A-A_{nn})e + \left(B_{nn} \rho(x_d)- B_{nn} \rho(x_d-e)  \right) \leq \; c_1(x_d,u,t) e^\top K e, \\
& \forall x_d \in C^1_{n_d}(\mathbb{R}_{\geq 0}), \;
u \in C^1_m(\mathbb{R}_{\geq 0}),  \; t \geq 0,
\end{align*}

under the conditions that $\rho$ is a Lipschitz function, the norm of $\gamma(t)$ at $t \geq 0$ is finite, and $f(s(t),\ell(s(t)) $ $< +\infty$ holds for $\| s \|_{\mathbb{R}} < +\infty$. Therefore, we then impose the second assumption on the nonlinearity $\rho$ and the operators $\gamma, f$ as follows:

The nonlinearity $\rho$  is a Lipschitz function, the operator $\gamma$ satisfies $\| \gamma \|_{\mathbb{R}} < +\infty$, and $f(s(t),\ell(s(t)) < +\infty$ holds for $\| s \|_{\mathbb{R}} < +\infty$.

The condition on the norm $\| \cdot \|_{\mathbb{R}}$ is more restrictive than $\| \cdot \|_{[r^1,r^2]}$, where $r^1:= \inf_{t \geq 0}  r(t), \; r^2:= \sup_{t \geq 0}   r(t)$, where $r(t)$ is a scalar-valued function. The latter case can also be considered for further relaxation for the system runs on the finite time interval $[r^1,r^2]$.

We are now ready to give the third assumption:

Assume that there exist matrices $P = P^\top \succ 0, W = W^\top \succ 0$ such that for the system (<ref>), the following matrix inequality:

\[
A^\top P+ PA + P L^{-1} P + c^1 K  +  W \preceq 0,
\]

has a solution $(P,W)$, where $L,c^1,K$ can be found in Assumption <ref>.

Assumptions <ref>-<ref> are used in the proof of the following proposition:

Let Assumptions <ref>-<ref> be satisfied. For the NDAE system (<ref>) and the DNN (<ref>), the error $e(t)$ satisfies the following property:

\begin{align*}
& \lim_{t \to +\infty} \sup \| e(t) \| \leq \sqrt{\frac{c^0}{\lambda_{\min}(P) \lambda_{\min}(P^{-\frac{1}{2}} W P^{-\frac{1}{2}}) }}.
\end{align*}

where $P^{-\frac{1}{2}}$ denotes the principal square root of $P^{-1} \succ 0$.

Consider the Lyapunov function

\begin{align*}
& V(e) = e^\top P e, \; \text{where} \; P = P^\top \succ 0.
\end{align*}

Taking the time derivative of $V(e)$, we have

\begin{align*}
& \dot{V}(e) = 2 e^\top P \dot{e} = 2 e^\top P \Big( Ae + \phi(e,x_d,u,t) \Big) = e^\top (A^\top P + PA)e + 2 e^\top P \phi .
\end{align*}

Here, for the term $2 e^\top P \phi$, by Matrix Young Inequalities <cit.>, we derive that for any $L = L^\top \succ 0$ (w.l.o.g., we select the matrix $L$ in Assumption <ref>), the following equality:

\begin{align*}
& 2 e^\top P \phi = (e^\top P) \phi + \phi^\top (P e) \leq \phi^\top L \phi+ e^\top P L^{-1} P e,
\end{align*}

holds. Applying Assumption <ref>, we get

\begin{align*}
& 2 e^\top P \phi \leq c_0(x_d,u,t)  + c_1(x_d,u,t) e^\top K e + e^\top P L^{-1} P e.
\end{align*}

Therefore, it can be deduced that

\begin{align*}
\dot{V}(e) \leq ~& e^\top (A^\top P + PA)e + c_0(x_d,u,t)  + c_1(x_d,u,t) e^\top K e + e^\top P L^{-1} P e \\
\leq ~& e^\top \Big( A^\top P + PA + P L^{-1} P + c^1 K \Big) e +  c^0.
\end{align*}

Then, by Assumption <ref>, we have the following inequality:

\begin{align*}
\dot{V}(e)  \leq & - e^\top W e + c^0
\leq  - e^\top P^{\frac{1}{2}} \Big( P^{-\frac{1}{2}} W P^{-\frac{1}{2}} \Big) P^{\frac{1}{2}} e + c^0
\leq  -\lambda_{\min}(P^{-\frac{1}{2}} W P^{-\frac{1}{2}})  V(e) + c^0.
\end{align*}

Let $\tilde{W}:= P^{-\frac{1}{2}} W P^{-\frac{1}{2}}$. We solve the above inequality and obtain

\begin{align} \label{eq:final_ineq}
\lambda_{\min}(P) \| e(t) \|^2  \leq ~& V(0) e^{-\lambda_{\min}(\tilde{W}) t} + \frac{c^0}{\lambda_{\min}(\tilde{W})} e^{-\lambda_{\min}(\tilde{W}) t} \Big( e^{\lambda_{\min}(\tilde{W}) t} -1 \Big)
\leq \frac{c^0}{\lambda_{\min}(\tilde{W})},
\end{align}

due to $0< e^{-\lambda_{\min}(\tilde{W}) t} \Big( e^{\lambda_{\min}(\tilde{W}) t} -1 \Big) \leq 1$ (observe that $\lambda_{\min}(\tilde{W}) \ge 0$ holds as $-\tilde{W} \preceq 0$ is satisfied.). The inequality (<ref>) results in the property in Proposition <ref>. This completes the proof.

It is noteworthy that in <cit.>, the identification error analysis of the same DNN is considered. However, it imposes more restrictive conditions (or assumptions) on the DNN; for instance, in the current paper, we assume that $ x_d \in C^1_{n_d}(\mathbb{R}_{\geq 0})$ holds, so that the condition: $f(x_d, \ell({x_d}))$ $< +\infty$ is milder since $f(x_d, \ell({x_d})) < +\infty$ for all $x_d \in \mathbb{R}^{n_d}$ greatly limits the scope of the operator $f$, while $ x_d \in C^1_{n_d}(\mathbb{R}_{\geq 0})$ is also more consistent with the physical properties of the DNN, compared with $x_d \in \mathbb{R}^{n_d}$.

§ CASE STUDY

In this section, we conduct system identification simulations for the power system modeled by the NDAE (<ref>). In particular, we consider the IEEE 9-bus test system which is most commonly used in various power system studies <cit.>. To ensure accurate and efficient identifications, we use the numerical simulation environment composed of Windows $11$, MATLAB R$2023$a, and a $12$th Gen Intel Core i$9$-$12900$H processor with $64$GB RAM.

In the simulations, we use the MATLAB index $1$ DAE solver  to simulate NDAE models, which is configured with the following parameter values: relative tolerance = $10^{-5}$, absolute tolerance = $10^{-6}$, and a maximum step size of $10^{-3}$.

[The trajectory of relative error $e_a^r(t)$.]   [The trajectory of relative error $e_d^r(t)$.]

[The comparison of $x_{d,4}$ and $\hat{x}_{d,4}$ (in the case of minimum identification error).]
[The comparison of $x_{d,5}$ and $\hat{x}_{d,5}$ (in the case of maximum identification error).]
The identification results of the proposed NN-based system identification technique.

In the model of (<ref>), the variables $x_d,x_a$ are in the spaces $\mathbb{R}^{12}, \mathbb{R}^{24}$, respectively. We randomly select the initial conditions for (<ref>) and measure the output $y^1$ as defined in (<ref>), under a given input $u$ (here, we assume that none of the parameters of (<ref>) is known).  From $y^1$, we then obtain the data set of $\{x_d^{n,i}, x_a^{n,i}\}_{i=0}^\eta$, where the pairs $\{x_d^{n,i},x_a^{n,i} \}$ are again randomly selected and close to each other.

First, an FNN is utilized (the relevant MATLAB built-in function is utilized) to approximate the following algebraic expression:

\begin{align*}
x_a & = \hat{\ell}(x_d),
\end{align*}

where the explicit form of the function $\hat{\ell}$ is to be estimated. The corresponding simulation results of the identification are depicted by Fig. <ref>, in which the relative error $e_a^r(t)$ is defined as follows:

\begin{align*}
e_a^r(t)~(\%) := 100 \times \frac{\| x_a(t) - \hat{x}_a(t) \| }{ \| x_a(t) \|} \Big(\text{respectively for $x_d$:} \; e_d^r(t)~(\%) := 100 \times \frac{\| x_d(t) - \hat{x}_d(t) \| }{ \| x_d(t) \|} \Big),
\end{align*}

where $\hat{x}_a$ represents the predicted value of $x_a$ in the identified model of (<ref>).

The second step is to identify the dynamic Eq. (<ref>). Similarly, $\hat{x}_d$ represents the predicted value of $x_d$ in the identified model of (<ref>). In this experiment, we have the following dynamic equation:

\begin{align*}
\dot{x}_d &= A_d x_d + C_d f(x_d,\hat{\ell}(x_d)) +  B u + hw_0,
\end{align*}

to be identified. To that end, we use the proposed NN structure (see Fig. <ref>, in which the output of the FNN is delivered to the inputs of the DNN). We implement the identification of (<ref>) and obtain the identification errors. We visualize the trajectory of relative error $e^r_d$ in Fig. <ref>. It illustrates that the identification error remains small. Furthermore, for brevity, we only visualize the minimum and maximum identification errors, which are $4^{\text{th}}$ and $5^{\text{th}}$ variables of $x_d$ (denoted by $x_{d,4}$ and $x_{d,5}$), respectively. The trajectories of $x_{d,4}$, $\hat{x}_{d,4}$ are shown in Fig. <ref>, by which we see that the error between the true output $x_{d,4}$ and its predicted value $\hat{x}_{d,4}$ is also small. Moreover, the trajectories of $x_{d,5}$, $\hat{x}_{d,5}$ are presented in Fig. <ref>, and one can realize that although the corresponding error is larger, it still illustrates the good performance of the proposed DNN. The above simulation results corroborate that our proposed NN scheme has satisfactory performance in the identification of the NDAE-modeled power system (<ref>).

§ CONCLUSION

This paper introduces an NN framework designed to effectively identify and simulate solution trajectories of NDAEs, such as those characterizing the dynamics of power networks known for exhibiting infinite stiffness. The proposed framework capitalizes on the synergy between IRK time-stepping schemes tailored for NDAEs and NNs (with a DNN). The framework compels an NN to adhere to the algebraic equation of NDAEs as hard constraints and its DNN component is suitable for simulating the ODE-modeled dynamic equation of NDAEs, using a penalty-based method. The paper showcases the efficacy of the proposed NN via the identification and simulation of solution trajectories for the considered DAE-modeled power system. Future work will consider the effect of disturbances that
are present in models, identifying power systems with partially unknown parameters and extending the proposed results to advanced power networks.

This work is supported by the National Science Foundation under
Grants 2152450 and 2151571.


Charged particle to photon ratio in heavy-ion collisions at RHIC: a different
point of view

Mohammad Asif Bhat111Corresponding Author<EMAIL_ADDRESS><EMAIL_ADDRESS>and Supriya Das

Department of Physics and Centre for Astroparticle Physics and Space Science
(CAPSS), Bose Institute, EN-80, Sector V, Kolkata 700091, India

Abstract

We report an alternative explanation of the variation of charged particle to
photon ratio in Cu+Cu collisions and Au+Au collisions from 1.4 $\pm$ 0.1 to
1.2 $\pm$ 0.1 at $\sqrt{\it{s}_{\rm NN}}$ = 62.4 GeV and 200 GeV respectively
in STAR experiment at RHIC. Based on the theoretical predictions and
experimental results we argue that the additional contribution of direct
photons at $\sqrt{\it{s}_{\rm NN}}$ = 200 GeV in both Cu+Cu and Au+Au
collisions are improving the ratio $\frac{N_{ch}}{N_{\gamma}}$ to 1.2 $\pm$
0.1 from 1.4 $\pm$ 0.1.

##  1 Introduction

High energy experiments are being performed for last so many years to study
the Quark Gluon Plasma (QGP) which is a strongly interacting hot and dense
medium of quarks and gluons having very short life time. Measurement of the
multiplicity of charged particles provides the information about the
underlying physics processes occuring during the medium formation such as
fragmentation [1, 2] and hadronization [3, 4]. Measurement of photon
multiplicity provides the complementary information to that of charged
particles thus helps in understanding of QGP medium [5, 6].

The measurement of the particle density in pseudorapidity $\eta$ is the
conventional way of describing particle production in heavy ion collisions.
The information of energy density, initial temperature, and velocity of sound
in the medium formed in the collisions is provided by the particle density in
pseudorapidity within the framework of certain model assumptions [7]. The
longitudinal flow and rescattering effects [8] strongly modifies the widths of
the pseudorapidity distributions. The relative importance of soft versus hard
processes in particle production can be studied by studying the variation of
particle density in $\eta$ with collision centrality, expressed in terms of
the number of participating nucleons $N_{part}$ and/or the number of binary
collisions $N_{coll}$. The test ground for various particle production models,
such as those based on ideas of parton saturation [9] and semiclassical QCD,
also known as the color glass condensate (CGC) [10] is provided by the
particle density in pseudorapidity $\eta$. The particle production mechanism
could be different in different regions of pseudorapidity at Relativistic
Heavy-ion Collider (RHIC). A significant increase in charged particle
production normalized to the number of participating nucleons has been
observed from peripheral to central Au+Au collisions [11] at midrapidity. This
is because of the onset of hard scattering processes, which scales with the
number of binary collisions. However, the total charged particle multiplicity
per participant pair, integrated over the whole pseudorapidity range, is
independent of centrality in Au+Au collisions. The centrality dependence of
particle production at midrapidity reflects the increase of gluon density due
to the decrease in the effective strong coupling constant in the framework of
the color glass condensate picture of particle production.

The detailed study of the increase in particle production at midrapidity with
increase in the center-of-mass energy has been done at Relativistic Heavy-ion
Collider (RHIC) [12]. The multiplicity of hadrons and their energy, centrality
and rapidity dependence from experimental data so far have been consistent
with the approach based on ideas of parton saturation. It has been argued that
this onset of saturation occurs somewhere in the center-of-mass energy
$\sqrt{\it{s}_{\rm NN}}$ = 17 GeV to 130 GeV [13].

##  2 Photon and charged particle production at forward rapidities

The number of particles (both charged particles and photons) produced per
participant pair as a function of $\eta-y_{beam}$, where $y_{beam}$ is the
beam rapidity, is observed to be independent of beam energy at forward
rapidities. This phenomenon is known as limiting fragmentation [14], which is
the longitudinal scaling of particle production at forward rapidties. The
centrality independent limiting fragmentation behavior is observed to be
followed by the inclusive photon production (primarily from decay of $\pi^{0}$
) at $\sqrt{\it{s}_{\rm NN}}$ = 62.4 GeV. The centrality dependent behavior of
limiting fragmentation is also observed to be followed by the inclusive
charged particles at 19.6 GeV and 200 GeV. The centrality dependent limiting
fragmentation for inclusive charged particles may be due to the baryons coming
from nuclear remnants and baryon transport [15], both of which change with
centrality.

The ratios of the number of charged particles to photons in the pseudorapidity
range $-3.7<\eta<-2.3$ are found to be $1.4\pm 0.1$ and $1.2\pm 0.1$ for
$\sqrt{\it{s}_{\rm NN}}$ = 62.4 GeV and 200 GeV, respectively in Cu+Cu and
Au+Au collisions [16]. Theoretically the ratio should be one in any system at
any energy if we consider only decay photons as majority of the particles
produced in any collision are pions and neutral pions decay into two photons.
The ratio here is varying as a function of collision energy but not with the
system size. It has been conjuctured that baryons contribute to charged
particles [17, 18], while mesons are the dominant contributors to photon
production in a given system at given energy so the ratio is more than one.
Also the contribution of baryons to charged particles is suppressed as we go
to the higher energies so the ratio is approaching to unity.

In this paper we will try to convince that the energy dependence of the
charged particle to photon ratio may be due to the varying contribution of the
direct photons keeping the baryon contribution fixed based on the theoretical
model predictions and experimental results.

The paper is organized in the following way:
Section 1 is the introduction of the paper and section 2 discusses the photon
and charged particle production at forward rapidities. Section 3 contains the
results discussing the evidences supporting our argument that the energy
dependence of the charged particle to photon ratio may be due to the varying
contribution of the direct photons keeping the baryon contribution fixed. In
section 4 we have summarized the content of the paper.

##  3 Results supporting our argument

Decay photons are those photons which are produced by the decay of hadrons
while direct photons are those which are produced not by the decay of hadrons.
Direct photons are divided into two catogories:
(1) Prompt direct photons: prompt direct photons are produced by the initial
hard scattering of partons in a collision such as Compton scattering of
partons [19].
(2) Thermal direct photons: thermal direct photons are emitted by the hot
Quark Gluon Plasma in the form of radiation or by quark-antiquark annihilation
[20].

The results supporting our argument that the energy dependence of the charged
particle to photon ratio may be due to the varying contribution of the direct
photons keeping the baryon contribution fixed are the following:

###  3.1 Increase in the direct photon cross section with the center of mass
energy

It has been observed that the direct photon cross section increases with
increase in the center of mass energy. The experimental result which is
available to us is direct photon production in S+Au collisions at
$\sqrt{\it{s}_{\rm NN}}$ = 19 GeV at WA80 experiment taking S+Au same as Au+Au
and Au+Au collisions at $\sqrt{\it{s}_{\rm NN}}$ = 200 GeV measured by PHENIX
experiment at RHIC.

Figure 1: Upper limits on direct photon production in central S+Au collisions
at $\sqrt{\it{s}_{\rm NN}}$ = 19 GeV as measured by the WA80 experiment [21].
The bars at the foot of each arrow mark the 90% C.L. upper limit on the
differential yield per event (the lengths of the arrows have no significance).
Figure 2: Direct photon invariant yields as a function of transverse momentum
for different centrality selections and minimum bias Au+Au collisions at
$\sqrt{\it{s}_{\rm NN}}$ = 200 GeV [22]. The vertical error bar on each point
indicates the total error. Arrows indicate measurements consistent with zero
yield with the tail of the arrow indicating the 90% confidence level upper
limit. The solid curves are pQCD predictions described in the text. Figure 3:
The ratio of number of charged particles to the number of photons as a
function of average number of participating nucleons $<N_{part}>$ for Au+Au
and Cu+Cu collisions at $\sqrt{\it{s}_{\rm NN}}$ = 200 GeV and 62.4 GeV for
$-3.7<\eta<-2.3$.

The comparison between the two results in Fig 1 and Fig 2 shows that both
multiplicity and energy of the produced direct photons increase with the
center of mass energy.

The ratios of the number of charged particles to photons in Fig 3 varying from
1.4 to 1.2 while going from $\sqrt{\it{s}_{\rm NN}}$ = 62.4 GeV to 200 GeV in
both Cu+Cu and Au+Au collisions may be due to the reason that direct photons
at $\sqrt{\it{s}_{\rm NN}}$ = 62.4 GeV in both Cu+Cu and Au+Au collisions are
produced in small quantity and with less energy so that they are unable to
reach the forward region $(-3.7<\eta<-2.3)$ where the Photon Multiplicity
Detector is placed, therefore only decay photons are contributing to the
photon multiplicity. On the other hand direct photons produced at
$\sqrt{\it{s}_{\rm NN}}$ = 200 GeV in both Cu+Cu and Au+Au collisions are
produced in large quantity and with high energy so that they are able to reach
Photon Multiplicity Detector and are contributing to the photon multiplicity
along with decay photons. This additional contribution of direct photons at
$\sqrt{\it{s}_{\rm NN}}$ = 200 GeV in both Cu+Cu and Au+Au collisions may be
the reason of improving the ratio $\frac{N_{ch}}{N_{\gamma}}$ to 1.2 $\pm$ 0.1
from 1.4 $\pm$ 0.1 keeping the energy to charged particle production fraction
same.

###  3.2 Existence and enhancement of direct photons in the forward region

Direct photons are radiated by the electric current of the projectile quarks,
which mostly stay in the fragmentation region of the beam, and tend to form a
peak at forward rapidities as shown in Fig 4. But at very large $p_{T}$ and
$\eta$ the kinematic limit pushes photon radiation to more central rapidities
and the peak at forward rapidities will be replaced by a kind of plateau at
central rapidities. Also gluons are radiated via nonabelian mechanisms by the
color current across the whole rapidity interval and tend to form a plateau at
midrapidity. The result shown below in Fig 5 is obtained from the Golec-
Biernat and W$\ddot{u}$sthoff (GBW) model which is a popular parametrization
for the q$\bar{q}$ dipole cross-section on a nucleon target [23, 24]. The Fig
5

Figure 4: Invariant cross-section for direct photon production in pp
collisions at LHC as a function of rapidity $\eta$ calculated with the Golec-
Biernat and W$\ddot{u}$sthoff (GBW) color dipole model for various fixed
$p_{T}$ [25]. Figure 5: Direct photon spectra obtained from the GBW dipole
model at the RHIC and CDF energies for pp collisions [25]. Experimental data
are from the PHENIX experiment [26] at $\eta$ = 0, and from the CDF experiment
[27, 28] at $\eta<0.9$. The error bars are the linear sum of the statistical
and systematic uncertainties.

shows the direct photon yield predicted by the (GBW) color dipole model is
almost in good agreement with the RHIC data in pp collisions at
$\sqrt{\it{s}}$ = 200 GeV. The experimental result shown in Fig 6 is from
ATLAS experiment showing the enhancement of direct photon yield in the more
forward region $1.52\leq|\eta|<2.37$ in the centrality class 10-20% in Pb+Pb
collisions at $\sqrt{\it{s}_{\rm NN}}$ = 2.76 TeV.

Figure 6: The normalized yields of prompt photons in Pb+Pb collisions at
$\sqrt{\it{s}_{\rm NN}}$ = 2.76 TeV divided by JETPHOX predictions for p+p
collisions as a function of $p_{T}$ in central (top panels, $|\eta|<1.37$) and
forward (bottom panels, $1.52<|\eta|<2.37$) rapidities in different centrality
selections from peripheral (left panels) to central (right panels). The
combined scale and PDF uncertainty on the JETPHOX calculation is shown by the
yellow area. Two other scenarios of JETPHOX calculation divided to the same
JETPHOX p+p is shown, including the isospin effect from Pb+Pb collisions
(black curve), and nuclear modification of the PDFs (blue area) [29].

The theoretical prediction and experimental results showing that direct
photons not only exist in the forward region but are enhanced in the more
forward region supports our argument that the additional contribution of
direct photons at $\sqrt{\it{s}_{\rm NN}}$ = 200 GeV in both Cu+Cu and Au+Au
collisions.

###  3.3 Estimation of direct photon contribution

It was observed in WA80 experiment that for O+Au collisions direct photons
constituted no more than 15% of the inclusive photon yield at each $p_{T}$ in
the range $0.4<p_{T}<2.4$ GeV$/c$.

Figure 7: Ratio of inclusive photon cross section to $\pi^{0}$ cross section
as a function of transverse momentum for 200 GeV p+C, 200 GeV p+Au , 200
A$\cdot$GeV $O^{16}$+C, and 60 A$\cdot$GeV $O^{16}$+Au events under minimum-
bias trigger conditions. The squares indicate the data, and the shaded
histograms are the Monte Carlo results representing hadronic decays. [31]

Fig 7 shows the collision systems 200 GeV p+C, 200 GeV p+Au, 200 A$\cdot$GeV
$O^{16}$+C, and 60 A$\cdot$GeV $O^{16}$+Au under minimum-bias trigger
conditions. Experimental data in the Fig 7 is indicated by the squares and
Monte Carlo results are indicated by the shaded histograms based on measured
$\pi^{0}$ spectra and represent the hadronic background. The direct photon
contribution can be attributed to the difference between data and the
histograms. The $p_{T}$ dependence of the data is reproduced by the Monte
Carlo simulations for all four systems. The data for the $O^{16}$ induced
reactions are approximately 10% above the prediction, but this finding is
still within the systematic uncertainty [30] of the data. After taking all the
systematic uncertainties into consideration the contribution of direct photons
is approximately about 15% of the neutral meson production [31].

Based on this observation we assume that at $\sqrt{\it{s}_{\rm NN}}$ = 62.4
GeV the direct photons do not reach the forward rapidities where Photon
Multiplicity Detector is situated, thus keeping the ratio of
$\frac{N_{ch}}{N_{\gamma}}$ high. But at $\sqrt{\it{s}_{\rm NN}}$ = 200 GeV if
15% additional contribution from direct photons is added to the yield of decay
photons, the experimentally observed ratio $\frac{N_{ch}}{N_{\gamma}}$ = 1.2
for both Cu+Cu and Au+Au could be explained.

##  4 Summary

In this work we report an alternative explanation of the variation of charged
particle to photon ratio in Cu-Cu collisions and Au-Au collisions from 1.4
$\pm$ 0.1 to 1.2 $\pm$ 0.1 at $\sqrt{\it{s}_{\rm NN}}$ = 62.4 GeV and 200 GeV
respectively in STAR experiment at RHIC. The experimental results such as
increase in the direct photon cross section with the center of mass energy,
existence and enhancement of direct photons in the forward region and 15%
contribution of direct photons in the total neutral meson production validates
our claim of additional direct photon contribution. The theoretical
predictions from the Golec-Biernat and W$\ddot{u}$sthoff (GBW) model which
predicted the direct photon yield in good agreement with the RHIC data in pp
collisions at $\sqrt{\it{s}}$ = 200 GeV, predicts the existence and
enhancement of direct photons in the forward region as discussed in the
section 3. Based on these experimental results and theoretical predictions we
argue that the additional contribution of direct photons at $\sqrt{\it{s}_{\rm
NN}}$ = 200 GeV in both Cu+Cu and Au+Au collisions are improving the ratio
$\frac{N_{ch}}{N_{\gamma}}$ to 1.2 $\pm$ 0.1 from 1.4 $\pm$ 0.1 keeping the
energy to charged particle production fraction same.

##  5 Acknowledgement

This work is supported by the Department of Science and Technology (DST),
Govt. of India (DST/INSPIRE Fellowship/2018/IF180058). We would like to thank
high energy physics group of Department of Physics and Centre for
Astroparticle Physics and Space Science, Bose Institute for fruitful
discussions.

## References

  * [1] W. Ke and X. N. Wang,“QGP modification to single inclusive jets in a calibrated transport model,” JHEP 05, 041 (2021)
doi:10.1007/JHEP05(2021)041 [arXiv:2010.13680 [hep-ph]].

  * [2] K. J. Gonçalves, A. V. Giannini, D. D. Chinellato and G. Torrieri, “Limiting fragmentation as an initial state probe in heavy ion collisions,” Phys. Rev. C 100, no.5, 054901 (2019)
doi:10.1103/PhysRevC.100.054901 [arXiv:1906.00947 [nucl-th]].

  * [3] M. He and R. Rapp, “Hadronization and Charm-Hadron Ratios in Heavy-Ion Collisions,” Phys. Rev. Lett. 124, no.4, 042301 (2020)
doi:10.1103/PhysRevLett.124.042301 [arXiv:1905.09216].

  * [4] S. Cao, G. Y. Qin and S. A. Bass, “Energy loss, hadronization and hadronic interactions of heavy flavors in relativistic heavy-ion collisions,” Phys. Rev. C 92, no.2, 024907 (2015)
doi:10.1103/PhysRevC.92.024907 [arXiv:1505.01413 [nucl-th]].

  * [5] J. Casalderrey-Solana, G. Milhano, D. Pablos and K. Rajagopal, “Modification of Jet Substructure in Heavy Ion Collisions as a Probe of the Resolution Length of Quark-Gluon Plasma,” JHEP 01, 044 (2020)
doi:10.1007/JHEP01(2020)044 [arXiv:1907.11248 [hep-ph]].

  * [6] C. Aidala et al. [PHENIX], “Creation of quark–gluon plasma droplets with three distinct geometries,” Nature Phys. 15, no.3, 214-220 (2019)
doi:10.1038/s41567-018-0360-0 [arXiv:1805.02973 [nucl-ex]].

  * [7] J. D. Bjorken, “Highly Relativistic Nucleus-Nucleus Collisions: The Central Rapidity Region,” Phys. Rev. D 27, 140-151 (1983)
doi:10.1103/PhysRevD.27.140

  * [8] J. L. Klay et al. [E895], “Longitudinal flow from 2-A-GeV to 8-A-GeV Au+Au collisions at the Brookhaven AGS,” Phys. Rev. Lett. 88, 102301 (2002)
doi:10.1103/PhysRevLett.88.102301 [arXiv:nucl-ex/0111006 [nucl-ex]].

  * [9] L. V. Gribov, E. M. Levin and M. G. Ryskin, “Semihard Processes in QCD,” Phys. Rept. 100, 1-150 (1983)
doi:10.1016/0370-1573(83)90022-4

  * [10] L. D. McLerran and R. Venugopalan, “Computing quark and gluon distribution functions for very large nuclei,” Phys. Rev. D 49, 2233-2241 (1994)
doi:10.1103/PhysRevD.49.2233 [arXiv:hep-ph/9309289 [hep-ph]].

  * [11] T. Abbott, L. Ahle, Y. Akiba, D. Alburger, D. Beavis, L. Birstein, M. A. Bloomer, P. D. Bond, H. C. Britt and B. Budick, et al. “Further observations on midrapidity et distributions with aperture corrected scale,” Phys. Rev. C 68, 034908 (2003)
doi:10.1103/PhysRevC.68.034908

  * [12] B. B. Back et al. [PHOBOS], “Charged particle pseudorapidity density distributions from Au+Au collisions at $\sqrt{s_{NN}}$ = 130-GeV,” Phys. Rev. Lett. 87, 102303 (2001)
doi:10.1103/PhysRevLett.87.102303 [arXiv:nucl-ex/0106006 [nucl-ex]].

  * [13] D. Kharzeev, E. Levin and M. Nardi, “The Onset of classical QCD dynamics in relativistic heavy ion collisions,” Phys. Rev. C 71, 054903 (2005)
doi:10.1103/PhysRevC.71.054903 [arXiv:hep-ph/0111315 [hep-ph]].

  * [14] R. Beckmann, S. Raha, N. Stelte and R. M. Weiner, “Limiting Fragmentation in High-energy Heavy Ion Reactions and Preequilibrium,” Phys. Lett. B 105, 411-416 (1981)
doi:10.1016/0370-2693(81)91194-1

  * [15] J. Adams et al. [STAR], “Multiplicity and pseudorapidity distributions of charged particles and photons at forward pseudorapidity in Au + Au collisions at s(NN)**(1/2) = 62.4-GeV,” Phys. Rev. C 73, 034906 (2006)
doi:10.1103/PhysRevC.73.034906 [arXiv:nucl-ex/0511026].

  * [16] B. I. Abelev et al. [STAR], “Center of mass energy and system-size dependence of photon production at forward rapidity at RHIC,” Nucl. Phys. A 832, 134-147 (2010)
doi:10.1016/j.nuclphysa.2009.11.011 [arXiv:0906.2260 [nucl-ex]].

  * [17] I. G. Bearden et al. [BRAHMS], “Nuclear stopping in Au + Au collisions at s(NN)**(1/2) = 200-GeV,” Phys. Rev. Lett. 93, 102301 (2004)
doi:10.1103/PhysRevLett.93.102301 [arXiv:nucl-ex/0312023 [nucl-ex]].

  * [18] D. Kharzeev, “Can gluons trace baryon number?,” Phys. Lett. B 378, 238-246 (1996)
doi:10.1016/0370-2693(96)00435-2 [arXiv:nucl-th/9602027].

  * [19] J. F. Owens, “Large Momentum Transfer Production of Direct Photons, Jets, and Particles,” Rev. Mod. Phys. 59, 465 (1987)
doi:10.1103/RevModPhys.59.465

  * [20] O. Linnyk, V. P. Konchakovski, W. Cassing and E. L. Bratkovskaya, “Photon elliptic flow in relativistic heavy-ion collisions: hadronic versus partonic sources, ” Phys. Rev. C 88, 034904 (2013) doi:10.1103/PhysRevC.88.034904 [arXiv:1304.7030 [nucl-th]].
  * [21] R. Albrecht et al. [WA80], “Limits on the production of direct photons in 200-A/GeV S-32 + Au collisions,” Phys. Rev. Lett. 76, 3506-3509 (1996)
doi:10.1103/PhysRevLett.76.3506

  * [22] S. S. Adler et al. [PHENIX], “Centrality dependence of direct photon production in s(NN)**(1/ 2) = 200-GeV Au + Au collisions,” Phys. Rev. Lett. 94, 232301 (2005)
doi:10.1103/PhysRevLett.94.232301 [arXiv:nucl-ex/0503003 [nucl-ex]].

  * [23] J. Jalilian-Marian, “Rapidity dependence of the photon to pion production ratio in high energy collisions,” Nucl. Phys. A 806, 305-311 (2008)
doi:10.1016/j.nuclphysa.2008.02.300 [arXiv:nucl-th/0703069 [nucl-th]].

  * [24] S. Jeon, J. Jalilian-Marian and I. Sarcevic, “Prompt photon and inclusive pi0 production at RHIC and LHC,” Nucl. Phys. A 715, 795 (2003)
doi:10.1016/S0375-9474(02)01491-4 [arXiv:nucl-th/0211084 [nucl-th]].

  * [25] A. H. Rezaeian and A. Schafer, “Hadrons and direct photon in pp and pA collisions at LHC and saturation effects,” Phys. Rev. D 81, 114032 (2010)
doi:10.1103/PhysRevD.81.114032 [arXiv:0908.3695 [hep-ph]].

  * [26] S. S. Adler et al. [PHENIX], “Measurement of direct photon production in p+p collisions at s**(1/2) = 200-GeV,” Phys. Rev. Lett. 98, 012002 (2007)
doi:10.1103/PhysRevLett.98.012002 [arXiv:hep-ex/0609031].

  * [27] D. Acosta et al. [CDF], “Direct photon cross section with conversions at CDF,” Phys. Rev. D 70, 074008 (2004)
doi:10.1103/PhysRevD.70.074008 [arXiv:hep-ex/0404022].

  * [28] F. Abe et al. [CDF], “A Precision measurement of the prompt photon cross-section in $p\bar{p}$ collisions at $\sqrt{s}=1.8$ TeV,” Phys. Rev. Lett. 73, 2662-2666 (1994) [erratum: Phys. Rev. Lett. 74, 1891-1893 (1995)]
doi:10.1103/PhysRevLett.73.2662.

  * [29] G. Aad et al. [ATLAS], “Centrality, rapidity and transverse momentum dependence of isolated prompt photon production in lead-lead collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV measured with the ATLAS detector,” Phys. Rev. C 93, no.3, 034914 (2016)
doi:10.1103/PhysRevC.93.034914 [arXiv:1506.08552].

  * [30] T. Åkesson et al. [Helios], “Inclusive Negative Particle $p(T$) Spectra in $P$ \- Nucleus and Nucleus-nucleus Collisions at 200-GeV Per Nucleon,” Z. Phys. C 46, 361-368 (1990)
doi:10.1007/BF01621023

  * [31] R. Albrecht et al. [WA80], “Upper limit for thermal direct photon production in heavy ion collisions at 60-A/GeV and 200-A/GeV,” Z. Phys. C 51, 1-10 (1991)
doi:10.1007/BF01579555