Data-driven discovery of multiscale chemical reactions

governed by the law of mass action

Juntao Huang [Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA. E-mail<EMAIL_ADDRESS> Yizhou Zhou [School of Mathematical Sciences, Peking University, Beijing, China. E-mail<EMAIL_ADDRESS>Wen-An Yong [Department of Mathematical Sciences, Tsinghua University, Beijing, China. E-mail<EMAIL_ADDRESS>
In this paper, we propose a data-driven method to discover multiscale chemical reactions governed by the law of mass action.
First, we use a single matrix to represent the stoichiometric coefficients for both the reactants and products in a system without catalysis reactions.
The negative entries in the matrix denote the stoichiometric coefficients for the reactants and the positive ones for the products.
Second, we find that the conventional optimization methods usually get stuck in the local minima and could not find the true solution in learning the multiscale chemical reactions. To overcome this difficulty, we propose a partial-parameters-freezing (PPF) technique to progressively determine the network parameters by using the fact that the stoichiometric coefficients are integers. With such a technique, the dimension of the searching space is gradually reduced in the training process and the global mimina can be eventually obtained. Several numerical experiments including the classical Michaelis–Menten kinetics, the hydrogen oxidation reactions and the simplified GRI-3.0 mechanism verify the good performance of our algorithm in learning the multiscale chemical reactions. The code is available at <https://github.com/JuntaoHuang/multiscale-chemical-reaction>.

Key Words:
Chemical Reactions; Multiscale; Machine Learning; Nonlinear Regression; Ordinary Differential Equations

§ INTRODUCTION

Chemical reactions are fundamental in many scientific fields including biology, material science, chemical engineering and so on.
To identity the reactions from experimental data, the traditional methods are mainly based on some empirical laws and expert knowledge [10]. Recently, thanks to the rapid development of machine learning [19] and data-driven modeling [34, 22, 6, 3, 29, 20, 30, 14], it is desirable to develop a data-driven method of discovering the underlying chemical reactions from massive data automatically.

Consider a reaction system with $n_s$ species participating in $n_r$ reactions:

\nu_{i1}'\mathcal{S}_1 + \nu_{i2}'\mathcal{S}_2 + \cdots + \nu_{in_s}'\mathcal{S}_{n_s} \ch{<=>[$k\sb{if}$][$k\sb{ir}$]}
\nu_{i1}''\mathcal{S}_1 + \nu_{i2}''\mathcal{S}_2 + \cdots + \nu_{in_s}''\mathcal{S}_{n_s}

for $i = 1,2, \cdots, n_r$. Here $\mathcal{S}_k$ is the chemical symbol for the $k$-th species, the nonnegative integers $\nu_{ik}'$ and $\nu_{ik}''$ are the stoichiometric coefficients of the $k$-th species in the $i$-th reaction, and $k_{if}$ and $k_{ir}$ are the direct and reverse reaction rates of the $i$-th reaction. The reaction is reversible if both $k_{if}$ and $k_{if}$ are positive. Strictly speaking, all elementary chemical reactions are reversible due to microscopic reversibility. However, in real applications, some of the rate constants are  negligible, thus the corresponding reactions can be omitted and the retained ones can be considered as irreversible.

Denote by $u_k = u_k(t)$  the concentration of the $k$-th species at time $t$ for $k=1,2,\cdots,n_s$.
According to the law of mass action [38], the evolution of $u_k$ obeys the ordinary differential equations (ODEs) [27]

\begin{equation}\label{eq:reaction-ODE}
\frac{du_k}{dt} = \sum_{i=1}^{n_r}(\nu_{ik}''-\nu_{ik}')\left(k_{if}\prod_{j=1}^{n_s}u_j^{\nu_{ij}'}-k_{ir}\prod_{j=1}^{n_s}u_j^{\nu_{ij}''} \right),
\end{equation}

for $k=1,2,\cdots,n_s$.
Given the concentration time series data $\{u_k(t_{n}), \ k=1,\cdots,n_s, \ n=1,\cdots,N \}$, our goal is to learn the stoichiometric coefficients $\nu_{ik}'$, $\nu_{ik}''$ and reaction rates $k_{if}$ and $k_{ir}$.

In the literature there are already some works on this topic.
In [5], the authors applied linear regressions to infer the chemical reactions, with the assumption that the reactions are at most the result of bimolecular collisions and the total reaction order is not greater than two. In [40], the linear regression was utilized with an L1 objective, which transforms the problem into a mixed-integer linear programming (MILP). This approach suffers from the same restrictive assumptions as in [5]. In [18], the authors presented an approach to infer the stoichiometric subspace of a chemical reaction network from steady-state concentration data profiles, which is then cast as a series of MILP.
In [25], some chemically reasonable requirements were considered such as the mass conservation and the principle of detailed balance.
The deep neural networks (DNNs) were applied to extract the chemical reaction rate information in [31, 32], but the weights are difficult to interpret physically.
In [13], the authors adapted the sparse identification of nonlinear dynamics (SINDy) method [4, 9] to the present problem. However, the approach relies on expert knowledge, which precludes the application in a new reaction system with unknown reaction pathways.
Within the framework of SINDy, other works are [2, 1, 24]. In order to improve the performance of SINDy, two additional steps including least-squares regression and stepwise regression in the identification were introduced in [2], which are based on the traditional statistical methods. In [1], SINDy was combined with the DNNs to adaptively model and control the process dynamics. An implicit-SINDy was proposed and applied to infer the Michaelis-Menten enzyme kinetics in [24].
Additionally, a statistical learning framework was proposed based on group-sparse regression which leverage prior knowledge from physical principles in [23]. For example, the mass conservation is enforced in the JAK-STAT reaction pathway for signal transduction in [23].

Our work is mainly motivated by [15], where the authors proposed a Chemical Reaction Neural Network (CRNN) by resorting to the feature of the equations in (<ref>). The discovery of chemical reactions usually involves two steps: the identification of the reaction pathways (i.e., the stoichiometric coefficients) and the determination of the reaction rates. For complex reaction processes, one could not even identify the reaction pathways and has to infer both the stoichiometric coefficients and the rate constants from data.
The work in [15] presents a neural network approach for discovering unknown reaction pathways from concentration data. The parameters in CRNN correspond to the stoichiometric coefficients and reaction rates and the network has only one hidden layer with the exponential activation functions.

Different from CRNN in [15], we use a single matrix of order $n_r\times n_s$ to represent the stoichiometric coefficients for both the forward and reverse reactions by assuming no catalysis reactions. The negative entries in the matrix denote the stoichiometric coefficients for the reactants and the positive for the products.

On the other hand, the reaction rates often differ in a wide range of magnitudes, which causes a lot of troubles in learning the multiscale chemical reactions.
To provide some insights into this difficulty, we design a nonlinear regression problem to fit a polynomial with two terms, see (<ref>) in Section <ref>. The given coefficients of the polynomial differ in several orders of magnitudes and the polynomial degree is to be determined. We find numerically that the conventional optimization algorithm usually gets stuck in the local minima and could not find the true solution. Another observation in the numerical experiment is that the learned polynomial degree of the terms with larger coefficient is close to the true solution. Inspired by this observation, we propose a partial-parameters-freezing (PPF) technique to escape from the local minima. Specifically, we perform a round operation on the learned polynomial degree which are close to integer in the optimization process if the loss function does not decrease. The revised algorithm works well for this problem. Some theoretical analysis is also provided to explain the numerical phenomenon.

We then generalize the PPF technique to learn the multiscale chemical reactions. Notice that the stoichiometric coefficients are integers. In the training process, if the loss function stops to decrease, the stoichiometric coefficients which are close to integers are rounded and then frozen afterwards. With such a treatment, the stoichiometric coefficients are gradually determined, the dimension of the searching space is  reduced in the training process, and eventually the global mimina can be obtained. Several numerical experiments including the classical Michaelis–Menten kinetics, the hydrogen oxidation reactions and the simplified GRI-3.0 mechanism verify that our method performs much better in learning the mutiscale chemical reactions.

This paper is organized as follows. In Section <ref>, we investigate a multiscale nonlinear regression problem numerically and theoretically. Our algorithm for learning the multiscale chemical reactions is presented in Section <ref>. In Section <ref>, the performance of the algorithm is validated through several numerical examples. Finally, conclusions and the outlook of future work are presented in Section <ref>.

§ MULTISCALE NONLINEAR REGRESSION PROBLEM

To provide some insights into the difficulties in learning the multiscale chemical reactions, we consider a nonlinear regression problem to fit the following function:

\begin{equation}\label{eq:regression-function}
y = f(x;\theta_1, \theta_2) = c_1 x^{\theta_1} + c_2 x^{\theta_2}.
\end{equation}

Here $c_1$ and $c_2$ are two given constants satisfying $\abs{c_1}\ll\abs{c_2}$, and $\theta_1, \theta_2$ are two integers to be determined. This simple toy model captures two key features of the multiscale chemical reactions. The first feature is that the right-hand side of the chemical reaction ODEs (<ref>) is polynomials and the stoichiometric coefficients are integers. The second one is that the multiscale chemical reactions often have reaction rates which differ in several orders of magnitudes.

Given the dataset $\{(x_i, y_i): ~ i=1,\cdots,N \}$, we define the loss function to be the mean squared error (MSE):

\begin{equation}\label{eq:regression-loss}
\mathcal{L}(\theta_1, \theta_2) = \frac{1}{N}\sum_{i=1}^N(f(x_i; \theta_1, \theta_2) - y_i)^2,
\end{equation}

to estimate the parameters $\theta_1$ and $\theta_2$.
Next, conventional optimization methods can be used to obtain the estimation of $\theta_1$ and $\theta_2$.

In the numerical experiment, we take $c_1=1$ and $c_2=100$. The ground truth solutions are $\theta_1=1$ and $\theta_2=2$. The data $x_i$ for $i=1,\cdots,N$ are randomly sampled from a uniform distribution in $(0, 1)$ with the number of data $N=1000$, and $y_i=c_1x_i+ c_2x_i^2$. The Adam optimization method [17] is applied with the full batch gradient decent. The learning rate is taken to be $10^{-4}$. The initial guess of $\theta_1$ and $\theta_2$ is randomly chosen in $(-1, 1)$.

For this toy model, we numerically find that the naive implementation will get stuck in the local minima $(\theta_1, \theta_2) = (3.8286, 1.9745)$ and could not find the true solution. The history of the loss function and the parameters $\theta_1$ and $\theta_2$ in the training is presented in Figure <ref>, see the dashed lines.

Although the naive optimization could not find the global minima, we notice that $\theta_2=1.9745$ in this local minima is close to the true solution $\theta_2=2$. Inspired by this observation, we propose a partial-parameters-freezing (PPF) technique to escape from the local minima. To be more specific, we keep track of the loss function in the training. If the loss does not decrease, we check the parameters $\theta_1$ and $\theta_2$: if any of these is close to its nearest integer with a given threshold, we round it to the integer and do not update it in the afterwards optimization process.

For comparison, we also plot the history of the loss function and the parameters with the PPF technique in Figure <ref>, see the solid lines. The threshold is taken to be 0.05 in this test. The loss stops decreasing with the epoch around 7000. Then $\theta_2$ is rounded to 2 and only $\theta_1$ is updated afterwards. The true solution is eventually obtained when the epoch is around 10000.

[loss vs. epoch]

[parameters $\theta_1$ and $\theta_2$ vs. epoch]

Multiscale nonlinear regression problem: the history of loss function in (<ref>) and the parameters $\theta_1$ and $\theta_2$ in the training process. Solid lines: the method with the PPF technique; dashed lines: the method without the PPF technique.

To better understand why it is easy to get stuck in the local minima without the PPF treatment, we investigate the landscape of the loss function. In Figure <ref>, we plot the 3D surface and the contour map for the loss as a function of $(\theta_1,\theta_2)$. In Figure <ref> (a), it is observed that the loss function has several local minima in which $\theta_2$ is close to 2. Moreover, the local minima $(\theta_1, \theta_2) = (3.8286, 1.9745)$ in the naive implementation is also labeled in Figure <ref> (b).

[loss function surface plot]

[loss function contour map]

Multiscale nonlinear regression problem: the landscape of the loss function in (<ref>). Left: loss function surface plot (in log scale); right: loss function contour map (in log scale), local minima $(\theta_1, \theta_2) = (3.8286, 1.9745)$.

We also plot the profiles of the loss function with fixed $\theta_2=1.99$, 2 and 2.01 in Figure <ref>. It is observed that slight perturbations in $\theta_2$ have a considerable impact on the minima of the loss function. Moreover, the loss as a 1D function with fixed $\theta_2=2$ is well-behaved. This explains why our algorithm is easy to find the global minima after freezing the integer parameter $\theta_2$.

Multiscale nonlinear regression problem: loss function in (<ref>) with fixed parameters $\theta_2=1.99$, 2 and 2.01.

We mention that we also test other cases with different coefficients $c_1$ and $c_2$ satisfying $\abs{{c_2}/{c_1}}=10^3, 10^4, 10^5$ and different integers $\theta_1$ and $\theta_2$. The results are similar and thus omitted here.

We conclude this section with some theoretical analysis to  explain the local minima phenomenon observed above. By taking gradient of the loss function in (<ref>), we have

\begin{equation}\label{eq:regression-grad-loss}
\frac{\partial\mathcal{L}}{\partial\theta_j} =  \frac{2c_jc_2 }{N}\sum_{i=1}^N  \brac{\frac{c_1}{c_2} (x_i^{\theta_1} - x_i^{\theta_1^{\textrm{e}}}) +  (x_i^{\theta_2} - x_i^{\theta_2^{\textrm{e}}})} x_i^{\theta_j} \ln x_i, \quad j=1,2.
\end{equation}

Here $\theta_i^{\textrm{e}}$ denotes the true solution of the parameter $\theta_i$ for $i=1,2$.
From the expression (<ref>), we can provide some insights on the phenomenon that the local minina $\theta_2$ is close to the true solution $\theta_2^{\textrm{e}}$.
To reach the local minima, the gradient should be zero.
Refer to the expression (<ref>).
Since $\abs{c_1/c_2}\ll 1$, whether or not the gradient is close to zero depends mainly on the fact that $\theta_2$, instead of $\theta_1$, is close to the ground truth.

§ ALGORITHM

In this section, we present our algorithm for learning the multiscale chemical reactions. First, we use a single matrix to represent the stoichiometric coefficients for both the reactants and products. Each row of the matrix represents one reaction, where the negative entries denote the stoichiometric coefficients for the reactants and the positive ones for the products. This setup is valid for systems without catalysis reactions. In addition, we adapt the PPF technique for the multiscale nonlinear regression problem proposed in Section <ref> to learn the multiscale chemical reactions (<ref>).

We assume that the data are given in the form of the concentrations and the time derivatives in different time snapshots $\{(u_k(t_{n}), u_k'(t_{n})), \ k=1,\cdots,n_s, \ n=1,\cdots,N \}$, our goal is to learn the stoichiometric coefficients and the reaction rates. Realistically, often only $u_k(t_{n})$ is available, and the time derivatives $u_k'(t_{n})$ could be approximated using numerical differentiations [33, 7].

To better illustrate the algorithm, we firstly introduce some vector notations. We denote the forward and reverse reaction rates in (<ref>) by $\bm{k}_f = (k_{1f},k_{2f},...,k_{n_rf})$ and $\bm{k}_r = (k_{1r},k_{2r},...,k_{n_rr})$. The stoichiometric coefficients in (<ref>) are collected in two matrices:

\begin{equation}
\bm{V}' =
\begin{pmatrix}
\nu_{11}' & \nu_{12}' & \cdots & \nu_{1n_s}'\\
\nu_{21}' & \nu_{22}' & \cdots & \nu_{2n_s}'\\
\vdots & \vdots & \ddots & \vdots\\
\nu_{n_r1}' & \nu_{n_r2}' & \cdots & \nu_{n_rn_s}'
\end{pmatrix}
, \qquad
\bm{V}'' =
\begin{pmatrix}
\nu_{11}'' & \nu_{12}'' & \cdots & \nu_{1n_s}''\\
\nu_{21}'' & \nu_{22}'' & \cdots & \nu_{2n_s}''\\
\vdots & \vdots & \ddots & \vdots\\
\nu_{n_r1}'' & \nu_{n_r2}'' & \cdots & \nu_{n_rn_s}''
\end{pmatrix}.
\end{equation}

Assume that there is no catalysis reactions. Therefore, only one of $\nu_{ik}'$ and $\nu_{ik}''$ can be non-zero for any $(i,k)$.
In this case, the matrix $\bm{V}=(\nu_{ik}):=\bm{V}''-\bm{V}'$ satisfies

\nu_{ik}=\left\{ {\begin{array}{*{20}c}
\vspace{1.5mm} \nu_{ik}'',\qquad & \textrm{if} \quad \nu_{ik}\ge0,\\[2mm]
-\nu_{ik}',\qquad & \textrm{if} \quad \nu_{ik}<0.
\end{array}} \right.

According to this property, we only need to pin down the matrix $\bm{V}$. Then $\bm{V}'$ and $\bm{V}''$ can be recovered by $\nu_{ik}''=\max(0,\nu_{ik})$ and $\nu_{ik}'=-\min(0,\nu_{ik})$, respectively.

Next, we define the neural network $\mathcal{N}=\mathcal{N}(u_1,\cdots,u_{n_s}):\mathbb{R}^{n_s}\rightarrow\mathbb{R}^{n_s}$ which has the input $\bm{u}:=(u_1,\cdots,u_{n_s})$ and the parameters $\bm{l}_f = (l_{1f}, l_{2f},..., l_{n_rf})$, $\bm{l}_r = (l_{1r}, l_{2r},..., l_{n_rr})$ and $\bm{V}$:

\begin{equation*}
\mathcal{N}(u_1,\cdots,u_{n_s})_k = \sum_{i=1}^{n_r}\nu_{ik}\left(\exp({l_{if}})\prod_{j=1}^{n_s}u_j^{-\min(0,\nu_{ik})} - \exp({l_{ir}})\prod_{j=1}^{n_s}u_j^{\max(0,\nu_{ik})} \right)
\end{equation*}

for $k=1,\cdots,n_s$.
Here the parameters $l_{if}$ and $l_{ir}$ denote the logarithms of the reaction rates $k_{if}$ and $k_{ir}$ [15]. This change of variables technique has two advantages. The first one is that the positivity of the reaction rates is guaranteed automatically. The second one is that the reaction rates for the multiscale chemical reactions usually differ in several orders of magnitudes. The slight changes of $l_{if}$ and $l_{ir}$ will make $k_{if}$ and $k_{ir}$ change a lot, which could potentially make the neural network to be more robust in the training process.

The loss function is defined as the mean squared error (MSE) between the data for the time derivatives and the output of the neural network:

\begin{equation}\label{3.2}
\mathcal{L} =  \frac{1}{N}\sum_{n=1}^N\sum_{k=1}^{n_s}\brac{\mathcal{N}(u_1(t_{n}),\cdots,u_{n_s}(t_{n}))_k - u_k'(t_{n})}^2 + \lambda  \mathcal{L}_r.
\end{equation}

Here $\lambda  \mathcal{L}_r$ is a regularization term with $\lambda>0$ the regularization constant and

\begin{equation}
\mathcal{L}_r = \sum_{i=1}^{n_r}\sum_{k=1}^{n_s} \abs{\nu_{ik}} + \sum_{i=1}^{n_r}\sum_{k=1}^{n_s} \nu_{ik}^2 + \sum_{i=1}^{n_r}(\abs{l_{if}} + \abs{l_{ir}}) + \sum_{i=1}^{n_r}(l_{if}^2 + l_{ir}^2).
\end{equation}

Here both $L_1$ and $L_2$ regularization terms are included.

This neural network works quite well for non-stiff chemical reactions. However, for stiff reactions, we observe that the optimization usually gets stuck in the local minima in the training process and could not find the true solution. The common techniques such as the mini-batch and reducing the learning rate do not work in such a situation. To attack this problem, we adapt the PPF technique proposed in the previous section.

The training procedure is split into two parts. The first part is to learn the matrix $\bm{V}$. To better illustrate the algorithm, we introduce some notation. Denote the vector in the $j$-th row of $\bm{V}$ by $\bm{v}_j$ for $j=1,\cdots,n_r$. Define the distance to the nearest integer for any vector $\bm{v}\in\mathbb{R}^{n_s}$ as

\begin{equation}\label{3.4}
d_{\textrm{int}}(\bm{v}) := \norm{\bm{v} - \nint{\bm{v}}}_{\infty} = \max_{i\in\{1,\cdots,n_s\}} \abs{v_i - \nint{v_i}},
\end{equation}

where $\nint{}$ denotes the function rounding an arbitrary real number to its nearest integer and it is defined to work element-wise on vectors. We keep track of the loss function in the training process. If the loss function stops decreasing, we check if any row of $\bm{V}$ is close to the nearest integers, i.e., $d_{\textrm{int}}(\bm{v}_j)\le\epsilon$. Here, $\epsilon>0$ is a hyperparameter and we take $\epsilon=0.05$ in all the numerical examples in Section <ref>. If the $j$-th row of $\bm{V}$ satisfies the condition $d_{\textrm{int}}(\bm{v}_j)\le\epsilon$, then we round $\bm{v}_j$ to $\nint{\bm{v}_j}$ and do not update it in the afterwards training. In addition, to help the optimization algorithm escape from the local minima, we randomly reinitialize other non-integer entries in $\bm{V}$ when the loss stops decreasing. After all the entries in $\bm{V}$ reach integer, we freeze them and then learn the parameters $\bm{l}_f$ and $\bm{l}_r$ related to the reaction rates. We remark that the SINDy algorithms [4, 13] can also be applied in learning the reaction rates when the stoichiometric coefficients $\bm{V}$ are known. The algorithm is summarized in Algorithm <ref>.

Here we assume that all the reactions are reversible. However, the algorithm can be also applied to irreversible reactions without any modification. The expected result is that the learned reverse reaction rates for the irreversible reactions will be close to zero. This will be demonstrated numerically in Example <ref> in Section <ref>.

The number of reactions can be learned by repeatedly executing the algorithm with different $n_r$. The ground truth of $n_r$ can be inferred from the best one. This will be shown in the numerical examples in the next section.

In many chemical reaction systems, the rate constants usually depend on the temperature. For example, the Arrhenius law can describe such a dependence:

\begin{equation}
k = A \exp\brac{-\frac{E_a}{RT}},
\end{equation}

where $k$ is the reaction rate, $A$ is the pre-exponential factor, $E_a$ is the activation energy and $R$ is the gas constant. In this case, the unknown parameters will include the pre-exponential factor, the activation energy and the stoichiometric coefficients. Our PPF technique can be directly applied without much modification. The performance will be verified numerically in the test in the next section.

time series data $\{(u_k(t_{n}), u_k'(t_{n})), \ k=1,\cdots,n_s, \ n=1,\cdots,N \}$
stoichiometric coefficient matrix $\bm{V}$, chemical reaction rates $\bm{k}_f$ and $\bm{k}_r$

initialize hyperparameters: number of reactions $n_r$, total number of epoch $M$, learning rate $lr$, regularization coefficient $\lambda$, integer threshold $\epsilon$

initialize parameters: $\bm{V}$, $\bm{l}_f$ and $\bm{l}_r$

step 1: learning ${V}$
$L_{\textrm{rec}}$ = np.zeros($M$);  *[f]record loss function in each epoch

Compute loss $\mathcal{L}$

Compute $\frac{\partial \mathcal{L}}{\partial\theta}$ by backpropagation

Update parameters (excluding the integer entries in $\bm{V}$) by Adam method

if loss increase, then check if any row of ${V}$ is close to integer
$L_{rec}[i] \ge L_{rec}[i-1]$

$\bm{v}_j \leftarrow \nint{\bm{v}_j} $

j $\notin$ $S_{\textrm{int}}$

$\bm{v}_j \leftarrow \textrm{rand}(-2, 2)$;		*[f]random reinitialize non-integer entries in $V$

if all the entries in $V$ are integers, then stop learning $V$
$S_{\textrm{int}} = \{ 1,\cdots,n_r \}$

step 2: learning $k_f$ and $k_r$

Compute loss $\mathcal{L}$

Compute $\frac{\partial \mathcal{L}}{\partial\theta}$ by backpropagation

Update parameters $\theta$ (excluding $\bm{V}$) by Adam method

$k_{if} \leftarrow \exp({l_{if}})$
$k_{ir} \leftarrow \exp({l_{ir}})$

Algorithm for learning chemical reactions

§ NUMERICAL RESULTS

Here the performance of our algorithm will be shown with five examples. The first example is an artificial reaction mechanism with two reactions [21]. The second one is the well-known Michaelis-Menten kinetics [16] in biochemistry. The third one is the hydrogen oxidation reactions [11, 8]. The fourth one is the extended Zeldovich mechanism, a typical chemical mechanism describing the oxidation of nitrogen and NOx formation [41]. The last one is the simplified GRI-3.0 mechanism, a chemical mechanism describing the methane oxidation [37].

In each numerical example, we randomly take 100 different initial conditions to generate the data. For each initial condition, we take uniform time snapshots at $t_n=n\Delta t$ with $n=0,\dots,10$ and $\Delta t=0.1$. The data is generated by solving the governing ODEs numerically using implicit Runge-Kutta method of the Radau IIA family of the fifth order [39] with small enough tolerance. The datasets are randomly split into the training datasets and the validation datasets by a ratio of 4:1. It is worthy to note that here we do not take $\Delta t$ to be too small so that the datasets could be potentially replaced by the experiment data in the future. The algorithm is implemented with PyTorch [28].

Now we present some details of the training and validation for the following four numerical tests. In the training process, all the parameters in the neural network are first randomly initialized from the uniform distribution in the interval $(-0.5, 0.5)$. Then, we update the parameters by minimizing the loss in (<ref>) using the standard Adams algorithm [17]. The learning rate is taken to be $10^{-3}$ and the regularization coefficient $\lambda$ in (<ref>) is $10^{-8}$. Recall the training method following (<ref>), we take the integer threshold to be 0.05. Besides, the total epoch number is $10^6$ and the mini-batch gradient descent is applied with the batch size 10.
For the validation, we use the following relative $L^2$ error:

\begin{equation*}
E =  \sqrt{\frac{\sum_{n=1}^N\sum_{k=1}^{n_s}\abs{\mathcal{N}(u_1(t_{n}),\cdots,u_{n_s}(t_{n}))_k - u_k'(t_{n})}^2}{\sum_{n=1}^N\sum_{k=1}^{n_s}\abs{ u_k'(t_{n})}^2}}.
\end{equation*}

Here the $(u_k(t_{n}), u_k'(t_{n}))$'s
come from the validation dataset.
For the other details, we refer the interested readers to our code in <https://github.com/JuntaoHuang/multiscale-chemical-reaction>.

[hypothetical stiff reaction network]
The first test case is an artificial reaction network with two reactions, taken from [21]:

\begin{align}
\ch{F <=>[$k\sb{1}\sp{+}$][$k\sb{1}\sp{-}$] R} \label{eq:artificial-reaction-1} \\
\ch{R <=>[$k\sb{2}\sp{+}$][$k\sb{2}\sp{-}$] P} \label{eq:artificial-reaction-2}
\end{align}

Here F, R and P indicate the fuel, radical and product in combustions, respectively. The reaction rates are taken to be $k_1^+ = k_1^- = 1$ and $k_2^+ = k_2^- = 10^3$. The two reversible reactions in (<ref>) have dramatically different reaction rates. Thus, the second reaction (<ref>) will quickly approach to equilibrium after a transient period, after which the first one (<ref>) becomes rate-limiting. This simple model is chosen to test the correctness of our code for stiff reactions.

The corresponding ODE system for (<ref>) is linear. The eigenvalues of the coefficient matrix are $\lambda_1=-2000$, $\lambda_2=-1.5$ and $\lambda_3=0$, which differ in several orders of magnitudes. This indicates that the ODE system is stiff [39].

To illustrate the advantage of the PPF technique, we compare the performance of the algorithm with and without this technique. The history of the training and validation errors is shown in Figure <ref>. The relative error stays around $10^{-3}$ without this technique, and decreases to $10^{-6}$ after applying this technique. The learned parameters are listed in Table <ref>. The upper part of the table is the learned parameters with the PPF technique, which agrees well with the ground truth in (<ref>). By contrast, the algorithm without imposing this technique could not generate the correct result. Moreover, it is interesting to see that, without using the technique, the learned stoichiometric coefficients in the first reaction and the second one has the opposite sign. We also notice that the summation of the forward rate $k_f$ of the first reaction and the reverse rate $k_r$ of the second one is close to the true reaction rate $10^3$. The same holds true for the reverse rate of the first reaction and the forward rate of the second one. This indicates that the effect of these two learned reactions is identical to the fast reaction (<ref>) and the slow reaction (<ref>) is not captured here. This is similar to the phenomenon we observed in the multiscale nonlinear regression problem in Section <ref>.

Example <ref>: the history of the relative error for the training data and the verification data. Solid line: the method with the PPF technique; dashed line: the method without the PPF technique.

freezing     $x_1$     $x_2$     $x_3$     $k_f$     $k_r$

$1$ 	    $0.000$      $1.000$     $-1.000$      1.000e+03     1.000e+03

$2$ 	    $-1.000$      $1.000$      $0.000$     1.000e+00     1.000e+00

no freezing     $x_1$     $x_2$     $x_3$     $k_f$     $k_r$

$1$ 	    $-0.001$      $0.999$     $-0.999$      7.448e+02     5.731e+02

$2$ 	    $0.000$      $-0.999$      $0.999$      4.277e+02     2.559e+02

Example <ref>: learned parameters. Upper part: with the PPF technique; lower part: without the PPF technique. Here $(x_1, x_2, x_3)$ denotes the row vector of the matrix $\bm{V}$.

Next, we test the algorithm with different number of chemical reactions. We take the number of reactions ranging from 1 to 4. The relative errors in the training data and the validation data are shown in Figure <ref>. The relative error decreases by three magnitudes when increasing the number of proposed reactions from one to two and reaches a plateau after that. Moreover, it is observed from Table <ref> that some of the learned stoichiometric coefficients or reaction rates are close to zero if the number of reactions are larger than two. It then can be inferred that the kinetics could be well described with two reactions.

reaction num 1     $x_1$     $x_2$     $x_3$     $k_f$     $k_r$

$1$ 	    $0.000$      $-1.000$     $1.000$      1.000e+03     1.000e+03

reaction num 2     $x_1$     $x_2$     $x_3$     $k_f$     $k_r$

$1$ 	    $0.000$      $1.000$     $-1.000$      1.000e+03     1.000e+03

$2$ 	    $-1.000$      $1.000$      $0.000$     1.000e+00     1.000e+00

reaction num 3     $x_1$     $x_2$     $x_3$     $k_f$     $k_r$

$1$ 	    $0.000$      $1.000$     $-1.000$      1.000e+03     1.000e+03

$2$ 	    $-1.000$      $1.000$      $0.000$     9.217e+02     1.218e+02

$3$ 	    $0.750$      $-0.101$      $0.384$     7.887e$-$04     1.813e$-$03

reaction num 4     $x_1$     $x_2$     $x_3$     $k_f$     $k_r$

$1$ 	    $0.000$      $-1.000$     $1.000$      1.000e+03     1.000e+03

$2$ 	    $0.000$      $0.000$      $0.000$     5.931e+01     5.929e+01

$3$ 	    $0.000$      $0.000$      $0.000$     2.526e+01     2.524e+01

$4$ 	    $1.000$      $-1.000$      $0.000$     1.000e+00     1.000e+00

Example <ref>: learned parameters with different number of reactions. Here $(x_1, x_2, x_3)$ denotes the row vector of the matrix $\bm{V}$.

Example <ref>: relative error for the training data and the validation data with different number of reactions.

[enzyme kinetics]

In this example, we consider the Michaelis–Menten kinetics [16], one of the best-known models of enzyme kinetics in biochemistry. It involves an enzyme E, binding to a substrate S, to form a complex ES, which in turn releases a product P, regenerating the original enzyme. This can be represented schematically as [16]

\begin{equation}\label{eq:enzyme-reaction}
\ch{E + S {<=>[$k\sb{f}$][$k\sb{r}$]} ES {->[$k\sb{cat}$]} E + P}
\end{equation}

Here $k_f$ denotes the forward rate constant, $k_r$ the reverse rate constant, and $k_{cat}$ the catalytic rate constant. This model is used in a variety of biochemical situations other than enzyme-substrate interaction, including antigen–antibody binding, DNA-DNA hybridization, and protein–protein interaction [26]. Moreover, the reaction rates vary widely between different enzymes. In our test case, we follow [36] and take $k_f=10^6$, $k_r=10^3$ and $k_{cat}=10$.

Note that the second reaction in (<ref>) is not reversible. Here, we show that the exactly same algorithm can be applied to this situation. The results with and without the PPF technique are listed in Table <ref>. In the upper part of the table, the reverse rate for the second reaction is $1.949\times10^{-4}$. It then can be inferred that the system could be well described using two reactions with the second one to be irreversible. Again, the algorithm without this treatment could only get the correct result for the first faster reaction in (<ref>). The evolution of the loss function is similar to that in Example <ref> and thus omitted here.

freezing     $x_1$     $x_2$     $x_3$     $x_4$     $k_f$     $k_r$

$1$ 	    $-1.000$      $-1.000$     $1.000$     $0.000$     1.000e+06     1.000e+03

$2$ 	    $1.000$      $0.000$      $-1.000$     $1.000$     1.000e+01     1.949e-04

no freezing     $x_1$     $x_2$     $x_3$     $x_4$     $k_f$     $k_r$

$1$ 	    $-0.999$      $-0.999$     $0.989$     $0.000$     9.921e+05     9.929e+02

$2$ 	    $-1.001$      $-1.000$     $2.385$     $0.000$     7.956e+03     1.392e+01

Example <ref>: learned parameters. Upper part: with the PPF technique; lower part: without the PPF technique. Here $(x_1, x_2, x_3, x_4)$ denotes the row vector of the matrix $\bm{V}$.

Next, we test the performance of the algorithm when the reaction rates depend on temperature. We assume that the rate constants in (<ref>) satisfy the Arrhenius law:

\begin{equation}
k_f = A_f \exp\brac{-\frac{E_{a,f}}{R T}}, \quad k_r = A_r \exp\brac{-\frac{E_{a,r}}{R T}}, \quad k_{cat} = A_{cat} \exp\brac{-\frac{E_{a,cat}}{R T}}
\end{equation}

where the pre-exponential factors are given by

\begin{equation}
A_f = 1, \quad A_r = 4, \quad A_{cat} = 10^3
\end{equation}

and the activation energy are

\begin{equation}
E_f = 1600, \quad E_r = 3680, \quad E_{cat} = 2240
\end{equation}

and the gas constant $R = 8.3145$. The temperature is randomly taken in a uniform distribution in the interval $[200, 400]$. In this case, the unknown parameters will include the pre-exponential factor, the activation energy in the Arrhenius law, and the stoichiometric coefficients. Our PPF technique can be directly applied without much modification.

We compare the performance of the algorithm with and without the PPF technique. The history of the relative error for the training data and the verification data with variable temperature is shown in Figure <ref>. We see clearly that the errors with the PPF technique are much smaller than those without the technique. We also show the learned parameters in Table <ref>. The upper part of the table is the learned parameters with the PPF technique, which agrees well with the ground truth. By contrast, the algorithm without imposing this technique could not generate the correct result.

Example <ref>: the history of the relative error for the training data and the verification data with variable temperature. Solid line: the method with the PPF technique; dashed line: the method without the PPF technique.

freezing     $x_1$     $x_2$     $x_3$     $x_4$     $A_f$     $A_r$     $E_f$     $E_r$

$1$     $1.000$      $0.000$      $-1.000$     $1.000$     1.000e+03     1.801e-05     2.240e+03     5.203e+03

$2$     $-1.000$     $-1.000$     $1.000$      $0.000$     1.000e+00     4.000e+00     1.600e+03     3.680e+03

no freezing     $x_1$     $x_2$     $x_3$     $x_4$     $A_f$     $A_r$     $E_f$     $E_r$

$1$     $1.061$      $0.001$     $-1.000$     $2.887$     3.391e+02     1.665e-06     2.240e+03     2.244e+03

$2$     $0.969$      $0.002$     $-1.000$     $0.032$     6.640e+02     2.689e-01     6.727e+01     6.684e+02

Example <ref>: learned parameters with variable temperatures. Upper part: with the PPF technique; lower part: without the PPF technique. Here $(x_1, x_2, x_3, x_4)$ denotes the row vector of the matrix $\bm{V}$.

[hydrogen oxidation reaction]
In this example, we consider a model for hydrogen oxidation reaction where six species H2 (hydrogen), O2 (oxygen), H2O (water), H, O, OH (radicals) are involved in six steps in a closed system under constant volume and temperature [11, 8]:

\begin{align}
\ch{H2 		&<=>[$k\sb{1}\sp{+}$][$k\sb{1}\sp{-}$] 2 H} \\
\ch{O2 		&<=>[$k\sb{2}\sp{+}$][$k\sb{2}\sp{-}$] 2 O} \\
\ch{H2O 	&<=>[$k\sb{3}\sp{+}$][$k\sb{3}\sp{-}$] H + OH} \\
\ch{H2 + O 	&<=>[$k\sb{4}\sp{+}$][$k\sb{4}\sp{-}$] H + OH} \\
\ch{O2 + H 	&<=>[$k\sb{5}\sp{+}$][$k\sb{5}\sp{-}$] O + OH} \\
\ch{H2 + O 	&<=>[$k\sb{6}\sp{+}$][$k\sb{6}\sp{-}$] H2O}
\end{align}

with the reaction rates $k_1^+=2$, $k_2^+=k_3^+=1$, $k_4^+=k_5^+=1\times10^3$, $k_1^+=1\times10^2$, $k_1^- = 2.16\times10^2$, $k_2^- = 3.375\times10^2$, $k_3^- = 1.4\times10^3$, $k_4^- = 1.08\times10^4$, $k_5^- = 3.375\times10^4$, $k_6^- = 7.714285714285716\times10^{-1}$.
The system (<ref>) corresponds to the simplified picture of this chemical process and the reaction rates reflect only orders of magnitude for relevant real-word systems.
The magnitude of the reaction rates vary from $10^{-1}$ to $10^4$, which leads to the multiscale phenomena. This reaction network has much more reactions and is more realistic than the first two test cases.

We first compare the performance of our algorithm with and without the PPF treatment. The history of the training and the validation error is shown in Figure <ref>. Again, we observe that this technique greatly reduces the training and validation errors. The learned parameters are listed in Table <ref>. The algorithm can generate the correct result with this technique. On the other hand, without using this technique, the phenomenon of the opposite signs observed in Table <ref> also appears.

Example <ref>: the history of the relative error for the training data and the verification data. Solid line: the method with the PPF technique; dashed line: the method without the PPF technique.

freezing     $x_1$     $x_2$     $x_3$     $x_4$     $x_5$     $x_6$     $k_f$     $k_r$

$1$ 	    $0.000$      $1.000$     $0.000$     $1.000$     $-1.000$     $-1.000$
3.375e+04     1.000e+03

$2$ 	    $1.000$      $0.000$      $0.000$     $-1.000$     $1.000$     $-1.000$
1.080e+04     1.000e+03

$3$ 	    $0.000$      $0.000$      $1.000$     $-1.000$     $0.000$     $-1.000$
1.400e+03     1.000e+00

$4$ 	    $0.000$      $1.000$      $0.000$     $0.000$     $-2.000$     $0.000$
3.375e+02     1.000e+00

$5$ 	    $1.000$      $0.000$      $0.000$     $-2.000$     $0.000$     $0.000$
2.160e+02     2.000e+00

$6$ 	    $-1.000$      $0.000$      $1.000$     $0.000$     $-1.000$     $0.000$
1.000e+02     7.714e$-$01

no freezing     $x_1$     $x_2$     $x_3$     $x_4$     $x_5$     $x_6$     $k_f$     $k_r$

1      $0.000 $       $0.938$       $0.000$       $ 0.923$       $-1.000$       $-1.001$
1.971e+04     5.512e+02

2      $0.000 $       $1.087$       $0.000$       $ 1.108$       $-0.999$       $-0.999$
1.407e+04     4.488e+02

3      $0.885 $       $0.000$       $0.115$       $-1.000$       $ 0.885$       $-1.000$
1.217e+04     2.112e+01

4     $-1.004$        $0.000$       $0.087$       $ 0.910$       $-1.008$       $ 0.913$
1.077e+03     2.926e+01

5      $0.000 $       $0.996$       $0.008$       $ 0.000$       $-1.997$       $ 0.000$
3.379e+02     8.912e$-$01

6      $0.987 $       $0.000$       $0.007$       $-1.984$       $ 0.004$       $ 0.004$
2.170e+02     3.377e+00

Example <ref>: learned parameters. Upper part: with the PPF technique; lower part: without the PPF technique. Here $(x_1, x_2, x_3, x_4, x_5, x_6)$ denotes the row vector of the matrix $\bm{V}$.

We also test the performance of the algorithm with Gaussian noise. The algorithm can get the correct prediction of the stoichiometric coefficients with the noise level $10^{-4}$ and $10^{-3}$. The learned reaction rates with noise are shown in Table <ref>. The relative errors for reaction rates are typically less than the order of $10^{-2}$ for $10^{-3}$ noise and $10^{-3}$ for $10^{-4}$ noise.

noise $10^{-3}$     $k_f$     relative error     $k_r$     relative error

1      3.375e+04       5.706e$-$05       1.002e+03       2.097e$-$03

2      1.080e+04       2.789e$-$04       1.001e+03       1.374e$-$03

3      1.399e+03       4.413e$-$04       9.631e$-$01       3.836e$-$02

4      3.399e+02       7.103e$-$03       9.235e$-$01       8.278e$-$02

5      2.161e+02       6.927e$-$04       2.130e+00       6.107e$-$02

6      9.764e+01       2.417e$-$02       8.047e$-$01       4.137e$-$02

noise $10^{-4}$     $k_f$     relative error     $k_r$     relative error

1      3.375e+04          6.482e$-$06       1.000e+03         2.099e$-$04

2      1.080e+04          2.803e$-$05       1.000e+03         1.379e$-$04

3      1.400e+03          4.491e$-$05       9.963e$-$01         3.707e$-$03

4      3.377e+02          7.161e$-$04       9.923e$-$01         7.717e$-$03

5      2.160e+02          6.965e$-$05       2.013e+00         6.469e$-$03

6      9.976e+01          2.366e$-$03       7.748e$-$01         4.294e$-$03

Example <ref>: learned reaction rates with noise.

Moreover, we plot the evolution of the concentrations of the six species with the noise level $10^{-3}$ in Figure <ref>. We observe a good agreement of the solution generated by our learned model and the exact solution.  We also measure the prediction errors of the learned model at 100 uniformly points in the time interval $[0, 10]$. The prediction errors are $1.953\times10^{-6}$ with zero noise, $9.152\times10^{-4}$ with noise level $10^{-4}$ and $8.710\times10^{-4}$ with noise level $10^{-3}$.

Example <ref>: the evolution of the concentration of the 6 species in the hydrogen oxidation reaction problem obtained by solving the original ODEs (<ref>) and our learned ODEs. noise level $10^{-3}$.

[extended Zeldovich mechanism]
In this example, we test our algorithm on the extended Zeldovich mechanism, which is a chemical mechanism describing the oxidation of nitrogen and NOx formation [41]. Similar to Example <ref>, this is another realistic test case. The reaction mechanisms read as

\begin{align}
\ch{N2 + O  &<=>[$k\sb{1}\sp{+}$][$k\sb{1}\sp{-}$] NO + N} \\
\ch{N + O2 	&<=>[$k\sb{2}\sp{+}$][$k\sb{2}\sp{-}$] NO + O} \\
\ch{N + OH 	&<=>[$k\sb{3}\sp{+}$][$k\sb{3}\sp{-}$] NO + H}
\end{align}

and the reaction rates are given by the Arrhenius law [12]:

\begin{equation}
\begin{aligned}
k_1^+ &= 1.8\times10^{11} \exp(-38370/T), \quad k_1^- = 3.8\times10^{10} \exp(-425/T), \\
k_2^+ &= 1.8\times10^{7} \exp(-4680/T), \quad k_2^- = 3.8\times10^{6} \exp(-20820/T), \\
k_3^- &= 7.1\times10^{10} \exp(-450/T), \quad k_3^- = 1.7\times10^{11} \exp(-24560/T),
\end{aligned}
\end{equation}

with $T$ the temperature.

In the numerical test, we fix the temperature to be $T=3000$, which is a reasonable temperature in real applications [12]. At this temperature, the reaction rates are

\begin{equation}
\begin{aligned}
& k_1^+ = 5.019\times10^5, \quad k_2^+ = 3.782\times10^6, \quad k_3^+ = 6.111\times10^{10}, \\
& k_1^- = 3.298\times10^{10}, \quad k_2^- = 3.679\times10^3, \quad k_3^- = 4.732\times10^7.
\end{aligned}
\end{equation}

Then, we follow the same procedure in the previous examples to generate the data and execute the algorithm to discover the stoichiometric coefficients and the reaction rates. Again, the algorithm with the PPF treatment can predict the correct result, which is shown in Table <ref>. We observe that the accurate reaction rates are obtained.

freezing     $x_1$     $x_2$     $x_3$     $x_4$     $x_5$     $x_6$     $x_7$     $k_f$     $k_r$

$1$ 	    $0.000$      $0.000$     $1.000$     $-1.000$     $0.000$     $-1.000$     $1.000$
6.111e+10     4.732e+07

$2$ 	    $1.000$      $1.000$      $-1.000$     $-1.000$     $0.000$     $0.000$     $0.000$
3.298e+10     5.019e+05

$3$ 	    $0.000$      $1.000$      $1.000$     $-1.000$     $-1.000$     $0.000$     $0.000$
3.782+06     3.931e+03

Example <ref>: learned parameters with the PPF technique. Here $(x_1, x_2, x_3, x_4, x_5, x_6, x_7)$ denotes the row vector of the matrix $\bm{V}$.

[simplified GRI-3.0 mechanism]
In this example, we test our algorithm on the simplified GRI-3.0 mechanism, which is a chemical mechanism describing the methane oxidation [37]. This is the most complicated reaction system tested in the paper. The mechanism includes 16 species with 12 reactions and reads as

\begin{align}
\ch{CH_4 + H & ->[$k\sb{1}$] CH_3 + H_2} \\
\ch{CH_2O + H_2 & ->[$k\sb{2}$] CH_3 + OH} \\
\ch{CH_2O & ->[$k\sb{3}$] CO + H_2} \\
\ch{C_2H_6 & ->[$k\sb{4}$] C_2H_4 + H_2} \\
\ch{C_2H_4 + OH & ->[$k\sb{5}$] CH_3 + CO + H_2} \\
\ch{2 CO + H_2 & ->[$k\sb{6}$] C_2H_2 + O_2} \\
\ch{CO + OH + H & ->[$k\sb{7}$] CO_2 + H_2} \\
\ch{H + OH & ->[$k\sb{8}$] H_2O} \\
\ch{2 H + 2 OH & ->[$k\sb{9}$] 2 H_2 + O_2} \\
\ch{H_2 & ->[$k\sb{10}$] 2 H} \\
\ch{H_2 + O_2 & ->[$k\sb{11}$] HO_2 + H} \\
\ch{H_2O_2 + H & ->[$k\sb{12}$] H_2 + HO_2}
\end{align}

The reaction rates are given in [37], which are derived from the reaction rates of the standard GRI-3.0 Mech [35]. We compute the reaction rates with the temperature $T=3000$ and list them in Table 4.8. Here, the reaction rates are normalized such that the smallest one is of order 1.

$k_1$     $k_2$     $k_3$     $k_4$     $k_5$     $k_6$

5.088e+00     1.891e+00     2.607e+00     6.268e+00     5.446e+00     1.283e+01

$k_7$     $k_8$     $k_9$     $k_{10}$     $k_{11}$     $k_{12}$

1.349e+00     5.264e+03     3.268e+01     4.873e+03     2.978e+02     5.227e+03

Example <ref>: reaction rates in simplified GRI-3.0 Mech

Note that all the reactions in (<ref>) are not reversible. Here, we apply exactly the same algorithm to this situation, similar to Example <ref>.
To illustrate the advantage of the PPF technique, we first compare the performance of the algorithm with and without this technique. The history of the training and validation errors is shown in Figure <ref>. The relative error stays around $10^{-3}$ without this technique, and decreases to $10^{-6}$ after applying this technique.

Example <ref>: the history of the relative error for the training data and the verification data. Solid line: the method with the PPF technique; dashed line: the method without the PPF technique.

$k^+_1$     $k^+_2$     $k^+_3$     $k^+_4$     $k^+_5$     $k^+_6$

5.088e+00     1.891e+00     2.607e+00     6.268e+00     5.446e+00     1.283e+01

$k^+_7$     $k^+_8$     $k^+_9$     $k^+_{10}$     $k^+_{11}$     $k^+_{12}$

1.349e+00     5.264e+03     3.268e+01     4.873e+03     2.978e+02     5.227e+03

$k^-_1$     $k^-_2$     $k^-_3$     $k^-_4$     $k^-_5$     $k^-_6$

2.546e-04     1.695e-04     1.091e-04     1.751e-04     1.103e-04     9.052e-06

$k^-_7$     $k^-_8$     $k^-_9$     $k^-_{10}$     $k^-_{11}$     $k^-_{12}$

8.462e-05     1.146e-05     5.472e-04     2.625e-07     8.566e-07     3.863e-04

Example <ref>: learned reaction rates in simplified GRI-3.0 Mech. Upper part: reaction rates in the forward reaction; lower part: reaction rates in the reverse reaction.

We also list the learned parameters with the PPF technique in Table 4.9.
Here, the learned stoichiometric coefficients are the same with the true coefficients in (<ref>) and they are omitted here. The upper part of the table is the learned rates in the forward reactions with the PPF technique, which agrees well with the ground truth in Table 4.8.
The learned rates in the reverse reactions are in the magnitude of $10^{-7}$ to $10^{-4}$. It then can be inferred that the system can be well described using only forward reactions. By contrast, the algorithm without imposing this technique could not generate the correct result and we omit the results here.

§ CONCLUSION

In this paper, we propose a data-driven method to discover multiscale chemical reactions governed by the law of mass action. The method mainly contains two novel points.
First, we use a single matrix to represent the stoichiometric coefficients for both the reactants and products in a system without catalysis reactions.
The negative entries in the matrix denote the stoichiometric coefficients for the reactants and the positive ones for the products.
Second, by considering a multiscale nonlinear regression problem, we find that the conventional optimization methods usually get stuck in the local minima and could not find the true solution. To escape from the local minima, we propose a PPF technique. Notice that the stoichiometric coefficients are integers. In the training process, if the loss function stops to decrease, the stoichiometric coefficients which are close to integers are rounded and then frozen afterwards. With such a treatment, the stoichiometric coefficients are gradually determined, the dimension of the searching space is reduced in the training process, and eventually the global mimina can be obtained. Several numerical experiments including the classical Michaelis–Menten kinetics, the hydrogen oxidation reactions and simplified GRI-3.0 mechanism verify the validity of our algorithm in learning the multiscale chemical reactions.

There are still some problems to be addressed in order to develop a robust and general framework for discovering multiscale chemical reactions from data. We shall highlight some of the challenges that could guide future advances. First, it is interesting to generalize the PPF technique to the catalysis reactions.

Second, the number of species $n_s$ cannot be determined from our algorithm. In principle, to infer the unknown chemical reaction systems, we should have the concentration time series data for all the species. Our algorithm cannot treat the problem when the concentrations for some partial species are unknown. This difficulty may be overcomed by combining the current algorithm with the Neural ODE approach in [15].
The third challenge is that for very complex reaction networks with large number of reactions (hundreds or thousands), our algorithm may not always find out the correct solution. New ideas are needed at this point.

[1]
B. Bhadriraju, M. S. F. Bangi, A. Narasingam, and S. I. Kwon.
Operable adaptive sparse identification of systems (OASIS):
application to chemical processes.
AIChE Journal, 66(11):e16980, 2020.

[2]
B. Bhadriraju, A. Narasingam, and J. S.-I. Kwon.
Machine learning-based adaptive model identification of systems:
Application to a chemical process.
Chemical Engineering Research and Design, 152:372–383, 2019.

[3]
S. L. Brunton, B. R. Noack, and P. Koumoutsakos.
Machine learning for fluid mechanics.
Annual Review of Fluid Mechanics, 52:477–508, 2020.

[4]
S. L. Brunton, J. L. Proctor, and J. N. Kutz.
Discovering governing equations from data by sparse identification of
nonlinear dynamical systems.
Proceedings of the National Academy of Sciences,
113(15):3932–3937, 2016.

[5]
S. C. Burnham, D. P. Searson, M. J. Willis, and A. R. Wright.
Inference of chemical reaction networks.
Chemical Engineering Science, 63(4):862–873, 2008.

[6]
K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton.
Data-driven discovery of coordinates and governing equations.
Proceedings of the National Academy of Sciences,
116(45):22445–22451, 2019.

[7]
R. Chartrand.
Numerical differentiation of noisy, nonsmooth data.
ISRN Applied Mathematics, 2011, 2011.

[8]
E. Chiavazzo and I. V. Karlin.
Quasi-equilibrium grid algorithm: Geometric construction for model
Journal of Computational Physics, 227(11):5535–5560, 2008.

[9]
B. de Silva, K. Champion, M. Quade, J.-C. Loiseau, J. N. Kutz, and S. Brunton.
PySINDy: A Python package for the sparse identification of
nonlinear dynamical systems from data.
Journal of Open Source Software, 5(49):1–4, 2020.

[10]
C. W. Gao, J. W. Allen, W. H. Green, and R. H. West.
Reaction mechanism generator: Automatic construction of chemical
kinetic mechanisms.
Computer Physics Communications, 203:212–225, 2016.

[11]
A. N. Gorban and I. V. Karlin.
Invariant manifolds for physical and chemical kinetics, volume
Springer Science & Business Media, 2005.

[12]
R. K. Hanson and S. Salimian.
Survey of rate constants in the n/h/o system.
In Combustion chemistry, pages 361–421. Springer, 1984.

[13]
M. Hoffmann, C. Fröhner, and F. Noé.
Reactive SINDy: Discovering governing reactions from concentration
The Journal of Chemical Physics, 150(2):025101, 2019.

[14]
J. Huang, Z. Ma, Y. Zhou, and W.-A. Yong.
Learning thermodynamically stable and galilean invariant partial
differential equations for non-equilibrium flows.
Journal of Non-Equilibrium Thermodynamics, 2021.

[15]
W. Ji and S. Deng.
Autonomous discovery of unknown reaction pathways from data by
chemical reaction neural network.
The Journal of Physical Chemistry A, 125(4):1082–1092, 2021.

[16]
J. P. Keener and J. Sneyd.
Mathematical Physiology, volume 1.
Springer, 1998.

[17]
D. P. Kingma and J. Ba.
Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980, 2014.

[18]
D. Langary and Z. Nikoloski.
Inference of chemical reaction networks based on concentration
profiles using an optimization framework.
Chaos: An Interdisciplinary Journal of Nonlinear Science,
29(11):113121, 2019.

[19]
Y. LeCun, Y. Bengio, and G. Hinton.
Deep learning.
Nature, 521(7553):436–444, 2015.

[20]
L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis.
DeepXDE: A deep learning library for solving differential
arXiv preprint arXiv:1907.04502, 2019.

[21]
T. Lu and C. K. Law.
On the applicability of directed relation graphs to the reduction of
reaction mechanisms.
Combustion and Flame, 146(3):472–483, 2006.

[22]
B. Lusch, J. N. Kutz, and S. L. Brunton.
Deep learning for universal linear embeddings of nonlinear dynamics.
Nature Communications, 9(1):1–10, 2018.

[23]
S. Maddu, B. L. Cheeseman, C. L. Müller, and I. F. Sbalzarini.
Learning physically consistent mathematical models from data using
group sparsity.
arXiv preprint arXiv:2012.06391, 2020.

[24]
N. M. Mangan, S. L. Brunton, J. L. Proctor, and J. N. Kutz.
Inferring biological networks by sparse identification of nonlinear
IEEE Transactions on Molecular, Biological and Multi-Scale
Communications, 2(1):52–63, 2016.

[25]
T. Nagy, J. Tóth, and T. Ladics.
Automatic kinetic model generation and selection based on
concentration versus time curves.
International Journal of Chemical Kinetics, 52(2):109–123,

[26]
D. L. Nelson, A. L. Lehninger, and M. M. Cox.
Lehninger principles of biochemistry.
Macmillan, 2008.

[27]
H. G. Othmer.
Analysis of complex reaction networks.
Lecture Notes, School of Mathematics, University of Minnesota,

[28]
A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen,
Z. Lin, N. Gimelshein, L. Antiga, et al.
Pytorch: An imperative style, high-performance deep learning library.
Advances in Neural Information Processing Systems,
32:8026–8037, 2019.

[29]
M. Raissi, P. Perdikaris, and G. E. Karniadakis.
Physics-informed neural networks: A deep learning framework for
solving forward and inverse problems involving nonlinear partial differential
Journal of Computational Physics, 378:686–707, 2019.

[30]
M. Raissi, A. Yazdani, and G. E. Karniadakis.
Hidden fluid mechanics: Learning velocity and pressure fields from
flow visualizations.
Science, 367(6481):1026–1030, 2020.

[31]
R. Ranade, S. Alqahtani, A. Farooq, and T. Echekki.
An ANN based hybrid chemistry framework for complex fuels.
Fuel, 241:625–636, 2019.

[32]
R. Ranade, S. Alqahtani, A. Farooq, and T. Echekki.
An extended hybrid chemistry framework for complex hydrocarbon fuels.
Fuel, 251:276–284, 2019.

[33]
L. I. Rudin, S. Osher, and E. Fatemi.
Nonlinear total variation based noise removal algorithms.
Physica D: Nonlinear Phenomena, 60(1-4):259–268, 1992.

[34]
S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz.
Data-driven discovery of partial differential equations.
Science Advances, 3(4):e1602614, 2017.

[35]
G. P. Smith, D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eiteneer,
M. Goldenberg, C. T. Bowman, R. K. Hanson, S. Song, W. C. J. Gardiner, V. V.
Lissianski, and Z. Qin.
Gri-mech home page, http://www.me.berkeley.edu/gri_mech/.

[36]
V. Srinivasan and R. Aiken.
Stage-wise parameter estimation for stiff differential equations.
AIChE journal, 32(2):195–199, 1986.

[37]
C. Sung, C. Law, and J.-Y. Chen.
Augmented reduced mechanisms for no emission in methane oxidation.
Combustion and Flame, 125(1-2):906–919, 2001.

[38]
E. O. Voit, H. A. Martens, and S. W. Omholt.
150 years of the mass action law.
PLoS Comput Biol, 11(1):e1004012, 2015.

[39]
G. Wanner and E. Hairer.
Solving Ordinary Differential Equations II.
Springer Berlin Heidelberg, 1996.

[40]
M. J. Willis and M. von Stosch.
Inference of chemical reaction networks using mixed integer linear
Computers & Chemical Engineering, 90:31–43, 2016.

[41]
I. Zeldovich, G. I. Barenblatt, V. Librovich, and G. Makhviladze.
Mathematical theory of combustion and explosions.


# The parabolic Sturmian-function basis representation of the six-dimensional
Coulomb Green’s function

S. A. Zaytsev<EMAIL_ADDRESS>Pacific National University,
Khabarovsk, 680035, Russia

###### Abstract

The square integrable basis set representation of the resolvent of the
asymptotic three-body Coulomb wave operator in parabolic coordinates is
obtained. The resulting six-dimensional Green’s function matrix is expressed
as a convolution integral over separation constants.

###### pacs:

03.65.Nk

##  I Introduction

It is well known that the Schrödinger equation for a three-body Coulomb system
is asymptotically separable in terms of parabolic coordinates Klar

$\xi_{j}=r_{ls}+\hat{\bf k}_{ls}\cdot{\bf
r}_{ls},\quad\eta_{j}=r_{ls}-\hat{\bf k}_{ls}\cdot{\bf r}_{ls},$ (1)

where ${\bf r}_{ls}$ and ${\bf k}_{ls}$ are the relative coordinate and
momentum vectors between the particles $l$ and $s$. Here $j$, $l$, $s$ is a
cyclic permutation of 1, 2, 3. The long-ranged six-dimensional operator, which
provides a three-body continuum wave function, 3C-function C31 ; C32 , that
satisfies exact asymptotic boundary conditions for Coulomb systems in the
region where the distances between the particles are large, reads as the sum
of two-dimensional operators Klar :

$\sum_{j=1}^{3}\frac{1}{\mu_{ls}(\xi_{j}+\eta_{j})}\left[\hat{h}_{\xi_{j}}+\hat{h}_{\eta_{j}}+2k_{ls}t_{ls}\right],$
(2)

where $t_{ls}=\frac{Z_{l}Z_{s}\mu_{ls}}{k_{ls}}$,
$\mu_{ls}=\frac{m_{l}m_{s}}{m_{l}+m_{s}}$. Here the one-dimensional operators
$\hat{h}_{\xi_{j}}$ and $\hat{h}_{\eta_{j}}$ are defined by

$\begin{array}[]{ccc}\hat{h}_{\xi_{j}}&=&-2\left(\frac{\partial}{\partial\xi_{j}}\xi_{j}\frac{\partial}{\partial\xi_{j}}+ik_{ls}\xi_{j}\frac{\partial}{\partial\xi_{j}}\right),\\\
\hat{h}_{\eta_{j}}&=&-2\left(\frac{\partial}{\partial\eta_{j}}\eta_{j}\frac{\partial}{\partial\eta_{j}}-ik_{ls}\eta_{j}\frac{\partial}{\partial\eta_{j}}\right).\end{array}$
(3)

In the previous paper previous  we introduce the six-dimensional operator
$\mathfrak{h}$ which is obtained by multiplying (2) on the left by
$\prod_{j=1}^{3}\mu_{ls}(\xi_{j}+\eta_{j})$:

$\hat{\mathfrak{h}}=\mu_{13}(\xi_{2}+\eta_{2})\,\mu_{12}(\xi_{3}+\eta_{3})\hat{\mathfrak{h}}_{1}+\mu_{23}(\xi_{1}+\eta_{1})\,\mu_{12}(\xi_{3}+\eta_{3})\hat{\mathfrak{h}}_{2}+\mu_{23}(\xi_{1}+\eta_{1})\,\mu_{13}(\xi_{2}+\eta_{2})\hat{\mathfrak{h}}_{3},$
(4)

where

$\hat{\mathfrak{h}}_{j}=\hat{h}_{\xi_{j}}+\hat{h}_{\eta_{j}}+2k_{ls}t_{ls}.$
(5)

The resolvent of the operator (4) can be used in the corresponding Lippmann-
Schwinger equation for the three-body Coulomb wave function.

It has been suggested previous  to treat the operator (4) within the context
of $L^{2}$ parabolic Sturmian basis set Ojha1

$\left|\mathcal{N}\right>=\prod_{j=1}^{3}\phi_{n_{j}\,m_{j}}(\xi_{j},\,\eta_{j}),$
(6)
$\phi_{n_{j}\,m_{j}}(\xi_{j},\,\eta_{j})=\psi_{n_{j}}(\xi_{j})\,\psi_{m_{j}}(\eta_{j}),$
(7) $\psi_{n}(x)=\sqrt{2b}e^{-bx}L_{n}(2bx),$ (8)

where $b$ is the scaling parameter. In particular, a matrix representation
${\bf G}_{j}$ of the resolvent for the two-dimensional operator

$\hat{\mathfrak{h}}_{j}+\mu_{ls}C_{j}(\xi_{j}+\eta_{j})$ (9)

has been obtained which is formally the matrix inverse to the infinite matrix
$\left[{\bf h}_{j}+\mu_{ls}C_{j}{\bf Q}_{j}\right]$ of the operator (9):

$\left[{\bf h}_{j}+\mu_{ls}C_{j}{\bf Q}_{j}\right]{\bf
G}_{j}(t_{ls};\;\mu_{ls}C_{j})={\bf I}_{j}.$ (10)

Here

${\bf h}_{j}={\bf h}_{\xi_{j}}\otimes{\bf I}_{\eta_{j}}+{\bf
I}_{\xi_{j}}\otimes{\bf h}_{\eta_{j}}+2k_{ls}t_{ls}{\bf I}_{j}$ (11)

is the matrix of the operator $\hat{\mathfrak{h}}_{j}$ (5) in the basis (7),
${\bf I}_{\xi_{j}}$, ${\bf I}_{\eta_{j}}$ and ${\bf I}_{j}={\bf
I}_{\xi_{j}}\otimes{\bf I}_{\eta_{j}}$ are the unit matrices. ${\bf
Q}_{j}={\bf Q}_{\xi_{j}}\otimes{\bf I}_{\eta_{j}}+{\bf I}_{\xi_{j}}\otimes{\bf
Q}_{\eta_{j}}$, where ${\bf Q}_{\xi_{j}}$ and ${\bf Q}_{\eta_{j}}$ are the
matrices of $\xi_{j}$ and $\eta_{j}$ in basis (8), respectively.

In the previous paper previous  it has been suggested that the matrices ${\bf
G}_{j}$ of the two-dimensional Green’s functions be used to construct the
matrix $\underline{\mathfrak{G}}$ which is inverse to the operator (4) matrix

$\underline{\mathfrak{h}}=\mu_{13}\mu_{12}{\bf h}_{1}\otimes{\bf
Q}_{2}\otimes{\bf Q}_{3}+\mu_{23}\mu_{12}{\bf Q}_{1}\otimes{\bf
h}_{2}\otimes{\bf Q}_{3}+\mu_{23}\mu_{13}{\bf Q}_{1}\otimes{\bf
Q}_{2}\otimes{\bf h}_{3}.$ (12)

Namely, we proposed to express the six-dimensional Green’s function matrix
$\underline{\mathfrak{G}}$ as the convolution integral

$\underline{\mathfrak{G}}=\aleph\int\int dC_{1}dC_{2}\,{\bf
G}_{1}(t_{23};\;\mu_{23}C_{1})\otimes{\bf
G}_{2}(t_{13};\;\mu_{13}C_{2})\otimes{\bf
G}_{3}(t_{12};\;-\mu_{12}(C_{1}+C_{2})),$ (13)

where $\aleph$ is a normalizing factor. Thus, our problem is now to determine
the pathes of integration over the separation constants $C_{1}$ and $C_{2}$ in
(13) and to find the corresponding normalizing factor $\aleph$ such that the
condition

$\underline{\mathfrak{h}}\underline{\mathfrak{G}}={\bf I}_{1}\otimes{\bf
I}_{2}\otimes{\bf I}_{3}$ (14)

holds. For this purpose consider the product
$\underline{\mathfrak{h}}\;\underline{\mathfrak{G}}$. From the relation (10)
we obtain

$\begin{array}[]{c}\underline{\mathfrak{h}}\;\underline{\mathfrak{G}}=\aleph\left\\{\int\int
dC_{1}dC_{2}\left[{\bf I}_{1}-\mu_{23}C_{1}{\bf Q}_{1}{\bf
G}_{1}(t_{23};\;\mu_{23}C_{1})\right]\otimes\mu_{13}{\bf Q}_{2}{\bf
G}_{2}(t_{13};\;\mu_{13}C_{2})\right.\\\\[8.53581pt] \otimes\mu_{12}{\bf
Q}_{3}{\bf G}_{3}(t_{12};\;-\mu_{12}(C_{1}+C_{2}))\\\\[8.53581pt] \int\int
dC_{1}dC_{2}\,\mu_{23}{\bf Q}_{1}{\bf
G}_{1}(t_{23};\;\mu_{23}C_{1})\otimes\left[{\bf I}_{2}-\mu_{13}C_{2}{\bf
Q}_{2}{\bf G}_{2}(t_{13};\;\mu_{13}C_{2})\right]\\\\[8.53581pt]
\otimes\mu_{12}{\bf Q}_{3}{\bf
G}_{3}(t_{12};\;-\mu_{12}(C_{1}+C_{2}))\\\\[8.53581pt] \int\int
dC_{1}dC_{2}\,\mu_{23}{\bf Q}_{1}{\bf
G}_{1}(t_{23};\;\mu_{23}C_{1})\otimes\mu_{13}{\bf Q}_{2}{\bf
G}_{2}(t_{13};\;\mu_{13}C_{2})\\\\[8.53581pt] \left.\otimes\left[{\bf
I}_{3}+\mu_{12}(C_{1}+C_{2}){\bf Q}_{3}{\bf
G}_{3}(t_{12};\;-\mu_{12}(C_{1}+C_{2}))\right]\right\\},\\\ \end{array}$ (15)

and hence

$\begin{array}[]{c}\underline{\mathfrak{h}}\;\underline{\mathfrak{G}}=\aleph\left\\{\int\int
dC_{1}dC_{2}\left[{\bf I}_{1}\otimes\mu_{13}{\bf Q}_{2}{\bf
G}_{2}(t_{13};\;\mu_{13}C_{2})\otimes\mu_{12}{\bf Q}_{3}{\bf
G}_{3}(t_{12};\;-\mu_{12}(C_{1}+C_{2}))\right.\right.\\\\[8.53581pt]
\left.+\mu_{23}{\bf Q}_{1}{\bf G}_{1}(t_{23};\;\mu_{23}C_{1})\otimes{\bf
I}_{2}\otimes\mu_{12}{\bf Q}_{3}{\bf
G}_{3}(t_{12};\;-\mu_{12}(C_{1}+C_{2}))\right]\\\\[8.53581pt]
\left.+\left[\mu_{23}{\bf Q}_{1}\int dC_{1}\,{\bf
G}_{1}(t_{23};\;\mu_{23}C_{1})\right]\otimes\left[\mu_{13}{\bf Q}_{2}\int
dC_{2}\,{\bf G}_{2}(t_{13};\;\mu_{13}C_{2})\right]\otimes{\bf
I}_{3}\right\\}.\end{array}$ (16)

As a first step towards our goal we consider the integrals

$\int dC_{j}\,{\bf G}_{j}(t_{ls};\;\mu_{ls}C_{j})$ (17)

inside the figure brackets on the right-hand side of (16).

In Sec. II completeness of the eigenfunctions of the two-dimensional operator
(9) is considered. In particular, an integral representation of the matrix
${\bf A}$ which is inverse to the infinite matrix ${\bf Q}$ of the operator
$(\xi+\eta)$ in the basis (7) is obtained. In Sec. III it is demonstrated that
the integral (17) taken along an appropriate contour is proportional to the
matrix ${\bf A}$ obtained in the previous section. Finally, Sec. IV presents a
convolution integral representation of the six-dimensional Coulomb Green’
function matrix.

##  II Completeness relations

###  II.1 The continuous spectrum

Of particular interest are the regular solutions

$f(\gamma,\,\tau,\,\xi,\,\eta)=u(\gamma,\,\tau,\,\xi)\,v(\gamma,\,\tau,\,\eta)$
(18)

of the system

$\left[\hat{h}_{\xi}+2kt+\mu\,C\xi\right]u(\gamma,\,\tau,\,\xi)=0,$ (19)
$\left[\hat{h}_{\eta}+2k(t_{0}-t)+\mu\,C\eta\right]v(\gamma,\,\tau,\,\xi)=0.\\\
$ (20)

Obviously, the regular solutions $u$ and $v$ are proportional to confluent
hypergeometric functions Abramowitz :

$u(\gamma,\,\tau,\,\xi)=e^{\frac{i}{2}(\gamma-k)\xi}\,{{}_{1}F_{1}}\left(\frac{1}{2}+i\tau,\,1,\,-i\gamma\xi\right)$
(21)

and

$\begin{array}[]{l}v(\gamma,\,\tau,\,\eta)=e^{-\frac{i}{2}(\gamma-k)\eta}\,{{}_{1}F_{1}}(\frac{1}{2}+i(\tau-\tau_{0}),\,1,\,i\gamma\eta)\hfill\\\\[8.53581pt]
\hfill=e^{\frac{i}{2}(\gamma+k)\eta}\,{{}_{1}F_{1}}(\frac{1}{2}+i(\tau_{0}-\tau),\,1,\,-i\gamma\eta),\\\
\end{array}$ (22)

where

$\mu\,C=\frac{k^{2}}{2}-\frac{\gamma^{2}}{2},\quad\tau=\frac{k}{\gamma}\left(t+\frac{i}{2}\right),\quad\tau_{0}=\frac{k}{\gamma}t_{0}.$
(23)

It should be noted that since the representation of the two- and six-
dimensional Coulomb Green’s functions matrix elements involves an integration
over $\tau$ from $-\infty$ to $\infty$, in the subsequent discussion we assume
that $\tau$ is real. With this assumption it is readily seen that the
solutions (21) and (22) coincide, except for normalization and the phase
factors $e^{-\frac{i}{2}k\xi}$ and $e^{\frac{i}{2}k\eta}$, with parabolic
Coulomb Sturmians treated in Ref. Gasaneo3 . In this case $\gamma$ plays the
role of the momentum and $\mathcal{E}=\frac{\gamma^{2}}{2}$ is the energy.
From this we conclude that for $\gamma>0$ the solutions (21) and (22)
correspond to the continuous spectrum of $\mathcal{E}$.

It is readily verified that the solutions $u(\gamma,\,\tau,\,\xi)$ and
$v(\gamma,\,\tau,\,\eta)$ are expressed by basis set (8) expansions

$u(\gamma,\,\tau,\,\xi)=\frac{2\sqrt{2b}}{2b-i(\gamma-k)}\left(\frac{2b-i(\gamma-k)}{2b+i(\gamma+k)}\right)^{i\tau+\frac{1}{2}}\sum\limits_{n=0}^{\infty}\theta^{n}\,p_{n}\left(\tau;\;\zeta\right)\,\psi_{n}(\xi),$
(24)
$v(\gamma,\,\tau,\,\eta)=\frac{2\sqrt{2b}}{2b-i(\gamma+k)}\left(\frac{2b-i(\gamma+k)}{2b+i(\gamma-k)}\right)^{i(\tau_{0}-\tau)+\frac{1}{2}}\sum\limits_{n=0}^{\infty}\lambda^{-n}\,p_{n}\left(\tau_{0}-\tau;\;\zeta\right)\,\psi_{n}(\eta),$
(25)

where

$\theta=\frac{2b+i(\gamma-k)}{2b-i(\gamma-k)},\qquad\lambda=\frac{2b-i(\gamma+k)}{2b+i(\gamma+k)},\qquad\zeta=\frac{\lambda}{\theta}.$
(26)

The expansion (24) and (25) coefficients contain the polynomials previous

$p_{n}\left(\tau;\;\zeta\right)=\frac{(-1)^{n}}{n!}\frac{\Gamma\left(n+\frac{1}{2}-i\tau\right)}{\Gamma\left(\frac{1}{2}-i\tau\right)}\;{{}_{2}F_{1}}{\left(-n,\,\frac{1}{2}+i\tau;\;-n+\frac{1}{2}+i\tau;\;\zeta\right)}.$
(27)

The basis set (8) representation of the equation (19) is the three-term
recursion relation previous

$a_{n}\,y_{n-1}+b_{n}\,y_{n}+d_{n+1}\,y_{n+1}=0,\quad n\geq 1$ (28)

where

$\begin{array}[]{c}b_{n}=(b+\frac{\mu C}{2b}+ik)+2(b+\frac{\mu
C}{2b})n+2kt,\\\\[5.69054pt] a_{n}=(b-\frac{\mu C}{2b}-ik)n,\quad
d_{n}=(b-\frac{\mu C}{2b}+ik)n.\\\ \end{array}$ (29)

The functions

$s_{n}(t;\;\mu C)=\theta^{n}\,p_{n}(\tau;\;\zeta)$ (30)

are the “regular” solutions of (28) with the initial conditions: $s_{0}\equiv
1$, $s_{-1}\equiv 0$. The polynomials $p_{n}$ (27) of degree $n$ in $\tau$ are
orthogonal with respect to the weight function previous

$\rho(\tau;\zeta)=\frac{\Gamma\left(\frac{1}{2}-i\tau\right)\Gamma\left(\frac{1}{2}+i\tau\right)}{2\pi
i}\,(-\zeta)^{i\tau+\frac{1}{2}},$ (31)

where it is considered that $\left|\arg(-\zeta)\right|<\pi$. The corresponding
orthogonality relation reads

$\frac{i}{\zeta^{m}}\left(\frac{\zeta-1}{\zeta}\right)\int\limits_{-\infty}^{\infty}d\tau\,\rho(\tau;\;\zeta)\,p_{n}\left(\tau;\;\zeta\right)\,p_{m}\left(\tau;\;\zeta\right)=\delta_{n\,m}.$
(32)

###  II.2 The discrete spectrum

For $t_{0}<0$ the eigenfunctions of (19) corresponding to the discrete
spectrum: $\mathcal{E}^{(\ell)}=\frac{\gamma_{\ell}^{2}}{2}$,
$\gamma_{\ell}=i\kappa_{\ell}$, $\kappa_{\ell}=-\frac{kt_{0}}{\ell}$,
$\ell=1,2,\,\ldots,\,\infty$, are Landau

$f_{\ell,\,m}(\xi,\,\eta)=u_{\ell,\,m}(\xi)\,v_{\ell,\,\ell-m-1}(\eta),\qquad
m=0,\,1,\ldots\,,\,\ell-1,$ (33)

where

$u_{\ell,\,m}(\xi)=e^{-\frac{i}{2}k\xi}e^{-\frac{\kappa_{\ell}\xi}{2}}\,{{}_{1}F_{1}}(-m,\,1;\;\kappa_{\ell}\xi)=e^{-\frac{i}{2}k\xi}e^{-\frac{\kappa_{\ell}\xi}{2}}\,L_{m}(\kappa_{\ell}\xi)$
(34)

and

$v_{\ell,\,m}(\eta)=e^{\frac{i}{2}k\eta}e^{-\frac{\kappa_{\ell}\eta}{2}}\,{{}_{1}F_{1}}(-m,\,1;\;\kappa_{\ell}\eta)=e^{\frac{i}{2}k\eta}e^{-\frac{\kappa_{\ell}\eta}{2}}\,L_{m}(\kappa_{\ell}\eta).$
(35)

The solutions $f_{\ell,\,m}$ meet the orthogonality relation

$\frac{\kappa_{\ell}^{3}}{2\ell}\int\limits_{0}^{\infty}\int\limits_{0}^{\infty}(\xi+\eta)d\xi
d\eta\,f_{\ell,\,m}(\xi,\,\eta)\,\left[f_{{\ell\,}^{\prime},\,m^{\prime}}(\xi,\,\eta)\right]^{*}=\delta_{\ell,\,{\ell\,^{\prime}}}\,\delta_{m,\,m^{\prime}}.$
(36)

It is readily verified that the expansions of $u_{\ell,\,m}(\xi)$ and
$v_{\ell,\,m}(\eta)$ in a basis function (8) series are

$u_{\ell,\,m}(\xi)=\sum\limits_{n=0}^{\infty}\mathcal{S}^{(\ell,\,m)}_{n}\,\psi_{n}(\xi)$
(37)

and

$v_{\ell,\,m}(\eta)=\sum\limits_{n=0}^{\infty}\left[\mathcal{S}^{(\ell,\,m)}_{n}\right]^{*}\psi_{n}(\eta),$
(38)

where the coefficients are given by

$\begin{array}[]{l}\mathcal{S}^{(\ell,\,m)}_{n}=2\sqrt{2b}\,(-)^{n}\,\frac{(m+1)_{n}}{n!}\frac{(2b-\kappa_{\ell}-ik)^{n}(2b-\kappa_{\ell}+ik)^{m}}{(2b+\kappa_{\ell}+ik)^{n+m+1}}\hfill\\\\[5.69054pt]
\hfill\times{{}_{2}F_{1}}\left(-n,\,-m,\,-n-m;\;\frac{(2b+\kappa_{\ell})^{2}+k^{2}}{(2b-\kappa_{\ell})^{2}+k^{2}}\right).\\\
\end{array}$ (39)

### One-dimensional completeness relations

The eigensolution $u(\gamma,\,\tau,\,\xi)$ of (19) corresponding to the
continuous spectrum ($\gamma>0$) for large $\xi$ behaves as

$u(\gamma,\,\tau,\,\xi)\mathop{\sim}\limits_{\xi\rightarrow\infty}\,\frac{e^{-\frac{i}{2}k\xi}\,e^{\frac{\pi\tau}{2}}}{\left|\Gamma(\frac{1}{2}+i\tau)\right|}\,\frac{2i}{\sqrt{\gamma\xi}}\,\sin\left(\frac{\gamma\xi}{2}-\tau\ln(\gamma\xi)+\frac{\pi}{4}+\sigma\right),$
(40)

where $\sigma=\arg\Gamma\left(\frac{1}{2}+i\tau\right)$. Therefore (see e. g.
Michel ),

$\frac{\gamma}{4}\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}\,e^{-\pi\tau}\,\int\limits_{0}^{\infty}\xi\,d\xi\,u(\gamma,\,\tau,\,\xi)\,\left[u(\gamma^{\prime},\,\tau,\,\xi)\right]^{*}=\pi\delta(\gamma-\gamma^{\prime}),$
(41)

and for $\tau>0$

$\xi\int\limits_{0}^{\infty}\gamma\,d\gamma\,\frac{\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}}{4\pi}\,e^{-\pi\tau}\,u(\gamma,\,\tau,\,\xi)\,\left[u(\gamma,\,\tau,\,\xi^{\prime})\right]^{*}=\delta(\xi-\xi^{\prime}).$
(42)

On the other hand, if the functions $u(\gamma,\,\tau,\,\xi)$ are regarded as
charge Sturmians Gasaneo3 , i. e. the parameter $\tau$ is considered as
eigenvalue of the problem, whereas the momentum $\gamma$ remains constant, the
corresponding orthogonality and completeness relations are given by Gasaneo3

$\gamma\,\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}\,e^{-\pi\tau}\int\limits_{0}^{\infty}d\xi\,u(\gamma,\,\tau,\,\xi)\left[u(\gamma,\,\tau^{\prime},\,\xi)\right]^{*}=2\pi\delta(\tau-\tau^{\prime})$
(43)

and

$\gamma\,\int\limits_{-\infty}^{\infty}d\tau\,\frac{\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}}{2\pi}\,e^{-\pi\tau}\,u(\gamma,\,\tau,\,\xi)\left[u(\gamma,\,\tau,\,\xi^{\prime})\right]^{*}=\delta(\xi-\xi^{\prime}).$
(44)

Taking matrix elements of the completeness relation (44) we find

$\frac{8\,b\,\gamma\,\theta^{n-m}}{\sqrt{\left[4b^{2}+(\gamma+k)^{2}\right]\left[4b^{2}+(\gamma-k)^{2}\right]}}\,\int\limits_{-\infty}^{\infty}d\tau\,\frac{\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}}{2\pi}\,(-\zeta)^{i\tau}p_{n}\left(\tau;\;\zeta\right)\left[p_{n}\left(\tau;\;\zeta\right)\right]^{*}=\delta_{n,\,m}.$
(45)

It may be noted that (45) is closely related to (32). To see this, let
$\gamma>0$. Then, it follows from the definitions (26) that

$\frac{(\zeta-1)}{\zeta}(-\zeta)^{\frac{1}{2}}=\frac{8\,b\,\gamma}{\sqrt{\left[4b^{2}+(\gamma+k)^{2}\right]\left[4b^{2}+(\gamma-k)^{2}\right]}}.$
(46)

Further, the regular solution $\theta^{n}p_{n}\left(\tau;\;\zeta\right)$ of
the three-term recursion relation (28) is an even function of $\gamma$, since
the coefficients $a_{n}$, $b_{n}$ and $d_{n}$ depend on $\gamma$ only through
$\mu C=\frac{1}{2}(k^{2}-\gamma^{2})$. Thus, replacing $\gamma$ by $-\gamma$,
and hence $\tau$ by $-\tau$ ($\theta\rightarrow\lambda$ and
$\lambda\rightarrow\theta$) and $\zeta$ by $\zeta^{-1}$ in Eq. (27) gives

$\theta^{n}\,p_{n}\left(\tau;\;\zeta\right)=\lambda^{n}\frac{(-1)^{n}}{n!}\frac{\Gamma\left(n+\frac{1}{2}+i\tau\right)}{\Gamma\left(\frac{1}{2}+i\tau\right)}\;{{}_{2}F_{1}}{\left(-n,\,\frac{1}{2}-i\tau;\;-n+\frac{1}{2}-i\tau;\;\zeta^{-1}\right)}.$
(47)

Comparing Eqs. (27) and (47) then yields the relation

$\theta^{n}\,p_{n}\left(\tau;\;\zeta\right)=\lambda^{n}\left[p_{n}\left(\tau;\;\zeta\right)\right]^{*},$
(48)

and hence

$\left[p_{n}\left(\tau;\;\zeta\right)\right]^{*}=\zeta^{-n}\,p_{n}\left(\tau;\;\zeta\right).$
(49)

From the argument above, we conclude that for $\gamma>0$ the orthogonality
relation (32) reduces to (45).

### The two-dimensional completeness relation

It follows from the relations (41) and (43) and analogous relations for
$v(\gamma,\,\tau,\,\xi)$ that the two-dimensional orthogonality relation for
$\gamma>0$ is given by

$\begin{array}[]{c}e^{-\pi\tau_{0}}\,\frac{\gamma^{2}}{4}\,\frac{\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}}{2\pi}\frac{\left|\Gamma\left(\frac{1}{2}+i(\tau_{0}-\tau)\right)\right|^{2}}{2\pi}\,\int\limits_{0}^{\infty}\int\limits_{0}^{\infty}(\xi+\eta)\,d\xi
d\eta\left\\{f(\gamma,\,\tau,\,\xi,\,\eta)\right.\\\\[8.53581pt]
\qquad\qquad\qquad\left.\times\left[f(\gamma^{\prime},\,\tau^{\prime},\,\xi,\,\eta)\right]^{*}\right\\}=\delta(\gamma-\gamma^{\prime})\,\delta(\tau-\tau^{\prime}).\\\
\end{array}$ (50)

In turn, in the case $t_{0}>0$ (where there are no bound states) it would
appear reasonable that the two-dimensional completeness relation would be
given by

$\begin{array}[]{l}(\xi+\eta)\left\\{\alpha\int\limits_{0}^{\infty}d\gamma\,\gamma^{2}e^{-\pi\tau_{0}}\int\limits_{-\infty}^{\infty}d\tau\frac{\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}}{2\pi}\frac{\left|\Gamma\left(\frac{1}{2}+i(\tau_{0}-\tau)\right)\right|^{2}}{2\pi}\right.\hfill\\\\[8.53581pt]
\hfill\left.\times
f(\gamma,\,\tau,\,\xi,\,\eta)\left[f(\gamma,\,\tau,\,\xi^{\prime},\,\eta^{\prime})\right]^{*}\vphantom{\int\limits_{-\infty}^{\infty}}\right\\}=\delta(\xi-\xi^{\prime})\,\delta(\eta-\eta^{\prime}).\\\
\end{array}$ (51)

The integration over $\tau$ in (51) is performed on the assumption that $\tau$
is independent of $\gamma$. To test this hypothesis and determine the
normalizing factor $\alpha$, we carried out the following numerical
experiments. First with some parameters $t_{0}>0$, $k>0$ and $b$ we calculate
the matrix elements $A_{n_{1},\,n_{2};\;m_{1},\,m_{2}}$ for the expression in
the figure braces on the right-hand side of (51) in the basis (7):

$\begin{array}[]{l}A_{n_{1},\,n_{2};\;m_{1},\,m_{2}}=\alpha\int\limits_{0}^{\infty}d\gamma\,\frac{64\,b^{2}\,\gamma^{2}\,(-\zeta)^{i\tau_{0}}}{\left[4b^{2}+(\gamma+k)^{2}\right]\left[4b^{2}+(\gamma-k)^{2}\right]}\,\theta^{n_{1}-m_{1}}\lambda^{m_{2}-n_{2}}\hfill\\\\[8.53581pt]
\times\int\limits_{-\infty}^{\infty}d\tau\frac{\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}}{2\pi}\frac{\left|\Gamma\left(\frac{1}{2}+i(\tau_{0}-\tau)\right)\right|^{2}}{2\pi}\,p_{n_{1}}\left(\tau;\;\zeta\right)p_{n_{2}}\left(\tau_{0}-\tau;\;\zeta\right)\left[p_{m_{1}}\left(\tau;\;\zeta\right)p_{m_{2}}\left(\tau_{0}-\tau;\;\zeta\right)\right]^{*}.\\\
\end{array}$ (52)

It should be noted that the value of $(-\zeta)^{i\tau_{0}}$ in this formula is
determined by the condition $\left|\arg(-\zeta)\right|<\pi$. Then the
resulting matrix ${\bf A}$ is multiplied by the matrix

${\bf Q}={\bf Q}_{\xi}\otimes{\bf I}_{\eta}+{\bf I}_{\xi}\otimes{\bf
Q}_{\eta}$ (53)

of the operator $(\xi+\eta)$. Finally, using the condition

${\bf Q}\,{\bf A}={\bf I}_{\xi}\otimes{\bf I}_{\eta}.$ (54)

we have obtained that $\alpha=\frac{1}{2}$. Notice that the infinite symmetric
matrices ${\bf Q}_{\xi}$ and ${\bf Q}_{\eta}$ are tridiagonal previous ,
therefore the equations on the left hand-side of the linear system (54) each
contain no more than six terms.

For for $t_{0}<0$ the completeness relation (51) transforms into

$\begin{array}[]{c}(\xi+\eta)\left\\{\frac{1}{2}\int\limits_{0}^{\infty}d\gamma\,\gamma^{2}e^{-\pi\tau_{0}}\int\limits_{-\infty}^{\infty}d\tau\frac{\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}}{2\pi}\frac{\left|\Gamma\left(\frac{1}{2}+i(\tau_{0}-\tau)\right)\right|^{2}}{2\pi}\,f(\gamma,\,\tau,\,\xi,\,\eta)\left[f(\gamma,\,\tau,\,\xi^{\prime},\,\eta^{\prime})\right]^{*}\right.\\\\[8.53581pt]
\left.\hfill+\sum\limits_{\ell=1}^{\infty}\frac{\kappa_{\ell}^{3}}{2\ell}\,\sum\limits_{m=0}^{\ell-1}f_{\ell,\,m,\,\ell-m-1}(\xi,\,\eta)\left[f_{\ell,\,m,\,\ell-m-1}(\xi^{\prime},\,\eta^{\prime})\right]^{*}\right\\}=\delta(\xi-\xi^{\prime})\,\delta(\eta-\eta^{\prime}).\\\
\end{array}$ (55)

