# The use of the Born-Oppenheimer factorization in the phase-space
representation of the time-independent Schrödinger equation for bilinearly
coupled harmonic oscillators.

Carlos A. Arango<EMAIL_ADDRESS>Department of Chemical Sciences,
Universidad Icesi, Cali, Colombia

###### Abstract

A system of two bilinearly coupled harmonic oscillators has been solved
analytically by using the Born-Oppenheimer (BO) product wavefunction ansatz
and the phase-space bound trajectory approach [J. S. Molano et al., Chem.
Phys. Lett. 76(12), 138171 (2021)]. The bilinearly coupled oscillator system
allows to obtain the analytical expression of the quantum system, facilitating
comparison with the results of using the BO ansatz product. The analytical and
BO wavefunctions are obtained as a product of parabolic cylinder functions.
The arguments of the parabolic cylinder functions of the exact and the BO
wavefunctions are related to the physical coordinates by linear
transformations, represented by matrices $\mathsf{G}$ and $\mathsf{\tilde{G}}$
respectively. A $QU$ decomposition of the matrix $\mathsf{G}$ outputs an upper
triangular matrix $\mathsf{U}$ that is closely related to the BO
$\mathsf{\tilde{G}}$ matrix. The effect of the BO non-adiabatic coupling is
analyzed on the stability of the phase-space equations of motion. The
eigenvalues and eigenfunctions obtained by the Born-Oppenheimer phase-space
trajectory approach show excellent agreement with the analytical solutions.

stability analysis, phase-space, Born-Oppenheimer, bilinearly coupled harmonic
oscillators, non-adiabatic coupling, $QR$ decomposition

##  I Introduction

Coupled oscillators models have been extensively used in chemical physics (1,
2). Intramolecular vibrational energy transfer (1, 2), vibrational
spectroscopy of polyatomic molecules (3, 4), and quantum entanglement (5, 6)
are among the successful applications in chemical physics. Bilinearly coupled
harmonic oscillators have been vastly studied and used to represent physical
systems due to the possibility of having analytical solutions of the
Schrödinger equation by using algebraic decoupling approaches and generalized
transformations among other mathematical techniques (7, 8, 9, 10, 11).
Approximate solution methods based on the self consistent field (SCF) theory
and the adiabatic approximation (AA) have been used to solve the Schrödinger
equation for coupled harmonic oscillators (12, 13, 14, 15, 16, 17, 18). These
method are based on product factorizations of the wavefunction. Self
consistent field is based on Hartree products, meanwhile the adiabatic
approximation method is based on Born-Oppenheimer products. The SCF is
corrected by allowing state interaction (19). The AA is improved by including
the non-adiabatic corrections by means of perturbation theory (20, 21, 18).

In this work a phase-space trajectory method, proposed previously (22), for
solving the time independent Schrödinger equation (TISE) has been extended to
systems of two degrees of freedom. This phase-space method is based on the use
of bound trajectories in phase-space as eigenstates of the TISE. The extension
of the phase-space method to systems with two degrees of freedom has been done
by using the Born-Oppenheimer product ansatz for the wave function. The
equations developed by the use of the phase-space trajectory approach are
applied to a system of two harmonic oscillators with bilinear coupling. The
resulting differential equations can be written in matrix form allowing to
study the dynamical structure and stability of the system (23). Exact and
Born-Oppenheimer wavefunctions are written in terms of parabolic cylinder
functions $D_{\nu}$ (24). The arguments of the $D_{\nu}$ functions for the
exact and BO solutions are given as linear combination of the oscillators
coordinates, with well defined transformation matrices. The analysis of the
elements of these transformation matrices allows to relate the error in the BO
approximation with the coupling between the oscillators.

##  II Theory and Methods

A system of two interacting particles, labeled 1 and 2, with coordinates $x$
and $y$, respectively, is described by the time independent Schrödinger
equation (TISE) $\hat{H}\psi(x,y)=\mathcal{E}\psi(x,y)$. The Hamiltonian for
this system can be written as

$\hat{H}=\hat{T}_{1}(\hat{p}_{x})+\hat{T}_{2}(\hat{p}_{y})+\hat{V}_{12},$ (1)

with kinetic energy operators
$\hat{T}_{1}=-\tfrac{1}{2m_{1}}\partial^{2}_{x}$, and
$\hat{T}_{2}=-\tfrac{1}{2m_{2}}\partial^{2}_{y}$, and potential energy
operator $\hat{V}_{12}=\hat{V}_{12}(x,y)$. This Hamiltonian can be separated
as $\hat{H}=\hat{T}_{2}+\hat{h}$, with $\hat{h}=\hat{h}(\hat{p}_{x},x,y)$
given by

$\hat{h}=\hat{T}_{1}+\hat{V}_{12}.$ (2)

The Born-Oppenheimer approximation resorts to a product wave function,
$\psi=\chi\varphi$, of the form

$\psi(x,y)=\chi(x,y)\varphi(y).$ (3)

The $\chi(x,y)$ functions are obtained as the bound solutions of the
$y$-parametrized TISE

$\hat{h}\chi_{n}(x,y)=\varepsilon_{n}(y)\chi_{n}(x,y).$ (4)

The eigenvalues $\varepsilon_{n}(y)$ are the potential energy curves of the
$y$ degree of freedom. The eigenfunctions $\chi_{n}(x,y)$ obey the
orthonormality condition,

$\braket{\chi_{l}}{\chi_{n}}=\int_{-\infty}^{\infty}{\chi_{l}^{*}(x,y)\chi_{n}(x,y)dx}=\delta_{ln}(y),$
(5)

with $\delta_{ln}(y)$ as the $y$ dependent Kronecker delta.

The action of the Hamiltonian operator on the wave function $\psi$,
$\hat{H}\psi=\mathcal{E}\psi$, produces

$-\tfrac{1}{2m_{1}}\varphi\partial^{2}_{x}\chi-\tfrac{1}{2m_{2}}\left(\varphi\partial^{2}_{y}\chi+\chi\varphi^{\prime\prime}+2\varphi^{\prime}\partial_{y}\chi\right)+\hat{V}\chi\varphi=\mathcal{E}\chi\varphi,$
(6)

with $\mathcal{E}$ as the total energy of the full system. This result can be
rearranged to,

$-\tfrac{1}{2m_{2}}\chi\varphi^{\prime\prime}=\left(\tfrac{1}{2m_{1}}\partial^{2}_{x}\chi-(\hat{V}-\mathcal{E})\chi\right)\varphi+\tfrac{1}{2m_{2}}\left(2\varphi^{\prime}\partial_{y}\chi+\varphi\partial^{2}_{y}\chi\right).$
(7)

The last term of this equation can be written compactly as

$\frac{1}{2m_{2}}\left(2\varphi^{\prime}\partial_{y}\chi+\varphi\partial^{2}_{y}\chi\right)=\frac{\partial_{y}(\varphi^{2}\partial_{y}\chi)}{2m_{2}\varphi},$
(8)

which can be used, together with the definition of $\hat{h}$, in equation (7)
to get

$-\tfrac{1}{2m_{2}}\chi\varphi^{\prime\prime}=(\mathcal{E}-\hat{h})\chi\varphi+\frac{\partial_{y}(\varphi^{2}\partial_{y}\chi)}{2m_{2}\varphi}.$
(9)

The use of the eigenfunctions $\chi_{n}$ in equation (9) gives

$-\tfrac{1}{2m_{2}}\chi_{n}\varphi_{n}^{\prime\prime}=\left(\mathcal{E}-\varepsilon_{n}(y)\right)\chi_{n}\varphi_{n}+\frac{\partial_{y}(\varphi_{n}^{2}\partial_{y}\chi_{n})}{2m_{2}\varphi_{n}}.$
(10)

A simpler differential equation for $\varphi_{n}(y)$ is obtained if the last
term on the right of (10) is zero, which is possible only if
$\varphi_{n}^{2}\partial_{y}\chi_{n}=f_{n}(x)$, with $f_{n}(x)$ as a function
of $x$ only. Integrating this last result produces

$\chi_{n}=g_{n}(x)+f_{n}(x)\int{\frac{dy}{\varphi_{n}^{2}}}.$ (11)

with $g(x)$ as an arbitrary function of $x$. The result (11) indicates that
the second term of equation (10) can be left out in the case of complete
separation of variables.

The use of the definition $\vartheta_{n}=\varphi^{\prime}_{n}$, allows to
write equation (10) as a system of equations, for given $x$ and $n$,

$\begin{split}\varphi_{n}^{\prime}&=\vartheta_{n},\\\
\chi_{n}\vartheta_{n}^{\prime}&=\left(-2m_{2}\left(\mathcal{E}-\varepsilon_{n}(y)\right)\chi_{n}-\partial^{2}_{y}\chi_{n}\right)\varphi_{n}-(2\partial_{y}\chi_{n})\vartheta_{n}.\end{split}$
(12)

These are the equations for the phase-space of $\varphi_{n}$ for a given value
of the coordinate $x$.

Multiplying the second equation (12) on the left by $\chi_{l}^{*}$ and
integrating over $x$ gives

$\delta_{ln}\vartheta_{n}^{\prime}=\left(-2m_{2}(\mathcal{E}-\varepsilon_{n}(y))\delta_{ln}-\beta_{nl}\right)\varphi_{n}-2\alpha_{nl}\vartheta_{n},$
(13)

with $\alpha_{nl}=\braket{\chi_{l}}{\partial_{y}\chi_{n}}$ and
$\beta_{nl}=\braket{\chi_{l}}{\partial_{y}^{2}\chi_{n}}$ as the first and
second order non-adiabatic couplings, respectively.

The diagonal part of equation (13) gives

$\vartheta_{n}^{\prime}=\left(-2m_{2}(\mathcal{E}-\varepsilon_{n}(y))-\beta_{nn}\right)\varphi_{n}-2\alpha_{nn}\vartheta_{n},$
(14)

with $\alpha_{nn}=\braket{\chi_{n}}{\partial_{y}\chi_{n}}$ as the Berry’s
connection. Since the eigenfunctions $\chi_{n}$ are real and normalized, we
have $\alpha_{nn}=0$, and

$\vartheta_{n}^{\prime}=\left(-2m_{2}(\mathcal{E}-\varepsilon_{n}(y))-\beta_{nn}\right)\varphi_{n}.$
(15)

The off-diagonal elements of equation (13) give the equations

$-\beta_{nl}\varphi_{n}=2\alpha_{nl}\vartheta_{n}.$ (16)

By defining the vector $\mathbf{\Phi}_{n}=\\{\varphi_{n},\vartheta_{n}\\}$,
the diagonal equation (14) can be written in matrix form
$\mathbf{\Phi}^{\prime}_{n}=\mathsf{A}_{n}\mathbf{\Phi}_{n}$ with

$\mathsf{A}_{n}=\begin{pmatrix}0&1\\\
-2m_{2}\left(\mathcal{E}-\varepsilon_{n}(y)\right)-\beta_{nn}&-2\alpha_{nn}\\\
\end{pmatrix}.$ (17)

The system of equations
$\mathbf{\Phi}^{\prime}_{n}=\mathsf{A}_{n}\mathbf{\Phi}_{n}$ displays a fixed
point at $\mathbf{\Phi}_{n}=0$, with stability given by the eigenvalues of
matrix $\mathsf{A}_{n}$,

$\lambda_{1,2}=-\alpha_{nn}\mp\sqrt{\alpha_{nn}^{2}-2m_{2}\left(\mathcal{E}-\varepsilon_{n}(y)\right)-\beta_{nn}}.$
(18)

Elliptic phase-space trajectories around the fixed point are given by values
of $y$ where the condition
$\alpha_{nn}^{2}-2m_{2}\left(\mathcal{E}-\varepsilon_{n}(y)\right)-\beta_{nn}<0$
is fulfilled. Unlike 1-dimensional systems, classical turning points are not
boundaries between elliptic and hyperbolic behaviour of the phase-space
trajectories (22). Instead, 2-dimensional systems display non-adiabatic
turning points given by the $y$ solutions of
$\alpha_{nn}^{2}-2m_{2}\left(\mathcal{E}-\varepsilon_{n}(y)\right)-\beta_{nn}=0$.

##  III Bilinearly coupled 1-dimensional harmonic oscillators

Consider two 1-dimensional harmonic oscillators of frequencies $\omega$ and
$\Omega$, described by coordinates $x$ and $y$, and with masses $m$ and $M$,
respectively. The oscillators are coupled by a bilinear term, $cxy$. The
Hamiltonian operator is given by

$\hat{H}=-\frac{1}{2m}\frac{\partial^{2}}{\partial
x^{2}}-\frac{1}{2M}\frac{\partial^{2}}{\partial
y^{2}}+\frac{m\omega^{2}}{2}x^{2}+\frac{M\Omega^{2}}{2}y^{2}+cxy,$ (19)

It is convenient to describe the coupling by
$\delta=c\left(\sqrt{mM}\omega\Omega\right)^{-1}$, restricted to
$\left|\delta\right|<1$. The uncoupled case, $\delta=0$, gives as solution

