# Learning Potentials of Quantum Systems using Deep Neural Networks

Arijit Sehanobish,1 Hector H. Corzo,11footnotemark: 12 Onur Kara,3 David van
Dijk1
Equal Contribution

###### Abstract

Attempts to apply Neural Networks (NN) to a wide range of research problems
have been ubiquitous and plentiful in recent literature. Particularly, the use
of deep NNs for understanding complex physical and chemical phenomena has
opened a new niche of science where the analysis tools from Machine Learning
(ML) are combined with the computational concepts of the natural sciences.
Reports from this unification of ML have presented evidence that NNs can learn
classical Hamiltonian mechanics. This application of NNs to classical physics
and its results motivate the following question: Can NNs be endowed with
inductive biases through observation as means to provide insights into quantum
phenomena? In this work, this question is addressed by investigating possible
approximations for reconstructing the Hamiltonian of a quantum system in an
unsupervised manner by using only limited information obtained from the
system’s probability distribution.

##  1 Introduction

In the Machine Learning (ML) realm, Neural Networks (NNs) are among the most
used and exceptionally efficient models to learn and generalize information
from data. These data interpretative capabilities have provoked the widespread
use of NNs in Natural Language Processing (Torfi et al. 2020), Image
Classification (Kolesnikov et al. 2019), Video Captioning (Sun et al. 2019)
and Reinforcement Learning (Du and Narasimhan 2019; Higgins et al. 2016);
furthermore, recent works have shown the capabilities of NNs in symbolic
reasoning and mathematical problem solving (Lample and Charton 2019). In the
case of natural sciences, applying ML to physics is not new, several works
have been reported (Toth et al. 2019; Greydanus, Dzamba, and Yosinski 2019;
Cranmer et al. 2020; Tong et al. 2020) where different authors have combined
ML with Hamilton’s equations of motion to generate trajectories that obey
energy conservation principles and classical physical laws. In material
sciences, on the other hand, ML has proven to be an important interpretative
tool for the computational prediction of new materials (Schleder et al. 2019).
The encouraging results reported by the different applications of ML have
motivated the use of NNs as powerful tools to gain insight into the laws of
physics that govern the behavior of complicated natural classical phenomena.
Unlike classical physics, in quantum physics, objects have characteristics of
both particles and waves (wave–particle duality) for which the concept of
trajectory is no longer defined nor can their position and momentum, both, be
measured simultaneously (Sakurai and Commins 1995; Robinett 1997; Feynman,
Leighton, and Sands 1965; Robinett and Robinett 2006). Because of this
wave–particle duality, the state of a quantum mechanical system is fully
specified by its wave–function, which is typically obtained by solving the
Schrödinger equation (Sakurai and Commins 1995; Robinett 1997; Feynman,
Leighton, and Sands 1965; Robinett and Robinett 2006). In many cases, however,
not only solving this equation is difficult, but also its correct formulation
requires knowledge about the form of the potential energy operator, which
often may not be completely known. In contrast, the inverse form of the
Schrödinger equation (Nakatsuji 2002; Chadan and Sabatier 2012; Zakhariev and
Suzko 2012; Jensen and Wasserman 2018) presents an alternative for describing
quantum phenomena by reformulating the description of quantum mechanical
systems as solutions of inverse problems (Aster, Borchers, and Thurber 2018;
Groetsch and Groetsch 1993; Vogel 2002). Inverse problems are central to the
study of quantum mechanical systems, insomuch that much of what is known about
the electronic structure of matter has been mathematically characterized by
solutions of inverse problems (Beals and Greiner 2009; Zakhariev and Suzko
2012; Jensen and Wasserman 2018; Vogel 2002). Thus, numerical algorithms for
the inversion of the Schrödinger equation are important predictive tools for
the further development of approximate quantum mechanical methodologies such
as scattering approximated models (Zakhariev and Suzko 2012), Density
Functional Theory (DFT) (Jensen and Wasserman 2018), etc.

In this work, rather than handcrafting numerical solutions for the inverse
Schrödinger equation (Aster, Borchers, and Thurber 2018; Jensen and Wasserman
2018; Vogel 2002; Beals and Greiner 2009) to define a potential function and
describe a quantum phenomena, a neural network, termed Quantum Potential
Neural Network (QPNN), is instead designed to learn potential functions
directly from observables in an unsupervised manner. This proposed QPNN for
learning potential functions was developed based on the underlying formalism
for the inverse solution of the Schrödinger equation. Thus, the proposed QPNN
opens the possibility for generating simpler and succinct functions that can
be used to construct effective Hamiltonians for the description of a variety
of quantum systems using only a small portion of the available information
known about the system. Since these effective Hamiltonians can be generalized
to obtain other observables, QPNN may provide unique insights into complex
quantum phenomena were only a small amount of information is available.

##  2 Theory

The mathematical description of a quantum particle typically takes the form of
a complex function of spatial coordinates $\vec{x}$ and time coordinates $t$
called wave–function, $\Psi(\vec{x},t)$ (Sakurai and Commins 1995; Robinett
1997; Feynman, Leighton, and Sands 1965; Robinett and Robinett 2006).
$\Psi(\vec{x},t)$ is a complex–valued probability amplitude whose square
modulus ($|\Psi(\vec{x},t)|^{2}$) correspond to the probability of finding the
particle described by the wave–function at that given $\vec{x}$ and $t$. The
classically measured value of a physical observable, however, is not given
directly by $\Psi(\vec{x},t)$ but by the expectation values of the operators
that represent the desired measurement acting on $\Psi(\vec{x},t)$. In many
scenarios, wave–functions are obtained as direct solutions of the
time–dependent Schrödinger equation,

$i\hbar\frac{\partial\Psi(\vec{x},t)}{\partial t}=\hat{\rm H}\Psi(\vec{x},t),$
(1)

where $\hbar$ is Planck’s constant and $\hat{\rm H}$ is the Hamiltonian
operator of the system, which is an Hermitian operator acting on an infinite
dimensional space of $L^{2}$ functions. Thus, $\hat{\rm H}$ needs not be
compact and as much may not have any eigenvalues. When $\hat{\rm H}$ is
time–independent, equation 1 can be reduced to the following equation

$\hat{\rm H}\psi_{n}(\vec{x})=E\psi_{n}(\vec{x}),$ (2)

where $n$ indicates the quantum state of the system. In many cases, the
physical information contained in the time–independent wave–function
$\psi_{n}(\vec{x})$ may be enough for the characterization of the system under
study.

###  2.1 Hamiltonian

The Hamiltonian operator, $\hat{\rm H}$, is fundamental in many formulations
of quantum theory. This operator is expressed as the sum of the kinetic
($\hat{\rm T}$) and potential energy operators ($\hat{\rm V}$) for all
particles in the quantum system,

$\hat{\rm H}=\hat{\rm T}+\hat{\rm V}.$ (3)

Generally, the kinetic energy operator contained in $\hat{\rm H}$ only depends
on the second derivatives of the wave–function, with respect to its spatial
coordinates. The potential energy operator, however, depends on the physical
circumstances imposed onto the system, and varies from system to system. Thus
equation 3 may be expressed as

$\hat{\rm
H}=-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial\vec{x}^{2}}+\hat{\rm
V}(\vec{x},t)\equiv-\frac{\hbar^{2}}{2m}\nabla_{\vec{x}}^{2}+\hat{\rm
V}(\vec{x},t).$ (4)

In quantum mechanics, the problem of finding the $\hat{\rm H}$ that
characterizes a given phenomenon could be reduced to formulating the potential
operator that contains all the governing physical descriptors.

Figure 1: Graphical representation of our framework

###  2.2 Predicting potentials

The usual method for describing systems in quantum mechanics is by obtaining
the wave–function of the system as a solution of the Schrödinger equation.
These wave–functions strongly depend on the Hamiltonian, and in particular,
the definition of the potential used to describe the system. However, one
could also describe a quantum phenomena through the solution of the inverse
problem, i.e. by finding an effective potential or function that contains all
the important physical constrains that generated the observed outcomes.
Inverse problems like this one are common in quantum mechanics and electronic
structure, for example, the field of DFT (Jensen and Wasserman 2018; Parr and
Yang 1995; Burke, Werschnik, and Gross 2005) has, at its core, this type of
inverse problems. As discussed before, solutions for the full Schrödinger’s
equation, and thus wave–functions, are difficult to obtain except for some
simple models; the probability density $|\psi(\vec{x})|^{2}$, on the other
hand, may be experimentally inferred for several quantum systems. Thus, an
approximated wave–function, $|\psi|=\sqrt{|\psi(\vec{x})|^{2}}$, may be
defined for the construction of the effective potential.

##  3 Quantum Potential Neural Networks

In this section the proposed NN and the loss function that were used to
compute the effective potentials of various quantum systems are described.
This section is divided into two parts: (i) Time–independent systems and (ii)
Time–dependent systems

###  3.1 Time–independent Systems

In order to obtain the effective potential, a new parametric function
$U_{\theta}$ is learned in an unsupervised manner. This parametric function
corresponds to the effective potential of the quantum system and was obtained
by implementing a loss function that obeys the time–independent Schrödinger
(TISE) equation (Eq. 2),

$L_{TISE}(\theta)=\bigg{\lvert}\bigg{\rvert}D\left(-\frac{\hbar^{2}}{2m}\frac{\Delta_{\vec{x}}|\psi|}{|\psi|}+U_{\theta}(\vec{x})\right)\bigg{\lvert}\bigg{\rvert}_{2}^{2},$
(5)

where $D$ is the total derivative operator acting on multi-variate function
$-\frac{\hbar^{2}}{2m}\frac{\Delta_{\vec{x}}|\psi|}{|\psi|}+U_{\theta}(\vec{x})$
and $\lvert\rvert\cdot\lvert\rvert_{2}$ is the Frobenius norm. Because of the
definition of this loss function, energy conservation is effectively demanded
for time–independent systems. Since the $U_{\theta}$ function is given by a
differential equation, an initial condition was added to ensure that a unique
function is learned. The initial conditions used for the different quantum
systems are based on the inherent nature of each of the systems, more
information and explanations about some of these conditions can be found in
the literature (Romanowski 2007). Finally, the loss function after the
consideration of the initial condition reads

$L(\theta)=L_{TISE}(\theta)+(U_{\theta}(\vec{x})-y)^{2}$ (6)

where $\vec{x}$ is some point in the domain of the function and $y$ is the
expected ground truth value of the true potential at that point.
The main observation here is that using $|\psi|$ (instead of $\psi$) to solve
for the potential leads to the correct potential except, possibly, at finitely
many points where $\psi$ changes signs. However, this does not create any
difficulties for training the proposed model.

###  3.2 Time–dependent Systems

For time–dependent systems, the formulation for the time–dependent Schrödinger
equation (TDSE) loss reads

$L_{TDSE}(\theta)=\bigg{|}\bigg{|}Re\bigg{(}\frac{i\frac{\partial\psi}{\partial
t}+\frac{\hbar^{2}}{2m}\frac{\partial^{2}\psi}{\partial\vec{x}^{2}}}{\psi}\bigg{)}-U_{\theta}\bigg{|}\bigg{|}_{2}^{2}.$
(7)

It is important to mention that complex numbers may be more common to appear
in the time–dependent solution of the Schrödinger equation; however, for the
current study, the probability density of the considered systems are described
with an Hermitian Hamiltonian, and thus only real observables were considered
to avoid the handling of complex values. For the time–dependent results
presented in this report, the QPNN was trained with the full wave–function
instead of just the probability density. In a future work, the
density–to–potential results for time–dependent systems will be explored and
discussed in detail.

###  3.3 Model Architecture

For the construction of the NN, a $4$ layer feedforward network with hidden
sizes of $32,128$ and $128$ with a residual connection between second and
third layers was used. The input to the model are the $\vec{x}$ spatial
coordinates for time independent systems and ($\vec{x},t$) for time dependent
systems. For the network training, $3000$ of these coordinates were randomly
selected from the domain of definition of each particular system. The model
was trained for $500$ epochs with Adam optimizer (Kingma and Ba 2017).

Table 1: A quantitative analysis for the QPNN

System |  | RMSE between True and Learned
---
Potentials (QPNN)
| RMSE between True and Learned
---
Potentials (using RK4)
| RMSE between True and Learned
---
Energies (QPNN)
| RMSE between True and Learned
---
Energies (RK4)
Harmonic Oscillator | $$1\text{\times}{10}^{-2}$\pm$5\text{\times}{10}^{-3}$$ | $$9\text{\times}{10}^{-3}$\pm$4\text{\times}{10}^{-4}$$ | $$1\text{\times}{10}^{-2}$\pm$2\text{\times}{10}^{-3}$$ | $$5\text{\times}{10}^{-2}$\pm$7\text{\times}{10}^{-3}$$
Pöschl–Teller potential | $$1\text{\times}{10}^{-4}$\pm$6\text{\times}{10}^{-5}$$ | $$2\text{\times}{10}^{-4}$\pm$3\text{\times}{10}^{-5}$$ | $$8\text{\times}{10}^{-4}$\pm$6\text{\times}{10}^{-5}$$ | $$7\text{\times}{10}^{-3}$\pm$8\text{\times}{10}^{-4}$$
$\rm H_{2}$ molecule | $$2\text{\times}{10}^{-3}$\pm$4\text{\times}{10}^{-4}$$ | $$3\text{\times}{10}^{-3}$\pm$2\text{\times}{10}^{-4}$$ | $$9\text{\times}{10}^{-3}$\pm$7\text{\times}{10}^{-4}$$ | $$4\text{\times}{10}^{-3}$\pm$2\text{\times}{10}^{-4}$$
Soliton | $$3\text{\times}{10}^{-2}$\pm$4\text{\times}{10}^{-3}$$ | - | - | -

##  4 Related Work