In this case the matrix ${\bf A}$ with elements

$\begin{array}[]{l}A_{n_{1},\,n_{2};\;m_{1},\,m_{2}}=\int\limits_{0}^{\infty}d\gamma\,\frac{32\,b^{2}\,\gamma^{2}\,(-\zeta)^{i\tau_{0}}}{\left[4b^{2}+(\gamma+k)^{2}\right]\left[4b^{2}+(\gamma-k)^{2}\right]}\,\theta^{n_{1}-m_{1}}\lambda^{m_{2}-n_{2}}\hfill\\\\[8.53581pt]
\times\int\limits_{-\infty}^{\infty}d\tau\frac{\left|\Gamma\left(\frac{1}{2}+i\tau\right)\right|^{2}}{2\pi}\frac{\left|\Gamma\left(\frac{1}{2}+i(\tau_{0}-\tau)\right)\right|^{2}}{2\pi}\,p_{n_{1}}\left(\tau;\;\zeta\right)p_{n_{2}}\left(\tau_{0}-\tau;\;\zeta\right)\left[p_{m_{1}}\left(\tau;\;\zeta\right)p_{m_{2}}\left(\tau_{0}-\tau;\;\zeta\right)\right]^{*}\\\\[8.53581pt]
\hfill+\sum\limits_{\ell=1}^{\infty}\frac{\kappa_{\ell}^{3}}{2\ell}\,\sum\limits_{m=0}^{\ell-1}\mathcal{S}^{(\ell,\,m)}_{n_{1}}\,\mathcal{S}^{(\ell,\,\ell-m-1)}_{m_{2}}\left[\mathcal{S}^{(\ell,\,\ell-m-1)}_{n_{2}}\,\mathcal{S}^{(\ell,\,m)}_{m_{1}}\right]^{*},\\\
\end{array}$ (56)

