# Observing how deep neural networks understand physics through the energy
spectrum of one-dimensional quantum mechanics

Kenzo Ogure<EMAIL_ADDRESS>Department of Nuclear Engineering,
Kyoto University, Kyoto daigaku-katsura, Nishikyo-ku, Kyoto, 615-8540, Japan

###### Abstract

We investigated how neural networks(NNs) understand physics using one-
dimensional quantum mechanics. After training an NN to accurately predict
energy eigenvalues from potentials, we used it to confirm the NN’s
understanding of physics from four different aspects. The trained NN could
predict energy eigenvalues of different kinds of potentials than the ones
learned, predict the probability distribution of the existence of particles
not used during training, reproduce untrained physical phenomena, and predict
the energy eigenvalues of potentials with unknown matter effect. These results
show that NNs can learn the physical laws from experimental data, predict the
results of experiments under conditions different from those used for
training, and predict physical quantities of types not provided during
training. Because NNs understand physics in a different way than humans, they
will be a powerful tool for advancing physics by complementing the human way
of understanding.

A44 Neural networks, A2 Computational physics

##  1 Introduction

In recent times, deep neural networks (DNNs) have made remarkable progress in
image recognition, natural language processing, voice recognition, and anomaly
detection through numerous technological breakthroughs. Among these
breakthroughs, the residual connection prevents gradient loss, even when
neural networks (NNs) have deep layers, contributing to the very high image
recognition capabilityhe2016deep . Another example is the attention mechanism,
which has succeeded in connecting neurons in distant locations, a shortcoming
of convolutional neural networks (CNNs), and has significantly advanced fields
such as translation, where relationships between distant words are
essentialvaswani2017attention ; devlin2018bert . In addition to accuracy
improvement, NNs have begun generating new images or sentences by themselves,
and their performance has been rapidly improving kingma2013auto ;
goodfellow2014generative . In a slightly different field, the combination of
DNNs and reinforcement learning has rendered humans incapable of competing
with computers in table games such as Go and Shogi, where human intuition was
previously superiorsilverMastering2016 ; silverMastering2017 .

The significant difference between NNs and previous artificial intelligence is
that NNs do not require human intuition at the design stage. Earlier types of
artificial intelligence relied on humans to predetermine features of objects
and perform learning. In contrast, NNs do not require such human intuition and
find features in the learning process by themselves. In return, NNs comprise
numerous parameters and require high computing power for learning. In
addition, it is difficult for humans to understand how trained NNs operate
since they do not have the predetermined features humans want to observe to
understand the operation.

Another advantage of NNs is that the input and output can be exchanged in some
cases, enabling their application to inverse problems. If the inputs and
outputs of the problem are bijections, they can be swapped and trained to
create a NN that solves the inverse problem.

The properties of not requiring prior intuition and applicability to inverse
problems make NNs promising in natural sciences, where the goal is to reveal
unknown problems, and NNs have been used in many fields, including
physicscarleoMachine2019 . In this paper, the subject of this study is the
Schrödinger equation in quantum mechanics, and even if we focus only on the
surrounding area, we can see several recent developments related to NNs.
Partial differential equations, including the Schrödinger equation, have been
solved using NNsraissiPhysicsinformed2019 , potentials have been estimated
inversely from wave functionshongPredicting2021 ; sehanobishLearning2021 , and
soliton solutions to the nonlinear Schrödinger equation have been investigated
using NNspuSoliton2020 . A new family of NNs inspired by the Schrödinger
equation has also been proposednakajimaNeural2021 .

In these studies, the Schrödinger equation was used, and in this sense, human
understanding of physics was already assumed. However, since we consider this
study preparation for applying NNs to unsolved physical systems, we do not
directly solve the Schrödinger equation using NNs. From this viewpoint, the
research conducted in a context most similar to this paper is
Ref.millsDeep2017 . This pioneering study provided NNs with two-dimensional
(2D) potentials and the corresponding energy eigenvalues of the ground states.
This prior work and our study differ in technical aspects, such as the method
for generating the potentials and NN’s architecture, but the most significant
difference is that we are most interested in how NNs understand physics. In
this sense, despite the differences in physical systems and methods, our study
shares a common interest with Ref.itenDiscovering2020 .

The application of NNs to concrete physics problems seems to be progressing
well. However, how an NN learns physics is embedded in numerous parameters and
is difficult to find. In particular, whether the NN is merely learning its
output pattern or understanding the underlying physics is essential in
considering the future use of NNs. If the NN understands physical laws, it
should be able to predict physical quantities of various ranges and types that
were not used in training. It should also be possible to reproduce physical
phenomena and estimate physical quantities even if the laws governing that
physical system were not yet known. This study sheds light on these points
using one-dimensional (1D) quantum mechanics as a subject. Since our research
assumes that we will use NNs to solve unknown physics problems in the future,
we only provide NN potentials and a few energy eigenvalues, which are physical
quantities that can be measured experimentally.

This paper is organized as follows. First, in Sec.2, after a brief
introduction to NNs for physicists, we demonstrate how to generate datasets
and training results in this paper. Then, in Sec.3, where we state our main
results, we investigate four aspects of how NNs understand physics and the
potential for NNs to be used in physics research. Next, section 4 summarizes
the results of this paper and provides some insights into the future potential
of NNs in physics. Finally, the NN architecture used in this paper is
summarized in the appendix.

##  2 Setup and training of neural networks

This section first presents a brief introduction to NNs for physicists and
then explains how NNs are trained in this paper on 1D quantum mechanics
problems in Sec.2.1. Next, we explain how to prepare the datasets in Sec.2.2,
and finally, we present the training results in Sec.2.3.

###  2.1 Application of Neural Networks to one-dimensional Quantum Mechanics

Regression using an NN is a kind of variational method if we use a term
familiar to physicists. An NN is a function that contains parameters, which
are determined to reproduce the best ground truth outputs. This function,
inspired by the brain’s structure, has several characteristics:

  * •

While the variational trial function in physics is typically assumed to be an
elementary function or a superposition of elementary functions by some
physical intuition, NNs do not require such intuition. In return, NNs have far
more parameters than trial functions in traditional variational methods in
physics to be sufficiently flexible.

  * •

Since NNs are developed by a composite function of linear transformations and
simple functions called activation functions, the back-propagation technique
can quickly obtain their derivatives.

  * •

As a technical matter, NNs are beginning to be firmly recognized for their
usefulness, and supportive libraries such as Pytorch and Tensorflow have been
developed, as well as technologies such as Cuda for faster computation on
GPUs. This has enabled experimenting with NN technology to solve different
problems easily.

Here, we explain the training of NNs in more detail. Training an NN requires a
dataset recorded through experience, experimentation, or artificial data
generation. We write this as $\\{(v_{m}^{(k)},E_{i}^{(k)})\\}_{k=1,\cdots,N}$,
where $v_{m}^{(k)}\ (m=0,\cdots,M)$ is real-valued input vectors of dimension
$M+1$, and $E_{i}^{(k)}\ (i=0,\cdots,i_{max}-1)$ is real-valued output vectors
of dimension $i_{max}$. The index $k$ is the sample’s label, which can take
values from $1$ to $N$. The dataset size $N$ should be as large as possible to
train NNs.

An NN can be written as $f_{i}(v_{m}^{(k)};\\{w_{\alpha}\\})$. In addition to
taking $v_{m}^{(k)}$ as an argument, it has real valued parameter set
$\\{w_{\alpha}\\}$ with index $\alpha$, which can take values from 1 to the
total number of the NN’s parameters. The majority of $\\{w_{\alpha}\\}$ is the
weights and bias of the linear combination used to create the composite
function, but other parameters can also be included. Training an NN means
optimizing $\\{w_{\alpha}\\}$ so that $f_{i}(v_{m}^{(k)};\\{w_{\alpha}\\})$ is
as close to $E_{i}^{(k)}$ as possible for most of $v_{m}^{(k)}$.

The flow of this optimization is depicted in Fig.1. First, when the input data
$v_{m}^{(k)}$ is fed to the NN, the output
$f_{i}(v_{m}^{(k)};\\{w_{\alpha}\\})$ is obtained. Then, the loss resulting
from the error between the output and the ground truth data, $E_{i}^{(k)}$, is
calculated. In addition, the first-order derivative of the loss function to
each parameter, $w_{\alpha}$, can be calculated almost simultaneously using a
method called back-propagation. The value of $w_{\alpha}$ is updated to
decrease the loss using the derivative. Repeating this process until the loss
does not significantly decrease is called NN training.

In actual training, instead of adding up all samples in the loss function, a
few samples are selected, which is called mini-batch learning. Mini-batch
learning saves computer memory, and it is empirically known that mini-batch
learning is faster than taking all data. Furthermore, the final results are
often better when the batch size is smaller. At least, we do not know the
general theory for determining the optimal batch size and have empirically
fixed the size to 100 by trial and error in this paper.

Figure 1: Training a neural network(NN). An NN is a function, which takes a
multidimensional real-valued input and produces a multidimensional real-valued
output. Training an NN requires many input and ground truth output datasets.
NN learning is the optimization of real valued parameter set,
$\\{w_{\alpha}\\}$, to reduce the loss function, $MSE(\\{w_{\alpha}\\})$,
linked to the difference between the correct data and the prediction by the
output of the NN. In this optimization, the back-propagation method, which can
quickly calculate the derivative of the loss function to parameters, is
indispensable. In this paper, the input is 1D potentials’ values on mesh
points, and the output is the corresponding energy spectra. We describe the
architecture of our NN and the specific method for optimizing its parameters
in the appendix.

Next, we explain how we apply NN learning to a physical problem in this paper.
We deal with quantum mechanics in one dimension. A particle is confined to the
interval $-1<x<1$ by infinitely high walls. Then, we generate $N$ various
potentials between $-1<x<1$ and take their values at equally spaced points as
input data $v_{m}^{(k)}$. For the output data $E_{i}^{(k)}$, we take ten of
the smallest energy eigenvalues corresponding to the potentials
($i_{max}=10$). The following subsection explains how to generate the
potentials and calculate the corresponding eigenvalues. In this study,
$E_{i}^{(k)}$ is obtained by solving the Schrödinger equation numerically.
However, because the primary goal of this paper is to see if NN can be trained
to understand the physical laws using observed physical quantities,
$E_{i}^{(k)}$ is positioned as simulated experimental data.

We use a 1D CNN. In addition, we also use the residual connection and the
self-attention layer, which are effective when used in conjunction with CNNs
in the natural language processing and image generation fieldshe2016deep ;
zhang2019self . We adopt the mean squared error(MSE) between the ground truth
energy eigenvalues and the outputs from the NN for the loss function. We
describe the architecture of our NN and the specific method for optimizing its
parameters in the appendix.

###  2.2 How to prepare the datasets

We take the energy spectrum of a particle confined to an interval in one
dimension as our output data. Hence, we consider the following potential,