The use of deep learning for understating physical phenomena has been an
active field of development. In particular, there is a considerable amount of
literature where authors have endowed neural networks with classical
Hamiltonian mechanics (Toth et al. 2019; Greydanus, Dzamba, and Yosinski 2019;
Tong et al. 2020; Iten et al. 2020; Bondesan and Lamacraft 2019; Zhong, Dey,
and Chakraborty 2019; Chmiela et al. 2017); conservation of energy and
irreversibility in time are the key features of such networks. There are
recent reports extending these results in cases of damped pendula, i.e.,
systems where there is dissipation of energy (Zhong, Dey, and Chakraborty
2020). In computational quantum mechanics, deep neural networks have been
implemented to learn representations and extract the necessary features to
predict desired properties from raw unprocessed data (Goh, Hodas, and Vishnu
2017). Recently, two methods for estimating the density matrix for a quantum
system, the Quantum Maximum Likelihood (QML) and Quantum Variational Inference
(QVI) method, were introduced (Cranmer, Golkar, and Pappadopulo 2019). For
these methods, the authors used a flow based method (Toth et al. 2019; Rezende
and Mohamed 2015) to increase the expressivity of their variational family of
density matrices. The applicability of these methods, however, has been only
validated for the harmonic and anharmonic quantum oscillator models.
Applications of deep learning to quantum mechanics is still in its early
stages (Torfi et al. 2020; Raissi, Perdikaris, and Karniadakis 2017a, b;
Raissi, Perdikaris, and Karniadakis 2019; Dai et al. 2020; Carleo et al. 2019;
Amabilino et al. 2019; Unke and Meuwly 2019; Schmitz, Godtliebsen, and
Christiansen 2019; Schmidt et al. 2017; Hibat-Allah et al. 2020; Nakajima,
Tanaka, and Hashimoto 2020; Pu, Li, and Chen 2020; Mills, Spanner, and Tamblyn
2017; Manzhos 2020). Most of the deep learning quantum mechanic frameworks
introduced so far are focused on either solving the Schrödinger equation or
predicting the trends of specific observables such as the system’s energy.
Concerning inverse problems, Raissi et al, introduced the physics–informed
neural network for solving forward and inverse problems involving nonlinear
partial differential equations (Raissi, Perdikaris, and Karniadakis 2019).
Although the impressive results reported in this work, in terms of inverse
problem solutions, the deep learning framework reported by Raissi et al is
focused only on the solution of a partial differential equation for the
prediction of pressure profiles in an classical system. On the other hand, in
quantum mechanics observables may be inferred when a valid effective potential
is known for a given quantum system; thus, the solution of the
density–to–potential inversion problem to predict effective potential
functions (Jensen and Wasserman 2018) play an important role in the
understanding of the quantum phenomena, and in particular in the elucidation
of the electronic structure of molecules from a density functional theory
perspective. In this regard, Nagai and coworkers have proposed the
Neural–network Kohn–Sham exchange–correlation potential (Nagai et al. 2018),
which propose a supervised training scheme that uses information from well
defined potentials and probability densities to train a NN. However, to the
best of our knowledge, there are no reported works that uses deep learning to
solve the density–to–potential inverse problem to systematically estimate
potentials from observations in a completely unsupervised manner.

##  5 Experiments

The performance of the proposed Quantum Potential Neural Network is validated
on four different quantum systems, one of these systems describes the temporal
evolution of a quantum wave whereas the other three are examples of
time–independent systems. Among the time–independent systems, exact analytical
solutions for the time–independent Schrödinger equation can be obtained only
for the harmonic oscillator and the Pöschl–Teller (PT) potentials whereas for
the $\text{H}_{2}$ molecule, only approximate solutions are attained. Details
about the solutions, physical implications and interpretations of these
systems can be found in any standard book on quantum mechanics (Sakurai and
Commins 1995; Robinett 1997; Pronchik and Williams 2003). Finally, the
potentials learned by the QPNN were used to compute the total energy of each
of the systems. The quantitative results for all the reported experiments are
summarized in table 1. In order to compare the results obtained (in the form
of solutions to a differential equation) via the NN techniques to those
obtained through well established approaches, the differential equations were
solved numerically using the well–known and standard Runge–Kutta 4th Order (or
RK4) integrator implemented in the standard python libraries. The RK4
algorithm provides means of solving various kinds of differential equations
and is generally considered as a robust workhorse to bench mark new
computational techniques (Landau, Paez, and Bordeianu 2015). The differences
in accuracy of the values obtained by both the proposed QPNN method and the
RK4 numerical integrator were quantified through their root mean square error
(RMSE) values.

###  5.1 Density–to–Potential Experiments

For the density–to–potential experiments, the exact wave–functions,
$\psi(\vec{x})$, for the quantum Harmonic oscillator and the PT potential were
obtained by solving the time–independent Schrödinger equation. These
wave–functions were later used to define the probability distribution,
$|\psi(\vec{x})|^{2}$, for each of the systems. In the case of the Hydrogen
molecule, the probability density was defined according to the one–electron
$1s$ orbital function that delineate the approximated electronic density for
the Hydrogen molecule in a Born–Oppenheimer approximation. These probability
densities were used to define, for each of the quantum systems, an
approximated probability amplitude function,
$|\psi|=\sqrt{|\psi(\vec{x})|^{2}}$. The effective potential function for each
of the systems was obtained by training the QPNN with information provided by
randomly selected coordinates evaluated onto the approximated probability
amplitudes and into the loss function.

#### Quantum Harmonic Oscillator (QHO):

The motion of the the one–dimension QHO is, perhaps, the simplest quantum
mechanical system whose motion follows a linear differential equation with
constant coefficients. In the QPNN framework, for the prediction of the QHO
potential, the coordinate variable $\vec{x}$, was randomly sampled from
$[-5,5]$, and was used as the input to the model. In the analytical solution
of the time–independent Schrödinger equation, the wave–functions for the
different states of the QHO are given by Hermite polynomials
$H_{n},n=0,1,\cdots$, whereas the energies corresponding to these states
depends on the force constant $w$ and are given by $E_{n}=\hbar w(n+1/2)$. The
analytical wave–functions for the different states of the QHO defined the
probability densities used to train the QPNN. In this case, the initial
condition imposed is the fact that at the zero point of the reference
coordinates, all the energy in the system is kinetic, and thus the potential
energy at this point is zero, i.e. $V(0)=0$. Fig 2 shows the used probability
density, the learned potential and the energy computed by the QPNN.

Figure 2: Harmonic Oscillator system. (A) Wave–function and probability
Density used in the QPNN, (B) Ground and Learned Potential, (C) Total Energy

#### Pöschl-Teller potential:

The PT potential is a special class of anharmonic potentials for which the
one–dimensional Schrödinger equation can be solved in terms of special
functions. This potential may be used to model vibrational molecular
potentials with out–of–plane bending vibrations (Senn 1986; Jia, Zhang, and
Peng 2017) and observables of diatomic potentials (Pronchik and Williams
2003). For the PT potential, the wave–functions used to define the
approximated probability amplitude function are the Legendre functions
$P^{\mu}_{\lambda}$ (Riley 1974) with energy eigenvalues $E_{\mu}$ and
potential depth $V_{0}$(Hernández de la Peña 2018),

$\Big{\\{}P^{\mu}_{\lambda}(\tanh
x)\;\Big{|}\;E_{\mu}=\frac{-\mu^{2}}{2},V_{0}=\frac{-\lambda(\lambda+1)}{2},\;\begin{subarray}{c}\lambda=1,2,\cdots\\\
\;\mu=1,2,\cdots,\lambda\end{subarray}\Big{\\}}.$ (8)

Details about several suitable boundary terms and initial conditions for this
type of potentials are formulated and reviewed in the literature (Agboola
2010; Hernández de la Peña 2018). The input for this experiment is the spatial
coordinate $x$ randomly sampled from $[-3,3]$. For this experiment, the
wave–function defined by the Legendre function with $\lambda=2$ and $\mu=1$
was employed. From the density, the initial wave–function takes the form;

$\psi^{1}_{2}\left(x\right)=|\text{tanh }x|\text{ sech }x$ (9)

Fig 3 shows the probability density used to train the system, as well as the
learned potential and energy computed for the system.

Figure 3: Pöschl Teller system. (A) Wave–function and Density, (B) Ground and
Learned Potential, (C) Total Energy Figure 4: $\text{H}_{2}$ molecule system.
(A) Density, (B) Ground Truth Potential, (C,D) Learned Potential, (E) Computed
Energy

#### Hydrogen molecule:

The $\rm H_{2}$ molecule is the first multi-electronic system with
approximated probability density considered. For the training of the QPNN, an
ab initio electron density for the $\rm H_{2}$ molecule with an equilibrium
bond length ($x_{r_{e}}$) of 1.346Å and total energy of -0.958470046928 a.u.
was approximated by using a fast and systematic self–consistent field method
(Helgaker, Jorgensen, and Olsen 2014). This density was computed using three
Gaussian primitive functions for each $\rm H$ atom where the $\vec{x}$
coordinate defined on $[-3,3]$ was chosen as the reference internal
coordinate. The initial conditions for the system were defined following the
same lines as in (Rafi et al. 1995); specifically, the fact that
$V\left(x_{r_{e}}\right)=0$. Fig 4 shows the probability density used to train
the QPNN, as well as the learned potential and the energy computed using this
potential.

###  5.2 Exploration of a Time–Dependent System

In order to explore the behavior of the QPNN in systems with dependence on
time, the effective potential for a soliton model was computed. Solitons
represent solitary waves propagating without any temporal evolution in shape
or size when viewed in the reference frame moving with the group velocity of
the waves (Wazwaz 2009). This type of solitary waves are particularly
important in the Bose–Einstein condensation theory and arise in many contexts
such as the elevation of the surface of water and the intensity of light in
optical fibers. Solitons form a special class of solutions of model equations,
including the Korteweg de–Vries (KdV) and the Nonlinear Schrödinger (NLS)
equations. In this particular experiment, the one–dimensional soliton
satisfies the following differential equation:

$i\frac{\partial\psi}{\partial t}+\frac{\partial^{2}\psi}{\partial
x^{2}}+U(x,t)\psi=0.$ (10)

Thus, for this system, the loss function used to train the QPNN is given by
equation 7 where $\psi=2\text{sech}(\sqrt{2}(x-2t))e^{i(x+t)}$ and $U(x,t)$ is
$|\psi|^{2}$. For this experiment, the coordinates for the QPNN input were
defined on $\vec{x},\vec{t}\in[0,1]$. Fig 5 shows the potential obtained by
the QPNN for this system.

Figure 5: Left: Soliton Ground Truth Potential, Right: Soliton Learned
Potential

##  6 Discussion

Table 1 reports the RMSE values for the different models studied in this work.
In the case of the time–independent models, for the QPNN, the harmonic
oscillator presents the largest RMSE value between the true and learned
potential whereas the PT potential has the lowest RMSE among all the
time–independent systems. For all the systems but the harmoninc oscillator,
the RMSE values for the learned potentials are comparable in magnitude with
those obtained with the RK4 method. In the case of the energies, when compared
against the exact energy, the RMSE values for both the QPNN and the RK4 method
are around the same magnitude for all systems but the PT potential, where the
RMSE for the QPNN is one order of magnitude lower than for the RK4 method. In
terms of the energy values computed for each of the systems, when referenced
to the exact energy (blue line) a similar trend can be observed for the
energies computed with the RK4 method and the QPNN. If the RK4 method is
regarded as a more robust and mathematically superior method for the
calculation of the effective potential function, this trend may be interpreted
as an indicator of the reliability of the QPNN. In the case of the soliton
model, although solutions using the RK4 were not feasible due to the nature of
the system, the magnitude of RMSE values between the true potential and the
learned potential suggests the same qualitative behaviour as the one obtained
with the time–independent models.

##  7 Conclusion

In this work, the QPNN for learning the effective potential functions of
different quantum systems was presented. This new neural network is capable of
learning the effective potential functions of a variety of systems in a
completely unsupervised manner. The results obtained for the different studied
systems suggest that the QPNN has an accuracy comparable to the RK4
integrator. The potentials learned with this new QPNN can be used to calculate
observables like the energy of the system. The use of the QPNN was also
extended for time–dependent systems where the full wave–function was used. In
the time–dependent case, $|\psi|$ can not simply be taken as a proxy for a
wave–function as one needs to take into account some phase information, i.e.,
the wave–function in that case can not be a real valued function. In future
work, the density–to–potential problem will be analysed by incorporating phase
information to create a suitable proxy wave–function.