is also inverse to the matrix ${\bf Q}$ (53). The expression (56) can be
rewritten, in view of (46) and (49), as

$\begin{array}[]{l}A_{n_{1},\,n_{2};\;m_{1},\,m_{2}}=-\frac{1}{2}\int\limits_{0}^{\infty}d\gamma\left\\{\left(\frac{\zeta-1}{\zeta}\right)^{2}\,\frac{\theta^{n_{1}+m_{2}}}{\lambda^{n_{2}+m_{1}}}\,\int\limits_{-\infty}^{\infty}d\tau\rho(\tau;\;\zeta)\,\rho(\tau_{0}-\tau;\;\zeta)\right.\\\\[8.53581pt]
\hfill\left.\times
p_{n_{1}}(\tau;\;\zeta)\,p_{n_{2}}(\tau_{0}-\tau;\;\zeta)\,p_{m_{1}}(\tau;\;\zeta)\,p_{m_{2}}(\tau_{0}-\tau;\;\zeta)\vphantom{\left(\frac{\zeta-1}{\zeta}\right)^{2}}\right\\}\\\\[8.53581pt]
\hfill+\sum\limits_{\ell=1}^{\infty}\frac{\kappa_{\ell}^{3}}{2\ell}\,\sum\limits_{m=0}^{\ell-1}\mathcal{S}^{(\ell,\,m)}_{n_{1}}\,\mathcal{S}^{(\ell,\,\ell-m-1)}_{m_{2}}\left[\mathcal{S}^{(\ell,\,\ell-m-1)}_{n_{2}}\,\mathcal{S}^{(\ell,\,m)}_{m_{1}}\right]^{*}.\\\
\end{array}$ (57)