$\Psi(x,y)=\frac{\left(mM\omega\Omega\right)^{1/4}}{\sqrt{\pi
n_{1}!n_{2}!}}D_{n_{1}}(\sqrt{2m\omega}x)D_{n_{2}}(\sqrt{2M\Omega}y),$ (20)

where $D_{\nu}(\xi)$ are the parabolic cylinder (Weber) functions (24), and
$\nu=n_{1},n_{2}$, as non negative integers.

###  III.1 Analytical solution

It is convenient to use the matrix notation for the Hamiltonian (19),

$2\hat{H}=\mathbf{z}^{\mathsf{t}}\mathsf{H}\mathbf{z},$ (21)

with phase-space operator $\mathbf{z}$ as the column vector obtained by the
concatenation of $\mathbf{r}=\begin{pmatrix}x\\\ y\end{pmatrix}$ and
$\nabla_{\mathbf{r}}=\begin{pmatrix}\partial_{x}\\\
\partial_{y}\end{pmatrix}$. The matrix $\mathsf{H}$ is given by

$\mathsf{H}=\begin{pmatrix}m\omega^{2}&c&0&0\\\ c&M\Omega^{2}&0&0\\\
0&0&-\frac{1}{m}&0\\\ 0&0&0&-\frac{1}{M}\end{pmatrix}.$ (22)

The use of mass-weighted coordinates $\bm{\rho}=\sqrt{\mathsf{M}}\mathbf{r}$,
with $\mathsf{M}=\mathop{\mathrm{diag}}(m,M)$, and the corresponding gradient
operator $\nabla_{\bm{\rho}}=\frac{1}{\sqrt{\mathsf{M}}}\nabla_{\mathbf{r}}$,
allows to define the mass-weighted phase-space operator $\bm{\zeta}$ as the
concatenation of $\bm{\rho}$ and $\nabla_{\bm{\rho}}$. The use of $\bm{\zeta}$
simplifies equation (21) to

$2\hat{H}=\bm{\zeta}^{\mathsf{t}}\mathsf{K}\bm{\zeta},$ (23)

with $\mathsf{K}=\mathop{\mathrm{diag}}\left(\mathsf{B},-\mathsf{I}\right)$ as
a block diagonal matrix, with

$\mathsf{B}=\begin{pmatrix}\omega^{2}&\delta\omega\Omega\\\
\delta\omega\Omega&\Omega^{2}\end{pmatrix},$ (24)

and $\mathsf{I}$ as the $2\times 2$ identity matrix.

Matrix $\mathsf{B}$ deals only with the potential energy part of the
Hamiltonian (19), $2V=\bm{\rho}^{\mathsf{t}}\mathsf{B}\bm{\rho}$. The
eigenvectors of $\mathsf{B}$ are rows of an orthogonal matrix $\mathsf{R}$
that diagonalizes $\mathsf{B}$. The transformation of coordinates
$\mathbf{q}=\mathsf{R}\bm{\rho}$, with $\mathbf{q}=\begin{pmatrix}q_{1}\\\
q_{2}\end{pmatrix}$, produces a diagonal matrix $\mathsf{L}$ similar to
$\mathsf{B}$, $\mathsf{L}=\mathsf{RBR^{\mathsf{t}}}$, with
$\mathsf{L}=\mathop{\mathrm{diag}}{\left(\lambda_{1},\lambda_{2}\right)}$ as
the matrix of eigenvalues of $\mathsf{B}$. The potential energy after this
transformation is given by $2V=\mathbf{q^{t}}\mathsf{L}\mathbf{q}$.

The eigenvalues of $\mathsf{B}$ are explicitly given by

$\lambda_{1,2}=\tfrac{1}{2}\left(\omega^{2}+\Omega^{2}\pm\sqrt{\left(\omega^{2}-\Omega^{2}\right)^{2}+4\delta^{2}\omega^{2}\Omega^{2}}\right),$
(25)

which gives $0<\lambda_{2}<\lambda_{1}$ for $\Omega<\omega$. The normalized
eigenvectors of $\mathsf{B}$ are the rows of the orthogonal matrix

$\mathsf{R}=\begin{pmatrix}\cos{\theta}&\sin{\theta}\\\
-\sin{\theta}&\cos{\theta}\end{pmatrix},$ (26)

with $\theta=\arctan{\frac{\delta\omega\Omega}{\lambda_{1}-\Omega^{2}}}$,
restricted to
$\theta\in\left(-\nicefrac{{\pi}}{{2}},\nicefrac{{\pi}}{{2}}\right)$. The
kinetic energy part of the Hamiltonian (19) requires the relationship between
the gradient operators,
$\nabla_{\mathbf{q}}=\mathsf{R}^{\mathsf{t}}\nabla_{\bm{\rho}}$. Since the
transformation matrix $\mathsf{R}$ is orthogonal, the kinetic energy operator
is finally given by $2\hat{T}=-\nabla_{\mathbf{q}}\cdot\nabla_{\mathbf{q}}$,
and he Hamiltonian simplifies to

$\displaystyle\hat{H}$
$\displaystyle=-\tfrac{1}{2}\nabla_{\mathbf{q}}\cdot\nabla_{\mathbf{q}}+\tfrac{1}{2}\mathbf{q^{t}}\mathsf{L}\mathbf{q}$
(27)
$\displaystyle=-\tfrac{1}{2}\left(\tfrac{\partial^{2}}{\partial{q_{1}^{2}}}+\tfrac{\partial^{2}}{\partial{q_{2}^{2}}}\right)+\tfrac{1}{2}\left(\omega_{1}^{2}q_{1}^{2}+\omega_{2}^{2}q_{2}^{2}\right),$
(28)

with $\omega_{1,2}=\sqrt{\lambda_{1,2}}$. The transformation of coordinates
from the original $\mathbf{r}$ to the $\mathbf{q}$ is given by
$\mathbf{q}=\mathsf{R}\sqrt{\mathsf{M}}\mathbf{r}$ with volume element
transforming as $dxdy=\frac{dq_{1}dq_{2}}{\sqrt{mM}}$.

The analytical energy eigenvalues are:

$\mathcal{E}_{n_{1},n_{2}}=\omega_{1}\left(n_{1}+\tfrac{1}{2}\right)+\omega_{2}\left(n_{2}+\tfrac{1}{2}\right),$
(29)

with $n_{1,2}$ non negative integers. A simple inspection of equations (25)
and (29), for $\omega>\Omega$, shows that for the uncoupled case $\delta=0$,

$\mathcal{E}_{n_{1},n_{2}}=\omega\left(n_{1}+\tfrac{1}{2}\right)+\Omega\left(n_{2}+\tfrac{1}{2}\right).$
(30)

The normalized analytical wave function is given by
$\Psi_{n_{1},n_{2}}(q_{1},q_{2})=\psi_{n_{1}}(q_{1})\psi_{n_{2}}(q_{2})$, with

$\psi_{n_{i}}(q_{i})=\frac{1}{\sqrt{n_{i}!}}\left(\frac{mM\omega_{i}}{\pi}\right)^{1/4}D_{n_{i}}(\sqrt{2\omega_{i}}q_{i}),$
(31)

where $i=1,2$, and $D_{n_{i}}$ are the parabolic cylinder functions. The
normalized wave function is given by

$\begin{split}\Psi(q_{1},q_{2})&=\frac{1}{\sqrt{n_{1}!n_{2}!}}\left(\frac{\omega_{1}\omega_{2}}{\pi^{2}}\right)^{1/4}D_{n_{1}}(\sqrt{2\omega_{1}}q_{1})D_{n_{2}}(\sqrt{2\omega_{2}}q_{2})\\\
&=\frac{\left(\omega\Omega\sqrt{1-\delta^{2}}\right)^{1/4}}{\sqrt{\pi
n_{1}!n_{2}!}}D_{n_{1}}(\sqrt{2\omega_{1}}q_{1})D_{n_{2}}(\sqrt{2\omega_{2}}q_{2}).\end{split}$
(32)

The arguments of the $D_{n_{1}}$ and $D_{n_{2}}$ functions define the vector
$\bm{\xi}=\mathsf{G}\mathbf{r}$, with the matrix
$\mathsf{G}=\sqrt{2}\mathsf{F}\mathsf{R}\mathsf{M}^{\nicefrac{{1}}{{2}}}$, and
$\mathsf{F}=\mathop{\mathrm{diag}}{(\sqrt{\omega_{1}},\sqrt{\omega_{2}})}$.
Explicitly, matrix $\mathsf{G}$ is given by

$\mathsf{G}=\begin{pmatrix}\sqrt{2m\omega_{1}}\cos{\theta}&\sqrt{2M\omega_{1}}\sin{\theta}\\\
-\sqrt{2m\omega_{2}}\sin{\theta}&\sqrt{2M\omega_{2}}\cos{\theta}\end{pmatrix}.$
(33)

The volume element in terms of the components of $\bm{\xi}$ is given by

$d\xi_{1}d\xi_{2}=2(1-\delta^{2})^{\nicefrac{{1}}{{4}}}\sqrt{mM\omega\Omega}\,dxdy.$
(34)

###  III.2 Born-Oppenheimer solution

####  III.2.1 The $\chi_{n}(x,y)$ function

The system $x$ has Hamiltonian $\hat{h}$

$\hat{h}=-\frac{1}{2m}\frac{\partial^{2}}{\partial
x^{2}}+\frac{m\omega^{2}}{2}x^{2}+\frac{M\Omega^{2}}{2}y^{2}+cxy.$ (35)

The $y$-parametrized TISE for $x$, $\hat{h}\chi=\varepsilon_{n}(y)\chi$, has
analytical solution

$\chi(x,y)=C(y)D_{\nu}(z)$ (36)

$C(y)$ is a constant function of $y$. $D_{\nu}(z)$ is a parabolic cylinder
(Weber) function with

$z(x,y)=\left(\sqrt{2m\omega}x+\sqrt{\frac{2M\Omega^{2}}{\omega}}\delta
y\right),$ (37)

and

$\nu(y,\varepsilon)=\frac{\varepsilon}{\omega}-\frac{1}{2}+\frac{(\delta^{2}-1)M\Omega^{2}}{2\omega}y^{2}.$
(38)

The linear form of $z(x,y)$ allows to write

$z(x,y)=(\partial_{x}z)x+(\partial_{y}z)y,$ (39)

with

$\displaystyle\partial_{x}z$ $\displaystyle=\sqrt{2m\omega},$ (40)
$\displaystyle\partial_{y}z$
$\displaystyle=\sqrt{\frac{2M\Omega^{2}}{\omega}}\delta,$ (41)
$\displaystyle\partial_{y}x$
$\displaystyle=-\sqrt{\frac{m}{M}}\frac{\Omega\delta}{\omega}.$ (42)

Bound solutions of (35) are obtained by the condition $\nu(y,\varepsilon)=n$,
with $n$ a non negative integer. It is easily shown that

$\varepsilon_{n}(y)=\omega\left(n+\frac{1}{2}\right)+\frac{(1-\delta^{2})M\Omega^{2}}{2}y^{2}.$
(43)

These eigenvalues are the potential energy curves of the $y$ degree of
freedom. Since $|\delta|<1$, the effect of the coupling can be seen as a
reduction of $\Omega$ to an effective frequency $\Omega\sqrt{1-\delta^{2}}$.

The inner product on $x$ of the $D_{n}(z)$ functions is obtained by using the
chain rule,

$\begin{split}\braket{D_{l}(z)}{D_{n}(z)}&=\int_{-\infty}^{\infty}D_{l}(z)D_{n}(z)dx\\\
&=\int_{-\infty}^{\infty}D_{l}D_{n}(\partial_{x}z)^{-1}dz\\\
&=\frac{\sqrt{2\pi}}{\partial_{x}z}n!\delta_{ln}.\end{split}$ (44)

The normalized function $\chi_{n}$ is

$\begin{split}\chi_{n}(z)&=\frac{1}{\sqrt{n!}}\frac{\sqrt{\partial_{x}z}}{(2\pi)^{1/4}}D_{n}(z)\\\
&=\frac{1}{\sqrt{n!}}\left(\frac{m\omega}{\pi}\right)^{1/4}D_{n}(z).\end{split}$
(45)

####  III.2.2 Phase-space equations of motion and stability analysis

The equations of motion of $\left\\{\varphi_{n},\vartheta_{n}(y)\right\\}$ are
given by

$\begin{split}\varphi_{n}^{\prime}&=\vartheta_{n},\\\
\vartheta_{n}^{\prime}&=\left(-2M(\mathcal{E}-\varepsilon_{n}(y))-\beta_{nn}\right)\varphi_{n},\end{split}$
(46)

with

$\displaystyle\varepsilon_{n}(y)$
$\displaystyle=\omega\left(n+\tfrac{1}{2}+\tfrac{M\Omega^{2}}{2\omega}\left(1-\delta^{2}\right)y^{2}\right),$
$\displaystyle\beta_{nn}$
$\displaystyle=-\frac{M\Omega^{2}\delta^{2}}{\omega}\left(n+\tfrac{1}{2}\right).$

A detailed calculation of the non-adiabatic coupling
$\braket{\chi_{n}}{\partial^{2}_{y}\chi_{n}}$ can be seen in the Appendix,
section VII. This system of two bilinearly coupled harmonic oscillators has a
constant non-adiabatic coupling that only depends on the quantum number $n$.