## References

  * Agboola (2010) Agboola, D. 2010.  Solutions to the Modified Pöschl–Teller Potential in D-Dimensions.  _Chinese Physics Letters_ 27(4): 040301.
  * Amabilino et al. (2019) Amabilino, S.; Bratholm, L. A.; Bennie, S. J.; Vaucher, A. C.; Reiher, M.; and Glowacki, D. R. 2019.  Training neural nets to learn reactive potential energy surfaces using interactive quantum chemistry in virtual reality.  _The Journal of Physical Chemistry A_ 123(20): 4486–4499.
  * Aster, Borchers, and Thurber (2018) Aster, R. C.; Borchers, B.; and Thurber, C. H. 2018.  _Parameter estimation and inverse problems_.  Elsevier.
  * Athanassoulis (2008) Athanassoulis, A. G. 2008.  Exact equations for smoothed Wigner transforms and homogenization of wave propagation.  _Applied and Computational Harmonic Analysis_ 24(3): 378–392.
  * Beals and Greiner (2009) Beals, R.; and Greiner, P. C. 2009.  Strings, waves, drums: spectra and inverse problems.  _Analysis and Applications_ 7(02): 131–183.
  * Bondesan and Lamacraft (2019) Bondesan, R.; and Lamacraft, A. 2019.  Learning Symmetries of Classical Integrable Systems.
  * Burke, Werschnik, and Gross (2005) Burke, K.; Werschnik, J.; and Gross, E. K. U. 2005.  Time-dependent density functional theory: Past, present, and future.  _The Journal of Chemical Physics_ 123(6): 062206.  ISSN 1089-7690.  doi:10.1063/1.1904586.  URL http://dx.doi.org/10.1063/1.1904586.
  * Carleo et al. (2019) Carleo, G.; Cirac, I.; Cranmer, K.; Daudet, L.; Schuld, M.; Tishby, N.; Vogt-Maranto, L.; and Zdeborová, L. 2019.  Machine learning and the physical sciences.  _Reviews of Modern Physics_ 91(4).  ISSN 1539-0756.  doi:10.1103/revmodphys.91.045002.  URL http://dx.doi.org/10.1103/RevModPhys.91.045002.
  * Case (2008a) Case, W. B. 2008a.  Wigner functions and Weyl transforms for pedestrians.  _American Journal of Physics_ 76(10): 937–946.  doi:10.1119/1.2957889.  URL https://doi.org/10.1119/1.2957889.
  * Case (2008b) Case, W. B. 2008b.  Wigner functions and Weyl transforms for pedestrians.  _American Journal of Physics_ 76(10): 937–946.
  * Chadan and Sabatier (2012) Chadan, K.; and Sabatier, P. C. 2012.  _Inverse problems in quantum scattering theory_.  Springer Science & Business Media.
  * Chen, Xiong, and Shao (2018) Chen, Z.; Xiong, Y.; and Shao, S. 2018.  Numerical Methods for the Wigner Equation with Unbounded Potential.  _Journal of Scientific Computing_ 79(1): 345–368.  ISSN 1573-7691.  doi:10.1007/s10915-018-0853-0.  URL http://dx.doi.org/10.1007/s10915-018-0853-0.
  * Chen, Xiong, and Shao (2019) Chen, Z.; Xiong, Y.; and Shao, S. 2019.  Numerical methods for the Wigner equation with unbounded potential.  _Journal of Scientific Computing_ 79(1): 345–368.
  * Chmiela et al. (2017) Chmiela, S.; Tkatchenko, A.; Sauceda, H. E.; Poltavsky, I.; Schütt, K. T.; and Müller, K.-R. 2017.  Machine learning of accurate energy-conserving molecular force fields.  _Science advances_ 3(5): e1603015.
  * Cranmer, Golkar, and Pappadopulo (2019) Cranmer, K.; Golkar, S.; and Pappadopulo, D. 2019.  Inferring the quantum density matrix with machine learning.
  * Cranmer et al. (2020) Cranmer, M.; Greydanus, S.; Hoyer, S.; Battaglia, P.; Spergel, D.; and Ho, S. 2020\.  Lagrangian Neural Networks.
  * Curtright, Fairlie, and Zachos (1998) Curtright, T.; Fairlie, D.; and Zachos, C. 1998.  Features of time-independent Wigner functions.  _Physical Review D_ 58(2): 025002.
  * Dai et al. (2020) Dai, J.; Jasinski, A.; Montaner, J.; Forrey, R.; Yang, B.; Stancil, P.; Balakrishnan, N.; Vargas-Hernández, R.; and Krems, R. 2020.  Machine-learning-corrected quantum dynamics calculations.  _Bulletin of the American Physical Society_ .
  * Du and Narasimhan (2019) Du, Y.; and Narasimhan, K. 2019.  Task-Agnostic Dynamics Priors for Deep Reinforcement Learning.
  * Feynman, Leighton, and Sands (1965) Feynman, R.; Leighton, R.; and Sands, M. 1965.  The Feynman Lectures on Physics Vol. III, chap. 21, sec. 21-9.
  * Galleani and Cohen (2002) Galleani, L.; and Cohen, L. 2002.  Approximation of the Wigner distribution for dynamical systems governed by differential equations.  _EURASIP Journal on Advances in Signal Processing_ 2002(1): 514609\.
  * Goh, Hodas, and Vishnu (2017) Goh, G. B.; Hodas, N. O.; and Vishnu, A. 2017.  Deep learning for computational chemistry.  _Journal of computational chemistry_ 38(16): 1291–1307.
  * Gomes and Silva (2008) Gomes, D. A.; and Silva, J. D. 2008.  On the Wigner transform of solutions to the Schrodinger equation.  _São Paulo Journal of Mathematical Sciences_ 2(1): 85–97.
  * Greydanus, Dzamba, and Yosinski (2019) Greydanus, S.; Dzamba, M.; and Yosinski, J. 2019.  Hamiltonian neural networks.  In _Advances in Neural Information Processing Systems_ , 15353–15363.
  * Groetsch and Groetsch (1993) Groetsch, C. W.; and Groetsch, C. 1993.  _Inverse problems in the mathematical sciences_ , volume 52.  Springer.
  * Helgaker, Jorgensen, and Olsen (2014) Helgaker, T.; Jorgensen, P.; and Olsen, J. 2014.  _Molecular electronic-structure theory_.  John Wiley & Sons.
  * Heller (1976) Heller, E. J. 1976.  Wigner phase space method: Analysis for semiclassical applications.  _The Journal of Chemical Physics_ 65(4): 1289–1298.
  * Hernández de la Peña (2018) Hernández de la Peña, L. 2018.  A Simplified Pöschl–Teller Potential: An Instructive Exercise for Introductory Quantum Mechanics.  _Journal of Chemical Education_ 95(11): 1989–1995.  doi:10.1021/acs.jchemed.8b00029.
  * Hibat-Allah et al. (2020) Hibat-Allah, M.; Ganahl, M.; Hayward, L. E.; Melko, R. G.; and Carrasquilla, J. 2020\.  Recurrent Neural Network Wavefunctions.
  * Higgins et al. (2016) Higgins, I.; Matthey, L.; Glorot, X.; Pal, A.; Uria, B.; Blundell, C.; Mohamed, S.; and Lerchner, A. 2016.  Early Visual Concept Learning with Unsupervised Deep Learning.
  * Iten et al. (2020) Iten, R.; Metger, T.; Wilming, H.; del Rio, L.; and Renner, R. 2020.  Discovering Physical Concepts with Neural Networks.  _Physical Review Letters_ 124(1).  ISSN 1079-7114.  doi:10.1103/physrevlett.124.010508.  URL http://dx.doi.org/10.1103/PhysRevLett.124.010508.
  * Jensen and Wasserman (2018) Jensen, D. S.; and Wasserman, A. 2018.  Numerical methods for the inverse problem of density functional theory.  _International Journal of Quantum Chemistry_ 118(1): e25425.
  * Jia, Zhang, and Peng (2017) Jia, C.-S.; Zhang, L.-H.; and Peng, X.-L. 2017.  Improved Pöschl–Teller potential energy model for diatomic molecules.  _International Journal of Quantum Chemistry_ 117(14): e25383.  doi:https://doi.org/10.1002/qua.25383.
  * Kingma and Ba (2017) Kingma, D. P.; and Ba, J. 2017.  Adam: A Method for Stochastic Optimization.
  * Klimov, Sainz, and Romero (2020) Klimov, A.; Sainz, I.; and Romero, J. 2020.  Truncated Wigner approximation as non-positive Kraus map.  _Physica Scripta_ 95(7): 074006.
  * Kolesnikov et al. (2019) Kolesnikov, A.; Beyer, L.; Zhai, X.; Puigcerver, J.; Yung, J.; Gelly, S.; and Houlsby, N. 2019.  Big Transfer (BiT): General Visual Representation Learning.
  * Lample and Charton (2019) Lample, G.; and Charton, F. 2019.  Deep Learning for Symbolic Mathematics.
  * Landau, Paez, and Bordeianu (2015) Landau, R.; Paez, M. J.; and Bordeianu, C. 2015.  _Computational Physics: Problem Solving with Python_.  Wiley, 3rd edition.
  * Manzhos (2020) Manzhos, S. 2020.  Machine learning for the solution of the Schrödinger equation.  _Machine Learning: Science and Technology_ 1(1): 013002.
  * Mills, Spanner, and Tamblyn (2017) Mills, K.; Spanner, M.; and Tamblyn, I. 2017.  Deep learning and the Schrödinger equation.  _Physical Review A_ 96(4): 042113.
  * Nagai et al. (2018) Nagai, R.; Akashi, R.; Sasaki, S.; and Tsuneyuki, S. 2018.  Neural-network Kohn-Sham exchange-correlation potential and its out-of-training transferability.  _The Journal of chemical physics_ 148(24): 241737.
  * Nakajima, Tanaka, and Hashimoto (2020) Nakajima, M.; Tanaka, K.; and Hashimoto, T. 2020.  Neural Schrödinger Equation: Physical Law as Neural Network.  _arXiv preprint arXiv:2006.13541_ .
  * Nakatsuji (2002) Nakatsuji, H. 2002.  Inverse Schrödinger equation and the exact wave function.  _Physical Review A_ 65(5): 052122.
  * Parr and Yang (1995) Parr, R. G.; and Yang, W. 1995.  Density-functional theory of the electronic structure of molecules.  _Annual Review of Physical Chemistry_ 46(1): 701–728.
  * Pronchik and Williams (2003) Pronchik, J. N.; and Williams, B. W. 2003.  Exactly Solvable Quantum Mechanical Potentials: An Alternative Approach.  _Journal of Chemical Education_ 80(8): 918.  doi:10.1021/ed080p918.
  * Pu, Li, and Chen (2020) Pu, J.; Li, J.; and Chen, Y. 2020.  Soliton, Breather and Rogue Wave Solutions for Solving the Nonlinear Schrödinger Equation Using a Deep Learning Method with Physical Constraints.  _arXiv preprint arXiv:2011.04949_ .
  * Rafi et al. (1995) Rafi, M.; et al. 1995.  An empirical potential function of diatomic molecules.  _Physics Letters A_ 205(5-6): 383–387.
  * Raissi, Perdikaris, and Karniadakis (2017a) Raissi, M.; Perdikaris, P.; and Karniadakis, G. E. 2017a.  Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations.  _arXiv preprint arXiv:1711.10561_ .
  * Raissi, Perdikaris, and Karniadakis (2017b) Raissi, M.; Perdikaris, P.; and Karniadakis, G. E. 2017b.  Physics Informed Deep Learning (Part II): Data-driven Discovery of Nonlinear Partial Differential Equations.  _arXiv preprint arXiv:1711.10566_ .
  * Raissi, Perdikaris, and Karniadakis (2019) Raissi, M.; Perdikaris, P.; and Karniadakis, G. E. 2019.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ 378: 686–707.
  * Rezende and Mohamed (2015) Rezende, D. J.; and Mohamed, S. 2015.  Variational Inference with Normalizing Flows.
  * Riley (1974) Riley, K. F. 1974.  _Mathematical Methods for the Physical Sciences: An Informal Treatment for Students of Physics and Engineering_.  Cambridge University Press.  doi:10.1017/CBO9781139167550.
  * Robinett and Robinett (2006) Robinett, R.; and Robinett, R. W. 2006.  _Quantum mechanics: Classical results, modern systems, and visualized examples_.  Oxford University Press.
  * Robinett (1997) Robinett, R. W. 1997.  _Quantum mechanics_.  Oxford University Press, New York.
  * Romanowski (2007) Romanowski, Z. 2007.  Numerical Solution Of Kohn-Sham Equation For Atom.  _Acta Physica Polonica B_ 38(10).
  * Sakurai and Commins (1995) Sakurai, J. J.; and Commins, E. D. 1995.  Modern Quantum Mechanics, Revised Edition.  _AAPT_ .
  * Schleder et al. (2019) Schleder, G. R.; Padilha, A. C.; Acosta, C. M.; Costa, M.; and Fazzio, A. 2019.  From DFT to machine learning: recent approaches to materials science–a review.  _Journal of Physics: Materials_ 2(3): 032001.
  * Schmidt et al. (2017) Schmidt, J.; Shi, J.; Borlido, P.; Chen, L.; Botti, S.; and Marques, M. A. 2017\.  Predicting the thermodynamic stability of solids combining density functional theory and machine learning.  _Chemistry of Materials_ 29(12): 5090–5103.
  * Schmitz, Godtliebsen, and Christiansen (2019) Schmitz, G.; Godtliebsen, I. H.; and Christiansen, O. 2019.  Machine learning for potential energy surfaces: An extensive database and assessment of methods.  _The Journal of chemical physics_ 150(24): 244113.
  * Senn (1986) Senn, P. 1986.  The modified Poschl-Teller Oscillator.  _Journal of Chemical Education_ 63(1): 75.
  * Sun et al. (2019) Sun, C.; Myers, A.; Vondrick, C.; Murphy, K.; and Schmid, C. 2019.  VideoBERT: A Joint Model for Video and Language Representation Learning.
  * Tong et al. (2020) Tong, Y.; Xiong, S.; He, X.; Pan, G.; and Zhu, B. 2020.  Symplectic Neural Networks in Taylor Series Form for Hamiltonian Systems.
  * Torfi et al. (2020) Torfi, A.; Shirvani, R. A.; Keneshloo, Y.; Tavvaf, N.; and Fox, E. A. 2020.  Natural Language Processing Advancements By Deep Learning: A Survey.
  * Toth et al. (2019) Toth, P.; Rezende, D. J.; Jaegle, A.; Racanière, S.; Botev, A.; and Higgins, I. 2019.  Hamiltonian Generative Networks.
  * Unke and Meuwly (2019) Unke, O. T.; and Meuwly, M. 2019.  Machine Learning Potential Energy Surfaces.  _arXiv: Chemical Physics_ .
  * Vogel (2002) Vogel, C. R. 2002.  _Computational methods for inverse problems_ , volume 23.  Siam.
  * Wazwaz (2009) Wazwaz, A. 2009.  _Partial Differential Equations and Solitary Waves Theory_.  Springer Berlin Heidelberg.  doi:10.1007/978-3-642-00251-9.
  * Zakhariev and Suzko (2012) Zakhariev, B. N.; and Suzko, A. A. 2012.  _Direct and inverse problems: potentials in quantum scattering_.  Springer Science & Business Media.
  * Zhong, Dey, and Chakraborty (2019) Zhong, Y. D.; Dey, B.; and Chakraborty, A. 2019.  Symplectic ODE-Net: Learning Hamiltonian Dynamics with Control.
  * Zhong, Dey, and Chakraborty (2020) Zhong, Y. D.; Dey, B.; and Chakraborty, A. 2020.  Dissipative SymODEN: Encoding Hamiltonian Dynamics with Dissipation and Control into Deep Learning.

##  Appendix A Additional experiments using the full wave–function

In order to further explore the capabilities and accuracy of our NN, we start
off by describing additional experiments using the QPNN and full wave–function
as well as Wigner’s functions. The modifications to our loss function is
described in detail below.

Table 2: Wave–functions, Energies, and Potentials of various time independent
systems