$\displaystyle V^{(k)}(x)=\left\\{\begin{array}[]{ll}v^{(k)}(x)&(|x|\leq 1)\\\
\infty&(|x|>1)\end{array}\right.$ (3)

where $v^{(k)}(x)$ are arbitrary continuous functions and k takes values from
1 to N, which is the size of datasets.

Next, we use the kernel method to generate the potential functions,
$v^{(k)}(x)$. Note that we use the kernel method only to generate the
potential functions, not to perform regression using the Gaussian process. We
set mesh points on $-1<x<1$ with $\Delta x\ (=2/M)$ width as $x_{l}=(2l-M)/M\
(l=0,1,2,\cdots M)$ to use the kernel method. Let the kernel function be
$k(x,x^{{}^{\prime}})$, and then, under our discretization, the kernel matrix
is the $(M+1)\times(M+1)$ matrix defined by $K_{lm}=k(x_{l},x_{m})$. By using
this kernel matrix as the covariance matrix, we can probabilistically generate
the potential function, as follows:

${\bm{v}}\sim{\mathcal{N}}({\bf 0},{\bf K})$ (4)

where $\mathcal{N}$ represents the normal distribution function, with the mean
being ${\bf 0}$ and the covariance matrix being ${\bf K}$. Furthermore,
${\bm{v}}$ is an $(M+1)$-dimensional vector, whose components are regarded as
potential values at $x_{l}$. We obtain ${v}^{(k)}(x_{l})$ by performing this
potential generation for $N$ times.

This study uses four kernel functions to generate the training data: Gaussian
kernel $k^{(RBF)}(x,x^{{}^{\prime}})$, Matérn5 kernel
$k^{(M5)}(x,x^{{}^{\prime}})$, Matérn3 kernel $k^{(M3)}(x,x^{{}^{\prime}})$,
and exponential kernel $k^{(Exp)}(x,x^{{}^{\prime}})$. These functions are
defined, as follows:

$\displaystyle k^{(RBF)}(x,x^{{}^{\prime}})$ $\displaystyle=$
$\displaystyle\sigma^{2}\exp{\left(-\frac{r^{2}}{2L^{2}}\right)}$ (5)
$\displaystyle k^{(M5)}(x,x^{{}^{\prime}})$ $\displaystyle=$
$\displaystyle\sigma^{2}\left(1+\frac{\sqrt{5}r}{L}+\frac{5r^{2}}{3L^{2}}\right)\exp{\left(-\frac{\sqrt{5}r}{L}\right)}$
(6) $\displaystyle k^{(M3)}(x,x^{{}^{\prime}})$ $\displaystyle=$
$\displaystyle\sigma^{2}\left(1+\frac{\sqrt{3}r}{L}\right)\exp{\left(-\frac{\sqrt{3}r}{L}\right)}$
(7) $\displaystyle k^{(Exp)}(x,x^{{}^{\prime}})$ $\displaystyle=$
$\displaystyle\sigma^{2}\exp{\left(-\frac{r}{L}\right)}$ (8)

where we defined $r=|x-x^{{}^{\prime}}|$. These functions have two parameters:
$\sigma$, which controls the magnitude of the generated potential, and $L$,
which controls the correlation distance of the potential. These kernel
functions can be written concisely using the gamma function $\Gamma(\nu)$, and
the Bessel function of the second kind $K_{\nu}$, as follows:

$\displaystyle k_{\nu}(x,x^{{}^{\prime}})$ $\displaystyle=$
$\displaystyle\frac{2^{1-\nu}\sigma^{2}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}r}{L}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}r}{L}\right)$
(9)

The Gaussian, Matérn5, Matérn3, and the exponential kernels correspond to the
cases of $\nu=\infty,\ \frac{5}{2},\ \frac{3}{2},\ \frac{1}{2}$, respectively.
Kernels with larger $\nu$ generate smoother potential functions, and the
generated potentials are $[\nu]$ times differentiable. Here $[\nu]$ is the
largest integer not greater than $\nu$. Figure 2 shows the potentials
generated from each kernel. The four panels in this figure show that the
smoothness of the potentials changes as $\nu$ increases.

Figure 2: Potentials generated using the kernel method. The kernel parameters
were taken to be the correlation distance $L=0.2$ and the amplitude
$\sigma=100$ for all kernels. The smoothness varies depending on the kernel
used. (a) Potentials generated from the Gaussian kernel, Eq.(5) or Eq.(9) at
$\nu=\infty$. These are infinite-order differentiable. (b) Potentials
generated from the Matérn5 kernel, Eq.(6) or Eq.(9) at $\nu=\frac{5}{2}$.
These are second-order differentiable. (c) Potentials generated from the
Matérn3 kernel, Eq.(7) or Eq.(9) at $\nu=\frac{3}{2}$. These are first-order
differentiable. (d) Potentials generated from the exponential kernels, Eq.(8)
or Eq.(9) at $\nu=\frac{1}{2}$. These are not differentiable.

Now that we have generated the potentials, which are the input data needed to
train the NN, we need to find the energy spectrum of those potentials, the
output data. For this purpose, we solve the time-independent 1D Schrödinger
equation,

$\left\\{-\frac{1}{2}\frac{d^{2}}{dx^{2}}+V^{(k)}(x)\right\\}\psi^{(k)}(x)=E^{(k)}\psi^{(k)}(x),$
(10)

where $\psi^{(k)}(x)$ is the wave function and we use Planck’s constant
$h=2\pi$, and the particle’s mass $m=1$. Since we are considering the
potentials in Eq.(3), the boundary conditions for the wave function are
$\psi^{(k)}(\pm 1)=0$. In the case of the square potential, $v^{(k)}(x)=0$,
the energy eigenvalues are,

$E_{i}^{(k)}=\frac{\pi^{2}(i+1)^{2}}{8}\ \ \ (i=0,1,2,\cdots),$ (11)

and the corresponding wave functions are,

$\psi^{(k)}_{i}(x)=\sin{\left[\frac{\pi}{2}(i+1)(x+1)\right]}.$ (12)

Since we use the ten lowest energy eigenvalues as the ground truth output
data, it should have a rich structure. Referring to Eq.(11), we take
$\sigma=100$ in Eqs.(5)-(8) throughout this paper expect the dataset for
obtaining Fig.5. Figure 3 shows that this setting works in the actual
numerical experiment results.

We use the matrix methodizaacComputational2019a  to solve this eigenvalue
problem for general $v^{(k)}(x)$. This method is appropriate for our setting
because it enables us to find multiple energy eigenvalues and corresponding
wave functions simultaneously. Furthermore, since we have already discretized
the problem to generate the values of the potentials on the mesh points, we
can directly use this method.

Applying the matrix method to Eq.(10) is straightforward: rewriting Eq.(10) by
replacing the second-order derivative with a difference results in the
following matrix eigenvalue problem,

$\displaystyle H^{(k)}\Psi^{(k)}=E^{(k)}\Psi^{(k)}$ (13)

which we defined, as follows ($v_{m}^{(k)}\equiv v^{(k)}(x_{m})$),

$\displaystyle H^{(k)}=-\frac{1}{2(\Delta
x)^{2}}\left(\begin{array}[]{cccccc}-2&1&0&\ldots&0&0\\\ 1&-2&1&\ldots&0&0\\\
0&1&-2&\ldots&0&0\\\ \vdots&\vdots&\vdots&\ddots&\ddots&\vdots\\\
0&0&0&\ddots&-2&1\\\
0&0&0&\ldots&1&-2\end{array}\right)+\left(\begin{array}[]{cccccc}v^{(k)}_{1}&0&0&\ldots&0&0\\\
0&v^{(k)}_{2}&0&\ldots&0&0\\\ 0&0&v^{(k)}_{3}&\ldots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\ddots&\vdots\\\ 0&0&0&\ddots&v^{(k)}_{M-2}&0\\\
0&0&0&\ldots&0&v^{(k)}_{M-1}\end{array}\right)$ (26)

Here, we have used the boundary conditions
$\psi^{(k)}(x_{0})=\psi^{(k)}(x_{M})=0$, so the matrix $H^{(k)}$ is an
$(M-1)\times(M-1)$ dimensional matrix and the eigenvector $\Psi^{(k)}$ is
$(M-1)$-dimensional vector. This eigenvalue problem is easy to solve because
$H^{(k)}$ is a tridiagonal real symmetric matrix. We denote the eigenvectors
and corresponding eigenvalues as $E_{i}^{(k)}$ and

$\displaystyle\Psi^{(k)}_{i}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}\psi^{(k)}_{i}(x_{1})\\\
\psi^{(k)}_{i}(x_{2})\\\ \psi^{(k)}_{i}(x_{3})\\\ \vdots\\\
\psi^{(k)}_{i}(x_{M-2})\\\ \psi^{(k)}_{i}(x_{M-1})\end{array}\right),\hskip
28.45274pt\Psi^{(k)T}_{i}\Psi^{(k)}_{i}=1$ (33)

, respectively. Note that this normalization of the wave function is natural
as a normalization for the matrix eigenvalue problem but deviates by
$\sqrt{\Delta x}$ from the ordinary normalization condition of the wave
function, $\int_{-1}^{1}|\Psi^{(k)}(x)|^{2}dx=1$. There are M-1 eigenvalues,
labeling them as $E_{i}^{(k)}<E_{j}^{(k)}$ for $i<j$. Since this is 1D quantum
mechanics, the energy eigenvalues have no degeneracy. We adopt ten
eigenvalues($i=0,1,2\cdots,9$) for the ground truth output data of the NN.

We have depicted one potential generated from the Gaussian kernel and the ten
lowest energy eigenvalues of the potential in Fig.3(c). We note that we have
confirmed that these are the ten lowest energy eigenvalues by counting the
nodes of the corresponding wave functions.

This subsection explains how to create a set of ($M+1$)-dimensional potential
vectors, $v_{m}^{(k)}$, as input data and 10-dimensional energy eigenvalues,
$E_{i}^{(k)}$, as output data for training an NN. Here we used the Schrödinger
equation to prepare these datasets, but we did not feed the equation to the
NN. We only feed the NN with the set of potentials $v_{m}^{(k)}$ and energy
eigenvalues $E_{i}^{(k)}$ as if we observe them in a physical experiment. Of
course, the Schrödinger equation is not required if the datasets are prepared
by conducting experiments to measure energy eigenvalues while varying the
potential in various ways. The primary goal of this paper is to demonstrate
that NN understands physical laws using only the dataset
$\\{(v_{m}^{(k)},E_{i}^{(k)})\\}_{k=1,\cdots,N}$ principally observed in the
experiment.

###  2.3 Learning setup and results

We generated three training datasets, N1, N2, and N5 with the sizes, $1\times
10^{5}$, $2\times 10^{5}$, and $5\times 10^{5}$, respectively, for all four
kernels of Eq.(5)-(8) to train the NN. The reason why we constructed datasets
of different sizes was to examine the degree of overfitting. In addition, we
generated different datasets of sizes $1\times 10^{4}$ for all kernels for
validation. The parameters of the kernels were $\sigma=100$ and $L=0.2$. Since
the batch size was fixed at 100, one epoch of training would update the NN
parameters $\\{w_{\alpha}\\}$, $1\times 10^{3}$, $2\times 10^{3}$, and
$5\times 10^{3}$ times for the datasets N1, N2, and N5, respectively. To keep
the total number of parameter updates constant and compare the learning status
fairly, we trained $100$, $50$, and $20$ epochs for the datasets N1, N2, and
N5, respectively. This setup resulted in the parameter set,
$\\{w_{\alpha}\\}$, being updated $1\times 10^{5}$ times in total for any
dataset. The loss function is the MSE. The design of the NN and the precise
optimization method of the parameter set, $\\{w_{\alpha}\\}$, are included in
the appendix.