Equations (46) are the equations for the tangent field of
$\\{\varphi_{n},\vartheta_{n}\\}$. It is straightforward to verify that the
fixed point of this field is $\\{\varphi_{n},\vartheta_{n}\\}^{*}=\\{0,0\\}$.
The eigenvalues of the equations for the tangent field are given by equation
(18),

$\lambda_{1,2}=\mp\sqrt{-2M(\mathcal{E}-\varepsilon_{n}(y))-\beta_{nn}}.$ (47)

The condition for elliptic behaviour of the phase-space trajectory is given
for values of $y$ such that
$\mathcal{E}>\varepsilon_{n}(y)-\frac{1}{2M}\beta_{nn}$, explicitly

$\mathcal{E}>\omega\left(n+\tfrac{1}{2}\right)\left(1+\frac{\Omega^{2}\delta^{2}}{2\omega^{2}}\right)+\frac{M\Omega^{2}}{2}\left(1-\delta^{2}\right)y^{2}.$
(48)

The effect of the non-adiabatic coupling is twofold. On one hand, the higher
frequency oscillator displays an effective frequency
$\tilde{\omega}_{1}=\omega\left(1+\frac{\Omega^{2}\delta^{2}}{2\omega^{2}}\right)$,
with $\omega<\tilde{\omega}_{1}$, on the other hand, the lower frequency
oscillator shows an effective frequency
$\tilde{\omega}_{2}=\Omega\sqrt{1-\delta^{2}}$, with
$\tilde{\omega}_{2}<\Omega$.

####  III.2.3 The $\varphi_{nl}(y)$ and $\Psi_{nl}(x,y)$ functions

Bound solutions of the system of equations (46) are given by
$\varphi_{nl}=D_{\mu}(\tilde{y})$, with

$\displaystyle\mu$
$\displaystyle=-\frac{1}{2}\pm\frac{\mathcal{E}-\left(n+\frac{1}{2}\right)\left(\omega+\frac{\delta^{2}\Omega^{2}}{2\omega}\right)}{\Omega\sqrt{1-\delta^{2}}},$
(49) $\displaystyle\tilde{y}$
$\displaystyle=y\sqrt{2M\Omega}\left(1-\delta^{2}\right)^{1/4}.$ (50)

The condition $\mu=l(n)$, with $l(n)=0,1,2,...$, gives the Born-Oppenheimer
energy

$\begin{split}\tilde{\mathcal{E}}_{nl}&=\omega\left(1+\frac{\delta^{2}\Omega^{2}}{2\omega^{2}}\right)\left(n+\tfrac{1}{2}\right)+\Omega\sqrt{1-\delta^{2}}\left(l+\tfrac{1}{2}\right)\\\
&=\tilde{\omega}_{1}\left(n+\tfrac{1}{2}\right)+\tilde{\omega}_{2}\left(l+\tfrac{1}{2}\right).\end{split}$
(51)

The second line of equation (51) defines the Born-Oppenheimer frequencies, or
energy differences between consecutive levels, $\tilde{\omega}_{1}$ and
$\tilde{\omega}_{2}$. It can be proved that
$\tilde{\omega}_{2}<\tilde{\omega}_{1}$ under the condition $\Omega<\omega$.
It is useful to define the matrix of BO square-root frequencies by
$\tilde{\mathsf{F}}=\mathop{\mathrm{diag}}{(\sqrt{\tilde{\omega}_{1}},\sqrt{\tilde{\omega}_{2}})}$.

For the uncoupled case, a simple inspection of equation (51) gives

$\tilde{\mathcal{E}}_{nl}=\omega\left(n+\tfrac{1}{2}\right)+\Omega\left(l+\tfrac{1}{2}\right).$
(52)

The normalized solution for the $y$ degree of freedom is given by

$\begin{split}\varphi_{nl}(\tilde{y})&=\frac{1}{\sqrt{l!}}\frac{\sqrt{\partial_{y}\tilde{y}}}{(2\pi)^{1/4}}D_{l}(\tilde{y})\\\
&=\frac{\left(4M^{2}\Omega^{2}\left(1-\delta^{2}\right)\right)^{1/8}}{\sqrt{l!}(2\pi)^{1/4}}D_{l}(\tilde{y}).\end{split}$
(53)

The full Born-Oppenheimer product wavefunction is then

$\begin{split}\chi_{n}(z)\varphi_{nl}(\tilde{y})&=\frac{\left(4M^{2}\Omega^{2}\left(1-\delta^{2}\right)\right)^{1/8}}{\sqrt{l!}\sqrt{n!}}\left(\frac{m\omega}{2\pi^{2}}\right)^{1/4}D_{n}(z)D_{l}(\tilde{y})\\\
&=\frac{\left(mM\omega\Omega\sqrt{1-\delta^{2}}\right)^{1/4}}{\sqrt{\pi
n!l!}}D_{n}(z)D_{l}(\tilde{y}).\end{split}$ (54)

The arguments of the parabolic cylinder functions $D_{n}$ and $D_{l}$, $z$ and
$\tilde{y}$ respectively, define the vector,
$\tilde{\bm{\xi}}=\begin{pmatrix}z\\\ \tilde{y}\end{pmatrix}$, which is
related to the position vector $\mathbf{r}$ by
$\tilde{\bm{\xi}}=\tilde{\mathsf{G}}\mathbf{r}$ with

$\tilde{\mathsf{G}}=\begin{pmatrix}\sqrt{2m\omega}&\sqrt{\frac{2M\Omega^{2}}{\omega}}\delta\\\
0&\sqrt{2M\Omega}\left(1-\delta^{2}\right)^{\nicefrac{{1}}{{4}}}\end{pmatrix}.$
(55)

The volume element of the BO coordinates $\tilde{\bm{\xi}}$ is given by
$dzd\tilde{y}=2(1-\delta^{2})^{\nicefrac{{1}}{{4}}}\sqrt{mM\omega\Omega}\,dxdy$.
The use of this volume element, and the product (54) gives the normalized BO
solution

$\tilde{\Psi}_{nl}(z,\tilde{y})=\frac{1}{\sqrt{2\pi
n!l!}}D_{n}(z)D_{l}(\tilde{y}).$ (56)

The vector of arguments of the exact and BO wavefunctions, $\bm{\xi}$ and
$\tilde{\bm{\xi}}$ respectively, are related by using the result
$\mathbf{r}=\tilde{\mathsf{G}}^{-1}\tilde{\bm{\xi}}$ in the definition
$\bm{\xi}=\mathsf{G}\mathbf{r}$, to obtain

$\bm{\xi}=\mathsf{G}\tilde{\mathsf{G}}^{-1}\tilde{\bm{\xi}}.$ (57)

It can be proved that the determinant of the product $\tilde{\mathsf{G}}^{-1}$
is one. In analogy to the exact solution (32), the matrix $\tilde{\mathsf{G}}$
can be written as the product
$\tilde{\mathsf{G}}=\sqrt{2}\tilde{\mathsf{F}}\tilde{\mathsf{R}}\mathsf{M}^{\nicefrac{{1}}{{2}}}$,
which allows to obtain $\tilde{\mathsf{R}}$ as

$\displaystyle\tilde{\mathsf{R}}$
$\displaystyle=\tfrac{1}{\sqrt{2}}\tilde{\mathsf{F}}^{-1}\tilde{\mathsf{G}}\mathsf{M}^{\nicefrac{{-1}}{{2}}}$
(58)
$\displaystyle=\begin{pmatrix}\frac{\omega}{\sqrt{\omega^{2}+\delta^{2}\Omega^{2}/2}}&\frac{\delta\Omega}{\sqrt{\omega^{2}+\delta^{2}\Omega^{2}/2}}\\\
0&1\end{pmatrix}.$ (59)

This upper triangular matrix is clearly non orthogonal as matrix $\mathsf{R}$.

###  III.3 Error analysis of the Born-Oppenheimer solutions

####  III.3.1 Error of the Born-Oppenheimer energies

The exact and BO energies for the uncoupled case, equations (30) and (52)
respectively, are identical if the replacements $n_{1}\to n$ and $n_{2}\to l$
is made on equation (30). In the general case, equation (29), the exact
energies are given by

$\mathcal{E}_{nl}=\omega_{1}\left(n+\tfrac{1}{2}\right)+\omega_{2}\left(l+\tfrac{1}{2}\right).$
(60)

The relative error in the BO energy levels is given by

$\begin{split}\epsilon_{nl}^{BO}&=\frac{\mathcal{E}_{nl}^{BO}-\mathcal{E}_{nl}}{\mathcal{E}_{nl}}\\\
&=\frac{\epsilon_{\omega}\left(n+\tfrac{1}{2}\right)+\epsilon_{\Omega}\left(l+\tfrac{1}{2}\right)}{\bar{\omega}_{1}\left(n+\tfrac{1}{2}\right)+\bar{\omega}_{2}\left(l+\tfrac{1}{2}\right)}.\end{split}$
(61)

with the definitions

$\displaystyle\epsilon_{\omega}$
$\displaystyle=\omega^{-1}\left(\tilde{\omega}_{1}-\omega_{1}\right)$ (62)
$\displaystyle=\left(1+\tfrac{\delta^{2}\bar{\Omega}^{2}}{2}\right)-\bar{\omega}_{1},$
(63) $\displaystyle\epsilon_{\Omega}$
$\displaystyle=\omega^{-1}\left(\tilde{\omega}_{2}-\omega_{2}\right)$ (64)
$\displaystyle=\bar{\Omega}\sqrt{1-\delta^{2}}-\bar{\omega}_{2},$ (65)

and

$\begin{split}\bar{\omega}_{1,2}&=\omega_{1,2}/\omega\\\
&=\tfrac{1}{\sqrt{2}}\left(1+\bar{\Omega}^{2}\pm\sqrt{\left(1-\bar{\Omega}^{2}\right)^{2}+4\delta^{2}\bar{\Omega}^{2}}\right)^{1/2},\end{split}$
(66)

in terms of $\bar{\Omega}=\Omega/\omega$. The relative error of the BO
energies, equation (61), defines a surface in terms of $n$ and $l$, for a
given value of $\delta$ and $\bar{\Omega}$. The magnitude and sign of the BO
error can be analyzed by a detailed examination of the terms of equation (61).
Since the denominator of $\epsilon_{nl}^{BO}$ is a positive quantity, the
analysis of the sign of the BO error depends entirely on the two terms of the
numerator of (61). For the second term of the numerator, the inequality
$\epsilon_{\Omega}>0$ leads to the result
$\delta^{2}\Omega^{2}(\delta^{2}-1)<0$, which holds for $0<|\delta|<1$ ,
indicating that the BO result overestimates the energy of the $\omega_{2}$
oscillator for any value of the parameters $\delta$ and $\bar{\Omega}$. For
the first term of the numerator, the quantity $\epsilon_{\omega}$ is greater,
or less, than zero in the region of the $\delta-\bar{\Omega}$ plane implicitly
defined by $\mathcal{B}(\delta,\bar{\Omega})>0$, or
$\mathcal{B}(\delta,\bar{\Omega})<0$ respectively, with
$\mathcal{B}(\delta,\bar{\Omega})=4(5\delta^{2}-4)+4(2\delta^{2}-1)\bar{\Omega}^{2}+\delta^{6}\bar{\Omega}^{4}$,
and $\bar{\Omega}<1$. Figure 1 displays regions of the $\delta-\bar{\Omega}$
plane where the BO energy overestimates $\epsilon_{\omega}>0$, or
underestimates $\epsilon_{\omega}<0$, the energy of the oscillator
$\omega_{1}$, it can be seen in this figure, that the BO energy overestimate
the exact energy $\omega_{1}$ for large values of the coupling $\delta$.

Figure 1: Regions of the $\delta-\bar{\Omega}$ plane where the Born-
Oppenheimer approximation overestimates, $\epsilon_{\omega}>0$ in light gray,
or underestimate, $\epsilon_{\omega}<0$ in dark gray, the exact energy
$\omega_{2}$. The BO energies overestimate the exact energy of the
$\omega_{2}$ oscillator at large values of the coupling $\delta$. This figure
is used also to show the regions of the $\delta-\bar{\Omega}$ plane where the
sign of the slope of line (67) is positive (dark gray) or negative (light
gray).

The boundary between regions of the $n-l$ plane with positive and negative
relative error is given by the condition $\epsilon^{BO}=0$. This condition
gives the equation of a straight line,

$l=-\frac{\epsilon_{\omega}}{\epsilon_{\Omega}}n-\frac{1}{2}\left(1+\frac{\epsilon_{\omega}}{\epsilon_{\Omega}}\right),$
(67)