System | Potential $V(x)$ | wave–function $\psi(x)$ | Energy
---|---|---|---
Harmonic Oscillator | $\frac{1}{2}kx^{2}$ | ${\frac{1}{\sqrt{2^{n}\,n!}}}\left({\frac{m\omega}{\pi\hbar}}\right)^{1/4}e^{-{\frac{m\omega x^{2}}{2\hbar}}}H_{n}\left({\sqrt{\frac{m\omega}{\hbar}}}x\right)$ | $\hbar\omega(n+\frac{1}{2})$
Pöchl–Teller potential | $-\frac{\lambda(\lambda+1)}{2}\text{sech}^{2}(x)$ | $P_{\lambda}^{\mu}(\text{tanh}(x))$ | $-\frac{\hbar^{2}}{2m}(\lambda-\mu)$
Radial Hydrogen atom | $\frac{l(l+1)}{2r^{2}}-\frac{1}{r}$ | $e^{-r/n}(\frac{2r}{na_{0}})^{l}L^{2l+1}_{n+l}(\frac{2r}{na_{0}})$ | $-\frac{R_{H}}{(n+l)^{2}}$
2D Harmonic Oscillator | $\frac{1}{2}k(x^{2}+y^{2})$ | $H_{n_{x}}\left({\sqrt{\frac{m\omega}{\hbar}}}x\right)H_{n_{y}}\left({\sqrt{\frac{m\omega}{\hbar}}}y\right)e^{-{\frac{m\omega(x^{2}+y^{2})}{2\hbar}}}$ | $\hbar\omega(n_{x}+n_{y}+1)$

###  A.1 Using the 1D time–independent Schrödinger equation

In this section we consider some simple one–dimensional time–independent
systems. Wave–functions, potential energy and energy levels can be found in
Table 2. We report our learned potentials in figure 6 and show that our models
obey energy conservation laws in figure 7. The quantitative results for the
experiments can be found in table 3. For the derivation of these
wave–functions and the general properties of these systems, please see
(Sakurai and Commins 1995; Feynman, Leighton, and Sands 1965; Robinett 1997).
As before, $\hbar$, $m$ and $\omega$ were set equal to $1$.

Figure 6: Position (x-axis) vs Potential (y-axis); Ground (Blue) and Learned
(Green Dots) Potentials: From left to right: Harmonic Oscillator, Hydrogen
Atom (2p case), Pöschl-Teller (1,1 case)

.

Quantum Harmonic Oscillator: The wave–functions in this case are given by
Hermite polynomials $H_{n},n=0,1,\cdots$. We choose $x\in[-5,5]$ as input to
our model. Since $U_{\theta}$ is given by a differential equation (equation 5)
one needs to impose an initial condition to get a unique solution. However,
constraining the output of $U_{\theta}\in[0,12.5]$ removes the need for the
initial condition. Fig 6 and fig 7 (left) shows the learned potential and
energy of the system.
The Hydrogen Atom (2p case) : The general radial wave–functions are given by
generalized Laguerre polynomials $L^{l}_{n},n=1,2,\cdots\ \text{and}\
l=0,1,2,\cdots,n-1$ but in this case simplifies to
$\psi(r)=\frac{1}{8\sqrt{\pi}}re^{\frac{-r}{2}}$. We used $r\in[0.5,10]$ as
input to our model, the initial condition $U(1)=0.$ and the loss function
$\mathcal{L}(\theta)=L(\theta)+U_{\theta}(1)^{2}$ where $L(\theta)$ is given
by equation 5. Fig 6 and fig 7 (middle) shows the learned potential and energy
of the system.
Pöschl-Teller potential : The wave–function $\psi$ generated by this potential
is defined by Legendre functions
$P_{\lambda}^{\mu}(\text{tanh}(x)),\lambda=1,2,3\cdot;\
\mu=1,2,\cdot,\lambda-1,\lambda$. For simplicity, let $\mu=1$. We choose
$x\in[-3,3]$ as input to our model. We imposed an initial condition
$U_{\theta}(0)=-\frac{\lambda(\lambda+1)}{2}$ and used a similar auxiliary
loss function as above. Fig 6 and fig 7 (right) shows the learned potential
and energy of the system.

Figure 7: Position (x-axis) vs Energy (y-axis); Ground (Blue) and Learned
(Green Dots) Energies: From left to right: Harmonic Oscillator, Hydrogen Atom
(2p case), Pöschl-Teller (1,1 case).

###  A.2 Particle in a box (perturbed by some external potential)

Now, we turn our attention to a quantum system where the Schrödinger equation
cannot be solved exactly, but can be formulated in an approximate manner using
perturbation theory. A particle with no spin, of mass m, was placed in a
square one dimensional box, $x\in[0,L]$, of length $L$. Later the particle was
presented with the perturbation $V(x)=10x^{2}$. The wave–function for the
perturbed system was approximated by considering first order corrections for
the unperturbed particle in a box wave–function,

$\psi_{n}=\psi^{0}_{n}+\sum_{n\neq
k}\dfrac{<\psi^{0}_{n}|V(x)|\psi^{0}_{k}>}{E^{0}_{n}-E^{0}_{k}}\psi^{0}_{k},\quad
n,k=1,2,3,\ldots,\quad$ (11)

where $\psi^{0}_{n}$, and $E^{0}_{n}$ are the unperturbed particle in a box
$n$th state wave–function and its energy, whereas,
$<\psi^{0}_{n}|V(x)|\psi^{0}_{k}>$ indicates the following integral

$<\psi^{0}_{n}|V(x)|\psi^{0}_{k}>=\int{(\psi^{0}_{n})^{*}V(x)\psi^{0}_{k}dx}.$
(12)

For our computations, the wave–function, $\psi^{0}_{n}$, obtained as solution
of the Schröedinger equation for the particle in a box model reads,

$\psi^{0}_{n}=\sqrt{\frac{2}{L}}\text{sin}(\frac{n\pi}{L})x\quad
n=1,2,3,\ldots,$ (13)

and the energy for the system is given by

$E^{0}_{n}=\frac{n^{2}\hbar^{2}\pi^{2}}{2mL^{2}}\quad n=1,2,3,\ldots.$ (14)

Here, we use the wave–function only corrected up to a first order for the
particle in a box and $x\in[0,1]$. In this experiment we not only learn the
potential but also learn the perturbed wave–function based only in the systems
initial conditions without perturbation. We use two neural networks, one to
learn the potential and the other to learn the perturbed wave–function. The
perturbed wave–function was learned in a supervised manner, whereas the
potential was learned in an unsupervised manner. If $W_{\theta}$ is the neural
network learning the perturbed wave–function $\psi_{\text{pert}}$, then our
auxiliary loss function becomes

$\mathcal{L}_{\theta}=\big{\lvert}\big{\rvert}W_{\theta}-\psi_{\text{pert}}\big{\lvert}\big{\rvert}_{2}^{2}+L_{\theta}$
(15)

where $L_{\theta}$ is the time–independent Schrödinger loss in the main text
used to learn the potential and is calculated using the perturbed
wave–function. Fig 8 shows our results on this system. It seems that energy is
not conserved for this system but that is merely due to our perturbation
approximations.

Figure 8: From left to right: Wave–function of the particle in a box,
Potential (y-axis) vs Position (x-axis) of the particle, Conservation of an
”energy” like object.

###  A.3 2D Harmonic Oscillator

Our work scales easily and quickly to $2$-dimensions as well. The
wave–function here is a product of two Hermite polynomials defined above. We
choose $x,y\in[0,1]$ as input to our model and constrained our output to
$[0,1]$. Thus our loss function is exactly as $1$D Harmonic Oscillator. Fig 9
shows our results for this system and the middle figure shows our learned
energy is a good approximation to the total energy (z scale chosen from
$[4.99,5.01]$).

Figure 9: Potential and Energies of the Harmonic Oscillator

##  Appendix B Motivation behind the time–independent Schrödinger Loss

We present a brief explanation for our time–independent Schrödinger Loss
function. The Hamilton $\hat{H}$ is the sum of kinetic energy $\hat{T}$ and
potential energy $\hat{V}$. The kinetic energy is given by the Laplacian
operator.

$\hat{H}=-\frac{\hbar^{2}}{2m}\nabla_{x}^{2}+\hat{V}(x).$ (16)

For the time independent case, the Schrödinger’s equation boils down to

$\hat{H}\psi=E\psi$ (17)

where E is the energy of the system. For simplicity, let $\hbar=m=1$. Using
equation 16, we can write

$(-\frac{1}{2}\nabla_{x}^{2}+\hat{V}(x))\psi=E\psi$ (18)

Dividing the above equation by $\psi$ throughout we get

$\frac{-\frac{1}{2}\nabla_{x}^{2}\psi}{\psi}+\hat{V}(x)=E$ (19)

Since the energy is constant, the derivative with respect to $x$ on the left
hand side of equation 19 is $0$ and this is our time–independent Schrödinger
loss function.

##  Appendix C Wigner Functions

An alternative formulation of quantum dynamics may be given by the Wigner
function (Curtright, Fairlie, and Zachos 1998; Chen, Xiong, and Shao 2019).
The Wigner function, $W(x,p,t)$, is a phase space distribution function which
behaves similarly to $|\psi\left(x\right)|^{2}$ and momentum
$|\psi\left(p\right)|^{2}$ distribution functions (Case 2008a). Unlike
wave–functions, Wigner functions are real valued and bounded. However,
contrary to probability distributions, $W(x,p,t)$ can take negative values.
Thus, the Wigner distribution is termed as a quasi–probability distribution
and so in a sense loses some of it’s classical appeal. Using the Schrödinger’s
equation (equation 1) and the Taylor expansion, the time evolution of the
Wigner function is given by an infinite order partial differential equation
called Wigner–Moyal equation (Case 2008a).

$\begin{split}\frac{\partial W(x,p,t)}{\partial
t}&={}-\frac{p}{m}\frac{\partial W(x,p,t)}{\partial x}\\\
{}&+\sum_{s=0}^{\infty}(-h^{2})^{s}\frac{1}{(2s+1!)}(\frac{1}{2})^{2s}\frac{\partial^{2s+1}U(x)}{\partial
x^{2s+1}}\frac{\partial^{2s+1}W(x,p,t)}{\partial p^{2s+1}}\end{split}$ (20)

###  C.1 Learning Potentials from Wigner Functions

In the case of the Wigner function, our Neural Network was trained by
implementing a truncated Wigner–Moyal loss,

$L_{Wigner}(\theta)=\bigg{\lvert}\bigg{\rvert}\frac{\partial
W(x,p,t)}{\partial t}+\frac{p}{m}\frac{\partial W(x,p,t)}{\partial
x}-\sum_{s=0}^{k}(-h^{2})^{s}\frac{1}{(2s+1!)}(1/2)^{2s}\frac{\partial^{2s+1}U_{\theta}(x)}{\partial
x^{2s+1}}\frac{\partial^{2s+1}W(x,p,t)}{\partial
p^{2s+1}}\bigg{\lvert}\bigg{\rvert}_{2}^{2}$ (21)

where for all our experiments $k=0,1$. The case where $k=0$ is known as the
Liouville equation. However, we note that equation 21 determines $U_{\theta}$
up to a constant. Thus, an initial condition depending on each individual
system was added.

###  C.2 Experiments with the Wigner functions

#### Harmonic Oscillator :

The Wigner function for the harmonic oscillator has the following form (Case
2008a):

$\displaystyle W\left(x,p,t\right)$
$\displaystyle=e^{-(x^{2}+p^{2})}\left(x^{2}+p^{2}+\sqrt{2}x\cos{t}-\sqrt{2}p\sin{t}\right)$

Since $\frac{\partial^{n}U}{\partial x^{n}}=0,\ \forall\ n\geq 3$, the
Moyal–Wigner equation in this case degenerates to the classical Liouville
equation. Let $x,p,t\in[0,1]$ and $x$ is the input to the model. The initial
condition is $U_{\theta}(0)=0$ and our loss function
$\mathcal{L}(\theta)=L_{\text{Wigner}}(\theta)+U_{\theta}(0)^{2}$ where
$L_{\text{Wigner}}$ is given by equation 21. Fig 10(left) shows the potential
learned by the model.

#### Pöschl-Teller potential :

The Wigner function in this case (Chen, Xiong, and Shao 2018) is given by:

$W_{2,1,0}(x,k,t):=\frac{3}{8}\int^{\infty}_{-\infty}\text{sech}^{2}(x+\frac{y}{2})\text{sech}^{2}(x-\frac{y}{2})\times\bigg{[}2\text{sinh}(x+\frac{y}{2})\text{sinh}(x-\frac{y}{2})+\sqrt{2}\text{sinh}(x-\frac{y}{2})e^{\frac{i3t}{2}}+\sqrt{2}\text{sinh}(x+\frac{y}{2})e^{\frac{-i3t}{2}}+1\bigg{]}e^{-iky}dy$
(22)

The Wigner function is a real–valued bounded function. Thus by breaking the
integral in equation 22 into real and complex parts, we only focus on the real
part. Using Euler’s formula, we get the following:

$g_{2,1,0}(x,k,t)=\frac{3}{8}\int_{-\infty}^{\infty}\text{sech}^{2}(x+\frac{y}{2})\text{sech}^{2}(x-\frac{y}{2})\times\bigg{[}2\text{sinh}(x+\frac{y}{2})\text{sinh}(x-\frac{y}{2})\text{cos}(-ky)+\sqrt{2}\text{sinh}(x-\frac{y}{2})\text{cos}(\frac{3t}{2}-ky)+\sqrt{2}\text{sinh}(x+\frac{y}{2})\text{cos}(-\frac{3t}{2}-ky)+\text{cos}(-ky)\bigg{]}$
(23)

Note that the integral in equation 23 is invariant under the change of
variable $y\rightarrow-y$. This implies in order to calculate
$g_{2,1,0}(x,k,t)$, we only have to integrate from $0$ to $\infty$ and
multiply that integral by 2. Our final simplification comes from studying the
decay properties of the Wigner functions. Using
$\text{sech}(x)=\frac{2}{e^{x}+e^{-x}}$ and
$\text{sinh}(x)=\frac{e^{x}-e^{-x}}{2}$, we found that the integrand in
equation 23 behaves like $O(e^{-y})$ (resp. $O(e^{y})$) as
$y\rightarrow\infty$ (resp. $y\rightarrow-\infty$). We picked a threshold of
$10^{-9}$ to truncate the integral from positive real axis to a bounded
interval which gives the following form :