##  III Contour integrals

Notice that expressing the resolvent of the one-dimensional operator
$\left[\hat{h}_{\xi}+2kt+\mu\,C\xi\right]$ requires two linearly independent
solutions of (19). Irregular solutions of (19) are expressed in terms the
confluent hypergeometric function Abramowitz

$w^{(\pm)}(\gamma,\,\tau,\,\xi)=e^{\frac{i}{2}(\pm\gamma-k)\xi}\,{U\left(\frac{1}{2}\pm
i\tau,\,1;\;\mp\gamma\xi\right)}.$ (58)

The corresponding solutions of the three-term recursion relation (28) are

$\begin{array}[]{c}c_{n}^{(+)}(t;\;\mu
C)=\theta^{n+1}\,q^{(+)}_{n}\left(\tau;\;\zeta\right),\\\\[8.53581pt]
c_{n}^{(-)}(t;\;\mu C)=\lambda^{n+1}\,q^{(-)}_{n}\left(\tau;\;\zeta\right),\\\
\end{array}$ (59)

where

$\begin{array}[]{c}q_{n}^{(+)}(\tau;\;\zeta)=(-)^{n}\,\frac{n!\,\Gamma\left(\frac{1}{2}+i\tau\right)}{\Gamma\left(n+\frac{3}{2}+i\tau\right)}\;{{}_{2}F_{1}}\left(\frac{1}{2}+i\tau,\,n+1;\;n+\frac{3}{2}+i\tau;\;\zeta^{-1}\right),\\\\[8.53581pt]
q_{n}^{(-)}(\tau;\;\zeta)=(-)^{n}\,\frac{n!\,\Gamma\left(\frac{1}{2}-i\tau\right)}{\Gamma\left(n+\frac{3}{2}-i\tau\right)}\;{{}_{2}F_{1}}\left(\frac{1}{2}-i\tau,\,n+1;\;n+\frac{3}{2}-i\tau;\;\zeta\right).\\\
\end{array}$ (60)