Figure 3(a) shows the learning progress. Here, we defined the normalized mean
squared error(NMSE) to evaluate the error as a dimensionless quantity, as
follows,

$NMSE(\\{w_{\alpha}\\})=\frac{\sum_{k=1}^{\cal
N}\sum_{i=0}^{i_{max}-1}(E_{i}^{(k)}-f_{i}(v_{m}^{(k)};\\{w_{\alpha}\\}))^{2}}{\sum_{k=1}^{\cal
N}\sum_{i=0}^{i_{max}-1}(E_{i}^{(k)})^{2}}.$ (34)

Here, ${\cal N}$ is the size of the used dataset. In addition, we defined the
progress rate as the ratio of the performed update times of $\\{w_{\alpha}\\}$
to the total one, $1\times 10^{5}$. This figure displays the NN’s learning
process using datasets generated from the Gaussian kernel. The dashed line
represents the error for the training data, and the solid line represents the
error for the validation data. This figure shows that overfitting is
suppressed as the size of the dataset increases, and the effect of overfitting
is minimal at N5. The NMSE is less than $10^{-4}$ at the end of the training,
and the error in the energy eigenvalue is less than $1\%$ on average for N5.
We have confirmed the datasets generated from other kernels show the same
trend even though some quantitative differences exist.

Figure 3(b) shows the difference in the NMSE after the completion of training
depending on the kernel that generates the training data. In the figure legend
hereafter, Gaussian, Mattern 5, Mattern 3, and the exponential kernel are
abbreviated as G, M5, M3, and Exp, respectively. Validation data were
generated from the same kernel that generated the training data. This figure
shows that the error is minor for the dataset with a larger size. In addition,
the smoother the potential, the smaller the error.

Figure 3(c) shows the prediction of the energy eigenvalues of one potential
generated from the Gaussian kernel by an NN trained on the N5 dataset
generated from the same kernel. This figure shows that the deviation of the
prediction of the energy eigenvalue is small.

Figure 3: Training NNs. (a) NMSE in Eq.(34) is plotted against the training
progress. The dataset used for training varies between $1\times 10^{5}$(N1),
$2\times 10^{5}$(N2), and $5\times 10^{5}$(N5) to study the difference in
learning progress. The total number of parameter updates is $1\times 10^{5}$
in all cases. The dashed line represents the error for the training data, and
the solid line represents the error for the validation data. For validation,
we used data generated from the same kernel that generated the data used for
training. This figure shows that the larger the dataset, the less overfitting
there is. (b) The final epoch errors are displayed for each kernel that
generates data. This figure shows that smooth data are easier to learn. (c)
For one sample potential generated from the Gaussian kernel, we compared the
energy eigenvalues predicted by the NN with the ground truth values calculated
from the matrix diagonalization in Eq.(13). As shown in (b), for the data
generated from the Gaussian kernel, the average error of the energy
eigenvalues is less than $1\%$, which can barely be seen in this figure.

##  3 Understanding Physics with Neural Networks

In Sec.2, we showed that NNs can now predict energy eigenvalues accurately,
but it is still unclear whether they understand the physics described by
Eq.(10). As described in the Appendix, our NN has $4.9\times 10^{6}$ real
number parameters, $\\{w_{\alpha}\\}$, and the real number size of the output
of the largest dataset, $\\{E_{i}^{(k)}\\}_{k=1,\cdots,N}$, is $5\times
10^{6}$. This fact may lead us to suspect that the NN does not understand
physics but memorizes the output.

In this section, which describes the main results of this paper, we will show
in four aspects that NNs do understand physics and that they are helpful for
the study of physics. First, in Sec.3.1, we show that NNs can predict energy
eigenvalues even for potentials of a different shape than the training
dataset. Next, in Sec.3.2, we successfully extract the eigenvectors’
magnitudes not used for the training from the trained NN. We then show that
NNs can predict unknown physical phenomena in Sec.3.3. Finally, we treat
potentials with matter effects unknown to the experimenters to demonstrate the
broad applicability of the method in Sec.3.4.

These results show that NNs can learn the physical laws from experimental
data, predict the results of experiments under conditions different from those
used for training, and predict physical quantities of types not provided
during training. Because NNs understand physics in a different way than
humans, they will be a powerful tool for advancing physics by complementing
the human way of understanding.

###  3.1 Expand the scope of applicable data

In Sec.2.3, the training performance of the NN was verified with validation
data generated from the same kernel that generated the training data. However,
if the NN understands the physics governed by Eq.(10), it should predict the
energy eigenvalues of a potential generated from a different kernel than the
one that generated the training data. Therefore, this subsection shows how the
NN predicts energy eigenvalues for the validation data generated from a
different kernel than the training data generated.

Figure 4(a) shows the training progress using the N5 dataset generated from
the Gaussian kernel. This figure shows that the error for the validation data
generated from the Gaussian kernel is the smallest, but the prediction is not
broken, even for datasets generated from other kernels.

We show the error at the final epoch in Fig.4(b) by the kernel that generated
the training data. This figure shows that NNs trained with data such as the
one in Fig.2(d), generated from the exponential kernel, show promising
results, regardless of the dataset type. On the other hand, the NN trained on
the dataset generated from the Gaussian kernel is not good at predicting
potentials that are not smooth.

We used the NN trained on the N5 dataset generated from the exponential kernel
in Fig.4(c) to predict the energy eigenvalues of the same potential as in
Fig.3(c). This figure demonstrates that, although there is a slightly
significant deviation in the energy eigenvalue of the seventh excited state,
the overall prediction is accurate. Though the NN has seen only the jagged
potential in Fig. 2(d), it correctly predicted the energy eigenvalues for a
smooth potential generated from Gaussian kernel. This result indicates that it
understood the physical laws governing this system.

Figure 4: (a) NMSE versus learning progress when using data generated from the
Gaussian kernel for training dataset. The predictions for validation datasets
generated from the other kernels are also good, and it seems that predictions
for smooth potentials are easier than rough ones. (b) The error at the final
epoch is shown for each N5 dataset used for training. The NN trained with the
exponential kernel shows a good prediction for any validation dataset. (c) The
NN trained by the N5 dataset generated by the exponential kernel was used to
predict energy eigenvalues for the same potential as in Fig.3(c). The energy
eigenvalues of the seventh excited state are slightly different, but
considering that the prediction is for the potential generated by the the
different kernel from the training one, the prediction is excellent.

In addition, to test whether the learning works even with highly undulating
potentials, we took $\sigma=1000$ and generated a new dataset, N5’, with
$5\times 10^{5}$ potentials for each kernel. The same figures as in Fig. 4 are
shown in Fig. 5 for this new data set. As can be seen from Fig. 5(c), the ten
lower energy eigenvalues taken are well below the highest point of the
potentials. Figure 5(a) and (b) show that although the learning is well
advanced, the error is larger than in the case of N5. This error increase is
most likely due to the differences in larger energy eigenvalues; in N5, some
of the larger energy eigenvalues were approaching the obvious ones determined
only by the boundary conditions, but in N5’, all ten energy eigenvalues are
nontrivial.

Figure 5: The same figures as in Fig. 4 for the dataset, N5’($\sigma=1000$).

Following that, we quantitatively examine the relationship between the
similarity and loss between training and validation potentials. The potentials
generated by different kernels appear to have different shapes in Fig.2(a)
through 2(d). However, some of the $5\times 10^{5}$ potentials in N5 may be
close by chance even for different kernels, and the NNs may simply remember
them.

To discuss this, we define a measure of similarity between a training dataset,
$\\{v^{(k^{{}^{\prime}})}_{m}\\}_{k^{{}^{\prime}}=1,\cdots,N}$ , and a
validation potential, $u^{(k)}_{m}$, as

$\displaystyle
d(u^{(k)}_{m},\\{v^{(k^{{}^{\prime}})}_{m}\\}_{k^{{}^{\prime}}=1,\cdots,N})$
$\displaystyle\equiv$
$\displaystyle\min_{k^{{}^{\prime}}}\min_{c_{k^{{}^{\prime}}}}\sqrt{\frac{1}{M+1}\sum_{m=0}^{M}(u^{(k)}_{m}-v^{(k^{{}^{\prime}})}_{m}-c_{k^{{}^{\prime}}})^{2}}$
(35) $\displaystyle=$
$\displaystyle\min_{k^{{}^{\prime}}}\sqrt{\frac{1}{M+1}\sum_{m=0}^{M}(u^{(k)}_{m}-v^{(k^{{}^{\prime}})}_{m})^{2}-{\bar{\Delta}}_{k^{{}^{\prime}}}^{2}}$
(36)

, where $c_{k^{{}^{\prime}}}$ is taken at

$\displaystyle{\bar{\Delta}}_{k^{{}^{\prime}}}=\frac{1}{M+1}\sum_{m=0}^{M}(u^{(k)}_{m}-v^{(k^{{}^{\prime}})}_{m})$
(37)

for the minimization.

The smaller this number, $d$, the closer the validation potential is to the
training dataset. In this definition, we find the potential with the smallest
$L^{2}$ norm to the validation potential in the training dataset and use the
norm as our measure. However, because we are concerned with the shape of the
potentials, we add a constant, $c_{k^{{}^{\prime}}}$, to the training dataset
potential to minimize the $L^{2}$. For the validation dataset, $10^{4}$
potentials were generated from the Gaussian kernel, of which the one that
minimizes $d$ is shown in Fig 6(a). We used the N5 dataset generated from the
Gaussian kernel for the training dataset. The potential that gives the minimum
$d$ value in the training dataset,
$v^{(k^{{}^{\prime}})}_{m}+{\bar{\Delta}}_{k^{{}^{\prime}}}$ is also shown. In
other words, this figure shows our measure’s closest pair of the $10^{4}$
validation potentials and the $5\times 10^{5}$ training potentials both
generated from Gaussian kernel. We can see from this figure that if some of
the validation and training data are similar in shape, we can successfully
find them with our measure. Similarly, Fig.6(b) depicts the best match between
the $10^{4}$ Gaussian kernel validation potentials and the N5 training dataset
generated by the exponential kernel. However, it can be seen that their
overall undulations are similar, but their shapes are not. The closest match
between the validation and N5 training potentials generated by the exponential
kernel is shown in Fig.6(c). This figure also demonstrates that while their
overall undulations are similar, their shapes are not.

To evaluate the similarity between the validation datasets and the training
datasets, the following metric is defined,

$\displaystyle
D(\\{u^{(k)}_{m}\\}_{k=1,\cdots,N_{v}},\\{v^{(k^{{}^{\prime}})}_{m}\\}_{k^{{}^{\prime}}=1,\cdots,N})=\frac{1}{N_{v}}\sum_{k=1}^{N_{v}}d(u^{(k)}_{m},\\{v^{(k^{{}^{\prime}})}_{m}\\}_{k^{{}^{\prime}}=1,\cdots,N}).$
(38)