$f_{2,1,0}(x,k,t)=\frac{3}{4}\int_{0}^{20}\text{sech}^{2}(x+\frac{y}{2})\text{sech}^{2}(x-\frac{y}{2})\times\bigg{[}2\text{sinh}(x+\frac{y}{2})\text{sinh}(x-\frac{y}{2})\text{cos}(-ky)+\sqrt{2}\text{sinh}(x-\frac{y}{2})\text{cos}(\frac{3t}{2}-k)+\sqrt{2}\text{sinh}(x+\frac{y}{2})\text{cos}(-\frac{3t}{2}-ky)+\text{cos}(-ky)\bigg{]}$
(24)

The Wigner method to study the time–frequency properties of dynamical systems
involves taking the partial derivatives with respect to time of the Wigner
function. These derivatives on the Wigner function yield what is known as the
Wigner–Moyal equation. The physical interpretations, numerical difficulties
and approximations of the Wigner–Moyal equation have been widely discussed in
the literature, thus for information about the mathematical challenges
associated with the Wigner–Moyal equation, we recommend readers to consult
these references (Chen, Xiong, and Shao 2018; Case 2008b; Galleani and Cohen
2002; Heller 1976; Curtright, Fairlie, and Zachos 1998; Klimov, Sainz, and
Romero 2020; Athanassoulis 2008; Gomes and Silva 2008).
The potential is $U(x)=-3\text{sech}(x)^{2}$ which is infinitely
differentiable. We choose $x\in[0,1]$ as input to our model. In this
experiment we will attempt to approximate the infinite order PDE (equation 20)
by equation 21. One then cannot assume that any non–steady state solution
predicted by the truncated Wigner function is immediately valid, as it can be
shown that higher order quantum corrections are responsible for quantum
mechanical phase space behavior. The $0$th order truncation matches the
potential in a small neighborhood of $0$. Figure 10 summarizes some of our
findings, for more details on approximation procedures and their challenges
for the infinite order Wigner–Moyal PDE see supplementary material.

Figure 10: Left to right: Position (x-axis) vs Potential (y-axis) of Harmonic
oscillator using the Wigner function; Our approximation of the Wigner function
of Pöschl–Teller potential; Various approximations of the potential using the
Wigner–Moyal equation. Table 3: A quantitative analysis of our model

System | RMSE between True and Learned Potentials | RMSE between True and Learned Energies
---|---|---
Harmonic Oscillator | $$1.1\text{\times}{10}^{-1}$\pm$5.0\text{\times}{10}^{-2}$$ | $$1.0\text{\times}{10}^{-1}$\pm$2.0\text{\times}{10}^{-2}$$
Pöchl–Teller potential | $$1.0\text{\times}{10}^{-4}$\pm$6.0\text{\times}{10}^{-5}$$ | $$8.0\text{\times}{10}^{-4}$\pm$6.0\text{\times}{10}^{-5}$$
Radial Hydrogen atom | $$3.0\text{\times}{10}^{-4}$\pm$8.0\text{\times}{10}^{-5}$$ | $$3.0\text{\times}{10}^{-4}$\pm$7.0\text{\times}{10}^{-5}$$
2D Harmonic Oscillator | $$3.0\text{\times}{10}^{-3}$\pm$9.0\text{\times}{10}^{-4}$$ | $$4.0\text{\times}{10}^{-3}$\pm$8.0\text{\times}{10}^{-4}$$
Particle in a Box | $$4.3\text{\times}{10}^{-1}$\pm$6.0\text{\times}{10}^{-2}$$ | $$5.5\text{\times}{10}^{-1}$\pm$8.0\text{\times}{10}^{-1}$$
Soliton | $$2.9\text{\times}{10}^{-1}$\pm$4.0\text{\times}{10}^{-2}$$ | -
Harmonic Oscillator from Wigner | $$4.0\text{\times}{10}^{-3}$\pm$8.0\text{\times}{10}^{-5}$$ | -

##  Appendix D Some training details and hyperparameters

Our Neural Network is a 4-layer feedforward network with a residual connection
between the second and the third layers. The activation and the scaling in the
final layers varied from experiment to experiment. Our main motivation for
scaling and using different activation is to show that an appropriate
architecture can perfectly learn the correct potential without an initial
condition. All the models are trained for $1000$ epochs. Table 4 shows the
activation, scaling and the size of the training data for each of the studied
systems. All the training data was randomly sampled from the appropriate
domains and trained in a minibatch fashion with batch size $32$.

Table 4: Some training details and model hyperparameters

System | Final Layer Activation | Final Layer Scaling | Size of training data
---|---|---|---
Harmonic Oscillator | Sigmoid | $12.5$ | $2500$
Pöchl–Teller potential | None | None | $2500$
Radial Hydrogen atom | None | None | $2500$
2D Harmonic Oscillator | Sigmoid | None | $5000$
Potential for Particle in a Box | Sigmoid | 10 | $4000$
Perturbation for Particle in a Box | None | None | $4000$
Soliton | None | None | $3000$
Harmonic Oscillator from Wigner | None | Sigmoid | $5000$
Pöschl–Teller from Wigner | None | None | $2000$

##  Appendix E Additional Figures

Some additional figures of the wave–functions and Wigner functions used in our
experiments.

Figure 11: Left: $2p$ Radial wave–function for the Hydrogen atom. Right: $2p$
Radial probability distribution for the Hydrogen atom  Figure 12:
Wave–function (Left) and Probability distribution (Right) for the first bound
state of the Pöschl–Teller potential.

Figure 13: Left: Wigner quasi–probability distribution for the ground state of
the Harmonic Oscillator, Right: Our approximation of the Wigner function for
the Pöschl–Teller

Figure 14: Left: 2D Harmonic Oscillator Wave–function, Right: 2D Harmonic
Oscillator Wave–function


Department of Physics, The Ohio State University, Columbus, USA; GSI Helmholtz
Centre for Heavy Ion Research, Darmstadt, Germany;

# Proton-proton and proton-cluster femtoscopy at the HADES experiment

M. Stefaniak for HADES Collaborationins:xins:y ins:xins:xins:yins:y

###### Abstract

This work explores femtoscopic correlations in proton-proton and proton-
cluster systems at the HADES experiment, GSI. Through high-precision
correlation functions, it reveals the impact of strong interactions between
protons and light nuclei (deuteron, triton, Helium-3). For proton-proton
interactions, measurements in two centralities and five $k_{T}$ bins were
analyzed. Extracted radii were compared with different Equations of State in
UrQMD, providing valuable insights into collision dynamics. This research
enhances our understanding of particle interactions and contributes to
refining theoretical models in nuclear and particle physics.

##  1 Introduction

The exploration of the Quantum Chromodynamics phase diagram stands as a
central and ambitious goal for numerous nuclear facilities. A recent surge in
interest has been directed towards a region characterized by finite net-baryon
densities and relatively low temperatures. This focus is driven by the belief
that the study of heavy-ion collisions within this specific domain can offer
valuable insights, not only for nuclear physics but also as a reference for
astrophysical investigations into Neutron Star Mergers (NSM). Such studies are
pivotal for unraveling the intricacies of the Equation of State (EoS)
governing nuclear matter.

One of the key challenges in this pursuit lies in comprehending the Strong
Interactions (SI) between hadrons and the excited or unbound states of nuclear
matter. These interactions remain inadequately described, posing a significant
gap in our understanding of the EoS. To address this, the HADES experiment
emerges as a crucial player, employing a specialized tool known as Femtoscopic
Correlations. This approach allows for a detailed examination of the complex
interaction of particles through precise correlation functions.

##  2 HADES experiment

The High Acceptance DiElectron Spectrometer (HADES) stands as a fixed-target
experiment situated at the SIS-18 accelerator within the GSI Helmholtz Centre
for Heavy Ion Research in Germany. This state-of-the-art experiment boasts an
almost complete azimuthal acceptance and covers polar angles ranging from
$18^{0}$ to $85^{0}$. This broad acceptance enables HADES to comprehensively
capture particle emissions from heavy-ion collisions, offering a unique
perspective on the intricate dynamics within these interactions.

The created systems in these collisions can be characterized by temperatures
($T$ = 60-80 MeV) and net baryon densities ($\rho$ $<$ 2-3 $\rho_{0}$),
closely resembling the thermal and density parameters observed in NSM
scenarios [1]. This selection of collision parameters enhances the
applicability of our findings to astrophysical studies, contributing valuable
insights to the determination of the EoS of nuclear matter.

##  3 Femtoscopic correlations

Femtoscopic correlations (FC) serve as a powerful tool for probing the
geometric and dynamic properties of the particle-emitting source [2]. By
delving into these correlations, we gain valuable insights into potential
excited states of particles, extracting crucial information about the SI
between the involved particles. The theoretical framework for femtoscopic
correlations is often encapsulated by the Koonin-Pratt formula [2]:

$C(k^{\star})=\int S(r^{\star})|\Psi(k^{\star},r^{\star})|^{2}d^{3}r^{\star}$
(1)

Here, $S(r^{\star})$ represents the source function, describing the
distribution of the relative positions ($r^{\star}$) of the two particles,
while $\Psi(k^{\star},r^{\star})$ is the two-particle wave function. Derived
from the Schrödinger equation for a given potential characterizing particle
interactions, the wave function is a key element encoding the effect of two-
particle interactions on femtoscopic correlations. Finally, $k^{\star}$
denotes the relative momentum of the particles in the pair rest frame.

The femtoscopic correlation function’s shape is influenced by the convolution
of three primary effects:

\- Coulomb interaction: This factor depends on the electric charge of the
particles and contributes to shaping the correlation function: attractive
unlike-sign particles, repulsive same-sign particles.

\- Quantum statistics: Governed by the Pauli exclusion law, quantum statistics
introduce a repulsion for fermions, influencing negatively overall form of the
correlation function.

\- Strong Interactions: The femtoscopic correlation function is profoundly
influenced by the potentials between correlated particles, encompassing both
attractive and repulsive forces. The presence of resonances further
contributes to the varied forms of the correlation function, adding complexity
to the analysis.

Experimentally, the femtoscopic correlation function is based on the relative
momenta of particle pairs and can be expressed as
$C(k^{\star})=A(k^{\star})/B(k^{\star})$. Here, $A(k^{\star})$ and
$B(k^{\star})$ represent the distributions of $k^{\star}$ for all possible
pairs within one event and for mixed events, respectively. The latter serves
as a representation of uncorrelated particles, and the division helps suppress
non-femtoscopic correlations arising from factors such as flow, non-uniform
detector efficiency, and other extraneous influences [2].

##  4 Results

In this section, we present the results derived from an analysis of data
obtained by the HADES experiment for Ag+Ag collisions at a center-of-mass
energy of $\sqrt{s_{NN}}$ = 2.55 GeV.

To address potential experimental effects, we implemented the following
measures. The challenges posed by two-track tracking detector effects,
particularly merging, were effectively mitigated by the application of a
$\Delta\theta>0.04$ and $\Delta\phi>0.25$ cut.

Flow effects, a common phenomenon in heavy-ion collisions, were systematically
suppressed by adopting mixing events with similar event planes. This is
achieved by requiring the respective event plane angles to the in the same
intervals of 30 degree width.

Additionally,an energy loss correction was applied to the reconstucted tracks
to account for variations in particle energies due to interactions with the
medium.

###  4.1 Proton-light nuclei correlation functions

We conducted a series of measurements involving the correlation functions of
proton-light nuclei systems. The sequence began with proton-proton (p-p)
correlations, followed by the addition of a neutron to form proton-deuteron
(p-d), and subsequently, the inclusion of an additional neutron to create
proton-triton (p-t), respectively of a proton to form proton-helium-3 (p-3He)
systems (see Figure 1). These measurements were carried out within the
$10\%-30\%$ centrality range in the transverse momentum ($k_{T}$) bin of
350-500 MeV/c and are presented in Fig. 2.

Figure 1: Sequence of analyzed proton-light nuclei systems. Figure 2: Measured
correlation functions of p-p (blue), p-d (red), p-t(green), p-3He (yellow) in
Ag+Ag collisions at $\sqrt{s_{NN}}$ = 2.55 GeV.

The proton-proton correlation, represented in blue, exhibits a positive
correlation attributed to SI. Additionally, it features a negative correlation
induced by the repulsive Coulomb interaction and Quantum Statistics, aligning
with theoretical expectations.

In the case of proton-deuteron correlations, depicted in red, the introduction
of a neutron to the correlated system results in repulsive interactions
between p-d pairs. Our correlation measurements are in good agreement with
previous FOPI results but also provide superior precision when compared to
other published results [7, 8].

Moving to proton-triton correlation, depicted in green, the inclusion of an
additional neutron in the correlated system (forming p-t) yields a strong
positive correlation. This correlation displays a sharp peak, possibly linked
to light nuclei decays, with greater prominence compared to FOPI results.

Finally, for p-3He correlation, shown in yellow, the addition of a proton to
the correlated system (p-3He) reveals a strong positive correlation. As there
is no comparison to FOPI shown in Fig. 2 this statement remains unclear and I
would suggest to remove it here.

Figure 3: Proton-triton (blue) and proton-3He (red) correlation functions.
Arrows correspond to positions of the possible decays of nuclear states.

In Figure 3, correlation functions for p-t and p-3He pairs are compared.
Despite triton and 3He having similar masses and consisting of three nucleons,
their correlation functions exhibit significant differences. The observed
variations in correlation function values are attributed to the decays of
excited nuclear states. The specific decays are summarized as follows:
${}^{4}He^{\star}\rightarrow p+t$:

  * •

E = 20.21 MeV, $J_{\pi}$ = $0_{+}$, $\Gamma$ = 0.5 MeV, $\Gamma_{p}/\Gamma$ =
1, $k_{1}^{\star}$ = 20 MeV/c

  * •

E = 21.01 MeV, $J_{\pi}$ = $0_{+}$, $\Gamma$ = 0.84 MeV, $\Gamma_{p}/\Gamma$ =
0.76, $k_{2}^{\star}$ =53.3 MeV/c

  * •

E = 21.84 MeV, $J_{\pi}$ = 2, $\Gamma$ =2.01 MeV, $\Gamma_{p}/\Gamma$ = 0.63,
$k_{3}^{\star}$ = 56.6 MeV/c

${}^{4}Li\rightarrow p+^{3}He$:

  * •