with slope $-\epsilon_{\omega}/\epsilon_{\Omega}$, and $l$-intercept
$-\frac{1}{2}\left(1+\frac{\epsilon_{\omega}}{\epsilon_{\Omega}}\right)$.
Since $\epsilon_{\Omega}>0$, the sign of the slope of line (67) depends only
on $\text{sgn}(-{\epsilon_{\omega}})$. Figure 1 displays the regions of the
$\delta-\bar{\Omega}$ plane where the error $\epsilon_{\omega}$ is positive or
negative, this figure also shows the regions of the $\delta-\bar{\Omega}$
plane where the slope of the line (67) is positive or negative. It can be seen
in Figure 1 that the slope of line (67) is positive only for large values of
$\delta$. It can be proved after some algebraic manipulations that the
inequality $\left(1+\frac{\epsilon_{\omega}}{\epsilon_{\Omega}}\right)>0$
holds for $0<\bar{\Omega}<1$ and $0<|\delta|<1$, which shows that the upper
limit of the slope is one, i.e.
$-\frac{\epsilon_{\omega}}{\epsilon_{\Omega}}<1$, and that the $l$-intercept
is always negative. These two results about the slope and the $l$-intercept of
line (67) show that the BO approximation overestimates the energy of the
ground state for linearly coupled oscillators. The overestimation of the BO
ground state energy contrasts with previous works: Bratsev (25), Epstein (26)
and McCoy (27), who obtained that the BO ground state energy is always smaller
than the true ground state energy. Figure 2 displays the contours of constant
relative error $\epsilon_{nl}^{BO}$ as functions of $n$ and $l$ for
$\delta=0.6$ and $\bar{\Omega}=0.2$, the zero contour, $\epsilon_{nl}^{BO}=0$,
gives the straight line (67). The relative error of the BO states are less
than one percent, as can be seen in figure 2 with the aid of the legend. In
figure 2 it is seen that BO states with $l=0$ display the smallest relative
error, while the quantum states with $n=0$ show the biggest relative errors.

Figure 2: Contours of constant relative error, equation (61), in the $n-l$
plane for $\delta=0.6$ and $\bar{\Omega}=0.2$. The lattice of quantum states
is displayed as black points. Colored regions, as indicated by the legend,
give the range of the error for each quantum state.

####  III.3.2 Wave functions: error analysis

Figure 3: The inner product $\sigma_{nl}$ as a function of $\delta$ for
$\bar{\Omega}=0.2$. In all the panels: the blue, orange, green, red, purple,
and brown curves represent the cases $l=0,1,2,3,4,5$, respectively. Panels a,
b, c, d, e, and f, correspond to $n=0,1,2,3,4$, and $5$, respectively. The
smallest error is given by the ground state, $\sigma_{00}$. For a given value
of $\delta$ the $\sigma_{nl}$ decreases as the quantum number $n$ and $l$
increase.

The exact and BO wave functions, equations (32) and (54), are products of
parabolic cylinder functions $D_{\nu}$, with $\nu$ non negative integer. Since
the $D_{\nu}$ are $L^{2}$ functions, the exact and the BO functions must be
$L^{2}$ functions, allowing to use the inner product
$\sigma_{nl}=\braket{\Psi_{nl}}{\tilde{\Psi}_{nl}}$ as a measure of error of
the BO function compared to the exact wavefunctions. Figure 3 displays the
value of $\sigma_{nl}$ as function of $\delta$ for $\bar{\Omega}=0.2$. In all
the panels of the figure, the value of $\sigma_{nl}$ is very close to one,
which reinforces the observed in figure 2, and shows that the BO approximation
works very well for linearly coupled harmonic oscillators. In all the panels
of this figure is clear that for a given value of $\delta$, the higher the
values of $n$ and $l$, the more the BO and exact wavefunctions differ. It is
interesting to compare the relative error of figure 2 with the wavefunction
overlap of figure 3 for the same values of $n$ and $l$. The relative error for
$n=0$ displays the largest values, meanwhile the wavefunction overlap displays
larger values as $n$ and $l$ grow larger.

Figure 4: Components of the matrices $\mathsf{G}$ (solid) and
$\tilde{\mathsf{G}}$ (dashed) as functions of the coupling parameter $\delta$
for $\bar{\Omega}=0.2$. Components $\mathsf{G}_{11}$, $\mathsf{G}_{12}$,
$\mathsf{G}_{21}$ and $\mathsf{G}_{22}$ are displayed in blue, orange, green
and red, respectively, with same color code for the elements of
$\tilde{\mathsf{G}}$.

All the panels of figure 3 show that $\lim_{\delta\to 0}\sigma_{nl}=1$, and
$\lim_{\delta\to 1}\sigma_{nl}=1$. This can be explained by the vectors of
arguments of the parabolic cylinder functions:
$\bm{\xi}=\mathsf{G}\mathbf{r}$, and
$\tilde{\bm{\xi}}=\tilde{\mathsf{G}}\mathbf{r}$. The components of the
transformation matrices $\mathsf{G}$ and $\tilde{\mathsf{G}}$ are displayed in
figure 4. Three of the four elements of the BO transformation matrix
$\tilde{\mathsf{G}}$ display excellent agreement with the elements of
$\mathsf{G}$. Equation (55) show that the BO $\tilde{\mathsf{G}}_{21}$ element
is always zero, as shown by the dashed green curve of figure 4. The green
solid curve shows that $\mathsf{G}_{21}$ is negative for all values of
$\delta$, being zero only at the endpoints $\delta=0$ and $\delta=1$. The
component $\mathsf{G}_{21}$ displays its lowest value for $\delta\approx 0.8$,
which agrees fairly well with the panels of figure 3. More precisely, the
matrix element $\mathsf{G}_{21}$ can be minimized with respect to $\delta$ to
obtain

$\delta^{*}=\frac{\sqrt{\left(1-\bar{\Omega}^{2}\right)\left(3+5\bar{\Omega}^{2}-\sqrt{9-2\bar{\Omega}^{2}-7\bar{\Omega}^{4}}\right)}}{2\sqrt{2}\bar{\Omega}},$
(68)

with $\bar{\Omega}=\Omega/\omega$. The case of figure 4, with
$\bar{\Omega}=0.2$, gives $\delta^{*}=\tfrac{2}{5}\sqrt{60-6\sqrt{87}}\approx
0.803565$, and $\mathsf{G}^{*}_{21}=-0.174359\sqrt{m}$.

###  III.4 $\mathsf{QU}$ decomposition of matrix $\mathsf{G}$

Matrices $\mathsf{G}$ and $\tilde{\mathsf{G}}$, equations (33) and (55)
respectively, are related by the $\mathsf{QU}$ decomposition of $\mathsf{G}$
(28),

$\mathsf{G}=\mathsf{QU},$ (69)

with $\mathsf{Q}$ an orthogonal matrix and $\mathsf{U}$ an upper triangular
matrix. Let define the matrix $\mathsf{\tilde{U}}=-\sigma_{3}\mathsf{U}$, with
$\sigma_{3}=\mathop{\mathrm{diag}}{(1,-1)}$ as the Pauli’s diagonal matrix. It
is observed that the BO matrix $\tilde{\mathsf{G}}$ is approximately given by
$\mathsf{\tilde{U}}$, $\tilde{\mathsf{G}}\approx\mathsf{\tilde{U}}$. Figure 5
displays the elements of matrices $\mathsf{\tilde{G}}$ and
$\mathsf{\tilde{U}}$, $\tilde{G}_{ij}$ and $\tilde{U}_{ij}$ respectively, as
function of $\delta$ for $\bar{\Omega}=0.2$. As the value of $\bar{\Omega}$
increases, the approximation $\tilde{\mathsf{G}}\approx\mathsf{\tilde{U}}$
becomes more accurate.

The $\mathsf{QU}$ decomposition can be used to improve the BO wavefunction.
The fact that $\tilde{\mathsf{G}}\approx\mathsf{\tilde{U}}$ can be used to
improve the BO matrix $\tilde{\mathsf{G}}$ by using the $\mathsf{QU}$
decomposition inversely. Multiplying the matrix
$\mathsf{U}=-\sigma_{3}\tilde{\mathsf{G}}$ by an orthogonal transformation
$\mathsf{Q}$, for example the rotation matrix $\mathsf{Q}(\phi)$, is possible
to improve the Born-Oppenheimer approximation. This approach is being analyzed
in further details by the author and collaborators.

Figure 5: Elements of matrices $\mathsf{\tilde{G}}$ and $\mathsf{\tilde{U}}$
as function of $\delta$ for $\bar{\Omega}=0.2$. Solid and dashed curves
represent $\mathsf{\tilde{G}}$ and $\mathsf{\tilde{U}}$ respectively.
$\tilde{G}_{11}$ and $\tilde{U}_{11}$ are displayed in red, $\tilde{G}_{12}$
and $\tilde{U}_{12}$ in green, $\tilde{G}_{21}$ and $\tilde{U}_{21}$ in blue,
and $\tilde{G}_{22}$ and $\tilde{U}_{22}$ in black.

##  IV Conclusions

It has been proposed a method for finding bound states, for systems with two
degrees of freedom, based on the BO factorization ansatz and the phase-space
bound trajectory approach. This method has been tested by using a system
consisting of two 1-dimensional harmonic oscillators bilinearly coupled. The
use of bilinearly coupled harmonic oscillators makes straightforward to solve
the phase-space differential equations for both oscillators, since the non-
adiabatic coupling $\beta_{nn}$ is a constant. The systems of bilinearly
coupled oscillators is analytically solvable in terms of a linear
transformation of coordinates, this is also the case for the solution obtained
by using the Born-Oppenheimer ansatz and the phase-space bound trajectory
method. The analytical non-adiabatic eigenvalues and eigenfunctions have been
compared with the analytical solutions showing an excellent accuracy of the
BO-phase-space trajectory approach for the system of two bilinearly coupled
harmonic oscillators. An extensive mathematical analysis of the energy
eigenvalues and their relative errors has been done. This analysis has been
useful identifying regions on the parameter space where the BO approximation
overestimate or underestimate the exact quantum energy. The relative error of
the BO approximation for the lattice of quantum states has been analyzed in
detail. It has been shown that the relative error is mostly positive and grows
larger as the quantum numbers increase. The Born-Oppenheimer wave functions
are compared with the exact quantum wavefunctions by using an overlap integral
as a measure of error. The error of the BO wavefunctions has been obtained as
a function of the coupling parameters $\delta$. The use of transformation
matrices between the physical coordinates and the arguments of the parabolic
cylinder functions has allowed to analyze in detail the error in the
wavefunction by comparing the quantum and BO matrix elements, matrices
$\mathsf{G}$ and $\mathsf{\tilde{G}}$ respectively. A $\mathsf{QU}$
decomposition of matrix $\mathsf{G}$ outputs a upper triangular matrix
$\mathsf{U}$ which is closely related to the BO transformation matrix
$\mathsf{\tilde{G}}$. The use of the $\mathsf{QU}$ decomposition to improve
the Born-Oppenheimer approximation has been outlined, this issue is currently
investigated by the author and collaborators.

##  V Acknowledgments

This project has been fully financed by the internal research grants of
Universidad Icesi.

##  VI Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  VII Appendix: Calculations with the $D_{\nu}$ functions

This appendix displays some of the mathematical facts, equations, and
procedures used in section III.2

###  VII.1 Recurrence relations for the $D_{\nu}$ functions

$\displaystyle D_{\nu+1}(z)-zD_{\nu}(z)+\nu D_{\nu-1}(z)$ $\displaystyle=0,$
(70) $\displaystyle D^{\prime}_{\nu}(z)+\tfrac{z}{2}D_{\nu}(z)-\nu
D_{\nu-1}(z)$ $\displaystyle=0,$ (71) $\displaystyle
D^{\prime}_{\nu}(z)-\tfrac{z}{2}D_{\nu}(z)+D_{\nu+1}(z)$ $\displaystyle=0.$
(72)

###  VII.2 Matrix elements of the $D_{\nu}$ functions

The matrix elements $\braket{D_{l}(z)}{xD_{n}(z)}$ and
$\braket{D_{l}(z)}{x^{2}D_{n}(z)}$ (integral over $x$) are obtained from the
recurrence relations (70), and the definition of $z$ given by equation (37)
and (40),

$\displaystyle\braket{D_{l}(z)}{xD_{n}(z)}$
$\displaystyle=\frac{1}{\partial_{x}z}\int_{-\infty}^{\infty}{xD_{l}D_{n}dz}$
(73)
$\displaystyle=\frac{1}{(\partial_{x}z)^{2}}\int_{-\infty}^{\infty}{\left(z-(\partial_{y}z)y\right)D_{l}D_{n}dz}$
(74)
$\displaystyle=\frac{1}{(\partial_{x}z)^{2}}\int_{-\infty}^{\infty}{zD_{l}D_{n}dz}-\frac{(\partial_{y}z)y}{(\partial_{x}z)^{2}}\int_{-\infty}^{\infty}{D_{l}D_{n}dz}$
(75)
$\displaystyle=\frac{(\partial_{y}x)y}{(\partial_{x}z)}\sqrt{2\pi}n!\delta_{ln},$
(76)

and