Here, we take $N_{v}=10^{4}$ and $N=5\times 10^{5}$.

Figure 6(d) depicts the similarity of each dataset using this metric. This
figure shows that the training dataset generated by the exponential kernel is
not even similar to the validation dataset generated by the exponential kernel
itself. However, this diversity of potentials, prevents over-fitting and is
most likely one of the reasons why the training dataset generated by the
exponential kernel can train the NN to make good predictions for any of the
validation data. Next, the relationship between the similarity and the loss is
shown in Fig.7. Figure 7(a) redraws Fig.4(b) with $D$ on the axis. Here, the
annotations show the used training datasets. This graph demonstrates that a
lower $D$ does not imply a lower loss. Instead, the validation datasets
generated by Gaussian, Matérn 5, and Matérn 3 kernels have the lowest loss in
the training datasets generated by the Matérn 5, Matérn 3, and exponential
kernels, all of which have larger $D$ than the original kernels.

Figure 6:  (a) The closest potentials between validation and training dataset
both generated from Gaussian kernel. (b) The potentials with the highest
similarity between the validation and training datasets generated by the
Gaussian and exponential kernels, respectively. (c) The exponential kernel
generated the closest potentials between the validation and training datasets.
(d) The resemblance of the validation and training datasets.

Figure 7(b) shows the scatter plots of MSE for the ten eigenvalues value and
$d$ for the individual validation potentials. The training dataset is N5,
which was generated by the exponential kernel. This graph also demonstrates
that a lower d does not always imply a lower loss. These findings suggest that
while NNs do not remember potential shapes and the corresponding energy
eigenvalues, they do understand the law, Eq.(10).

Figure 7:  (a) NMSE in Fig 4(b) are shown on the metric, $D$. The annotations
denote the datasets used for training. This figure shows that the validation
dataset can have a lower loss when trained by a different kernel with a larger
$D$. (b) Scatter plots of validation potentials on MSE for the ten energy
eigenvalues and the measure, $d$. For training, the N5 exponential kernel
dataset is used. This graph demonstrates that a lower $d$ does not always
imply a lower loss.

Finally, we look at how the residual connection and the self-attention
mechanism affect NN predictions. For this examination, we omit portions of the
NN layers and observe how the loss changes. Figure 17 summarizes our NN
architecture. To begin, to see the effect of deepening NN via the residual
connection, we examine how the loss changes as the depth of the NN changes.

As shown in Fig.17, our NN has five blocks with residual connections, which we
refer to as ”Res-SA blocks.” To see the difference when we reduce these blocks
from the higher level, we define the following quantity,

$\displaystyle R_{D}=\frac{\text{NMSE( Number of Res-SA blocks
})}{\text{NMSE(1)}}$ (39)

, where NMSE(1) is NMSE resulting from the NN design that directly leads from
the initial Res-SA block to the Flatten layer. The loss rapidly decreases as
the NN is deepened using the residual connection for any training dataset, as
shown in the left panel of Fig.8. This change shows that it is essential to
capture the overall shape of the potential by broadening the field of view of
neurons.

Next, we see the effect of the self-attention mechanism: the self-attention
layer is in the Res-SA block in Fig.17, and the contents of the Res-SA block
are shown in Fig 18. The ratio of the NMSE evaluated from the NN with and
without the self-attention layer is defined as follows,

$\displaystyle R_{SA}\equiv\frac{\text{NMSE( With self-attention
mechanism})}{\text{NMSE( Without self-attention mechanism )}}.$ (40)

The quantity, $R_{SA}$, represents how much the loss is reduced by the self-
attention mechanism. As shown in the right panel of Fig.8, the self-attention
mechanism reduces losses for all training and validation dataset combinations.
It is especially effective when smooth potentials are used for validation. We
examine the attention maps to determine what the self-attention mechanism is
focusing on. Figure 9 depicts the first and second self-attention maps for the
potential depicted in Fig.4(c). We used the NN trained with N5 dataset
generated from the exponential kernel. To make the figure easier to
understand, we aligned the potential shapes on the left and top sides of the
self-attention maps. This self-attention map represents the importance of the
relationship between two points on the potential learned by the NN. The figure
on the left is the self-attention map placed in the layer closest to the input
data and shows that the NN perceived the relationship between two minima of
the potential as essential. The figure on the right is a self-attention map
placed in the second self-attention layer. The NN perceived the relationship
of maxima locations as meaningful in addition to the relationship of minima
locations. Thus, attention to the relationship between distant feature points
seems to reduce losses. This observation also shows that the NN learns the
necessary features by itself through learning.

Figure 8:  The left figure shows the mean loss ratio, $R_{D}$, defined in
Eq.(39) on the number of used Res-SA blocks in the NN. It can be seen that by
deepening the NN layers, the loss can be significantly reduced. The right
figure depicts the loss ratio between NMSE with and without the self-attention
mechanism, $R_{SA}$, defined in Eq.(40) for all training and validation
dataset combinations. The self-attention mechanism is especially useful for
smooth validation data.  Figure 9: Attention maps for the potential in
Fig.3(c). The figure on the left shows the attention map of the first self-
attention layer in the NN. The top and left sides of the figure show the
corresponding potential shape. We can see which potential points the NN
considers relevant for each reference point on the left side. The NN is
intensely interested in the relationship between the minima of the potentials.
The figure on the right shows the attention map of the NN’s second self-
attention layer. In this attention map, we can see that the NN is interested
in the relationship between minima and maxima.

The findings in this subsection show that NNs can help physics experiments
produce more results with fewer practical resources. For example, even when
only materials with rough surfaces are available, it would be possible to
train an NN with the results of experiments using those materials and then use
the trained NN to predict the results of experiments using materials with
smooth surfaces. In addition, for a physical system with wide parameter space,
we can experiment with a part of it and use the observed results to train the
NN. Then, we may raise candidates for parameter regions where interesting
physical phenomena would be observed by using the trained NN to search for
unknown parameter regions. A schematic diagram is shown in Fig.16.

###  3.2 Elicit information that is not being fed to the neural network

When we trained the NN, we only gave it the potential and the corresponding
energy eigenvalues. However, if the NN understands the physics described by
Eq.(10) in this process, we may extract a physical quantity other than the
energy eigenvalues. This subsection attempts to find the probability of the
particle’s existence, which is not used in training.

The function $f_{0}(v_{m};\\{w_{\alpha}\\})$ of the learned NN predicts the
energy eigenvalues of the ground state using values at each point of the
potential. If $f_{0}(v_{m};\\{w_{\alpha}\\})$ sensitively changes when the
value of the potential at a point, $v_{m}$, is changed, there is a high
probability that there is a particle in the vicinity. Therefore, the first-
order derivative of $f_{0}(v_{m};\\{w_{\alpha}\\})$ to the potential value is
considered the existence probability of the particle for the ground state. The
same can be said for excited states.

This can also be explained, as follows: $f_{i}(v_{m};\\{w_{\alpha}\\})$ should
have been learned to approximate $E_{i}$’s representation. From the
discretized quantum mechanics, Eq.(13)-(33), the first-order derivative of the
energy eigenvalue to the potential value is the square of the eigenvector, as
shown below,

$\displaystyle\frac{dE_{i}}{dv_{m}}=\frac{d}{dv_{m}}\left(\Psi^{T}_{i}H\Psi_{i}\right)=\Psi^{T}_{i}\left(\frac{d}{dv_{m}}H\right)\Psi_{i}+E_{i}\frac{d}{dv_{m}}\left(\Psi^{T}_{i}\Psi_{i}\right)=\psi^{2}_{i}(x_{m}).$
(41)

This should result in the following relationship if the learning is
sufficient:

$\displaystyle\left.\frac{\partial f_{i}(v_{m};\\{w_{\alpha}\\})}{\partial
v_{m}}\right|_{\\{w_{\alpha}\\}={\rm learned\
values}}\simeq\psi^{2}_{i}(x_{m}).$ (42)

As a property of NNs, the first-order derivative of the argument can be easily
obtained using back-propagation. When training an NN, only the derivative to
the parameters $\\{w_{\alpha}\\}$ is used, not the derivative to the input
data $v_{m}$, but, here, we reverse the position and fix the parameters
$\\{w_{\alpha}\\}$ on the trained values and obtain the derivative to the
input data $v_{m}$.

The left and right sides of Eq.(42) are displayed in Fig.10. The NN was
trained with N5 dataset generated from the exponential kernel. The derivative
values were not normalized again, and the values were displayed as they were.
Although there are some unnatural parts, such as negative values, the general
agreement regarding the location of the peaks and the number of nodes is
excellent. This result of obtaining the probability distribution of the
particle’s existence is surprising, considering that the used NN is only given
the jagged potential such as Fig.2(d) and the corresponding energy
eigenvalues.

This subsection shows that the trained NN provides observable quantity other
than the output data used during training. As a technical aspect, we also
found that differentiation on input data using back-propagation may be helpful
in physics. The same technic called saliency map has been used in object
detectionsimonyan2013deep . For example, by differentiating the dog
probability in the outputs concerning the input image, one can get an idea of
where a dog is in an image. It is very similar to the calculation in this
subsection in that it focuses on the response to the input perturbation.

Figure 10:  A neural network estimated the probability distribution of
particle existence. For the ground state and the first excited state, both
sides of Eq.(42), the eigenvector squared (ES) and the predicted energy
derivative (ED), are shown on the left. The same quantities are displayed in
the figure on the right for an even higher excited state.

###  3.3 Predicting Physical Phenomena with Neural Networks

We reproduce a well-known physical phenomenon here, but the same procedure
would be able to predict unknown phenomena. 1D quantum mechanics has no
degeneracy of energy eigenvalues, as can be seen from the fact that even the
symmetric double-well potential has no degeneracy of energy eigenvalues. Two
energy eigenvalues do not intersect, and their eigenstates cross over when the
potential changes slowly and continuously. This phenomenon plays an essential
role in the Mikheyev-Smirnov-Wolfenstein effect of the solar neutrino
problemwolfensteinNeutrino1978 ; mikheyev1986resonant ; parkeNonadiabatic1986
; bethePossible1986 .

The following setup is used to reproduce this phenomenon. First, we generate
two potentials, $v^{(1)}(x)$ and $v^{(2)}(x)$, using the Gaussian kernel, and
subtract the energy eigenvalues of the respective ground states to create
potentials $\tilde{v}^{(1)}(x)$ and $\tilde{v}^{(2)}(x)$ adjusted so that the
energy eigenvalues of the ground states become zero. We then introduce a real
mixing parameter $\lambda\in[0,1]$ and create a potential $v(x)$ that varies
continuously from the linear combination, as follows:

$v(x)=(1-\lambda)\tilde{v}^{(1)}(x)+\lambda\tilde{v}^{(2)}(x).$ (43)

This continuously varying potential is shown in Fig.11(a).