$J_{\pi}$ = 2-, $\Gamma$ = 6.0 MeV,$\Gamma_{p}/\Gamma$ = 1, $k_{1}^{\star}$ 72
MeV/c

While the enhancements in the correlation functions are visibly apparent,
estimating the number of decaying light nuclei is non-trivial due to the
convolution of repulsive Coulomb interaction and attractive SI shaping the
correlation function’s profile.

###  4.2 Proton-proton correlation functions

We present experimental results obtained for two centralities ($20\%-30\%$ and
$0-5\%$) and across five $k_{T}$ bins, as illustrated in Figure 4. The high-
quality statistics in our dataset enables thorough investigations of the
dependence on $k_{T}$. As $k_{T}$ increases, the correlation strength
amplifies, implying a reduction of the size (R) of the particle-emitting
source. At lower $k_{T}$ values, the correlation signals are notably
influenced by the resolution of momentum measurements and a less effective
suppression of track merging, causing an increase of systematic uncertainties.

Figure 4: Proton-proton correlation functions at $20\%-30\%$ (left panel) and
$0-5\%$ (right panel) centralities.

####  4.2.1 Femtoscopy vs EoS

In these simulations, depicted in Figure 5, different EoS were implemented
using the hadronic transport code UrQMD [11, 12]. For further insights into
the specific characteristics of the employed EoS, readers are directed to the
details provided in [9, 10]. The EoS and simulations were prepared by Jan
Steinheimer [14].

Figure 5: Pressure dependence on net baryon density of five EoS implemented in
UrQMD. HRG correspond to cascade mode, Soft and Hard EoS are based on Skyrme
potential, while CMF is chiral magnetic field, and CMF-PT includes the phase
transition [14].

The primary aim of these simulation studies is not to definitively establish
the ultimate EoS. Instead, the focus is on examining how different EoS relate
to femtoscopic measurements of p-p pairs. To achieve this, extreme scenarios
were considered; for instance, a soft EoS is characterized by a remarkably low
stiffness ($\kappa$ = 200 MeV), while an EoS incorporating a phase transition
(CMF-PT) is accentuated to a high degree.

It’s worth noting that while analogous studies have been conducted for pion-
pion pairs [10], this investigation extends the scope to the realm of p-p
femtoscopy, offering novel insights into the interplay between EoS and
particle correlations.

The theoretical correlation functions were derived using the CorAl software
[6], a development of Scott Pratt. This software incorporates the Reid
potential from the referenced work [13] in the calculation of the pair wave
function. Additionally, the source function is constructed by utilizing the
output from UrQMD simulations. This dual integration of theoretical tools
ensures a comprehensive and accurate representation of the correlation
functions, taking into account both potential interactions and simulated
source characteristics.

Figure 6: Particle emitting source dependence on $k_{T}$ for experimental (no
systematic uncertainties included) and simulated data. Left panel correspond
to midcentral collisions ($20\%-30\%$) and right one to central ($0-5\%$).

The source radius is extracted by fitting correlation functions using CorAl,
assuming a Gaussian source. In Figure 6, the comparison between HADES
experimental data (black points) and UrQMD simulations with different EoS for
two centralities reveals consistently larger source radii in the experimental
data. The discrepancies between experimental and simulated source radii can be
attributed to the exclusion of light nuclei production in the simulations,
challenges in accurately estimating proton ”feed-down” contributions (e.g.,
from 5Li), and the assumption of a Gaussian source in the fitting procedure.

In midcentral collisions, where differences between EoS are relatively small.
They start to become more pronounced as we transition to central collisions
with higher net baryon densities. Significantly, the comparison reveals
notable differences between the various EoS. The Skyrme EoS stands out as the
furthest from the experimental data.

##  5 Conclusions

In conclusion, our study on proton-light nuclei interactions has yielded
precise measurements, revealing the intricate effects of SI between protons
and light nuclei. The investigation further uncovered and the influence of
decays of excited states, specifically 4He$\star$ and 4Li. Our investigation
into proton-proton interactions revealed a robust and pronounced dependence of
correlation functions on $k_{T}$. This observation underscores the
significance of the $k_{T}$ parameter in shaping the dynamics of p-p
collisions, providing valuable insights into the complexities of these
interactions. The comparative analysis of different EoS within simulation
revealed shared features. Intriguingly, while there were nearly no discernible
differences between EoS when considering variations in centrality within the
same collision energy, the transition to higher baryon densities in more
central events accentuated more pronounced distinctions. These findings
underscore the sensitivity of correlation functions to changes in EoS,
providing valuable insights for refining EoS models in the context of high-
energy proton-proton collisions.

###### Acknowledgements.

I express sincere gratitude to Scott Pratt for invaluable support and to Jan
Steinheimer for providing simulated data with EoS along with engaging in
constructive discussions. This work was supported by Humboldt Foundation grant
for postdoctoral fellows and U.S. Department of Energy grant DE-SC0020651.

## References

  * [1] J. Adamczewski-Musch et al.: Nature Physics volume 15, pages 1040–1045 (2019)
  * [2] M. Lisa, S. Pratt, R. Soltz, U. Wiedemann, Ann.Rev. Nucl.Part.Sci. 55 (2005) 357-402
  * [3] V.G.J. Stoks et al., Phys. Rev. C 49, 2950 (1994)
  * [4] T.C. Black et al., Phys.Lett.B 471 (1999) 103-107
  * [5] T. V. Daniels et al., PRC 82, 034002 (2010)
  * [6] CorAl by Scott Pratt: https://github.com/scottedwardpratt/coral
  * [7] FOPI Collaboration: Eur. Phys. J. A 6, 185–195 (1999)
  * [8] Ke Mi: Acta Phys.Polon.Supp. 16 1, 1-A91 (2023)
  * [9] J. Steinheimer, at al: Eur.Phys.J.C 82 911 (2022)
  * [10] P. Li, at al: Sci. China-Phys. Mech. Astron. 66, 232011  (2023)
  * [11] S.A. Bass, at al: Prog. Part. Nucl. Phys. 41 225-370 (1998)
  * [12] M. Bleicher, at al: J. Phys. G: Nucl. Part. Phys. 25 1859-1896 (1999)
  * [13] V.G.J. Stoks et al., Phys. Rev. C 49, 2950 (1994)
  * [14] Private cominucation with Jan Steinheimer


Supplemental materials to

# Impacts of Point Defects on Shallow Doping in Cubic Boron Arsenide:
A First Principles Study

Shuxiang Zhou1, Zilong Hua1, Kaustubh K. Bawane1, Hao Zhou2, and Tianli Feng2
1Idaho National Laboratory, Idaho Falls, Idaho 83415, USA
2Department of Mechanical Engineering, University of Utah, Salt Lake City,
Utah 84112, USA

##  I  Formation energy of dopants with the presence of point defects

The formation energy $E_{f}$ of a defect (or a dopant) $X$ in pristine BAs is
calculated by:

$E_{f}(X)=E(X)-E_{0}-\sum_{i}{n_{i}^{X}(E_{i}+\mu_{i})}$ (S1)

where $E_{0}$ is the total energy of pristine system, $E(X)$ is the total
energy of the defected system, $n_{i}^{X}$ represents the number of atoms of
species $i$ changed according to the defect $X$ (positive if atoms are added,
negatives if atoms are removed), $E_{i}$ is the energy per atom in the
elemental phase, and $\mu_{i}$ is the chemical potential of species $i$.
Therefore, the formation energy of a coupled complex of dopant and defect
$D\mathchar 45\relax X$ is calculated by:

$E_{f}(D\mathchar 45\relax X)=E(D\mathchar 45\relax
X)-E_{0}-\sum_{i}{n_{i}^{D\mathchar 45\relax X}(E_{i}+\mu_{i})}=E(D\mathchar
45\relax
X)-E_{0}-\sum_{i}{n_{i}^{D}(E_{i}+\mu_{i})}-\sum_{i}{n_{i}^{X}(E_{i}+\mu_{i})}$
(S2)

As for the formation energy of a dopant $D$ in presence of a defect $X$, the
reference state is the defected BAs rather than pristine BAs:

$E_{f}^{X}(D)=E(D\mathchar 45\relax
X)-E(X)-\sum_{i}{n_{i}^{D}(E_{i}+\mu_{i})}$ (S3)

By combining Eq. S1-S3, one can obtain

$E_{f}^{X}(D)=E_{f}(D\mathchar 45\relax X)-E_{f}(X)$ (S4)

##  II  Formation energy calculation of dopants in BAs

The formation energy of charged p-type dopants (Be${}_{\text{B}}$,
C${}_{\text{As}}$, Si${}_{\text{As}}$, and Ge${}_{\text{As}}$) and n-type
dopants (Se${}_{\text{As}}$, C${}_{\text{B}}$, Si${}_{\text{B}}$, and
Ge${}_{\text{B}}$) as a function of the Fermi level $E_{F}$ was computed using
PBEsol (see Fig. S1). For a dopant $D$ with charge $q$, the formation energy
$E_{f}(D^{q})$ given by ):

$E_{f}(D^{q})=E(D^{q})-E_{0}-\sum_{i}{n_{i}^{D}(E_{i}+\mu_{i})}+q(E_{F}+E_{VBM})+E_{corr}(D^{q})$
(S5)

The formation energy of a charged dopant depends on the Fermi level $E_{F}$,
which is referenced to the energy of the valence band maximum (VBM),
$E_{VBM}$. Specifically, $E_{F}$ spans from zero to the band gap $E_{gap}$,
wherein $E_{F}=0$ signifies the Fermi level resting on the VBM, and
$E_{F}=E_{gap}$ indicates that the Fermi level aligns with the conduction band
minimum (CBM). In DFT calculations involving non-zero charges, it becomes
necessary to simulate a neutralizing background to circumvent infinite charge
issues under periodic boundary conditions. Thereby, a correction energy
$E_{corr}$ is introduced to account for the nonphysical electrostatic
interaction between periodic charged-defect images, which is computed using
SXDEFECTALIGN package [1]. The charges $q$ computed in the systems varied from
-2 to 2.

Comparing with the HSE results from Ref. [2, 3, 4], beside the underpredicted
band gap, the ionization energy is overpredicted in our PBEsol calculations.
However, the formation energy is comparable between our PBEsol results and
previous HSE results, where Be${}_{\text{B}}$ and Si${}_{\text{As}}$ are the
two p-type dopants with the lowest formation energy, while Se${}_{\text{As}}$
and Si${}_{\text{B}}$ are the two n-type dopants with the lowest formation
energy.

Figure S1:  Calculated formation energies of (a) p-type dopants and (b) n-type
dopants, as function of Fermi energy, which the minimum and maximum values are
corresponding to VBM and CBM, respectively.

##  III  Configurations of coupled complex of dopants and point defects

For each coupled complex of dopants and point defects, we compute multiple
configurations, where the dopants and point defects occupy different sites. We
first tabulate the the fractional coordinates of each site (see Table
LABEL:table:sm_sites), then provide the total energy of each configuration
(see Table S2-S7). For each coupled complex, the lowest total energy is
highlighted in bold text, to indicate the specific configuration used in the
DOS and formation energy calculation. The dopants include Be${}_{\text{B}}$,
Si${}_{\text{As}}$, Se${}_{\text{As}}$, and Si${}_{\text{B}}$, and the point
defects include C${}_{\text{As}}$, O${}_{\text{B}}$O${}_{\text{As}}$,
Si${}_{\text{As}}$, C${}_{\text{As}}$Si${}_{\text{B}}$,
O${}_{\text{B}}$Si${}_{\text{As}}$, As${}_{\text{B}}$B${}_{\text{As}}$ (as
discussed in the main text), As${}_{\text{B}}$, B${}_{\text{As}}$, and
N${}_{\text{As}}$ (as discussed in the following section). Therefore, the
coupled complexes can be categorized into six groups:

  * •

D${}_{\text{B}}$-X${}_{\text{As}}$, where D = Be or Si, X = C, Si, B, or N

  * •

D${}_{\text{As}}$-X${}_{\text{As}}$, where D = Si or Se, X = C, Si, B, or N,
excluding that both D and X = Si;

  * •

D${}_{\text{B}}$-X${}_{\text{B}}$Y${}_{\text{As}}$, where D = Be or Si, (X, Y)
= (O, O), (Si, C), (O, Si), or (As, B);

  * •

D${}_{\text{As}}$-X${}_{\text{B}}$Y${}_{\text{As}}$, where D = Si or Se, (X,
Y) = (O, O), (Si, C), (O, Si), or (As, B);

  * •

D${}_{\text{B}}$-X${}_{\text{B}}$, where D = Be or Si, X = As;

  * •

D${}_{\text{As}}$-X${}_{\text{B}}$, where D = Si or Se, X = As.