$\displaystyle\braket{D_{l}(z)}{x^{2}D_{n}(z)}$
$\displaystyle=\frac{1}{\partial_{x}z}\int_{-\infty}^{\infty}{x^{2}D_{l}D_{n}dz}$
(77)
$\displaystyle=\frac{1}{(\partial_{x}z)^{2}}\int_{-\infty}^{\infty}{\left(z-(\partial_{y}z)y\right)^{2}D_{l}D_{n}dz}$
(78)
$\displaystyle=\frac{1}{(\partial_{x}z)^{2}}\int_{-\infty}^{\infty}{z^{2}D_{l}D_{n}dz}+\frac{(\partial_{y}z)^{2}y^{2}}{(\partial_{x}z)^{2}}\int_{-\infty}^{\infty}{D_{l}D_{n}dz}$
(79)
$\displaystyle=\left(\frac{y\partial_{y}z}{\partial_{x}z}\right)^{2}\sqrt{2\pi}n!\delta_{ln}+\frac{1}{(\partial_{x}z)^{2}}\int_{-\infty}^{\infty}{(zD_{l})(zD_{n})dz}$
(80)
$\displaystyle=(y\partial_{y}x)^{2}\sqrt{2\pi}n!\delta_{ln}+\frac{1}{(\partial_{x}z)^{2}}\int_{-\infty}^{\infty}{(D_{l+1}+lD_{l-1})(D_{n+1}+nD_{n-1})dz}$
(81) $\displaystyle=(y\partial_{y}x)^{2}\sqrt{2\pi}n!\delta_{ln}$ (82)
$\displaystyle+\frac{1}{(\partial_{x}z)^{2}}\int_{-\infty}^{\infty}{(D_{l+1}D_{n+1}+nD_{l+1}D_{n-1}+lD_{l-1}D_{n+1}+lnD_{l-1}D_{n-1})dz}$
(83) $\displaystyle=(y\partial_{y}x)^{2}\sqrt{2\pi}n!\delta_{ln}$ (84)
$\displaystyle+\frac{\sqrt{2\pi}}{(\partial_{x}z)^{2}}\left((n+1)!\delta_{l+1,n+1}+n!\delta_{l+1,n-1}+l(n+1)!\delta_{l-1,n+1}+ln!\delta_{l-1,n-1}\right)$
(85) $\displaystyle=(y\partial_{y}x)^{2}\sqrt{2\pi}n!\delta_{ln}$ (86)
$\displaystyle+\frac{\sqrt{2\pi}n!}{(\partial_{x}z)^{2}}\left((n+1)\delta_{l+1,n+1}+\delta_{l+1,n-1}+l(n+1)\delta_{l-1,n+1}+l\delta_{l-1,n-1}\right)$
(87)

###  VII.3 Derivatives of the $\chi(z)$ function

By using the chain rule, and the recurrence relations (70), we obtain the
derivatives of $\chi_{n}(z)$,

$\begin{split}\partial_{y}\chi_{n}(z)&=\frac{1}{\sqrt{n!}}\frac{\sqrt{\partial_{x}z}}{(2\pi)^{1/4}}D_{n}^{\prime}(\partial_{y}z)\\\
&=\frac{(\partial_{y}z)}{\sqrt{n!}}\frac{\sqrt{\partial_{x}z}}{(2\pi)^{1/4}}\left(\tfrac{1}{2}zD_{n}-D_{n+1}\right)\\\
&=\frac{(\partial_{y}z)}{2\sqrt{n!}}\frac{\sqrt{\partial_{x}z}}{(2\pi)^{1/4}}\left(nD_{n-1}-D_{n+1}\right)\\\
&=\frac{(\partial_{y}z)}{2}\left(\sqrt{n}\chi_{n-1}-\sqrt{n+1}\chi_{n+1}\right)\end{split}$
(88)

and

$\begin{split}\partial^{2}_{y}\chi_{n}(z)&=\frac{1}{\sqrt{n!}}\frac{\sqrt{\partial_{x}z}}{(2\pi)^{1/4}}\partial_{y}(D^{\prime}_{n}\partial_{y}z)\\\
&=\frac{1}{\sqrt{n!}}\frac{\sqrt{\partial_{x}z}}{(2\pi)^{1/4}}(D^{\prime\prime}_{n}(\partial_{y}z)^{2}+D^{\prime}_{n}\partial_{y}^{2}z)\\\
&=\frac{(\partial_{y}z)^{2}}{\sqrt{n!}}\frac{\sqrt{\partial_{x}z}}{(2\pi)^{1/4}}D^{\prime\prime}_{n}\\\
&=-\frac{\left(\partial_{y}z\right)^{2}}{\sqrt{n!}}\frac{\sqrt{\partial_{x}z}}{(2\pi)^{1/4}}\left((n+\tfrac{1}{2})-(\tfrac{z}{2})^{2}\right)D_{n}\\\
&=-\left((n+\tfrac{1}{2})-(\tfrac{z}{2})^{2}\right)\left(\partial_{y}z\right)^{2}\chi_{n}.\end{split}$
(89)

####  VII.3.1 Matrix elements of $\chi_{n}(x,y)$

The use of the chain rule, and the recurrence relations (70), allow to obtain
the derivatives of $\chi_{n}(z)$,

$\begin{split}\partial_{y}\chi_{n}(z)&=\frac{(\partial_{y}z)}{2}\left(\sqrt{n}\chi_{n-1}-\sqrt{n+1}\chi_{n+1}\right)\\\
\partial^{2}_{y}\chi_{n}(z)&=-\left((n+\tfrac{1}{2})-(\tfrac{z}{2})^{2}\right)\left(\partial_{y}z\right)^{2}\chi_{n}.\end{split}$
(90)

The matrix elements of equations (14) and (16) are

$\begin{split}\braket{\chi_{l}}{\partial_{y}\chi_{n}}&=\frac{(\partial_{y}z)}{2}\left(\sqrt{n}\braket{\chi_{l}}{\chi_{n-1}}-\sqrt{n+1}\braket{\chi_{l}}{\chi_{n+1}}\right)\\\
&=\frac{(\partial_{y}z)}{2}\left(\sqrt{n}\delta_{l,n-1}-\sqrt{n+1}\delta_{l,n+1}\right),\end{split}$
(91)

and

$\begin{split}\braket{\chi_{l}}{\partial^{2}_{y}\chi_{n}}&=-\left(\partial_{y}z\right)^{2}\left((n+\tfrac{1}{2})\braket{\chi_{l}}{\chi_{n}}-\tfrac{1}{4}\braket{\chi_{l}}{z^{2}\chi_{n}}\right)\\\
&=-\left(\partial_{y}z\right)^{2}\left((n+\tfrac{1}{2})\delta_{l,n}-\tfrac{1}{4}\braket{\chi_{l}}{z^{2}\chi_{n}}\right).\end{split}$
(92)

Explicitly, the last integral on the RHS of equation (92) is given by

$\begin{split}\braket{\chi_{l}}{z^{2}\chi_{n}}&=\frac{(\partial_{x}z)}{\sqrt{2\pi
l!n!}}\int_{-\infty}^{\infty}{D_{l}(z)z^{2}D_{n}(z)dx}\\\
&=\frac{1}{\sqrt{2\pi l!n!}}\int_{-\infty}^{\infty}{z^{2}D_{l}D_{n}dz}\\\
&=\frac{\sqrt{n!}}{\sqrt{l!}}\left((n+1)\delta_{l+1,n+1}+\delta_{l+1,n-1}+l(n+1)\delta_{l-1,n+1}+l\delta_{l-1,n-1}\right).\end{split}$
(93)

The diagonal elements of matrices (91) and (92) are
$\braket{\chi_{n}}{\partial_{y}\chi_{n}}=0$, and

$\braket{\chi_{n}}{\partial^{2}_{y}\chi_{n}}=-\left(\frac{\partial_{y}z}{2}\right)^{2}(2n+1).$
(94)

The nonzero off-diagonal elements of matrices (91) and (92) are

$\displaystyle\braket{\chi_{n-1}}{\partial_{y}\chi_{n}}$
$\displaystyle=\frac{\partial_{y}z}{2}\sqrt{n},$ (95)
$\displaystyle\braket{\chi_{n+1}}{\partial_{y}\chi_{n}}$
$\displaystyle=-\frac{\partial_{y}z}{2}\sqrt{n+1},$ (96)

and

$\displaystyle\braket{\chi_{n-2}}{\partial^{2}_{y}\chi_{n}}$
$\displaystyle=\left(\frac{\partial_{y}z}{2}\right)^{2}n(n-1),$ (97)
$\displaystyle\braket{\chi_{n+2}}{\partial^{2}_{y}\chi_{n}}$
$\displaystyle=\left(\frac{\partial_{y}z}{2}\right)^{2}.$ (98)

The use of these matrix elements in equation (16) produces

$\displaystyle-\braket{\chi_{n-1}}{\partial_{y}^{2}\chi_{n}}\varphi_{n}$
$\displaystyle=2\braket{\chi_{n-1}}{\partial_{y}\chi_{n}}\vartheta_{n}$ (99)
$\displaystyle 0$ $\displaystyle=(\partial_{y}z)\sqrt{n}\vartheta_{n}$ (100)
$\displaystyle\vartheta_{n}$ $\displaystyle=0.$ (101)

and

$\displaystyle-\braket{\chi_{n+1}}{\partial_{y}^{2}\chi_{n}}\varphi_{n}$
$\displaystyle=2\braket{\chi_{n+1}}{\partial_{y}\chi_{n}}\vartheta_{n}$ (102)
$\displaystyle 0$ $\displaystyle=-(\partial_{y}z)\sqrt{n+1}\vartheta_{n}$
(103) $\displaystyle\vartheta_{n}$ $\displaystyle=0.$ (104)

## References

  * Uzer and Miller (1991) T. Uzer and W. Miller, Physics Reports 199, 73 (1991).
  * Rice (1981) S. A. Rice, “An overview of the dynamics of intramolecular transfer of vibrational energy: Dynamics of intramolecular transfer of vibrational energy,” in  _Advances in Chemical Physics_, Vol. 47, edited by J. Jortner, R. D. Levine, and S. A. Rice (John Wiley & Sons, Ltd, 1981) Chap. 2, pp. 117–200.
  * Herzberg (1990) G. Herzberg,  _Molecular Spectra and Molecular Structure: Volume II: Infrared and Raman Spectra of Polyatomic Molecules_ , Vol. 2 (Krieger Publishing Company, 1990).
  * Sathyanarayana (2015) D. N. Sathyanarayana,  _Vibrational spectroscopy: theory and applications_ (New Age International, 2015).
  * Makarov (2018) D. N. Makarov, Phys. Rev. E 97, 042203 (2018).
  * Chakraborty and Sarma (2018) S. Chakraborty and A. K. Sarma, Phys. Rev. A 97, 022336 (2018).
  * Moshinsky and Smirnov (1996) M. Moshinsky and Y. F. Smirnov,  _The Harmonic Oscillator in Modern Physics_, revised ed., Contemporary Concepts in Physics (CRC Press, 1996).
  * Feldmann, Cohen, and Burrows (2000) T. Feldmann, M. Cohen, and B. L. Burrows, Journal of Mathematical Physics 41, 5897 (2000).
  * Burrows and Cohen (2000) B. L. Burrows and M. Cohen, Molecular Physics 98, 1083 (2000).
  * Burrows, Cohen, and Feldmann (2003) B. L. Burrows, M. Cohen, and T. Feldmann, International Journal of Quantum Chemistry 92, 345 (2003).
  * Zúñiga, Bastida, and Requena (2017) J. Zúñiga, A. Bastida, and A. Requena, J. Math. Chem. 55, 1964–1984 (2017).
  * Caswell and Danos (1970) R. S. Caswell and M. Danos, Journal of Mathematical Physics 11, 349 (1970).
  * Bowman (1978) J. M. Bowman, The Journal of Chemical Physics 68, 608 (1978).
  * Zhi-Ding  _et al._ (1982) Q. Zhi-Ding, Z. Xing-Guo, L. Xing-Wen, H. Kono, and S. Lin, Molecular Physics 47, 713 (1982).
  * Ezra (1983) G. S. Ezra, Chemical Physics Letters 101, 259 (1983).
  * Certain and Moiseyev (1987) P. R. Certain and N. Moiseyev, The Journal of Chemical Physics 86, 2146 (1987).
  * Makarewicz (1985) J. Makarewicz, Journal of Physics A: Mathematical and General 18, 557 (1985).
  * Craig and Thirunamachandran (1994) D. Craig and T. Thirunamachandran, Theoret. Chim. Acta 89, 123–136 (1994).
  * Bowman, Christoffel, and Tobin (1979) J. M. Bowman, K. Christoffel, and F. Tobin, The Journal of Physical Chemistry 83, 905 (1979).
  * Born and Oppenheimer (1927) M. Born and R. Oppenheimer, Annalen der Physik 389, 457 (1927).
  * Frey and Holdship (1988) J. G. Frey and S. J. Holdship, Molecular Physics 64, 191 (1988).
  * Molano and Arango (2021) J. S. Molano and C. A. Arango, Chemical Physics Letters 762, 138171 (2021).
  * Strogatz (1994) S. Strogatz,  _Nonlinear Dynamics and Chaos_ (Addison-Wesley, Reading, MA, 1994).
  * Abramowitz and Stegun (1964) M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ , ninth dover printing, tenth gpo printing ed. (Dover, New York City, 1964).
  * Bratsev (1965) V. F. Bratsev, Dokl. Akad. Nauk SSSR 160, 570 (1965).
  * Epstein (1966) S. T. Epstein, The Journal of Chemical Physics 44, 836 (1966).
  * McCoy (2001) A. B. McCoy, Journal of Chemical Education 78, 401 (2001).
  * Strang (2019) G. Strang,  _Linear Algebra and Learning from Data_ , 1st ed. (Wellesley - Cambridge Press, Wellesley, MA, USA, 2019).