In particular, the functions

$\begin{array}[]{l}\widetilde{w}^{(\pm)}(\gamma,\,\tau,\,\xi)=\frac{2i\sqrt{2b}}{2b-i(\gamma-k)}\left(-\frac{2b+i(\gamma-k)}{2b-i(\gamma+k)}\right)^{i\tau+\frac{1}{2}}\hfill\\\\[8.53581pt]
\hfill\times\frac{e^{-\pi\tau}}{\theta\,\Gamma\left(\frac{1}{2}\pm
i\tau\right)}\sum\limits_{n=0}^{\infty}c_{n}^{(\pm)}(t;\;\mu
C)\,\psi_{n}(\xi)\\\ \end{array}$ (61)

tend to $w^{(\pm)}(\gamma,\,\tau,\,\xi)$ as $\xi\rightarrow\infty$.

The matrix elements of the resolvent of
$\left[\hat{h}_{\xi}+2kt+\mu\,C\xi\right]$ can be written in the form previous

$g^{(+)}_{n,\,m}(t;\;\mu
C)=\frac{i}{2\gamma}\left(\frac{\zeta-1}{\zeta}\right)\frac{\theta^{n-m}}{\zeta^{m}}\,p_{n_{<}}\left(\tau;\;\zeta\right)\,q^{(+)}_{n_{>}}(\tau;\;\zeta)$
(62)