Table S1: The fractional coordinates of each B and As sites in the $3\times 3\times 3$ supercell of the conventional cubic cell of BAs.  Sites | Fractional coordinates | Sites | Fractional coordinates | Sites | Fractional coordinates | Sites | Fractional coordinates
---|---|---|---|---|---|---|---
As1 | (0, 0, 0) | As55 | (0.5, 0, 0.1667) | B1 | (0.0833, 0.0833, 0.0833) | B55 | (0.5833, 0.0833, 0.25)
As2 | (0, 0, 0.3333) | As56 | (0.5, 0, 0.5) | B2 | (0.0833, 0.0833, 0.4167) | B56 | (0.5833, 0.0833, 0.5833)
As3 | (0, 0, 0.6667) | As57 | (0.5, 0, 0.8333) | B3 | (0.0833, 0.0833, 0.75) | B57 | (0.5833, 0.0833, 0.9167)
As4 | (0, 0.1667, 0.1667) | As58 | (0.5, 0.1667, 0) | B4 | (0.0833, 0.25, 0.25) | B58 | (0.5833, 0.25, 0.0833)
As5 | (0, 0.1667, 0.5) | As59 | (0.5, 0.1667, 0.3333) | B5 | (0.0833, 0.25, 0.5833) | B59 | (0.5833, 0.25, 0.4167)
As6 | (0, 0.1667, 0.8333) | As60 | (0.5, 0.1667, 0.6667) | B6 | (0.0833, 0.25, 0.9167) | B60 | (0.5833, 0.25, 0.75)
As7 | (0, 0.3333, 0) | As61 | (0.5, 0.3333, 0.1667) | B7 | (0.0833, 0.4167, 0.0833) | B61 | (0.5833, 0.4167, 0.25)
As8 | (0, 0.3333, 0.3333) | As62 | (0.5, 0.3333, 0.5) | B8 | (0.0833, 0.4167, 0.4167) | B62 | (0.5833, 0.4167, 0.5833)
As9 | (0, 0.3333, 0.6667) | As63 | (0.5, 0.3333, 0.8333) | B9 | (0.0833, 0.4167, 0.75) | B63 | (0.5833, 0.4167, 0.9167)
As10 | (0, 0.5, 0.1667) | As64 | (0.5, 0.5, 0) | B10 | (0.0833, 0.5833, 0.25) | B64 | (0.5833, 0.5833, 0.0833)
As11 | (0, 0.5, 0.5) | As65 | (0.5, 0.5, 0.3333) | B11 | (0.0833, 0.5833, 0.5833) | B65 | (0.5833, 0.5833, 0.4167)
As12 | (0, 0.5, 0.8333) | As66 | (0.5, 0.5, 0.6667) | B12 | (0.0833, 0.5833, 0.9167) | B66 | (0.5833, 0.5833, 0.75)
As13 | (0, 0.6667, 0) | As67 | (0.5, 0.6667, 0.1667) | B13 | (0.0833, 0.75, 0.0833) | B67 | (0.5833, 0.75, 0.25)
As14 | (0, 0.6667, 0.3333) | As68 | (0.5, 0.6667, 0.5) | B14 | (0.0833, 0.75, 0.4167) | B68 | (0.5833, 0.75, 0.5833)
As15 | (0, 0.6667, 0.6667) | As69 | (0.5, 0.6667, 0.8333) | B15 | (0.0833, 0.75, 0.75) | B69 | (0.5833, 0.75, 0.9167)
As16 | (0, 0.8333, 0.1667) | As70 | (0.5, 0.8333, 0) | B16 | (0.0833, 0.9167, 0.25) | B70 | (0.5833, 0.9167, 0.0833)
As17 | (0, 0.8333, 0.5) | As71 | (0.5, 0.8333, 0.3333) | B17 | (0.0833, 0.9167, 0.5833) | B71 | (0.5833, 0.9167, 0.4167)
As18 | (0, 0.8333, 0.8333) | As72 | (0.5, 0.8333, 0.6667) | B18 | (0.0833, 0.9167, 0.9167) | B72 | (0.5833, 0.9167, 0.75)
As19 | (0.1667, 0, 0.1667) | As73 | (0.6667, 0, 0) | B19 | (0.25, 0.0833, 0.25) | B73 | (0.75, 0.0833, 0.0833)
As20 | (0.1667, 0, 0.5) | As74 | (0.6667, 0, 0.3333) | B20 | (0.25, 0.0833, 0.5833) | B74 | (0.75, 0.0833, 0.4167)
As21 | (0.1667, 0, 0.8333) | As75 | (0.6667, 0, 0.6667) | B21 | (0.25, 0.0833, 0.9167) | B75 | (0.75, 0.0833, 0.75)
As22 | (0.1667, 0.1667, 0) | As76 | (0.6667, 0.1667, 0.1667) | B22 | (0.25, 0.25, 0.0833) | B76 | (0.75, 0.25, 0.25)
As23 | (0.1667, 0.1667, 0.3333) | As77 | (0.6667, 0.1667, 0.5) | B23 | (0.25, 0.25, 0.4167) | B77 | (0.75, 0.25, 0.5833)
As24 | (0.1667, 0.1667, 0.6667) | As78 | (0.6667, 0.1667, 0.8333) | B24 | (0.25, 0.25, 0.75) | B78 | (0.75, 0.25, 0.9167)
As25 | (0.1667, 0.3333, 0.1667) | As79 | (0.6667, 0.3333, 0) | B25 | (0.25, 0.4167, 0.25) | B79 | (0.75, 0.4167, 0.0833)
As26 | (0.1667, 0.3333, 0.5) | As80 | (0.6667, 0.3333, 0.3333) | B26 | (0.25, 0.4167, 0.5833) | B80 | (0.75, 0.4167, 0.4167)
As27 | (0.1667, 0.3333, 0.8333) | As81 | (0.6667, 0.3333, 0.6667) | B27 | (0.25, 0.4167, 0.9167) | B81 | (0.75, 0.4167, 0.75)
As28 | (0.1667, 0.5, 0) | As82 | (0.6667, 0.5, 0.1667) | B28 | (0.25, 0.5833, 0.0833) | B82 | (0.75, 0.5833, 0.25)
As29 | (0.1667, 0.5, 0.3333) | As83 | (0.6667, 0.5, 0.5) | B29 | (0.25, 0.5833, 0.4167) | B83 | (0.75, 0.5833, 0.5833)
As30 | (0.1667, 0.5, 0.6667) | As84 | (0.6667, 0.5, 0.8333) | B30 | (0.25, 0.5833, 0.75) | B84 | (0.75, 0.5833, 0.9167)
As31 | (0.1667, 0.6667, 0.1667) | As85 | (0.6667, 0.6667, 0) | B31 | (0.25, 0.75, 0.25) | B85 | (0.75, 0.75, 0.0833)
As32 | (0.1667, 0.6667, 0.5) | As86 | (0.6667, 0.6667, 0.3333) | B32 | (0.25, 0.75, 0.5833) | B86 | (0.75, 0.75, 0.4167)
As33 | (0.1667, 0.6667, 0.8333) | As87 | (0.6667, 0.6667, 0.6667) | B33 | (0.25, 0.75, 0.9167) | B87 | (0.75, 0.75, 0.75)
As34 | (0.1667, 0.8333, 0) | As88 | (0.6667, 0.8333, 0.1667) | B34 | (0.25, 0.9167, 0.0833) | B88 | (0.75, 0.9167, 0.25)
As35 | (0.1667, 0.8333, 0.3333) | As89 | (0.6667, 0.8333, 0.5) | B35 | (0.25, 0.9167, 0.4167) | B89 | (0.75, 0.9167, 0.5833)
As36 | (0.1667, 0.8333, 0.6667) | As90 | (0.6667, 0.8333, 0.8333) | B36 | (0.25, 0.9167, 0.75) | B90 | (0.75, 0.9167, 0.9167)
As37 | (0.3333, 0, 0) | As91 | (0.8333, 0, 0.1667) | B37 | (0.4167, 0.0833, 0.0833) | B91 | (0.9167, 0.0833, 0.25)
As38 | (0.3333, 0, 0.3333) | As92 | (0.8333, 0, 0.5) | B38 | (0.4167, 0.0833, 0.4167) | B92 | (0.9167, 0.0833, 0.5833)
As39 | (0.3333, 0, 0.6667) | As93 | (0.8333, 0, 0.8333) | B39 | (0.4167, 0.0833, 0.75) | B93 | (0.9167, 0.0833, 0.9167)
As40 | (0.3333, 0.1667, 0.1667) | As94 | (0.8333, 0.1667, 0) | B40 | (0.4167, 0.25, 0.25) | B94 | (0.9167, 0.25, 0.0833)
As41 | (0.3333, 0.1667, 0.5) | As95 | (0.8333, 0.1667, 0.3333) | B41 | (0.4167, 0.25, 0.5833) | B95 | (0.9167, 0.25, 0.4167)
As42 | (0.3333, 0.1667, 0.8333) | As96 | (0.8333, 0.1667, 0.6667) | B42 | (0.4167, 0.25, 0.9167) | B96 | (0.9167, 0.25, 0.75)
As43 | (0.3333, 0.3333, 0) | As97 | (0.8333, 0.3333, 0.1667) | B43 | (0.4167, 0.4167, 0.0833) | B97 | (0.9167, 0.4167, 0.25)
As44 | (0.3333, 0.3333, 0.3333) | As98 | (0.8333, 0.3333, 0.5) | B44 | (0.4167, 0.4167, 0.4167) | B98 | (0.9167, 0.4167, 0.5833)
As45 | (0.3333, 0.3333, 0.6667) | As99 | (0.8333, 0.3333, 0.8333) | B45 | (0.4167, 0.4167, 0.75) | B99 | (0.9167, 0.4167, 0.9167)
As46 | (0.3333, 0.5, 0.1667) | As100 | (0.8333, 0.5, 0) | B46 | (0.4167, 0.5833, 0.25) | B100 | (0.9167, 0.5833, 0.0833)
As47 | (0.3333, 0.5, 0.5) | As101 | (0.8333, 0.5, 0.3333) | B47 | (0.4167, 0.5833, 0.5833) | B101 | (0.9167, 0.5833, 0.4167)
As48 | (0.3333, 0.5, 0.8333) | As102 | (0.8333, 0.5, 0.6667) | B48 | (0.4167, 0.5833, 0.9167) | B102 | (0.9167, 0.5833, 0.75)
As49 | (0.3333, 0.6667, 0) | As103 | (0.8333, 0.6667, 0.1667) | B49 | (0.4167, 0.75, 0.0833) | B103 | (0.9167, 0.75, 0.25)
As50 | (0.3333, 0.6667, 0.3333) | As104 | (0.8333, 0.6667, 0.5) | B50 | (0.4167, 0.75, 0.4167) | B104 | (0.9167, 0.75, 0.5833)
As51 | (0.3333, 0.6667, 0.6667) | As105 | (0.8333, 0.6667, 0.8333) | B51 | (0.4167, 0.75, 0.75) | B105 | (0.9167, 0.75, 0.9167)
As52 | (0.3333, 0.8333, 0.1667) | As106 | (0.8333, 0.8333, 0) | B52 | (0.4167, 0.9167, 0.25) | B106 | (0.9167, 0.9167, 0.0833)
As53 | (0.3333, 0.8333, 0.5) | As107 | (0.8333, 0.8333, 0.3333) | B53 | (0.4167, 0.9167, 0.5833) | B107 | (0.9167, 0.9167, 0.4167)
As54 | (0.3333, 0.8333, 0.8333) | As108 | (0.8333, 0.8333, 0.6667) | B54 | (0.4167, 0.9167, 0.9167) | B108 | (0.9167, 0.9167, 0.75)
Table S2:  For the coupled complex of dopant and defect D${}_{\text{B}}$-X${}_{\text{As}}$, where D = Be or Si, X = C, Si, B, or N, the total energy $E($D${}_{\text{B}}$-X${}_{\text{As}})$ and the distance $d$ between D${}_{\text{B}}$ and X${}_{\text{As}}$ as functions of the site of D${}_{\text{B}}$. The site of X${}_{\text{As}}$ is at As1. | Site of
---
D${}_{\text{B}}$
| $d$
---
(Å)
|  $E($Be${}_{\text{B}}$-C${}_{\text{As}})$
---
(eV)
|  $E($Be${}_{\text{B}}$-Si${}_{\text{As}})$
---
(eV)
|  $E($Be${}_{\text{B}}$-B${}_{\text{As}})$
---
(eV)
|  $E($Be${}_{\text{B}}$-N${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{B}}$-C${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{B}}$-Si${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{B}}$-B${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{B}}$-N${}_{\text{As}})$
---
(eV)
B1 | 2.07 | -1317.62 | -1314.11 | -1314.81 | -1318.44 | -1319.89 | -1316.28 | -1317.07 | -1318.34
B21 | 3.96 | -1317.33 | -1314.60 | -1313.99 | -1317.18 | -1318.77 | -1315.97 | -1315.70 | -1317.40
B22 | 5.21 | -1317.54 | -1314.62 | -1314.22 | -1317.39 | -1318.84 | -1315.92 | -1315.68 | -1317.57
B2 | 6.21 | -1317.47 | -1314.69 | -1314.14 | -1317.25 | -1318.70 | -1315.87 | -1315.53 | -1317.44
Table S3:  For the coupled complex of dopant and defect D${}_{\text{As}}$-X${}_{\text{As}}$, where D = Si or Se, X = C, Si, B, or N, the total energy $E($D${}_{\text{As}}$-X${}_{\text{As}})$ and the distance $d$ between D${}_{\text{As}}$ and X${}_{\text{As}}$ as functions of the site of D${}_{\text{As}}$. The site of X${}_{\text{As}}$ is at As1. | Site of
---
D${}_{\text{As}}$
| $d$
---
(Å)
|  $E($Si${}_{\text{As}}$-C${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{As}}$-B${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{As}}$-N${}_{\text{As}})$
---
(eV)
|  $E($Se${}_{\text{As}}$-C${}_{\text{As}})$
---
(eV)
|  $E($Se${}_{\text{As}}$-Si${}_{\text{As}})$
---
(eV)
|  $E($Se${}_{\text{As}}$-B${}_{\text{As}})$
---
(eV)
|  $E($Se${}_{\text{As}}$-N${}_{\text{As}})$
---
(eV)
As21 | 3.38 | -1320.74 | -1317.49 | -1320.72 | -1319.26 | -1316.38 | -1316.23 | -1317.82
As2 | 4.78 | -1320.90 | -1317.56 | -1320.72 | -1319.07 | -1316.24 | -1315.93 | -1317.78
As23 | 5.85 | -1320.93 | -1317.54 | -1320.72 | -1319.03 | -1316.21 | -1315.87 | -1317.77
As38 | 6.76 | -1320.96 | -1317.58 | -1320.74 | -1319.00 | -1316.20 | -1315.83 | -1317.75

Table S4:  For the coupled complex of dopant and defect
D${}_{\text{B}}$-X${}_{\text{B}}$Y${}_{\text{As}}$, where D = Be or Si, (X, Y)
= (O, O), (Si, C), (O, Si), or (As, B), the total energy
$E($D${}_{\text{B}}$-X${}_{\text{B}}$Y${}_{\text{As}})$, the distance $d_{1}$
between D${}_{\text{B}}$ and X${}_{\text{B}}$, and the distance $d_{2}$
between D${}_{\text{B}}$ and Y${}_{\text{As}}$, as functions of the site of
D${}_{\text{B}}$. The sites of X${}_{\text{B}}$ and Y${}_{\text{As}}$ are at
B1 and As1, respectively.