# Hyperorder net baryon number fluctuations in nuclear matter at low
temperature

Xin-ran Yang  School of Physics, Xi’an Jiaotong University, Xi’an, 710049,
China  MOE Key Laboratory for Nonequilibrium Synthesis and Modulation of
Condensed Matter, Xi’an Jiaotong University, Xi’an, 710049, China  Guo-yun
Shao<EMAIL_ADDRESS>School of Physics, Xi’an Jiaotong University,
Xi’an, 710049, China  MOE Key Laboratory for Nonequilibrium Synthesis and
Modulation of Condensed Matter, Xi’an Jiaotong University, Xi’an, 710049,
China  Wei-bo He  School of Physics, Peking University, Beijing, 100871, China

###### Abstract

We calculate the density fluctuations of net baryon number up to sixth order
induced by the interactions of nuclear matter, and explore their relationship
with the nuclear liquid-gas phase transition (LGPT), including the stable and
metastable phase as well as the region far from the phase transition. The
results show that dramatic density fluctuations exist in the vicinity of LGPT,
and the hyperorder density fluctuations are more sensitive than the lower
order ones to the interactions and structural properties of nuclear matter.
The study also indicates that, even far away from the critical region of LGPT,
the hadronic interactions can still lead to significant hyperorder density
fluctuations. In combination with the chemical freeze-out line fitted from the
experimental data, the derived results can be referred to investigate the
chiral phase transition, nuclear LGPT, as well as the analysis of related
experimental signals.

##  I introduction

Exploring the phase structure of strongly interacting matter is a crucial
topic in high energy nuclear physics. One of the tasks is to investigate the
types of phase transition from quark-gluon plasma to hadronic matter. Studies
based on the hadron resonance gas (HRG) model and lattice QCD show that the
phase transition at high temperature and low chemical potential is a smooth
crossover Aoki06 ; Bazavov19 ; Borsanyi13 ; Bazavov14 ; Bazavov17 ; Borsanyi14
; Borsanyi20 . However, the application of lattice QCD in the high-density
region is limited due to the sign problem associated with the fermion
determinant. Nevertheless, most studies with the effective quark models (e.g.,
Fukushima04 ; Ratti06 ; Costa10 ; Fu08 ; Sasaki12 ; Ferreira14 ; Shao2018 ;
Schaefer10 ; Skokov11 ; Liu2018 ; Chen20 ; Zhao23 ), the Dyson-Schwinger
equation approach Qin11 ; Gao16 ; Gao20 ; Fischer14 ; Shi14 , the functional
renormalization group theory Fu20 ; Rennecke17 ; Fu21 , suggest that the
chiral phase transition in high chemical potential region undergoes a first-
order phase transition.

It is proposed that the phase structure of QCD can be studied through the
fluctuations and correlations of conserved charges Stephanov . The net proton
(proxy of net baryon) cumulants have been measured at the Relativistic Heavy
Ion Collider (RHIC) Aggarwal10 ; Adamczyk14 , and a nonmonotonic behavior in
the fourth-order fluctuations as a function of collision energy was discovered
Luo2014 ; Luo2017 , which has sparked extensive discussions about the
existence of the QCD critical endpoint (CEP). Recent publication from STAR
collaboration further investigated the energy dependence of net proton number
fluctuations up to sixth order, showing not only a nonmonotonic behavior but
also more strongly hyperorder density fluctuations at lower collision energies
Aboona23 . Particularly, the fluctuations of net proton number in collisions
at center-of mass energy $\sqrt{s_{NN}}=3\,$GeV is markedly different from
those at $7.7\,$GeV and above. The study further suggests that in collisions
at $\sqrt{s_{NN}}=3\,$GeV, the fluctuation distribution of net proton number
is primarily dominated by the interaction among hadrons Aboona23 .

The new experimental data from RHIC have drawn further attention to the
critical behavior of QCD phase transition and motivated the study of density
fluctuations induced by the interaction among hadrons Huang23 ; ShaoJ23 ;
Marczenko23 . It is well-known that the interaction in nuclear matter can lead
to the LGPT in the low-temperature region Chomaz04 ; Pochodzalla95 ;
Borderie01 ; Botvina95 ; Agostino99 ; Srivastava02 ; Elliott02 . The critical
point of the nuclear LGPT is estimated to be at a temperature of approximately
15 MeV and a chemical potential of around 923 MeV, although the values
slightly depend on the nuclear models. By far various properties of nuclear
matter and the liquid-gas phase transition have been extensively studied in
the literature He23 ; Shao20 ; Chomaz04 ; Pochodzalla95 ; Borderie01 ;
Botvina95 ; Agostino99 ; Srivastava02 ; Elliott02 ; Xu2023 ; Ma1999 ; Deng2022
; Hempel13 ; Mukherjee17 ; Xu13 ; Savchuk20 ; Wang20 .

The presence of strong density fluctuations in collisions at
$\sqrt{s_{NN}}=3\,$GeV raises the question of how the haddronic interaction
and LGPT affect the hyperorder density fluctuations at lower-energy regimes.
The authors in Ref. Poberezhnyuk19 ; Shao20 ; Xu2023  calculated the density
fluctuations up to fourth order induced by hadronic interactions and analyzed
the relationship between the lower-order density fluctuations and nuclear
LGPT. In Ref. Vovchenko17 ; Poberezhnyuk21 , a van der Waals model was taken
to study the higher order fluctuation distributions of net baryon number in
both the pure and mixed phase. In Ref. Marczenko23 , the properties of the
second-order susceptibility of net baryon number density for positive- and
negative-parity nucleons were recently examined near the chiral and nuclear
liquid-gas phase transitions using the double parity model, in which both the
chiral phase transition and nuclear LGPT are effectively included. However, a
systematic study on the hyperorder density fluctuations in combination with
the recent experimental data is still absent.

Considering the hyperorder density fluctuations dominated by the interactions
among hadrons at lower energy regimes (below $3\,$GeV) and the plan of the
HADES collaboration at GSI Helmholtzzentrum für Schwerionenforschung to
measure higher-order net proton and net charge fluctuations in central Au + Au
reactions at collision energies ranging from $0.2A$ to $1.0A\,$GeV to probe
the LGPT region Bluhm20 , it is significant to investigate the hyperorder
density fluctuations in nuclear matter. In this study, we will calculate the
hyperorder net baryon distributions up to sixth order using the nonlinear
Walecka model. The density fluctuations in the the stable and metastable
phases of nucelar LGPT, as well as the area far from the critical region will
be all explored. We will discuss the relations between the nucleon-nucleon
interaction, the nucelar LGPT and the organization structure of the hyperorder
density fluctuations of net baryon number. This investigation has a
consequence for indicating the chiral phase transition, nuclear LGPT, as well
as the analysis of related experimental signals.

The paper is organized as follows. In Sec. II, we introduce the formulas to
describe fluctuations of conserved charges and the nonlinear Walecka model in
the mean-field approximation. In Sec. III, we illustrate the numerical results
of net-baryon number fluctuations, and discuss their relations with nucleon-
nucleon interactions and the LGPT in nuclear matter. A summary is finally
given in Sec. IV.

##  II  Fluctuations of conserved charges and the nonlinear Walecka model

For a thermodynamic system, the fluctuations of conserved charges are
sensitive observables for a phase transition, in particular, for a critical
phenomenon. The pressure of a system in the grand-canonical ensemble is
related to the logarithm of the partition function Karsch15 :

$\frac{P}{T^{4}}=\frac{1}{VT^{3}}\ln[Z(V,T,\mu_{B},\mu_{Q},\mu_{S})],$ (1)

where $\mu_{B},\mu_{Q},\mu_{S}$ are the chemical potentials of conserved
charges, the baryon number, electric charge and strangeness in strong
interaction, respectively. The generalized susceptibilities can be derived by
taking the partial derivatives of the pressure with respect to the
corresponding chemical potentials Luo2017

$\chi^{BQS}_{ijk}=\frac{\partial^{i+j+k}[P/T^{4}]}{\partial(\mu_{B}/T)^{i}\partial(\mu_{Q}/T)^{j}\partial(\mu_{S}/T)^{k}}.$
(2)

The cumulants of multiplicity distributions of the conserved charges are
connected with the generalized susceptibilities by

$\\!C^{BQS}_{ijk}\\!=\\!\frac{\partial^{i\\!+\\!j\\!+\,k}\ln[Z(V,T,\mu_{B},\mu_{Q},\mu_{S})]}{\partial(\mu_{B}/T)^{i}\partial(\mu_{Q}/T)^{j}\partial(\mu_{S}/T)^{k}}\\!=\\!V\\!T^{3}\chi^{BQS}_{ijk}\,.$
(3)

To eliminate the volume dependence in heavy-ion collision experiments,
observables are usually constructed by the ratio of cumulants, and then can be
compared with theoretical calculations of the generalized susceptibilities
with

$\frac{C^{BQS}_{ijk}}{C^{BQS}_{lmn}}=\frac{\chi^{BQS}_{ijk}}{\chi^{BQS}_{lmn}}.$
(4)

The nonlinear Walecka model is taken to calculate fluctuations of net baryon
number in nuclear matter. This model is generally used to describe the
properties of finite nuclei and the equation of state of nuclear matter. The
study in Ref. Fukushima15  also indicates the approximate equivalence of this
model to the hadron resonance gas model at low temperature and small density.
The Lagrangian density for nucleons-meson system in the nonlinear Walecka
model Glendenning97  is

$\displaystyle\cal{L}\\!$ $\displaystyle\\!=\\!$
$\displaystyle\sum_{N}\bar{\psi}_{N}\\!\big{[}i\gamma_{\mu}\partial^{\mu}\\!-\\!(\\!m_{N}\\!-\\!g_{\sigma}\sigma)\\!\\!-\\!g_{\omega}\gamma_{\mu}\omega^{\mu}\\!-\\!g_{\rho}\\!\gamma_{\mu}\boldsymbol{\tau}\\!\cdot\\!\boldsymbol{\rho}^{\mu}\big{]}\\!\psi_{N}$
(5) $\displaystyle+\frac{1}{2}\left(\partial_{\mu}\sigma\partial^{\mu}\sigma-
m_{\sigma}^{2}\sigma^{2}\right)\\!-\\!\frac{1}{3}bm_{N}\,(g_{\sigma}\sigma)^{3}-\frac{1}{4}c\,(g_{\sigma}\sigma)^{4}$
$\displaystyle+\frac{1}{2}m^{2}_{\omega}\omega_{\mu}\omega^{\mu}-\frac{1}{4}\omega_{\mu\nu}\omega^{\mu\nu}$
$\displaystyle+\frac{1}{2}m^{2}_{\rho}\boldsymbol{\rho}_{\mu}\\!\cdot\\!\boldsymbol{\rho}^{\mu}\\!-\\!\frac{1}{4}\boldsymbol{\rho}_{\mu\nu}\\!\cdot\\!\boldsymbol{\rho}^{\mu\nu},$

where $\omega_{\mu\nu}=\partial_{\mu}\omega_{\nu}-\partial_{\nu}\omega_{\mu}$,
$\rho_{\mu\nu}=\partial_{\mu}\boldsymbol{\rho}_{\nu}-\partial_{\nu}\boldsymbol{\rho}_{\mu}$.
The interactions between baryons are mediated by $\sigma,\,\omega,\,\rho$
mesons.

The thermodynamical potential of the nucleons-meson system can be derived in
the mean-field approximation as

$\displaystyle\Omega\\!$ $\displaystyle\\!=$
$\displaystyle\\!-\\!\beta^{-1}\sum_{N}2\\!\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}\\!\bigg{[}\ln\\!\left(1+e^{-\beta\left(E_{N}^{*}(k)-\mu_{N}^{*}\\!\right)}\right)\\!$
(6)
$\displaystyle+\\!\ln\\!\left(1\\!+\\!e^{-\beta\left(E_{N}^{*}(k)+\mu_{N}^{*}\right)}\right)\\!\bigg{]}\\!+\\!\frac{1}{2}m_{\sigma}^{2}\sigma^{2}\\!+\\!\frac{1}{3}bm_{N}\left(g_{\sigma}\sigma\right)^{3}$
$\displaystyle+\frac{1}{4}c\left(g_{\sigma}\sigma\right)^{4}-\frac{1}{2}m_{\omega}^{2}\omega^{2}-\frac{1}{2}m_{\rho}^{2}\rho_{3}^{2},$

where $\beta=1/T$, $E_{N}^{*}=\sqrt{k^{2}+m_{N}^{*2}}$, and $\rho_{3}$ is the
third component of $\rho$ meson field. The effective chemical potential
$\mu_{N}^{*}$ is defined as
$\mu_{N}^{*}=\mu_{N}-g_{\omega}\omega-\tau_{3N}g_{\rho}\rho_{3}$
($\tau_{3N}=1/2$ for proton, $-1/2$ for neutron).