and

$g^{(-)}_{n,\,m}(t;\;\mu
C)=\frac{i}{2\gamma}\left(\frac{\zeta-1}{\zeta}\right)\frac{\theta^{n-m}}{\zeta^{m}}\,p_{n_{<}}\left(\tau;\;\zeta\right)\,\zeta^{n_{>}+1}\,q^{(-)}_{n_{>}}(\tau;\;\zeta),$
(63)

where $n_{>}$ and $n_{<}$ are the greater and lesser of $n$ and $m$. Notice
that $c_{n}^{(+)}$ $\left(c_{n}^{(-)}\right)$ are defined in the upper (lower)
half of the complex $\gamma$-plane where $|\zeta|\geq 1$ $\left(|\zeta|\leq
1\right)$. To analytically continue $c_{n}^{(+)}$ onto the lower half of the
$\gamma$-plane the relation previous

$c_{n}^{(+)}(t;\;\mu C)=c_{n}^{(-)}(t;\;\mu C)+2\pi
i\rho(\tau;\;\zeta)\,\theta^{n+1}p_{n}\left(\tau;\;\zeta\right)$ (64)

can be used.

In previous  we obtained the basis set (7) representation of the resolvent for
the two-dimensional operator
$\left[\hat{h}_{\xi}+2kt+\mu\,C\xi\right]+\left[\hat{h}_{\eta}+2k(t_{0}-t)+\mu\,C\eta\right]$.
In particular, the matrix elements of the two-dimensional Green’s function can
be expressed as the convolution integral