| Site of
---
D${}_{\text{B}}$
| $d_{1}$
---
(Å)
| $d_{2}$
---
(Å)
|  $E($Be${}_{\text{B}}$-O${}_{\text{B}}$O${}_{\text{As}})$
---
(eV)
|  $E($Be${}_{\text{B}}$-Si${}_{\text{B}}$C${}_{\text{As}})$
---
(eV)
|  $E($Be${}_{\text{B}}$-O${}_{\text{B}}$Si${}_{\text{As}})$
---
(eV)
|  $E($Be${}_{\text{B}}$-As${}_{\text{B}}$B${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{B}}$-O${}_{\text{B}}$O${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{B}}$-Si${}_{\text{B}}$C${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{B}}$-O${}_{\text{B}}$Si${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{B}}$-As${}_{\text{B}}$B${}_{\text{As}})$
---
(eV)
B106 | 3.38 | 2.07 | -1315.45 | -1316.95 | -1312.81 | -1312.29 | -1316.14 | -1316.68 | -1314.12 | -1313.84
B94 | 3.38 | 3.96 | -1314.11 | -1316.66 | -1313.34 | -1312.81 | -1314.40 | -1316.02 | -1314.63 | -1312.81
B3 | 4.78 | 3.96 | -1313.77 | -1316.52 | -1313.21 | -1312.57 | -1313.83 | -1316.12 | -1313.98 | -1312.99
B4 | 3.38 | 5.21 | -1313.87 | -1316.83 | -1313.58 | -1312.96 | -1314.37 | -1315.99 | -1314.19 | -1312.70
B108 | 5.85 | 3.96 | -1313.74 | -1316.50 | -1313.18 | -1312.54 | -1313.80 | -1316.13 | -1314.18 | -1313.02
B88 | 5.85 | 5.21 | -1314.05 | -1316.70 | -1313.24 | -1312.69 | -1314.19 | -1316.12 | -1314.23 | -1313.03
B2 | 4.78 | 6.21 | -1313.85 | -1316.67 | -1313.24 | -1312.78 | -1313.92 | -1316.04 | -1314.17 | -1312.85

Table S5:  For the coupled complex of dopant and defect
D${}_{\text{As}}$-X${}_{\text{B}}$Y${}_{\text{As}}$, where D = Si or Se, (X,
Y) = (O, O), (Si, C), (O, Si), or (As, B), the total energy
$E($D${}_{\text{As}}$-X${}_{\text{B}}$Y${}_{\text{As}})$, the distance $d_{1}$
between D${}_{\text{As}}$ and X${}_{\text{B}}$, and the distance $d_{2}$
between D${}_{\text{As}}$ and Y${}_{\text{As}}$, as functions of the site of
D${}_{\text{As}}$. The sites of X${}_{\text{B}}$ and Y${}_{\text{As}}$ are at
B1 and As1, respectively.

| Site of
---
D${}_{\text{As}}$
| $d_{1}$
---
(Å)
| $d_{2}$
---
(Å)
|  $E($Si${}_{\text{As}}$-O${}_{\text{B}}$O${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{As}}$-Si${}_{\text{B}}$C${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{As}}$-O${}_{\text{B}}$Si${}_{\text{As}})$
---
(eV)
|  $E($Si${}_{\text{As}}$-As${}_{\text{B}}$B${}_{\text{As}})$
---
(eV)
|  $E($Se${}_{\text{As}}$-O${}_{\text{B}}$O${}_{\text{As}})$
---
(eV)
|  $E($Se${}_{\text{As}}$-Si${}_{\text{B}}$C${}_{\text{As}})$
---
(eV)
|  $E($Se${}_{\text{As}}$-O${}_{\text{B}}$Si${}_{\text{As}})$
---
(eV)
|  $E($Se${}_{\text{As}}$-As${}_{\text{B}}$B${}_{\text{As}})$
---
(eV)
As19 | 2.07 | 3.38 | -1318.36 | -1320.33 | -1318.00 | -1316.81 | -1315.76 | -1317.08 | -1315.58 | -1313.28
As21 | 3.96 | 3.38 | -1317.40 | -1319.18 | -1316.66 | -1316.00 | -1314.40 | -1317.26 | -1314.67 | -1313.42
As2 | 3.96 | 4.78 | -1317.37 | -1319.35 | -1316.72 | -1316.27 | -1314.25 | -1317.09 | -1314.59 | -1313.16
As93 | 5.21 | 3.38 | -1317.28 | -1319.29 | -1316.65 | -1315.96 | -1314.28 | -1317.31 | -1314.34 | -1313.50
As23 | 3.96 | 5.85 | -1317.36 | -1319.39 | -1316.76 | -1316.35 | -1314.24 | -1317.05 | -1314.55 | -1313.10
As3 | 6.21 | 4.78 | -1317.29 | -1319.29 | -1316.70 | -1316.10 | -1314.22 | -1317.17 | -1314.21 | -1313.30
As42 | 5.21 | 5.85 | -1317.42 | -1319.47 | -1316.83 | -1316.25 | -1314.46 | -1317.13 | -1314.59 | -1313.21

Table S6:  For the coupled complex of dopant and defect D${}_{\text{B}}$-As${}_{\text{B}}$, where D = Be or Si, the total energy $E($D${}_{\text{B}}$-As${}_{\text{B}})$ and the distance $d$ between D${}_{\text{B}}$ and As${}_{\text{B}}$ as functions of the site of D${}_{\text{B}}$. The site of As${}_{\text{B}}$ is at B1. | Site of
---
D${}_{\text{B}}$
| $d$
---
(Å)
|  $E($Be${}_{\text{B}}$-As${}_{\text{B}})$
---
(eV)
|  $E($Si${}_{\text{B}}$-As${}_{\text{B}})$
---
(eV)
B22 | 3.38 | -1310.16 | -1309.24
B2 | 4.78 | -1309.97 | -1309.36
B23 | 5.85 | -1309.95 | -1309.37
B38 | 6.76 | -1309.88 | -1309.29
Table S7:  For the coupled complex of dopant and defect D${}_{\text{As}}$-As${}_{\text{B}}$, where D = Si or Se, the total energy $E($D${}_{\text{As}}$-As${}_{\text{B}})$ and the distance $d$ between D${}_{\text{As}}$ and As${}_{\text{B}}$ as functions of the site of D${}_{\text{As}}$. The site of As${}_{\text{B}}$ is at B1. | Site of
---
D${}_{\text{As}}$
| $d$
---
(Å)
|  $E($Si${}_{\text{As}}$-As${}_{\text{B}})$
---
(eV)
|  $E($Se${}_{\text{As}}$-As${}_{\text{B}})$
---
(eV)
As1 | 2.07 | -1314.18 | -1310.50
As2 | 3.96 | -1313.54 | -1309.64
As38 | 5.21 | -1313.50 | -1309.73
As44 | 6.21 | -1313.38 | -1309.69

##  IV  Influence of other defects on the dopants

In this section we present the results for defects As${}_{\text{B}}$,
B${}_{\text{As}}$, and N${}_{\text{As}}$. As${}_{\text{B}}$ and
B${}_{\text{As}}$ are common antisite defects, while N${}_{\text{As}}$ is an
charge-neutral impurity, as N is in the same group of As. In N single point
defects and the point defect pairs in the first nearest neighbors, the
formation energies per point defect are 2.56, 0.82, 0.90, 2.59, 1.00 eV for
N${}_{\text{B}}$, N${}_{\text{As}}$, N${}_{\text{B}}$N${}_{\text{As}}$,
N${}_{\text{B}}$N${}_{\text{B}}$, and N${}_{\text{As}}$N${}_{\text{As}}$,
respectively. Therefore, N${}_{\text{As}}$ is the most favorable among these
defects. The DOS and isosurface results are presented in Fig. S2-S5, and the
formation energy is in Fig. S6. As${}_{\text{B}}$ is harmful to p-type
dopants, and is beneficial to Se${}_{\text{As}}$; B${}_{\text{As}}$ is harmful
to n-type dopants, and is beneficial to Be${}_{\text{B}}$. While
N${}_{\text{As}}$ does not couple with As-site dopants Si${}_{\text{As}}$ and
Se${}_{\text{As}}$, it is beneficial to Be${}_{\text{B}}$, as a p-type B-site
dopant, and is harmful to Si${}_{\text{B}}$, as a n-type B-site dopant.

Figure S2:  Calculated electronic DOS of Be${}_{\text{B}}$-doped BAs in the
presence of (a) As${}_{\text{B}}$, (b) B${}_{\text{As}}$, and (c)
N${}_{\text{As}}$ point defect. The inset figures show the isosurfaces of
charge density near VBM ($<$0.1eV) or within the band gap.  Figure S3:
Calculated electronic DOS of Si${}_{\text{As}}$-doped BAs in the presence of
(a) As${}_{\text{B}}$, (b) B${}_{\text{As}}$, and (c) N${}_{\text{As}}$ point
defect. The inset figures show the isosurfaces of charge density near VBM
($<$0.1eV) or within the band gap.  Figure S4:  Calculated electronic DOS of
Si${}_{\text{B}}$-doped BAs in the presence of (a) As${}_{\text{B}}$, (b)
B${}_{\text{As}}$, and (c) N${}_{\text{As}}$ point defect. The inset figures
show the isosurfaces of charge density near CBM ($<$0.1eV) or within the band
gap.  Figure S5:  Calculated electronic DOS of Se${}_{\text{As}}$-doped BAs in
the presence of (a) As${}_{\text{B}}$, (b) B${}_{\text{As}}$, and (c)
N${}_{\text{As}}$ point defect. The inset figures show the isosurfaces of
charge density near CBM ($<$0.1eV) or within the band gap.  Figure S6:
Calculated formation energy of (a) Be${}_{\text{B}}$, (b) Si${}_{\text{As}}$,
(c) Si${}_{\text{B}}$, and (d) Se${}_{\text{As}}$ dopants in the presence of
selected point defects. The dashed lines represent formation energy of dopants
without point defects. The green color signifies that the shallow doping state
is preserved as the dopants in pristine BAs, whereas the red color indicates
otherwise.

##  V  Isosurfaces of charge density

In this section we provide enlarged figures of isosurfaces of charge density,
for all bands near CBM or VBM ($<$0.1 eV) and within the band gap. In each
case, from left to right, the energies of electronic bands are arranged from
small to large values, indicating that the band indexes also increase from
small to large values.

###  V.1 Dopants

Be${}_{\text{B}}$: A three-fold degenerated band at VBM

Si${}_{\text{As}}$: A three-fold degenerated band at VBM

Si${}_{\text{B}}$: Two bands near CBM

Se${}_{\text{As}}$: Two bands near CBM

###  V.2 Dopants with point defects

Be${}_{\text{B}}$-As${}_{\text{B}}$B${}_{\text{As}}$: Three bands near VBM

Be${}_{\text{B}}$-C${}_{\text{As}}$: One band near VBM:

Be${}_{\text{B}}$-O${}_{\text{B}}$O${}_{\text{As}}$: One band at VBM, and one
band at CBM

Be${}_{\text{B}}$-Si${}_{\text{As}}$: Three bands near VBM

Be${}_{\text{B}}$-C${}_{\text{As}}$Si${}_{\text{B}}$: Three bands near VBM

Be${}_{\text{B}}$-O${}_{\text{B}}$Si${}_{\text{As}}$: One band at VBM, and one
band in the band gap

Be${}_{\text{B}}$-B${}_{\text{As}}$: One two-fold degenerated band at VBM, and
one band in the band gap

Be${}_{\text{B}}$-N${}_{\text{As}}$: One band near VBM, and one two-fold
degenerated band at VBM

Si${}_{\text{As}}$-As${}_{\text{B}}$B${}_{\text{As}}$: Three bands near VBM

Si${}_{\text{As}}$-C${}_{\text{As}}$: Three bands near VBM

Si${}_{\text{As}}$-O${}_{\text{B}}$O${}_{\text{As}}$: One band in the band gap

Si${}_{\text{As}}$-C${}_{\text{As}}$Si${}_{\text{B}}$: Three bands near VBM

Si${}_{\text{As}}$-O${}_{\text{B}}$Si${}_{\text{As}}$: Two bands at VBM, and
one band in the band gap

Si${}_{\text{As}}$-B${}_{\text{As}}$: Three bands near VBM

Si${}_{\text{As}}$-N${}_{\text{As}}$: Three bands near VBM

Si${}_{\text{B}}$-As${}_{\text{B}}$B${}_{\text{As}}$: Two bands near CBM

Si${}_{\text{B}}$-O${}_{\text{B}}$O${}_{\text{As}}$: One band in the band gap

Si${}_{\text{B}}$-C${}_{\text{As}}$Si${}_{\text{B}}$: Two bands near CBM

Si${}_{\text{B}}$-O${}_{\text{B}}$Si${}_{\text{As}}$: One band in the band gap

Si${}_{\text{B}}$-As${}_{\text{B}}$: One band in the band gap, and two bands
near CBM

Si${}_{\text{B}}$-N${}_{\text{As}}$: One band at CBM

Se${}_{\text{As}}$-As${}_{\text{B}}$B${}_{\text{As}}$: Two bands near CBM

Se${}_{\text{As}}$-O${}_{\text{B}}$O${}_{\text{As}}$: One band in the band gap

Se${}_{\text{As}}$-C${}_{\text{As}}$Si${}_{\text{B}}$: Two bands near CBM

Se${}_{\text{As}}$-O${}_{\text{B}}$Si${}_{\text{As}}$: One band in the band
gap

Se${}_{\text{As}}$-As${}_{\text{B}}$: Two bands near CBM

Se${}_{\text{As}}$-N${}_{\text{As}}$: Two bands near CBM

## References

  * Freysoldt  _et al._ [2009] C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Physical Review Letters 102, 016402 (2009).
  * Chae  _et al._ [2018] S. Chae, K. Mengle, J. T. Heron, and E. Kioupakis, Applied Physics Letters 113, 212101 (2018).
  * Lyons  _et al._ [2018] J. L. Lyons, J. B. Varley, E. R. Glaser, J. A. Freitas, Jr., J. C. Culbertson, F. Tian, G. A. Gamage, H. Sun, H. Ziyaee, and Z. Ren, Applied Physics Letters 113, 251902 (2018).
  * Buckeridge and Scanlon [2019] J. Buckeridge and D. O. Scanlon, Physical Review Materials 3, 051601 (2019).