Figure 11: Continuously varying potential and corresponding energy eigenvalue
changes. (a) Two Gaussian kernel-generated potentials adjusted to zero ground
state energy were linearly combined with a mixing constant $\lambda$. (b) An
NN trained with the dataset generated from the exponential kernel was used to
predict the evolution of the energy eigenvalues of the potentials generated in
(a). The energy eigenvalues are changing without degenerating. Figure 12:
Crossing over between the ground state(GS) and the first excited state(1st)
for the changing potential energy(PE). (a) The probability of the existence of
a particle at $\lambda=0.75$. (b) Existence probability of the particle at
$\lambda=0.79$.(c) Existence probability of the particle at $\lambda=0.83$. We
used the left side of Eq.(42) for this estimation.

Figure 11(b) shows the predicted energy eigenvalues from the NN as a function
of $\lambda$. The energy eigenvalues are changing without crossing. The NN
trained with the N5 dataset generated from the exponential kernel was used.
The energy eigenvalues of the ground state and the energy eigenvalue of the
first excited state are comparable around $\lambda=0.8$. We investigate what
is happening in the vicinity of this value using the method of Sec.3.2. Figure
12(a) shows the existence probability distribution of the particle at
$\lambda=0.75$. The ground state is on the left and the first excited state is
on the right. Figure 12(b) shows the existence probability distribution of the
particle at $\lambda=0.79$. The existence probability distributions of the
ground state and the first excited state overlap. Figure 12(c) shows the
existence probability distribution of the particle at $\lambda=0.83$. The
ground state is on the right, and the first excited state is on the left. In
this way, we can see how the eigenstates continuously cross over.

This subsection showed that the trained NN could reproduce the well-known
physical phenomenon. Therefore, it is also expected to help humans to find new
physical phenomena.

###  3.4 Wide range of applications and experiments necessary for learning

So far, we’ve looked at the case where the experimenter is fully aware of the
potential. We also know that the physical system was governed by the
Schrödinger equation, so we can solve it numerically and obtain the energy
eigenvalues in this case. In this subsection, we consider simple toy models to
demonstrate that the NN method has a broader applicability range.

Consider the case of an experiment with a particle in the matter. The
experimenters can make tunable potentials, $V_{e}(x)$, to the particle by
employing an electric field. However, the potentials are expected to be
$V_{total}(x)$ in matter due to unknown matter effects to the experimenters.
Because the experimenters do not know how matter affects the potential in this
case, the Schrödinger equation cannot be used to determine the energy
eigenvalues. Even in this case, if they can observe the energy eigenvalues for
various $V_{e}(x)$ in the experiment, they can train NN with the results and
predict the energy eigenvalues using the trained NN without knowing the laws
governing the system.

The matter effects used in Sec.3.4.1 and Sec.3.4.2 are simple examples and do
not model specific physical systems. We employ simple matter effects of
comparable magnitude to $V_{e}(x)$. However, the method used here is
sufficiently general that it works independently of the choice of the matter
effect function.

####  3.4.1 Case 1

First, as a simple case, consider that a fixed potential is added as
follows($\sigma=100$),

$V_{total}(x)=V_{e}(x)+\sigma\sin{(\pi x)}.$ (44)

We assume the particle moves in $-1<x<1$, and the experimenters can vary
$V_{e}(x)$ in a variety of ways. The experimenters do not know the matter
effect, $\sigma\sin{(\pi x)}$, and then $V_{total}(x)$, but the energy
eigenvalues are assumed to be observable in the experiment. The experimenters
can train NN with their generated $V_{e}(x)$ as inputs and the observed energy
eigenvalues as outputs.

We here use the N5 dataset as the input potentials, $V_{e}(x)$, for the
numerical experiment. The energy eigenvalues used for the ground truth data
are calculated from $V_{total}(x)$ using Schrödinger equation to make
artificial experimental results datasets. These values correspond to the
observed energy eigenvalues without the experimental error.

Figure 13 shows the results of numerical experiments in this case. The left
figure shows the input potentials $V_{e}(x)$, $\sigma\sin{(\pi x)}$,
$V_{total}(x)$, the ground truth energy eigenvalues, and their predictions.
This figure shows that the energy eigenvalues are correctly predicted. The
NMSE is depicted in the figure on the right. The losses are comparable to
those shown in Fig.4(b), indicating that the learning has progressed
sufficiently. These figures demonstrate that NN comprehends the law composed
of Eqs. (10) and (44).

Figure 13:  (a) The energy eigenvalues predicted by NN for a modified
potential by the law, Eq.(44). The N5 exponential dataset is used for
training. (b) NMSE for the combinations of the training and validation
datasets in case 1.

####  3.4.2 Case 2

Next, consider more complex case where the matter effect depends on $V_{e}(x)$
as follows,

$V_{total}(x)=V_{e}(x)+\frac{1}{2\sigma}V_{e}(x)^{2}-\sigma.$ (45)

We also use the N5 dataset as the input potential, $V_{e}(x)$, and
$\sigma=100$.

Figure 14 shows the results of numerical experiments in this case. The left
figure shows the input potentials $V_{e}(x)$,
$\frac{1}{2\sigma}V_{e}(x)^{2}-\sigma$, $V_{total}(x)$, the ground truth
energy eigenvalues, and their predictions. This figure shows that the energy
eigenvalues are also correctly predicted in this case. This figure shows that
the energy eigenvalues are also correctly predicted in this case. The NMSE is
depicted in the figure on the right. In comparison to Fig.13(b) in Case 1, the
error is slightly larger, but learning has progressed sufficiently in this
case as well. This diagram shows that NN comprehends the law composed of
Eqs.(10) and (45).

Even though $V_{total}$ is not given, the NN can predict the ten energy
eigenvalues for given $V_{e}$ in both cases, 1 and 2. These numerical
experiments demonstrate that the NN method can be used even when the laws
governing the physical system are unknown, as long as sufficient experimental
results are available.

Figure 14:  (a) The energy eigenvalues predicted by NN for a modified
potential by the law, Eq.(45). The N5 exponential dataset is used for
training. (b) NMSE for the combinations of the training and validation dataset
in case 2.

####  3.4.3 Necessary experiments and NN’s versatility

As we have seen from the numerical experiments in cases 1 and 2, the NN can be
trained to make predictions even if the experimenters do not know the laws
governing the physical system. However, it is necessary to provide a variety
of inputs and obtain the corresponding experimental results as outputs. As a
result, an experimental setup should be prepared to automatically generate
various inputs while monitoring the results.

As shown in Fig.13(b) and Fig.14(b), the NN trained on the potentials
generated by the exponential kernel predicts well the potentials generated by
the other kernels. This finding suggests that if NNs are trained on a variety
of easily experimented-with inputs, they can predict the outcomes of
experiments even in difficult-to-experiment-with input domains. This
observation demonstrates that sufficiently versatile NNs can be created even
from experiments in a limited environment.

This NN-based method would also be appropriate when the laws governing a
physical system are known but difficult to solve numerically, due to high
dimension or multi-particle interactions, for example. Even in such cases, NN
can be trained using a variety of experimental results. Once the trained NN is
obtained, the results of experiments can be predicted.

##  4 Conclusion

This paper investigated how NNs understand physics using 1D quantum mechanics
as a subject. In the preparation stage in Sec.2, it was found that a DNN can
find energy eigenvalues only from the potentials with errors of approximately
1$\%$.

We showed from four aspects that NNs could understand physics and are helpful
for physics research in Sec.3.

  * •

The NN could predict energy eigenvalues, even for different kinds of
potentials from those used for training. For example, a NN trained on a
dataset created from jagged potentials could successfully predict the energy
eigenvalues of a smooth potential. We also found that capturing the overall
potential shape using the self-attention mechanism and deepening the NN
through residual-connection improved prediction accuracy.

  * •

The NN could predict the existence probability distribution of particles not
used during training. The results show that extracting physical quantities not
given during training from the NN is possible.

  * •

The NN predicted the crossover of states when the potential changes slowly and
continuously. The results show that NN can reproduce or predict physical
phenomena.

  * •

NNs are found to be trained with diverse inputs and the corresponding
experimental results, even when the laws governing the physical system are
unknown to the experimenters. Once the NN has been trained, it can be used to
predict the experimental results of the systems.

As shown in Fig.15, how NNs and humans understand physics seems different.
Humans understand physics using their intuition from experimental data,
discover unifying laws such as the Schrödinger equation that governs physical
systems, or build phenomenological approximation models. On the other hand,
NNs attempt to understand physical systems by adjusting numerous parameters,
taking advantage of their high computational power to explain a large amount
of experimental data. Although humans and NNs follow different paths to
understand physics, the underlying physics is the same, and, as a result, they
can make the same predictions.

Figure 15: Differences in the understanding of physics by humans and neural
networks. Humans look at the results of physical experiments and use their
intuition to derive concise physical laws and create models that describe
physical systems well. On the other hand, neural networks use significant
computational power to optimize numerous parameters from the results of
physical experiments. Both can make predictions about physics, but the process
is entirely different.

The understanding of physical systems by NNs has the following advantages.

  * •

It allows us to tackle problems without prior intuition. The NN used in this
paper is a combination of layers used in natural language processing and image
recognition, and no part of it has been specially devised for physics.

  * •

We can create a dataset in a region where experiments or numerical simulations
are simple to perform and use it to train NNs. As demonstrated in this paper,
this trained NN is versatile and can be used to make predictions in conditions
that are difficult to experiment with or numerically simulate. Figure 16
depicts a schematic diagram of the procedure.

On the other hand, it has the following drawbacks.

  * •

Massive amounts of experimental data or highly diverse numerical simulation
results are required to train a NN with high reproducibility while avoiding
overfitting. Experiments for these preparations must be automated. High
computer performance and clever algorithms are required when numerical
simulations, such as Monte Carlo methods, generate the training dataset.

  * •

It is possible to improve the approximation’s accuracy, but it does not appear
to be as simple to obtain the exact solution as the Schrödinger equation so
far. Humans still appear to be in charge of finding exact solutions or
approximating physical systems with a model with few parameters. However, the
ability of NNs to predict the outcomes of experiments under a variety of
conditions will aid humans in locating them.

These findings demonstrate that NNs can learn physical laws from experimental
data, predict experimental results under conditions different from those used
for training, and predict physical quantities of types not provided during
training. Because NNs understand physics in a different way than humans, they
will be a powerful tool for advancing physics by complementing the human way
of understanding.

Figure 16:  Using NN for physics research. Initially, a dataset is prepared in
a training area within a physical system that adheres to the same physical
laws. If the dataset is created through experimentation, this training area
will meet conditions such as temperature and material, making it easy to
conduct experiments. If the dataset was generated using a numerical
simulation, such as a Monte Carlo method, the training area corresponds to a
parameter region with fast convergence. NNs can be trained using data
generated by experiments or numerical simulations. As demonstrated in this
paper, NNs are versatile and can be used to predict experimental results under
difficult conditions to perform experiments or numerical simulations.

##  Appendix A Architecture of Neural Networks

This appendix describes the architecture of the NN we used in this paper. Our
NN consists of a combination of the Resnet and SAGAN structures, which have
recently been used with success in natural language processing and image
generationhe2016deep ; zhang2019self . The residual connection used in Resnet
mitigates gradient vanishing and increases the depth of the NN. On the other
hand, the self-attention mechanism used in SAGAN mitigates the shortcoming
that convolutional layers cannot relate to distant neurons and is very
powerful in natural language processing and image generation.