$\begin{array}[]{l}G^{(\pm)}_{n_{1},\,n_{2};\;m_{1},\,m_{2}}(t_{0};\;\mu
C)=i\left(\frac{\zeta-1}{\zeta}\right)\,\frac{\lambda^{m_{2}-n_{2}}}{\zeta^{m_{2}}}\int\limits_{-\infty}^{\infty}d\tau\,\rho(\tau_{0}-\tau;\;\zeta)\\\\[8.53581pt]
\hfill\times g^{(\pm)}_{n_{1},\,m_{1}}(t;\;\mu
C)\,p_{n_{2}}(\tau_{0}-\tau;\;\zeta)\,p_{m_{2}}(\tau_{0}-\tau;\;\zeta).\end{array}$
(65)

Notice that in this case only the regular solutions of (20) discrete analogues
$\lambda^{-n}p_{n}(\tau_{0}-\tau;\;\zeta)$ are used.

Let us consider the integral

$\mathcal{I}_{1}=\frac{1}{2\pi
i}\int\limits_{-\infty}^{\infty}\gamma\,d\gamma\,G^{(+)}_{n_{1},\,n_{2};\;m_{1},\,m_{2}}(t_{0};\;\mu
C).$ (66)

Notice that by replacing $\gamma\rightarrow-\gamma$ (and hence
$\theta\rightarrow\lambda$, $\lambda\rightarrow\theta$, $\zeta\rightarrow
1/\zeta$, $\tau\rightarrow-\tau$ and $\tau_{0}\rightarrow-\tau_{0}$) in Eq.
(65) $G^{(+)}_{n_{1},\,n_{2};\;m_{1},\,m_{2}}(t_{0};\;\mu C)$ is transformed
to $G^{(-)}_{n_{1},\,n_{2};\;m_{1},\,m_{2}}(t_{0};\;\mu C)$. Thus, for the
integral $\mathcal{I}_{1}$ we obtain