By minimizing the thermodynamical potential

$\frac{\partial\Omega}{\partial\sigma}=\frac{\partial\Omega}{\partial\omega}=\frac{\partial\Omega}{\partial\rho_{3}}=0,$
(7)

the meson field equations can be derived as

$g_{\sigma}\sigma\\!=\\!\left(\frac{g_{\sigma}}{m_{\sigma}}\right)^{2}\\!\left[\rho_{p}^{s}+\rho_{n}^{s}-bm_{N}\left(g_{\sigma}\sigma\right)^{2}-c\left(g_{\sigma}\sigma\right)^{3}\right],$
(8)
$g_{\omega}\omega=\left(\frac{g_{\omega}}{m_{\omega}}\right)^{2}\left(\rho_{p}+\rho_{n}\right),$
(9)
$g_{\rho}\rho_{3}=\frac{1}{2}\left(\frac{g_{\rho}}{m_{\rho}}\right)^{2}\left(\rho_{p}-\rho_{n}\right).$
(10)

In Eqs.(8)-(10), the nucleon number density

$\rho_{i}=2\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}[f\left(E_{i}^{*}-\mu_{i}^{*}\right)-\bar{f}\left(E_{i}^{*}+\mu_{i}^{*}\right)],$
(11)

and the scalar density

$\rho_{i}^{s}=2\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}\frac{m_{i}^{*}}{E_{i}^{*}}[f\left(E_{i}^{*}-\mu_{i}^{*}\right)+\bar{f}\left(E_{i}^{*}+\mu_{i}^{*}\right)],$
(12)

where $f(E_{i}^{*}-\mu_{i}^{*})$ and $\bar{f}(E_{i}^{*}+\mu_{i}^{*})$ are the
fermion and antifermion distribution functions.

For a given temperature and chemical potential (or baryon number density), the
meson field equations can be solved. In the present study, the symmetric
nuclear matter is considered to give an essential description of the density
fluctuations. The model parameters, $g_{\sigma},g_{\omega},g_{\rho},b$ and
$c$, are fitted with the compression modulus $K=240\,$MeV, the symmetric
energy $a_{sym}=31.3\,$MeV, the effective nucleon mass
$m^{*}_{N}=m_{N}-g_{\sigma}\sigma=0.75m_{N}$ ($m_{N}$ is the nucleon mass in
vacuum) and the binding energy $B/A=-16.0\,$MeV at nuclear saturation density
with $\rho_{0}=0.16\,fm^{-3}$.

##  III numerical results and discussions

We first plot in Fig. 1 the phase diagram of nuclear matter in the
temperature-chemical potential plane. The black solid line is the first-order
liquid-gas phase transition line with a critical endpoint at $T_{c}=15.9\,$MeV
and $\mu_{c}=909.8\,$MeV. The black dashed curve, marked with “Line A ”, is
derived with $\partial\sigma/\partial\mu_{B}$ taking the maximum value for a
given temperature, which is connected to the nuclear liquid-gas phase
transition. This plot is somewhat similar to the chiral crossover
transformation line and first-order phase transition in quark model. They both
correspond to the fastest change or jump of thermodynamic order parameter
($\sigma$ field in Walecka model and chiral condensate in quark model) related
to fermion mass. And the subsequent numerical results indeed indicate a
similar structure of density fluctuations to QCD phase transition. The
fluctuation distributions near the “line A” above the critical temperature are
most dramatic. These characteristics reveal essentially similar physical
properties of the two phase transitions because they belong to the same
universal class.

Figure 1:  Phase diagram of LGPT in nuclear matter, extracted chemical freeze-
out temperature versus baryon chemical potential from STAR experiments
Adamczyk17  and SIS data Cleymans99 ; Becattini01 ; Averbeck03 , as well as
the the fitted chemical freeze-out curve from Ref. Cleymans06 . “Line A ” is
derived with $\partial\sigma/\partial\mu_{B}$ taking the maximum value for a
given temperature. Figure 2:  Phase diagram of LGPT and the chemical freeze-
out line in the $T-\rho_{B}$ plane derived in the nonlinear Walecka model on
the basis of the chemical freeze-out line shown in Fig. 1 with the fitted
formula in Eq.(13).

In Fig. 1, we also plot the fitted chemical freeze-out curve and the data from
heavy-ion collision experiments Cleymans06  for the convenience of later
discussion of experimental signals. The fitted chemical freeze-out line can be
described with

$T\left(\mu_{B}\right)=a-b\mu_{B}^{2}-c\mu_{B}^{4},$ (13)

where $a=0.166\,\mathrm{GeV},b=0.139\,\mathrm{GeV}^{-1}$ and
$c=0.053\,\mathrm{GeV}^{-3}$.

We plot in Fig. 2 the phase diagram of liquid-gas transition in the
temperature-density plane. The spinodal line is given with a red dashed curve.
The area under the spinodal line is the unstable phase. The region between the
spinodal line and the first-order phase transition line is the metastable
phase. We also plot the chemical freeze-out curve, which will be referenced to
discuss the density fluctuations in the $T-\rho_{B}$ plane. The freeze-out
line in Fig. 2 is translated from the curve in Fig. 1 based on the calculation
in the nonlinear Walecka model. Using the fitted relation of $T$ and $\mu_{B}$
at freeze-out given by Eq.(13), we can calculate the baryon density by solving
Eqs.(8)-(10). A distinct feature is that the baryon density on the chemical
freeze-out line is nonmonotonic with the decrease of collision energy (or
temperature). This behavior has been confirmed by the experimental data in
combination with hadron resonance gas model Randrup06 ; Poberezhnyuk19 . A
similar result was also found in the $\sigma-\omega$ model in Ref. Fukushima15
.

Besides, Fig. 2 shows that the freeze-out baryon number density is close to
zero near $T=0$. This is consistent with the experimental data mentioned above
as well as the theoretic prediction, because the chemical freeze-out at low
temperature should occur in the gas phase of the first-order LGPT of nuclear
matter. From Fig. 2, it can be seen that the baryon density in the gas phase
is much smaller than that in the liquid phase. As a matter of fact, the baryon
number densities at chemical freeze-out derived from experiments with
$\sqrt{s_{NN}}=200-3\,$GeV are all smaller than the nuclear saturation density
$\rho_{0}$ Randrup06 .

We demonstrate in Fig. 3 the ratios of susceptibilities,
$\chi_{3}^{B}/\chi_{1}^{B}$, $\chi_{4}^{B}/\chi_{2}^{B}$,
$\chi_{5}^{B}/\chi_{1}^{B}$ and $\chi_{6}^{B}/\chi_{2}^{B}$, as functions of
baryon chemical potential for $T=75,50,25\,$MeV. The numerical results show
that all these ratios of susceptibilities change from one to zero as the
chemical potential increases. Such a behavior is analogous to the density
fluctuations of chiral phase transition from a chiral broken phase to a chiral
restored phase Shao2018 . The dynamical nucleon mass plays a similar role for
nuclear matter. The strong density fluctuation for each temperature are
roughly located in the region where $\partial\sigma/\partial\mu_{B}$ changes
rapidly. The collection of maxima of $\partial\sigma/\partial\mu_{B}$ at
different temperatures above the critical endpoint corresponds to the dashed
line A plotted in Fig. 1.

Figure 3:  Ratios of net baryon number susceptibilities as functions of
chemical potential for different temperatures. Panel (a):
$\chi_{3}^{B}/\chi_{1}^{B}$; Panel (b): $\chi_{4}^{B}/\chi_{2}^{B}$(kurtosis);
Panel (c): $\chi_{5}^{B}/\chi_{1}^{B}$(hyperskewness); Panel (d):
$\chi_{6}^{B}/\chi_{2}^{B}$(hyperkurtosis). The solid dots demonstrate the
density fluctuations at chemical freeze-out given in Fig. 1.

Figure 4: Ratios of net baryon number susceptibilities as functions of baryon
density for different temperatures. Panel (a): $\chi_{3}^{B}/\chi_{1}^{B}$;
Panel (b): $\chi_{4}^{B}/\chi_{2}^{B}$(kurtosis); Panel (c):
$\chi_{5}^{B}/\chi_{1}^{B}$(hyperskewness); Panel (d):
$\chi_{6}^{B}/\chi_{2}^{B}$(hyperkurtosis). The solid dots demonstrate the
density fluctuations at chemical freeze-out given in Fig. 1.

Fig. 3 also shows that the amplitude of density fluctuations is relatively
small at $T=75\,$MeV. As the temperature decreases and approaches the critical
region, the density fluctuations gradually increase and diverge at the CEP.
This indicates a close relationship between the density fluctuations of
conserved charges and the phase transition of nuclear matter. From Fig. 3, we
can see that the hyperorder density fluctuations (hyperskewness,
$\chi_{5}^{B}/\chi_{1}^{B}$ and hyperkurtosis, $\chi_{6}^{B}/\chi_{2}^{B}$)
are dramatic than those of lower orders ($\chi_{3}^{B}/\chi_{1}^{B}$ and
$\chi_{4}^{B}/\chi_{2}^{B}$), meaning that the higher-order fluctuations are
more sensitive to observe in experiments. Similar behaviors aroused by the
chiral phase transition have been found in the beam energy scan experiments at
RHIC STAR Aboona23 . Additionally, we can observe from Fig. 3 that
$\chi_{3}^{B}/\chi_{1}^{B}$ and $\chi_{4}^{B}/\chi_{2}^{B}$ have one maximum
and one minimum, while $\chi_{5}^{B}/\chi_{1}^{B}$
($\chi_{6}^{B}/\chi_{2}^{B}$) has two maximums and one minimum (two minimums).

The solid dots on the curves in Fig. 3 demonstrate the fluctuation
distributions at the points of chemical freeze-out, corresponding to the
chemical freeze-out line plotted in Fig. 1. It can be observed that
$\chi_{3}^{B}/\chi_{1}^{B}$ and $\chi_{4}^{B}/\chi_{2}^{B}$ are all positive
for $T=75,50$ and $25\,$MeV at chemical freeze-out. However,
$\chi_{5}^{B}/\chi_{1}^{B}$ and $\chi_{6}^{B}/\chi_{2}^{B}$ are negative at
$T=75$ and $50\,$MeV, and become positive at $T=25\,$MeV. Note that the
negative $\chi_{5}^{B}/\chi_{1}^{B}$ and $\chi_{6}^{B}/\chi_{2}^{B}$ also
appear at chemical freeze-out at higher temperatures for QCD phase transition
Aboona23 . Interestingly, it is found that the ratios of net proton number
distributions, $C_{4}/C_{2}$, $C_{5}/C_{1}$ and $C_{6}/C_{2}$ at
$\sqrt{s_{NN}}=3\,$GeV with the collision centrality $0\%-40\%$, are very
close to the numerical results of $\chi_{4}^{B}/\chi_{2}^{B}$,
$\chi_{5}^{B}/\chi_{1}^{B}$ and $\chi_{6}^{B}/\chi_{2}^{B}$ at $T=28\,$MeV and
$\mu_{B}=864\,$MeV in the nonlinear Walecka model. However, the physical
conditions, $T=28\,$MeV and $\mu_{B}=864\,$MeV, are different from the
chemical freeze-out condition derived in collisions at a center-of-mass energy
of $3\,$GeV. For the later case, $T\approx 70\,$MeV and $\mu_{B}\approx
750\,$MeV are extracted at the chemical freeze-out. We also notice that there
is a different argument based on the data of only the net proton cumulants
$C_{4}/C_{2}$ with collision centralities $0\%-5\%$ Xu2023 . Higher accuracy
data of hyperskewness and hyperkurtosis at lower collision energies are needed
to distinguish between the two scenarios.

We further plot the ratios of net baryon number susceptibilities in the
temperature-density plane in Fig. 4. Compared with the results in Fig. 3, the
distributions of density fluctuations appear to be broadened in the
$T-\rho_{B}$ panel. The solid dots in Fig. 4 indicate that the baryon density
at chemical freeze-out decreases as the temperature descends from $75$ to
$25\,$MeV. In contrast, Fig. 3 shows that the chemical potential at freeze-out
increases as the temperature decreases from $75$ to $25\,$MeV. These results
are consistent with the chemical freeze-out curve given in Figs. 1 and 2,
since the variation of baryon density along the chemical freeze-out line is
nonmonotonic in the temperature-density phase diagram.

Figure 5: (color online) Contour plots of the ratios of net baryon number
susceptibilities, $\chi_{3}^{B}/\chi_{1}^{B}$, $\chi_{4}^{B}/\chi_{2}^{B}$,
$\chi_{5}^{B}/\chi_{1}^{B}$ and $\chi_{6}^{B}/\chi_{2}^{B}$, in the
$T-\mu_{B}$ phase diagram.

Figure 6: (color online) Contour plots of the ratios of net baryon number
susceptibilities, $\chi_{3}^{B}/\chi_{1}^{B}$, $\chi_{4}^{B}/\chi_{2}^{B}$,
$\chi_{5}^{B}/\chi_{1}^{B}$ and $\chi_{6}^{B}/\chi_{2}^{B}$, in the
$T-\rho_{B}$ phase diagram.