We began with the overall structure, as shown in Fig.17. First, the input is a
101-dimensional vector corresponding to the potential value at each point.
Next, this input is stretched in the channel direction by the blocks
containing the convolutional layer, reducing it to 17 dimensions in the
spatial direction, whereas it becomes 256 dimensions in the channel direction.
Finally, it is flattened to a 1D vector and passed through full connection
layers, outputting a 10-dimensional vector. This 10-dimensional vector
corresponds to ten small energy eigenvalues.

Figure 18 shows the details of the Res-SA block used in Fig.17. The input
first passes through the convolutional layer and then passes through the
residual block three times. Finally, a self-attention layer is used to relate
the information about distant locations. We embedded the self-attention layer
in the same way as in the SAGAN studyzhang2019self . The only difference is
that SAGAN needed to flatten a 2D image to 1D to create the attention map, but
our 1D potential can skip this process.

In Fig.19, we show the structure of the residual block. This block consists of
three convolutional layers and one residual connection and has the feature
that the input and output have the same shape. This structure is the same as
the one used in Resnet.

We have employed the ReLU function for all activation functions, totaling 31
times. The total number of parameters used in this network is approximately
$4.9\times 10^{6}$.

Next, we explain the specific method of learning. We employed Adam as the
optimization method and varied the learning rate from $1\times 10^{-3}$
(initial value) to $1\times 10^{-5}$ (value at the last epoch) in a
geometrical progression. In addition, since learning is more efficient when
the input and output values are of order O(1), we multiplied the input and
output values by 0.01 and adjusted them to be of that magnitude. The batch
size was fixed to 100. Pytorch was the NN framework used, and the GPy library
provided the kernel functions used. The Numpy library generated the random
numbers from the covariance matrix. For diagonalization of the matrix, we used
the Scikit-learn library.

Figure 17: Design of the entire neural network. The rectangles show the shapes
of the tensor at that position. First, the input is the potential value at 101
points in a one-dimensional space. Then, passing through the Res-SA block five
times, the tensor becomes longer in the channel direction and shorter in the
spatial direction. Finally, it is transformed into a flattened 1-dimensional
tensor, which goes through full connection layers and becomes a 10-dimensional
vector. This 10-dimensional vector corresponds to the values of the ten lowest
energy eigenvalues. The Res-SA block in this neural network is depicted in
detail in Fig.18. Figure 18: Res-SA block design. The input tensor can change
the number of channels in the initial convolution layer. At that time, the
kernel size is 5, and the padding size is 0, so the tensor size in the spatial
direction is reduced by 4. If there is a pooling layer, the tensor size in the
spatial direction is halved there. For the other layers, the size of the
tensor does not change. The deep structure increases the flexibility of the
function, and the self-attention layer helps relate the information of points
far apart in position. For the self-attention layer, we used the same
structure as in the SAGAN paperzhang2019self . The details of the residual
block are depicted in Fig.19. Figure 19:  Design of the residual block. The
block consists of two pointwise convolutional layers corresponding to full
connection layers in the channel direction and a convolutional layer with a
kernel size of 3. The padding is adjusted so that the input and output tensors
have the same shape. This structure is the same as the one used in Resnet. The
presence of residual coupling allows learning to proceed without gradient
loss, even as the depth of the neural network increases.

## References

  * (1) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  Deep residual learning for image recognition,  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778 (2016).
  * (2) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin,  Attention is all you need,  In Proceedings of the 31st International Conference on Neural Information Processing Systems, pages 6000–6010 (2017).
  * (3) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova, arXiv:1810.04805 (2018).
  * (4) Diederik P Kingma and Max Welling, arXiv:1312.6114 (2013).
  * (5) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio,  Generative adversarial nets,  In Proceedings of the 27th International Conference on Neural Information Processing Systems, volume 2, pages 2672–2680 (2014).
  * (6) David Silver, Aja Huang, Chris J. Maddison, Arthur Guez, Laurent Sifre, George van den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, Sander Dieleman, Dominik Grewe, John Nham, Nal Kalchbrenner, Ilya Sutskever, Timothy Lillicrap, Madeleine Leach, Koray Kavukcuoglu, Thore Graepel, and Demis Hassabis, Nature, 529(7587), 484–489 (January 2016).
  * (7) David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, Yutian Chen, Timothy Lillicrap, Fan Hui, Laurent Sifre, George van den Driessche, Thore Graepel, and Demis Hassabis, Nature, 550(7676), 354–359 (October 2017).
  * (8) Giuseppe Carleo, Ignacio Cirac, Kyle Cranmer, Laurent Daudet, Maria Schuld, Naftali Tishby, Leslie Vogt-Maranto, and Lenka Zdeborová, Rev. Mod. Phys., 91(4) (December 2019).
  * (9) M. Raissi, P. Perdikaris, and G.E. Karniadakis, J.Comput.Phys., 378, 686–707 (February 2019).
  * (10) Rui Hong, Peng-Fei Zhou, Bin Xi, Jie Hu, An-Chun Ji, and Shi-Ju Ran, SciPost Phys. Core, 4(3), 022 (September 2021).
  * (11) Arijit Sehanobish, Hector H. Corzo, Onur Kara, and David van Dijk (January 2021), arXiv:2006.13297.
  * (12) Juncai Pu, Jun Li, and Yong Chen, Nonlinear Dynamic, 105, 1723 (2021), arXiv:2011.04949.
  * (13) Mitsumasa Nakajima, Kenji Tanaka, and Toshikazu Hashimoto, IEEE Trans. Neural Netw. Learning Syst., 33, 2686 (2022).
  * (14) Kyle Mills, Michael Spanner, and Isaac Tamblyn, Phys. Rev. A, 96(4), 042113 (October 2017).
  * (15) Raban Iten, Tony Metger, Henrik Wilming, Lidia del Rio, and Renato Renner, Phys. Rev. Lett., 124(1), 010508 (January 2020).
  * (16) Han Zhang, Ian Goodfellow, Dimitris Metaxas, and Augustus Odena,  Self-attention generative adversarial networks,  In International conference on machine learning, pages 7354–7363. PMLR (2019).
  * (17) Joshua Izaac and Jingbo Wang,  Computational Quantum Mechanics,  (Springer, New York, NY, 1st ed. 2018 edition edition, March 2019).
  * (18) Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman, arXiv:1312.6034 (2013).
  * (19) L. Wolfenstein, Phys. Rev. D, 17(9), 2369–2374 (1978).
  * (20) SP Mikheyev and A Yu Smirnov, Il Nuovo Cimento C, 9(1), 17–26 (1986).
  * (21) Stephen J. Parke, Phys. Rev. Lett., 57(10), 1275–1278 (1986).
  * (22) H. A. Bethe, Phys. Rev. Lett., 56(12), 1305–1308 (1986).


#  Resistive and magnetized accretion flows with convection

Kazem Faghei and Mobina Omidvand

###### Abstract

We considered the effects of convection on the radiatively inefficient
accretion flows (RIAF) in the presence of resistivity and toroidal magnetic
field. We discussed the effects of convection on transports of angular
momentum and energy. We established two cases for the resistive and magnetized
RIAFs with convection: assuming the convection parameter as a free parameter
and using mixing-length theory to calculate convection parameter. A self-
similar method was used to solve the integrated equations that govern the
behavior of the presented model. The solutions showed that the accretion and
rotational velocities decrease by adding the convection parameter, while the
sound speed increases. Moreover, by using mixing-length theory to calculate
convection parameter, we found that the convection can be important in RIAFs
with magnetic field and resistivity.

00footnotetext: School of Physics, Damghan University, Damghan, Iran
e-mail<EMAIL_ADDRESS>

Keywords accretion, accretion discs, convection, magnetohydrodynamics: MHD

##  1 Introduction

The existence of radiatively inefficient accretion flows (RIAFs) have been
confirmed in low-luminosity state of X-ray binaries and nuclei of galaxies
(Narayan et al. 1996; Esin et al. 1997; Di Matteo et al. 2003; Yuan et al.
2003). It was understood that RIAFs are likely to be convectively unstable in
the radial direction due to the inward increase of the entropy of accreting
gas (Narayan & Yi 1994). Moreover, hydrodynamical and magetohydrodynamical
simulations of low-viscosity RIAFs have confirmed these flows are convectively
unstable ( e. g. Igumenshchev et al. 1996; Stone et al. 1999; Machida et al.
2001; Hawley & Balbus 2002; McKinney & Gammie 2002; Igumenshchev et al. 2003).
Self-similar or global solutions for convection-dominated accretion flows
(CDAFs) were presented by several authors (e. g. Narayan et al. 2000; Quataert
& Gruzinov 2000; Abramowicz et al. 2002; Lu et al. 2004; Zhang & Dai 2008).

Igumenshchev et al. (2003) studied the resistive MHD simulations of RIAFs onto
black holes. They assumed two cases for the geometry of the injected magnetic
field: pure toroidal field and pure poloidal field. They found that in the
case of pure toroidal magnetic field, the accreting gas forms a nearly
axisymmetric, geometrically thick, turbulent accretion disc. Moreover, their
solutions represented that the flow resembles in many respects CDAFs found in
previous numerical and analytical investigations of viscous hydrodynamic
flows. Zhang & Dai (2008) investigated the effect of magnetic field on RIAFs
with convection by a semi-analytically method. By exploit of
$\alpha$-prescription for viscosity and convection, they used two methods to
study of magnetized flows with convection, i.e. they take the convective
coefficient $\alpha_{c}$ as a free parameter to discuss the effects of
convection for simplicity. They also established the $\alpha_{c}$-$\alpha$
relation for magnetized flows using the mixing-length theory and compare this
relation with the non-magnetized case. They found that the magnetic field
makes the $\alpha_{c}$-$\alpha$ relation be distinct from that of non-
magnetized flows.

Since the importance of toroidal magnetic field and resistivity in accretion
flows have been confirmed observationally (see Faghei 2011 and references
therein), Faghei (2011) considered the steady, radially self-similar solutions
of accretion flows in the presence of the toroidal magnetic field and the
resistivity. However, he ignored the effects of convection in his model.
Generally semi-analytical studies of magnetized CDAFs are related to non-
resistive magnetized CDAFs (e. g. Zhang & Dai 2008) and the resistive and
magnetized CDAF was studied in MHD simulations (e. g. Igumenshchev & Narayan
2002; Hawley & Balbus 2002; Igumenshchev et al. 2003). Thus, it will be
interesting to study the effects of resistivity on RIAFs with convection.
Here, we adopt the presented solutions by Narayan et al. (2000) and Faghei
(2011). Similar to Narayan et al. (2000), we will discuss the effects of
convection on angular momentum and energy equations. The paper is organized as
follow. In section 2, the basic equations of constructing a model for quasi-
spherical magnetized RIAFs with convection will be defined. In section 3, a
self-similar method for solving equations which govern the behavior of the
accreting gas was utilized. The summary of the model will appear in section 4.

##  2 Basic Equations