$\mathcal{I}_{1}=\frac{1}{2\pi
i}\int\limits_{0}^{\infty}\gamma\,d\gamma\left\\{G^{(+)}_{n_{1},\,n_{2};\;m_{1},\,m_{2}}(t_{0};\;\mu
C)-G^{(-)}_{n_{1},\,n_{2};\;m_{1},\,m_{2}}(t_{0};\;\mu C)\right\\}.$ (67)

Inserting Eqs. (65), (62) and (63) into Eq. (67), we find, in view of (64),
that $\mathcal{I}_{1}$ coincides with the integral on the right hand-side of
(57).

Now, we consider the integral

$\frac{1}{2\pi i}\int\limits_{\mathcal{C}}\,d\mathcal{E}\,{\bf
G}^{(+)}\left(t_{0};\;\frac{k^{2}}{2}-\mathcal{E}\right),$ (68)

taken along the contour in the complex $\mathcal{E}$-plane shown in Fig. 1.
The contour $\mathcal{C}$ passes in a negative direction (clockwise) round all
the points $\mathcal{E}^{(\ell)}=-\frac{\kappa_{\ell}^{2}}{2}$ (filled circles
in Fig. 1 which accumulate at the origin) and the cut along the right half of
the real axis and is closed at infinity (see e. g. Shakeshaft ). The
corresponding matrix element of the integral along the two sides of the cut is
equal to $\mathcal{I}_{1}$ (66) (this is circumstantial evidence that the
normalizing factor $\alpha$ in the completeness relation (51) is equal to
$\frac{1}{2}$). On the other hand, the integration along a contour enclosing
$\mathcal{E}^{(\ell)}$ reduces to $(-1)$ times the double sum of the residues
of the integrand at the points $\tau^{(m)}=i\left(m+\frac{1}{2}\right)$,
$m=0,\,1,\,\ldots$ and
${\mathcal{E}^{(\ell)}=-\frac{(kt_{0})^{2}}{2\ell^{2}}}$,
$\ell=1,\,2,\,\ldots\,$, which are the poles of the gamma functions
$\Gamma\left(\frac{1}{2}+i\tau\right)$ and
$\Gamma\left(\frac{1}{2}+i(\tau_{0}-\tau)\right)=\Gamma\left(m+1+\frac{kt_{0}}{\sqrt{-2\mathcal{E}^{(\ell)}}}\right)$,
respectively. It is readily shown that the matrix element of this part of the
integral (68) coincides with the double sum in (57). Thus, the integral (68)
is equal to the matrix ${\bf A}$. The contour $\mathcal{C}$ can be deformed,
for instance, into a straight line parallel to the real axis. The resulting
path $\mathcal{C}_{1}$ shown in Fig. 2 runs above the cut and the bound-state
poles of ${\bf G}^{(+)}\left(t_{0};\;\frac{k^{2}}{2}-\mathcal{E}\right)$.
Then, to make the integral amenable to numerical integration, rotate
$\mathcal{C}_{1}$ about some point $\mathcal{E}_{0}$ on the right half of the
real axis through a negative angle $\varphi$ Shakeshaft ; see the contour
$\mathcal{C}_{2}$ in Fig. 2. The part of $\mathcal{C}_{2}$ on the unphysical
sheet is depicted by the dashed line. Thus, we obtain the following
representation of the matrix ${\bf A}$ (which is inverse to the matrix ${\bf
Q}$ (53)):

${\bf A}=\frac{1}{2\pi i}\int\limits_{\mathcal{C}_{2}}\,d\mathcal{E}\,{\bf
G}^{(+)}\left(t_{0};\;\frac{k^{2}}{2}-\mathcal{E}\right)={\bf Q}^{-1}.$ (69)

##  IV Six-dimensional Green’s function matrix

Using the relation (69) we can rewrite the expression (13) for the six-
dimensional Coulomb Green’ function matrix as the contour integral

$\underline{\mathfrak{G}}=\frac{\aleph}{\mu_{23}\mu_{13}}\int\limits_{\mathcal{C}_{2}}d\mathcal{E}_{1}\int\limits_{\mathcal{C}_{2}}d\mathcal{E}_{2}\,{\bf
G}_{1}\left(t_{23};\;\frac{k_{23}^{2}}{2}-\mathcal{E}_{1}\right)\otimes{\bf
G}_{2}\left(t_{13};\;\frac{k_{13}^{2}}{2}-\mathcal{E}_{2}\right)\otimes{\bf
G}_{3}\left(t_{12};\;\frac{k_{12}^{2}}{2}-\mathcal{E}_{3}\right),$ (70)

where $\mathcal{E}_{j}=\frac{k_{ls}^{2}}{2}-\mu_{ls}C_{j}$ and
$C_{1}+C_{2}+C_{3}=0$. This also allows us to determine the normalizing factor
$\aleph$. Indeed, it follows from (69) that the third term inside the figure
brackets on the right-hand side of (16) is proportional to the unit matrix:

$\begin{array}[]{c}\left[{\bf
Q}_{1}\int\limits_{\mathcal{C}_{2}}d\mathcal{E}_{1}\,{\bf
G}_{1}\left(t_{23};\;\frac{k_{23}^{2}}{2}-\mathcal{E}_{1}\right)\right]\otimes\left[{\bf
Q}_{2}\int\limits_{\mathcal{C}_{2}}d\mathcal{E}_{2}\,{\bf
G}_{2}\left(t_{13};\;\frac{k_{13}^{2}}{2}-\mathcal{E}_{2}\right)\right]\otimes{\bf
I}_{3}\\\\[8.53581pt] =(2\pi i)^{2}\,{\bf I}_{1}\otimes{\bf I}_{2}\otimes{\bf
I}_{3}.\\\ \end{array}$ (71)

Consider the first two terms in the figure braces in (16). For the energy
$\mathcal{E}_{3}=\frac{k_{12}^{2}}{2}+\mu_{12}(C_{1}+C_{2})$ we have

$\mathcal{E}_{3}=\frac{k_{12}^{2}}{2}+\frac{\mu_{12}}{\mu_{23}}\left(\frac{k_{23}^{2}}{2}-\mathcal{E}_{1}\right)+\frac{\mu_{12}}{\mu_{13}}\left(\frac{k_{13}^{2}}{2}-\mathcal{E}_{2}\right).$
(72)

On the other hand, on the contours $\mathcal{C}_{2}$ the energy variables
$\mathcal{E}_{j}$, $j=1,\,2$ are given by

$\mathcal{E}_{j}=\mathcal{E}_{0j}+E_{j}e^{i\varphi},$ (73)

where $\varphi<0$, $\mathcal{E}_{0j}$ is an arbitrary positive parameter,
$E_{j}$ is real and runs from $-\infty$ to $\infty$. Hence, for the energy
$\mathcal{E}_{3}$ (72) we obtain

$\mathcal{E}_{3}=\left[\frac{k_{12}^{2}}{2}+\frac{\mu_{12}}{\mu_{23}}\left(\frac{k_{23}^{2}}{2}-\mathcal{E}_{01}\right)+\frac{\mu_{12}}{\mu_{13}}\left(\frac{k_{13}^{2}}{2}-\mathcal{E}_{02}\right)\right]+\left(-\frac{\mu_{12}}{\mu_{23}}\,E_{1}-\frac{\mu_{12}}{\mu_{13}}\,E_{2}\right)e^{i\varphi}.$
(74)

Thus, $\mathcal{E}_{3}$ can be expressed in the form

$\mathcal{E}_{3}=\mathcal{E}_{03}+E_{3}\,e^{i\varphi},$ (75)

where $\mathcal{E}_{03}$ and $E_{3}$ denote the term in the square braces and
the real factor in front of the exponent in (74), respectively. Since
$\mathcal{E}_{03}$ should be positive, the positive parameters
$\mathcal{E}_{01}$ and $\mathcal{E}_{02}$ have to satisfy the constraint

$\frac{\mu_{12}}{\mu_{23}}\,\mathcal{E}_{01}+\frac{\mu_{12}}{\mu_{13}}\,\mathcal{E}_{02}<\frac{k_{12}^{2}}{2}+\frac{\mu_{12}}{\mu_{23}}\frac{k_{23}^{2}}{2}+\frac{\mu_{12}}{\mu_{13}}\frac{k_{13}^{2}}{2}.$
(76)

Now, we consider the integral

$\mathcal{I}_{2}=\int\limits_{\mathcal{C}_{2}}d\mathcal{E}_{1}{\bf
G}_{3}\left(t_{12};\;\frac{k_{12}^{2}}{2}-\mathcal{E}_{3}\right).$ (77)

With fixed $\mathcal{E}_{2}$, in view (74), (75) and (69), we see that

$\begin{array}[]{c}\mathcal{I}_{2}=e^{i\varphi}\int\limits_{-\infty}^{\infty}dE_{1}{\bf
G}_{3}\left(t_{12};\;\frac{k_{12}^{2}}{2}-\mathcal{E}_{3}\right)=-\frac{\mu_{23}}{\mu_{12}}\,e^{i\varphi}\int\limits_{-\infty}^{\infty}dE_{3}{\bf
G}_{3}\left(t_{12};\;\frac{k_{12}^{2}}{2}-\mathcal{E}_{3}\right)\\\\[8.53581pt]
=-\frac{\mu_{23}}{\mu_{12}}\int\limits_{\mathcal{C}_{2}}d\mathcal{E}_{3}{\bf
G}_{3}\left(t_{12};\;\frac{k_{12}^{2}}{2}-\mathcal{E}_{3}\right)=-2\pi
i\,\frac{\mu_{23}}{\mu_{12}}\,{\bf Q}_{3}^{-1}.\end{array}$ (78)

Similarly, we obtain

$\int\limits_{\mathcal{C}_{2}}d\mathcal{E}_{2}{\bf
G}_{3}\left(t_{12};\;\frac{k_{12}^{2}}{2}-\mathcal{E}_{3}\right)=-2\pi
i\,\frac{\mu_{13}}{\mu_{12}}\,{\bf Q}_{3}^{-1}.$ (79)

Inserting (78), (79) and (71) into (16) then yields

$\underline{\mathfrak{h}}\;\underline{\mathfrak{G}}=\aleph\,4\pi^{2}\,{\bf
I}_{1}\otimes{\bf I}_{2}\otimes{\bf I}_{3}.$ (80)

Thus, from (80) we conclude that

$\aleph=\frac{1}{4\pi^{2}}.$ (81)

##  V Conclusion

The Sturmian basis-set representation of the resolvent for the asymptotic
three-body Coulomb wave operator is obtained, which can be used in the
discrete analog of the Lippmann-Schwinger equation for the three-body
continuum wave function. The six-dimensional Green’s function matrix is
expressed as a convolution integral over separation constants. The integrand
of this contour integral involves Green’s function matrices corresponding to
the two-dimensional operators which are constituents of the full six-
dimensional wave operator. The completeness relation of the eigenfunctions of
these two-dimensional operators is used to define the appropriate pathes of
integration of the convolution integral.

## References

  * (1) H. Klar, Z. Phys. D 16, 231 (1990).
  * (2) Dz. Belkic, J. Phys. B 11, 3529 (1978).
  * (3) C. R. Garibotti, J. E. Miraglia, Phys. Rev. A 21, 572 (1980).
  * (4) S. A. Zaytsev, J. Phys. A, to be published.
  * (5) P. C. Ojha, J. Math. Phys. 28, 392 (1987).
  * (6) M. Abramowitz and I. A. Stegun, Handbook of mathematical functions (New York: Dover), 1970\.
  * (7) G. Gasaneo and F. D. Colavecchia, J. Phys. A 36 8443 (2003).
  * (8) L. D. Landau and E. M. Lifshitz, Quantum Mechanics: Non-Relativistic Theory (Oxford: Pergamon), 1965.
  * (9) N. Michel, J. Math. Phys. 49 022109 (2008).
  * (10) R. Shakeshaft, Phys. Rev. A, 70, 042704 (2004).

Figure 1: The path of integration $\mathcal{C}$ in (68).

Figure 2: The path of integration $\mathcal{C}_{2}$ in (69).