In order to visualize the organization structure of density fluctuations in
the phase diagram, we present the contour plots of net baryon number
fluctuations in Figs. 5 and 6. The two figures explicitly demonstrate that the
ratios of hyperorder susceptibilities exhibit more pronounced density
fluctuations than those of lower-order ones. At the meanwhile, we can see that
the closer the density fluctuations are to the critical endpoint of the phase
transition, the more frequently the oscillatory behavior of the density
fluctuations appears. Additionally, in Fig. 6, we also present the net baryon
number density fluctuations in the metastable phase. It can be observed that
the ratios of net-baryon susceptibilities change continuously from the stable
phase to the metastable phase. On the boundary of first-order liquid-gas phase
transition (excluding the critical point), the density fluctuations have
finite values, while the density fluctuations tend to diverge on the spinodal
line (red line).

These behaviors of net baryon number fluctuations can be simply understood
according to the relationship between $P/T^{4}$ and $\mu_{B}/T$. For each
given temperature above the critical endpoint, $P/T^{4}$ increase monotonously
as a function of $\mu_{B}/T$. For a lower temperature, the slop of $P/T^{4}$
increases faster with $\mu_{B}/T$ near the region where
$\partial{\sigma}/\partial{\mu_{B}}$ changes rapidly, and the maximum slope is
also larger. Such a feature results in that the higher order susceptibility
oscillates much sharply in a relatively narrow range at lower temperature. The
center of oscillation moves roughly along Line A (shown in Fig. 1) towards the
CEP of liquid-gas phase transition.

To improve the theoretical calculation in the future, more degrees of freedom
of hadrons and interactions should be included at relatively higher
temperature. On the experimental side, we should keep in mind that the initial
volume fluctuation effects in experiments become significant due to lower
charged particle multiplicity at lower-energy collisions. Nevertheless, this
study qualitatively predicts the behavior of net baryon number fluctuations at
chemical freeze-out at low temperatures. With the increase of statistical data
and the improvement of data precision in the second phase of beam energy scan
at RHIC STAR, the deep combination of theory and experiment will bring
important opportunities for exploring the phase structure of strongly
interacting matter.

##  IV summary and conclusion

Fluctuations of conserved charges are sensitive probes to explore the phase
transition of strongly interacting matter. In this research, using the non-
linear Walecka model we investigated the fluctuations of net baryon number up
to sixth order caused by hadronic interactions in nuclear matter, and explored
the relationship between hyperorder baryon number fluctuations and nuclear
liquid-gas phase transition, including the stable and metastable phases as
well as the region far from the phase transition.

The calculation indicates that the fluctuations of net baryon number gradually
increase from the high-temperature region to critical region of nuclear LGPT.
In particular, the fluctuations are prominent near the region where the
$\sigma$ field or the nucleon mass changes rapidly for each temperature. It
exhibits similar behavior to the chiral phase transition of quark matter. This
can be further attributed to that the two phase transitions have the same
universal class.

Compared with the kurtosis and skewness of net baryon number fluctuations, the
values of hyperkurtosis ($\chi_{6}^{B}/\chi_{2}^{B}$) and hyperskewness
($\chi_{5}^{B}/\chi_{1}^{B}$) are more dramatic, which provide more sensitive
signals to detect the phase transformation. In combination with heavy-ion
collision experiments, we further extracted density fluctuations at chemical
freeze-out, and obtain the variation trend of baryon number fluctuations with
the decrease of temperature (collision energy). With the release of more
precise data in BES II program in the future, the results obtained in this
research can be referred to analyze the QCD phase transition and nuclear
liquid-gas phase transition. In addition, more physical mechanisms, such as
the initial volume fluctuation effects, the non-equilibrium evolution and more
hadronic degrees of freedom with complex interactions, need to be considered
to better understand the evolution of heavy-ion collision experiments at
relatively lower collision energy.

## Acknowledgements

This work is supported by the National Natural Science Foundation of China
under Grant No. 11875213.

## References

## References

  * (1) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, Nature 443, 675 (2006).
  * (2) A. Bazavov H.-T. Ding b, P. Hegde, O. Kaczmarek, F. Karsch, N. Karthiks, et al., (hotQCD Collaboration) Phys. Lett. B 795, 15 (2019).
  * (3) S. Borsányi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabó, Phys. Rev. Lett. 111, 062005 (2013).
  * (4) A. Bazavov, et al., (hotQCD Collaboration), Phys. Rev. D. 90, 094503 (2014).
  * (5) A. Bazavov, et al., (hotQCD Collaboration), Phys. Rev. D 96, 074510 (2017).
  * (6) S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Sabzó, Phys. Lett. B 730, 99 (2014).
  * (7) S. Borsányi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 125, 052001 (2020).
  * (8) K. Fukushima, Phys. Lett. B 591, 277 (2004).
  * (9) C. Ratti, M. A. Thaler, and W. Weise, Phys. Rev. D 73, 014019 (2006).
  * (10) P. Costa, M. C. Ruivo, C. A. de Sousa, and H. Hansen, Symmetry 2, 1338 (2010) .
  * (11) W. J. Fu, Z. Zhang, and Y. X. Liu, Phys. Rev. D 77, 014006 (2008).
  * (12) T. Sasaki, J. Takahashi, Y. Sakai, H. Kouno, and M. Yahiro, Phys. Rev. D 85, 056009 (2012) .
  * (13) M. Ferreira, P. Costa, and C. Providência, Phys. Rev. D 89, 036006 (2014).
  * (14) Y. Q. Zhao, S. He, D. F. Hou, L. Li, and Z. B. Li, JHEP 04, 115 (2023).
  * (15) G. Y. Shao, Z. D. Tang, X. Y. Gao, and W. B. He, Eur. Phys. J. C 78, 138 (2018).
  * (16) B. J. Schaefer, M. Wagner, and J. Wambach, Phys. Rev. D 81, 074013 (2010) .
  * (17) V. Skokov, B. Friman, and K. Redlich, Phys. Rev. C 83, 054904 (2011).
  * (18) L. M. Liu, J. Xu, and G. X. Peng, Phys. Rev. D 104, 076009 (2021).
  * (19) X. Chen, D. N. Li, D. F. Hou, and M. Huang, Journal of High Energy Physics 2020, 73 (2020).
  * (20) Y. D. Chen, D. N. Li, and M. Huang, Phys.Rev.D 106, 106002 (2022).
  * (21) S. X. Qin, L. Chang, H. Chen, Y. X. Liu, and C. D. Roberts, Phys. Rev. Lett. 106, 172301 (2011).
  * (22) F. Gao, J. Chen, Y. X. Liu, S. X. Qin, C. D. Roberts, and S. M. Schmidt, Phys. Rev. D 93, 094019 (2016).
  * (23) F. Gao and J. M. Pawlowski, Phys. Rev. D 102, 034027 (2020).
  * (24) C. S. Fischer, J. Luecker, and C. A. Welzbacher, Phys. Rev. D 90, 034022 (2014) .
  * (25) C. Shi, Y. L. Wang, Y. Jiang, Z. F. Cui, and H. S. Zong, JHEP 1407, 014 (2014).
  * (26) W. J. Fu, J. M. Pawlowski, F. Rennecke, Phys. Rev. D 101, 054032 (2020).
  * (27) F. Rennecke and B. J. Schaefer, Phys. Rev. D 96, 016009 (2017).
  * (28) W. J. Fu, X. F. Luo, J. M. Pawlowski, F. Rennecke, R. Wen, and S. Yin, Phys. Rev. D 104, 094047 (2021)
  * (29) M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009); Phys. Rev. Lett. 107, 052301 (2011).
  * (30) M. M. Aggarwal, et al., (STAR Collaboration), Phys. Rev. Lett. 105, 022302 (2010).
  * (31) L. Adamczyk, et al., (STAR Collaboration), Phys. Rev. Lett. 112, 032302 (2014).
  * (32) X. Luo (STAR Collaboration), Proc. Sci. CPOD2014 2015, 019.
  * (33) X. Luo and N. Xu, Nucl. Sci. Tech. 28, 112 (2017).
  * (34) B. E. Aboona, J. Adam, L. Adamczyk, J. R. Adams, I. Aggarwal, and M. M. Aggarwal, et al., (STAR Collaboration), 130, 082301 (2023).
  * (35) C. Huang, Y. Y. Tan, R. Wen, S. Yin, and W. J. Fu, arXiv:2303.10869v1.
  * (36) J. X. Shao, W. J. Fu, and Y. X. Liu, arXiv:2305.16763v2.
  * (37) M. Marczenko, K. Redlich, and C. Sasaki, Phys. Rev. D 107, 054046 (2023).
  * (38) P. Chomaz, M. Colonna, and J. Randrup, Phys. Rep. 389, 263 (2004).
  * (39) J. Pochodzalla, T. Mohlenkamp, T. Rubehn, A. Schuttauf, A. Worner, E. Zude, et al., (ALADIN Collaboration), Phys. Rev. Lett. 75, 1040 (1995).
  * (40) B. Borderie, G. Tabacaru, P. Chomaz, M. Colonna, A. Guarnera, M. Parlog, et al., (NDRA Collaboration), Phys. Rev. Lett. 86, 3252 (2001).
  * (41) A. S. Botvina, I.N. Mishustin, M. Begemann-Blaich, Hubele, G. Imme, I. Iori, et al., Nucl. Phys. A584, 737 (1995).
  * (42) M. D’Agostino, A.S. Botvina, l, M. Bruno, A. Bonasera, J. R Bondorf, R. Bougault, et al., Nucl. Phys. A650, 329 (1999).
  * (43) B. K. Srivastava, R. P. Scharenberg, S. Albergo, F. Bieser, F. P. Brady, Z. Caccia, et al., (EOS Collaboration), Phys. Rev. C 65, 054617 (2002).
  * (44) J. B. Elliott, et al., (ISiS Collaboration), Phys. Rev. Lett. 88, 042701 (2002).
  * (45) W. B. He, G. Y. Shao, C. L. Xie, Phys. Rev. C 107, 014903 (2023).
  * (46) G. Y. Shao, X. Y. Gao, and W. B. He, Phys. Rev. D 101, 074029 (2020).
  * (47) K. Xu, M. Huang, arXiv:2307.12600v1.
  * (48) Y. G. Ma, Phys. Rev. Lett 83, 3617 (1999).
  * (49) X. G. Deng, P. Danielewicz, Y. G. Ma, H. Lin, Y. X. Zhang, Phys. Rev. C 105, 064613 (2022).
  * (50) M. Hempel, V. Dexheimer, S. Schramm, and I. Iosilevskiy, Phys. Rev. C 88, 014906 (2013).
  * (51) A. Mukherjee, J. Steinheimer, and S. Schramm, Phys. Rev. C 96, 025205 (2017).
  * (52) J. Xu a, L. W. Chen, C. M. Ko d, B. A. Li, Y. G. Ma, Phys. Lett. B 727, 244 (2013).
  * (53) O. Savchuk, V. Vovchenko, R. V. Poberezhnyuk, M. I. Gorenstein, and H. Stoecker, Phys. Rev. C 101, 035205 (2020).
  * (54) R. Wang, Y. G. Ma, R. Wada, L. W. Chen, W. W. He, H. L. Liu and K. J. Sun, Phys. Rev. Res. 2, 043202 (2020).
  * (55) R. V. Poberezhnyuk, V. Vovchenko, A. Motornenko, M. I. Gorenstein, and H. Stöcker, Phys. Rev. C 100, 054904 (2019).
  * (56) V. Vovchenko, M. I. Gorenstein, and H. Stöecker, Phys. Rev. Lett 118, 182301 (2017).
  * (57) R. V. Poberezhnyuk, O. Savchuk, M. I. Gorenstein, V. Vovchenko, and H. Stöecker, Phys. Rev. C 103, 024912 (2021).
  * (58) M. Bluhm, A. Kalweit, M. Nahrgang, M. Arslandok, P. Braun-Munzinger, S. Floerchinger, Nucl. Phys. A1003 (2020) 122016.
  * (59) H. T. Ding, F. Karsch, and S. Mukherjee. Int. J. Mod. Phys. E 24, 530007 (2015).
  * (60) K. Fukushima, Phys. Rev. C 91, 044910 (2015).
  * (61) N. K. Glendenning, COMPACT STARS, Nuclear Physics, Particle Physics, and General Relativity (Springer-Verlag, New York, 1997).
  * (62) J. Cleymans, H. Oeschler, K. Redlich, and S. Wheaton, Phys. Rev. C 73, 034905 (2006).
  * (63) L. Adamczyk, et al., (STAR Collaboration), Phys. Rev. C 96, 044904 (2017).
  * (64) J. Cleymans, H. Oeschler, and K. Redlich, Phys. Rev. C 59, 1663 (1999).
  * (65) F. Becattini, J. Cleymans, A. Keränen, E. Suhonen, and K. Redlich, Phys. Rev. C 64 024901 (2001).
  * (66) R. Averbeck, R. Holzmann, V. Metag, and R. S. Simon, Phys. Rev. C 67, 024903 (2003).
  * (67) J. Randrup and J. Cleymans, Phys. Rev. C 74, 047901 (2006).