Analytical theory of CDAF is based on a self-similar solution of a simplified
set of equations describing RIAFs. We adopted the presented solutions by
Narayan et al. (2000) and Faghei (2011). By using spherical coordinate ($r$,
$\theta$, $\varphi$) centered on a accreting object, let us consider
stationary, axisymmetric, quasi-spherical equations describing an accretion
flow onto the black hole of mass $M$. For the sake of simplicity, the general-
relativistic effect has been neglected and the gravitational force on a fluid
is characterized by Newtonian potential of a point mass, $\psi=-GM/r$. As
magnetic fields, we consider only toroidal fields, $B_{\varphi}$.

Under these assumptions, the continuity equation with mass loss is

$\frac{1}{r^{2}}\frac{d}{dr}(r^{2}\rho v_{r})=\dot{\rho},$ (1)

where $\rho$, $v_{r}$ and $\dot{\rho}$ are the density, the accretion velocity
($v_{r}<0$), and the mass-loss per unit volume, respectively.

The radial momentum equation is

$v_{r}\frac{dv_{r}}{dr}=r\left(\Omega^{2}-\Omega^{2}_{K}\right)-\frac{1}{\rho}\frac{d}{dr}(\rho
c^{2}_{s})-\frac{c^{2}_{A}}{r}-\frac{1}{2\rho}\frac{d}{dr}(\rho c^{2}_{A}),$
(2)

where $c_{s}$ is sound speed, which is defined as $c_{s}^{2}\equiv
p_{gas}/\rho$, with being $p_{gas}$ as the gas pressure, $\Omega$ is the
angular velocity, $\Omega_{K}[=\left(GM/r^{3}\right)^{1/2}]$ is the Keplerian
angular velocity, and $c_{A}$ is the alfven speed, which is defined as
$c_{A}^{2}\equiv B_{\varphi}^{2}/4\pi\rho=2p_{mag}/\rho$, with being $p_{mag}$
as the magnetic pressure. The ram-pressure term $v_{r}dv_{r}/dr$ and last two
terms due to the magnetic field in this equation were ignored in the self-
similar CDAF model of Narayan et al. (2000), while we include them here in
order to consider their effects.

The angular momentum equations can be written in the form of the balance of
advection and diffusion transport terms (Narayan et al. 2000),

$\displaystyle\rho
v_{r}\frac{d}{dr}(r^{2}\Omega)=\frac{1}{r^{2}}\frac{d}{dr}\left[\nu\rho
r^{4}\frac{d\Omega}{\partial r}\right]+~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\left[\nu_{c}\rho
r^{(5+3g)/2}\frac{d}{dr}\left(\Omega r^{3(1-g)/2}\right)\right],$ (3)

where the two terms of right hand side represent the angular momentum
transport by viscosity and convection. Here, $\nu$ is the kinematic viscosity
coefficient, $\nu_{c}$ is the convective diffusion coefficient, and $g$ is the
parameter to determine the condition of convective angular momentum transport.
When $g=1$, the flux of angular momentum due to convection is

$\dot{J}_{c}=-\nu_{c}\,\rho\,r^{4}\,\frac{d\Omega}{dr}.$ (4)

The above equation implies that the convective angular momentum flux is
oriented down the angular velocity gradient. For a quasi-Keplerian angular
velocity, $\Omega\propto r^{-3/2}$, angular momentum is transported outward.
When $g=-1/3$, the convective angular momentum flux can be written as

$\dot{J}_{c}=-\nu_{c}\,\rho\,r^{2}\,\frac{d\left(\Omega r^{2}\right)}{dr}.$
(5)

This equation represents that the convective angular momentum flux is oriented
down the specific angular momentum gradient. For a quasi-Keplerian angular
velocity, $\Omega\propto r^{-3/2}$, angular momentum is transported inward.
Generally, convection transports angular momentum inward (or outward) for
$g<0$ (or $>0$), and the specific case $g=0$ corresponds to zero angular
momentum transport (Narayan et al. 2000).

In this paper, we assume the kinematic coefficient of viscosity and the
magnetic diffusivity due to turbulence in the accretion flow. So, we use these
parameters in analogy to the $\alpha$-prescription of Shakura & Sunyaev (1973)
for the turbulent,

$\nu=P_{m}\eta=\alpha\frac{c_{s}^{2}}{\Omega_{K}},$ (6)

where $P_{m}$ is the magnetic Prandtl number of the turbulence, which assumed
to be a constant less than unity, $\eta$ is the magnetic diffusivity, and
$\alpha$ is a free parameter less than unity. For the convective diffusion
coefficient, $\nu_{c}$, we adopt the assumptions of Narayan et al. (2000) and
Lu et al. (2004) that all transport phenomena due to convection have the same
diffusion coefficient, which is defined as

$\nu_{c}=\left(\frac{L_{M}^{2}}{4}\right)\sqrt{-N_{eff}^{2}},$ (7)

where $L_{M}$ is the characteristic mixing length and $N_{eff}$ is the
effective frequency of convective blobs. The characteristic mixing length
$L_{M}$ in terms of the pressure scale height, $H_{p}$, can be written as

$L_{M}=2^{-1/4}l_{M}H_{p},~{}~{}~{}~{}H_{p}=-\frac{dr}{d\ln p_{gas}},$ (8)

where $l_{M}$ is the dimensionless mixing-length parameter and its amount is
estimated to be equal to $\sqrt{2}$ in ADAFs (Narayan et al. 2000; Lu et al.
2004). the effective frequency of convective blobs, $N_{eff}$, is given by

$N_{eff}^{2}=N^{2}+\kappa^{2},$ (9)

where $N$ is Brunt-Väisälä frequency, which is defined as

$N^{2}=-{1\over\rho}{dp_{gas}\over dr}{d\over
dr}\ln\left({p_{gas}^{1/\gamma}\over\rho}\right),$ (10)

and $\kappa$ is epicyclic frequency, which is defined as

$\kappa^{2}=2\Omega^{2}\frac{d\ln(\Omega r^{2})}{d\ln r}.$ (11)

For a non-Keplerian flows $\kappa\neq\Omega$, while for a quasi-Keplerian
($\Omega\propto r^{-3/2}$), $\kappa=\Omega$ (Narayan et al. 2000; Lu et al.
2004). Convection appears in flows with $N_{eff}^{2}<0$. We also write the
convective diffusion coefficient in the form similar to usual viscosity of
Shakura & Sunyaev (1973),

$\nu_{c}=\alpha_{c}\frac{c_{s}^{2}}{\Omega_{K}}$ (12)

where $\alpha_{c}$ is a dimensionless coefficient that describes the strength
of convective diffusion. The $\alpha_{c}$ coefficient can be obtained by
equations (8) and (13)

$\alpha_{c}=\frac{\Omega_{K}}{c_{s}^{2}}\left(\frac{L_{M}^{2}}{4}\right)\sqrt{-N_{eff}^{2}}.$
(13)

The energy equation is

$\displaystyle\rho v_{r}T\frac{ds}{dr}\equiv\rho
v_{r}\left[\frac{1}{\gamma-1}\frac{dc_{s}^{2}}{dr}-\frac{c_{s}^{2}}{\rho}\frac{d\rho}{dr}\right]=$
$\displaystyle Q_{diss}+Q_{conv}-Q_{rad},$ (14)

where $T$ is the temperature, $s$ is the specific entropy, $\gamma$ is the
ratio of specific heats, $Q_{diss}$ is dissipative heating rate, $Q_{rad}$ is
the radiative cooling rate, and
$Q_{conv}=-\mathbf{\nabla}\cdot\mathbf{F}_{conv}$, with being
$F_{conv}[=-\rho\nu_{c}Tds/dr]$ as the outward energy flux due to convection.
For the right hand side of the energy equation, we can write

$\displaystyle
Q_{adv}=fQ_{diss}-\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}F_{conv}\right),$ (15)

where $Q_{adv}$ is the advective transport of energy, and
$f[=1-Q_{rad}/Q_{diss}]$ is the advection parameter. The parameter $f$
measures the degree to which the flow is advection-dominated (Narayan & Yi
1994). The dissipative heating rate can be written as

$Q_{diss}=(\nu+g\nu_{c})\rho r^{2}\left(\frac{\partial\Omega}{\partial
r}\right)^{2}+\frac{\eta}{4\pi}{\bf J}^{2},$ (16)

where the right-hand side terms are heating rate due to viscosity, convection,
and resistivity, respectively. In above equation,
${\mathbf{J}}[=\nabla\times{\mathbf{B}}]$ is the current density, with being
${\mathbf{B}}$ as the magnetic field.

Finally, the induction equation with creation/escape of magnetic field can be
written as

$\frac{1}{r}\frac{d}{dr}\left[rv_{r}B_{\varphi}-\eta\frac{d}{dr}(rB_{\varphi})\right]=\dot{B}_{\varphi}.$
(17)

where $B_{\varphi}$ is the toroidal component of magnetic field and
$\dot{B}_{\varphi}$ is the field escaping/creating rate due to a magnetic
instability or dynamo effect. This induction equation is rewritten as

$\displaystyle\dot{B}_{\varphi}=~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\frac{1}{r}\frac{d}{dr}\left[\sqrt{4\pi\rho
c^{2}_{A}}\left(rv_{r}-\frac{\alpha}{4\chi
P_{m}}\frac{1}{r\rho\Omega_{K}}\frac{d}{dr}(r^{2}\rho
c^{2}_{A})\right)\right],$ (18)

where $\chi$ is the ratio of the magnetic pressure to the gas pressure, which
is defined by

$\chi=\frac{p_{mag}}{p_{gas}}=\frac{1}{2}\left(\frac{c_{A}}{c_{s}}\right)^{2}.$
(19)

##  3 Self-Similar Solutions

We seek self-similar solutions in the following form (e.g. Narayan & Yi 1994;
Akizuki & Fukue 2006)

$v_{r}(r)=-c_{1}\alpha\sqrt{\frac{GM_{*}}{r}}$ (20)
$\Omega(r)=c_{2}\sqrt{\frac{GM_{*}}{r^{3}}}$ (21)
$c^{2}_{s}(r)=c_{3}\frac{GM_{*}}{r}$ (22)
$c^{2}_{A}(r)=\frac{B^{2}_{\varphi}}{4\pi\rho}=2\chi c_{3}\frac{GM_{*}}{r}$
(23)

where $c_{1}$, $c_{2}$, and $c_{3}$ are dimensionless constant to be
determined. We use a power-law relation for density

$\rho(r)=\rho_{0}r^{\lambda},$ (24)

where $\rho_{0}$ and $\lambda$ are constant. Using equations (20)-(24), the
mass-loss rate and the magnetic field escaping/creating rate can be written as

$\dot{\rho}(r)=\dot{\rho}_{0}r^{\lambda-3/2},$ (25)

$\dot{B}_{\varphi}(r)=\dot{B}_{0}r^{\frac{\lambda-4}{2}},$ (26)

where $\dot{\rho}_{0}$ and $\dot{B}_{0}$ are constant. Since we have not
applied the effects of wind in the momentum and energy equations, we will
assume a no wind case, $\dot{\rho}=0$ and $\lambda=-3/2$. In this case,
$\dot{B}_{\varphi}\propto r^{-11/4}$, which implies that creation/escape of
magnetic field increases with approaching to central object. This property is
qualitatively consistent with previous studies of accretion flows (Machida et
al. 2006; Oda et al. 2007; Faghei & Mollatayefeh 2012).

Using the self-similar solutions in the continuity, radial momentum, angular
momentum, convection parameter, energy, and induction equations [(1)-(3),
(13), (14), and (18)], we can obtain the following relations:

$\dot{\rho}_{0}=-\left(\lambda+\frac{3}{2}\right)\alpha\rho_{0}c_{1}\sqrt{GM_{*}},$
(27)
$-\frac{1}{2}c^{2}_{1}\alpha^{2}+1-c^{2}_{2}+c_{3}\left[\lambda-1+\chi(1+\lambda)\right]=0,$
(28) $\alpha c_{1}=3(\alpha+g\alpha_{c})(\lambda+2)c_{3},$ (29)
$\displaystyle\alpha
c_{1}\left[\frac{1}{\gamma-1}+\lambda\right]=~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\frac{9}{4}\alpha
f\left[(1+\frac{\alpha_{c}}{\alpha}g)c_{2}^{2}+\frac{2\chi}{9P_{m}}c_{3}(1+\lambda)^{2}\right]$
$\displaystyle-\alpha_{c}c_{3}(\lambda+\frac{1}{2})\left[\frac{1}{\gamma-1}+\lambda\right]$
(30)
$\alpha_{c}=\frac{l_{M}^{2}}{4\sqrt{2}c_{3}(\lambda-1)^{2}}\sqrt{\frac{c_{3}(\lambda-1)}{\gamma}[\lambda(1-\gamma)-1]-c_{2}^{2}},$
(31) $\dot{B}_{0}=-\frac{\alpha\lambda}{2}GM_{*}\sqrt{2\pi\rho_{0}\chi
c_{3}}\left[2c_{1}+\frac{c_{3}}{P_{m}}(1+\lambda)\right].$ (32)

Above equations express for $\lambda=-3/2$, there is no mass loss, while for
$\lambda>-3/2$ mass loss (wind) exists.

##  4 Results

Here, similar to Zhang & Dai (2008), we will study the presence of convection
in two cases: $\alpha_{c}$ as a free parameter and $\alpha_{c}$ as a variable.

Fig. 1 : Physical variables as functions of $\chi$ for several values of
convective viscosity. The input parameters are set to $\alpha=0.2$,
$\gamma=1.5$, $P_{m}=1/2$, $f=1$, $l=\sqrt{2}$, $g=-1/3$, and $\lambda=-3/2$.
The solid, dashed, and dotted lines represent $\alpha_{c}=0$, $0.05$, and
$0.1$, respectively.

Fig. 2 : Same as Figure 1, but $\alpha_{c}=0.1$, and the solid, dashed, and
dotted lines represent $P_{m}=\infty$, $1.0$, and $0.5$, respectively.

Fig. 3 : Physical variables as functions of $\chi$ for several values of
magnetic Prandtl number. The input parameters are set to $\alpha=0.5$,
$\gamma=1.5$, $f=1$, $l=\sqrt{2}$, $g=-1/3$, and $\lambda=-3/2$. The solid,
dashed, and dotted lines represent $P_{m}=\infty$, $1.0$, and $0.5$,
respectively.

###  4.1 Case 1: $\alpha_{c}$ as a free parameter

In this case, we take the convective coefficient $\alpha_{c}$ as a free
parameter to discuss the effects of convection for simplicity. Examples of
such solutions are presented in Figures 1 and 2.

In Figure 1, the self-similar coefficients $c_{1}$, $c_{2}$, and $c_{3}$ are
shown as functions of the parameter $\chi$. By adding the parameter $\chi$
which indicates the role of magnetic filed on the dynamics of accretion discs,
we see the coefficients of radial and rotational velocities and sound speed
decrease. This properties are qualitatively consistent with results of Faghei
(2011). In Figure 1, we also studied the effect of convection parameter
$\alpha_{c}$ on the physical variables. The value of $\alpha_{c}$ measures the
strength of convective viscosity and a larger $\alpha_{c}$ denotes a stronger
turbulence due to convection. Figure 1 implies that for non-zero $\alpha_{c}$,
the radial infall velocity is lower than the standard ADAF solution and for
larger $\alpha_{c}$ this reduction of radial infall velocity is more evident.
It can be due to decrease of efficiency of angular momentum transport by
adding the convection parameter $\alpha_{c}$ (see equation 29). The profiles
of angular velocity show that it decreases with the magnitude of $\alpha_{c}$,
while the sound speed increases. These properties are in accord with results
of Zhang & Dai (2008).

In Figure 2, the physical variables are shown as functions of parameter $\chi$
for several values of magnetic Prandtl number. Since inverse of magnetic
Prandtl number is proportional to magnetic diffusivity,
$P_{m}\propto\eta^{-1}$. Thus, reduce of magnetic Prandtl number denotes to
increase of resistivity of the fluid. The solutions in Figure 2 imply that the
accretion velocity and the sound speed both increase with the magnitude of
resistivity, while the rotational velocity decreases. These properties
qualitatively confirm the results of Faghei (2011).

###  4.2 Case 2: $\alpha_{c}$ as a variable

Here, we calculate the dimensionless coefficient $\alpha_{c}$ by using the
mixing-length theory. Because we used a steady self-similar method to derive
$\alpha_{c}$, it becomes a constant throughout of the accreting gas. However,
it is a function of position and time (e. g. Lu et al. 2004). The amount of
convection parameter $\alpha_{c}$ is calculated by equation (31). Using this
equation and equation (28)-(30), we can obtain the behavior of physical
quantities in the presence of convection. Such solutions are shown in Figure
3.

In Figure 3, the coefficients $c_{1}$, $c_{2}$, $c_{3}$, and convection
parameter $\alpha_{c}$ are shown as functions of the degree of magnetic
pressure. Similar to case 1, the accretion and rotational velocities, and
sound speed decrease by adding the parameter $\chi$. While, the convection
parameter $\alpha_{c}$ increases for stronger toroidal magnetic field. This
property is qualitatively consistent with result of Zhang & Dai (2008). In
Figure 3, the physical variables are also studied for several values of
magnetic Prandtl number. The profiles of convection parameter $\alpha_{c}$
imply that it increases by adding the magnetic diffusivity. As for non-zero
magnetic diffusivity, $\alpha_{c}$ is larger than the standard CDAF solution
and for larger magnetic diffusivity this increase of convection parameter
$\alpha_{c}$ is more evident.

##  5 Summary and Discussion

The observational features of low-luminosity state of X-ray binaries and
nuclei of galaxies can be successfully explained by the models of radiatively
inefficient accretion flow (RIAF). The importance of convection in RIAFs was
realized by semi-analytical and direct numerical simulation (e. g. Narayan et
al. 2000; Igumenshchev et al. 2003).

In this research, we considered the effects of convection on the presented
model of Faghei (2011). Similar to Narayan et al. (2000), we assumed the
convection affects on transports of angular momentum and energy. Using a
radially self-similar approach, we studied the effects of convection on the
model for several values of magnetic field and resistivity. The solutions
showed that the accretion and rotational velocities, and sound speed decrease
for stronger magnetic filed. Moreover, we found that the accretion velocity
and sound speed increase with the magnitude of the resistivity, while the
rotational velocity decreased. These properties are qualitatively consistent
with results of Faghei (2011). We studied the effects of convection on a
resistive and magnetized RIAF in two cases: assuming the convection parameter
as a free parameter and using mixing length theory to calculate the convection
parameter. In the first case, we found that by adding the convection
parameter, the radial and rotational velocities decrease and the sound speed
increases. In the second case, we found that the convection parameter
increases by adding the magnetic filed and resistivity. These properties are
in many aspects in accord with results of Zhang & Dai (2008).

The present model have some limitations that can be modified in the future
works. For example, the latitudinal dependence of physical variables have been
ignored in this paper. While, two-dimensional and three-dimensional MHD
simulations of RIAFs show that the disc geometry strongly depends on magnetic
field configuration (e. g. Igumenshchev et al. 2003). Thus, the study of
present model in two/three dimensions can be an interesting subject for future
research. Moreover, it has been understood the magnetic field can change the
criterion for convective instability (e. g. Balbus & Hawley 2002). While, we
igonred the effects of magnetic field on the instability criterion. Thus, the
presented criterion in this paper can be modified in the future research.

## Acknowledgements

I wish to thank the anonymous referee for very useful comments that helped us
to improve the initial version of the paper.

## References

  * (1) Abramowicz M. A., Igumenshchev I. V., Quataert E., Narayan R., 2002, ApJ, 565, 1101
  * (2)
  * (3) Akizuki C., Fukue J., 2006, PASJ, 58, 469
  * (4)
  * (5) Balbus S. A., Hawley J. F., 2002, ApJ, 573, 749
  * (6)
  * (7) Di Matteo T., Allen S. W., Fabian A. C., Wilson A. S., Young A. J., 2003, ApJ, 582, 133
  * (8)
  * (9) Esin A. A., McClintock J. E., Narayan R., 1997, ApJ, 489, 867
  * (10)
  * (11) Faghei K., 2011, JA&A, arXiv:1111.7302
  * (12)
  * (13) Faghei K., Mollatayefeh A., 2012, MNRAS,
  * (14) doi:10.1111/j.1365-2966.2012.20645.x
  * (15)
  * (16) Hawley J. F., Balbus S. A., 2002, ApJ, 573, 738
  * (17)
  * (18) Igumenshchev I. V., Chen X., Abramowicz M. A., 1996, MNRAS, 278, 236
  * (19)
  * (20) Igumenshchev I. V., Narayan R., 2002, ApJ, 566, 137
  * (21)
  * (22) Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ, 592, 1042
  * (23)
  * (24) Lu J.-F., Li S.-L., Gu W.-M., 2004, MNRAS, 352, 147
  * (25)
  * (26) Machida M., Matsumoto R., Mineshige S., 2001, PASJ, 53, L1
  * (27)
  * (28) Machida M., Nakamura K. E., Matsumoto R., 2006, PASJ, 58, 193
  * (29)
  * (30) McKinney J. C., Gammie C. F., 2002, ApJ, 573, 728
  * (31)
  * (32) Narayan R., Igumenshchev I. V., Abramowicz M. A., 2000, ApJ, 539, 798
  * (33)
  * (34) Narayan R., McClintock J. E., Yi I., 1996, ApJ, 457, 821
  * (35)
  * (36) Narayan R., Yi I., 1994, ApJ, 428, L13
  * (37)
  * (38) Oda H., Machida M., Nakamura K. E., Matsumoto R., 2007, PASJ, 59, 457
  * (39)
  * (40) Quataert E., Gruzinov A., 2000, ApJ, 539, 809
  * (41)
  * (42) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * (43)
  * (44) Stone J. M., Pringle J. E., Begelmen M. C., 1999, MNRAS, 310, 1002
  * (45)
  * (46) Yuan F., Quataert E., Narayan R., 2003, ApJ, 598, 301
  * (47)
  * (48) Zhang D., Dai Z. G., 2008, MNRAS, 388, 1409
  * (49)