Generalized Lennard-Jones Potentials, SUSYQM and Differential Galois Theory

Generalized Lennard-Jones Potentials,
SUSYQM and Differential Galois Theory

Manuel F. ACOSTA-HUMÁNEZ ${}^{{\dagger}^{1}}$, Primitivo B. ACOSTA-HUMÁNEZ
${}^{{\dagger}^{2}{\dagger}^{3}}$
and Erick TUIRÁN ${}^{{\dagger}^{4}}$

M.F. Acosta-Humánez, P.B. Acosta-Humánez and E. Tuirán

${}^{{\dagger}^{1}}$ Departamento de Física, Universidad Nacional de Colombia,
${}^{{\dagger}^{1}}$ Sede Bogotá, Ciudad Universitaria 111321, Bogotá,
Colombia<EMAIL_ADDRESS>
${}^{{\dagger}^{2}}$ Facultad de Ciencias Básicas y Biomédicas, Universidad
Simón Bolívar,
${}^{{\dagger}^{2}}$ Sede 3, Carrera 59 No. 58–135. Barranquilla, Colombia
<EMAIL_ADDRESS>http://www.intelectual.co/primi/

${}^{{\dagger}^{3}}$ Instituto Superior de Formación Docente Salomé Ureña -
ISFODOSU,
${}^{{\dagger}^{3}}$ Recinto Emilio Prud’Homme, Calle R. C. Tolentino #51,
esquina 16 de Agosto,
${}^{{\dagger}^{3}}$ Los Pepines, Santiago de los Caballeros, República
Dominicana

${}^{{\dagger}^{4}}$ Departamento de Física y Geociencias, Universidad del
Norte,
${}^{{\dagger}^{4}}$ Km 5 Vía a Puerto Colombia AA 1569, Barranquilla,
Colombia<EMAIL_ADDRESS>
Received May 01, 2018, in final form September 14, 2018; Published online
September 19, 2018

In this paper we start with proving that the Schrödinger equation (SE) with
the classical $12-6$ Lennard-Jones (L-J) potential is nonintegrable in the
sense of the differential Galois theory (DGT), for any value of energy; i.e.,
there are no solutions in closed form for such differential equation. We study
the $10-6$ potential through DGT and SUSYQM; being it one of the two partner
potentials built with a superpotential of the form $w(r)\propto 1/r^{5}$. We
also find that it is integrable in the sense of DGT for zero energy. A first
analysis of the applicability and physical consequences of the model is
carried out in terms of the so called De Boer principle of corresponding
states. A comparison of the second virial coefficient $B(T)$ for both
potentials shows a good agreement for low temperatures. As a consequence of
these results we propose the $10-6$ potential as an integrable alternative to
be applied in further studies instead of the original $12-6$ L-J potential.
Finally we study through DGT and SUSYQM the integrability of the SE with a
generalized $(2\nu-2)-\nu$ L-J potential. This analysis do not include the
study of square integrable wave functions, excited states and energies
different than zero for the generalization of L-J potentials.

Lennard-Jones potential; differential Galois theory; SUSYQM; De Boer principle
of corresponding states

12H05; 81V55; 81Q05

##  1 Introduction

The Lennard-Jones potential (L-J) was proposed in 1931 in order to model the
concurrence between the long-range attraction and the short-range repulsion in
radial interatomic interactions [34]. In a later work, the description of such
potential was employed in order to describe the equation of state of a gas in
terms of its interatomic forces [35], thus concluding and enhancing an
investigation started by Mie in 1903 [38]. The L-J potential is usually used,
at the level of classical statistical mechanics, to study the behavior of
fluid materials, ranging from simple molecules to polymers and proteins [24,
31, 37]. In theoretical quantum chemistry, among many applications, we point
out: its implementation in the theory of molecular orbitals, allowing to
compute the tendency of two electrons in the same space orbital to keep each
other apart because of the repulsive field between them [26]; the numerical
implementations in order to compute the transferable inter-molecular potential
functions (TIPS) in alcohols, ethers and water, that have given an
understanding of the interactions of these chemical compounds in solvents
[27]. Also a mathematical model that has been proposed for calculating the
isosteric heat of adsorption of simple fluids onto flat surfaces. On this
respect, theoretical and experimental results were compared in order to study
the influence of the choice of the intermolecular potential parameters [41].
Finally, a experimental methodology and theoretical calculations applying the
Lennard-Jones potential, for determining micropore-size distributions,
obtained from physical adsorption isotherm data, have provided valuable
microstructural information, which is still widely used today [25, 42, 48].

With the increase of numerical techniques, calculations with explicit
solutions in physical models don’t have in the present the same importance as
in past decades. Nevertheless exact solutions when available, have always
served as elucidating tools for finding general properties of the system,
which otherwise could remain hidden. The main motivation of this paper is the
application of supersymmetric quantum mechanics (SUSYQM) and differential
Galois theory (DGT) to obtain explicit solutions of the Schrödinger equation
(SE) with variants of the Lennard-Jones potential, as well as the set of
eigenvalues associated to each solution.

SUSYQM, introduced by E. Witten in 1981, is the simplest example where
supersymmetry can be dynamically broken [51]. In spite of its initial
character of a toy model; SUSYQM has earned importance in the recent decades,
because it served as a starting point to the development of attractive
theoretical features and concepts like shape invariance, isospectrality and
factorization, that give new perspectives to old problems in quantum
mechanics, like the integrability of the SE, see for example [1, 21, 17] and
the path integral formulation of classical mechanics [22]. On the other hand,
there is plenty of papers in mathematical physics wherein DGT has been
applied; see for example [2, 4, 5, 6, 7, 8] for applications to study the non-
integrability of Hamiltonian systems. For applications in the integrability of
the SE, see [1, 3, 9, 11, 12, 13, 14]. For applications of differential Galois
theory to other quantum integrable systems see [15, 46]. The main Galoisian
tools used in some of these papers are the Hamiltonian algebrization and the
Kovacic’s algorithm. These tools have led and still lead, to deduce exact
solutions in several areas of mathematical physics.

The structure of this paper is as follows. Section 2 is devoted to the
theoretical background necessary to understand the rest of the paper. It
summarizes topics such as the Schrödinger equation for central potentials,
Lennard-Jones potentials $12-6$, $10-6$ and $(2\nu-2)-\nu$, SUSYQM, the De
Boer principle of corresponding states, the virial equation and DGT. In
Section 3 we study the integrability of the SE with the usual $12-6$ Lennard-
Jones potential, as well as the alternative versions $10-6$ and
$(2\nu-2)-\nu$. Our contributions consist in the deduction of algebraic and
physical conditions over the parameters of such SE’s to get their
integrability in the sense of DGT and the superpotentials in the integrable
cases. A first study of physical consequences will also be detailed in this
section. In Section 4 some remarks concerning future works are established.

##  2 Preliminaries

### The Schrödinger equation for a central potential

We are interested in studying a physical model for a many-body system where
the main contribution of the interaction of its constituents is pairwise and
radial in nature. In addition, the physical conditions of the system
(temperature, density, etc) are such, that its quantum behavior is non-
negligible. In this section we set shortly the theoretical background, in
order to establish our physical model with a central potential, and also the
notation to be applied in the rest of the paper [16]. The Hamiltonian for a
system of two spinless particles with masses $m_{1}$ and $m_{2}$ interacting
via a radial potential $V\left(\left|\vec{r}_{1}-\vec{r}_{2}\right|\right)$ is
given by

$\displaystyle
H=T+V=\frac{p_{1}^{2}}{2m_{1}}+\frac{p_{2}^{2}}{2m_{2}}+V(|\vec{r}_{1}-\vec{r}_{2}|).$
(2.1)

It is an usual subject of textbooks in classical mechanics to show that (2.1)
can be separated into two parts, one related to the motion of the center of
mass $\vec{R}$ of the system and the other related to the relative motion of
the particles. The new coordinate system is given by the following
transformation rules

$\displaystyle\vec{R}=\frac{m_{1}\vec{r}_{1}+m_{2}\vec{r}_{2}}{m_{1}+m_{2}},\qquad\vec{r}=\vec{r}_{1}-\vec{r}_{2},\qquad\mu=\frac{m_{1}m_{2}}{m_{1}+m_{2}},\qquad
M=m_{1}+m_{2},$ $\displaystyle\vec{p}_{r}=\mu\frac{{\rm d}\vec{r}}{{\rm
d}t},\qquad\vec{p}_{R}=M\frac{{\rm d}\vec{R}}{{\rm d}t},$

where $M$ is the total mass of the system, $\mu$ is called the reduced mass.
The Hamiltonian in the new coordinates takes the form

$\displaystyle H=\frac{p_{r}^{2}}{2\mu}+\frac{p_{R}^{2}}{2M}+V(r),$ (2.2)

where $p_{r}$ and $p_{R}$ are the canonical momenta conjugated respectively to
the coordinates $r=|\vec{r}_{1}-\vec{r}_{2}|$ and $R=|\vec{R}|$. Since we are
not dealing with external forces, the motion of the center of mass is uniform
rectilinear. For several analysis it is suitable to work in a frame at rest
with the center of mass, which is still an inertial reference frame, in that
case the Hamiltonian (2.2) is reduced to

$\displaystyle H=\frac{p_{r}^{2}}{2\mu}+V(r).$ (2.3)

The Hamiltonian in (2.3) represents the energy of the relative motion of the
two particles; it describes the motion of a fictitious particle, the _relative
particle_ with a mass given by the reduced mass $\mu$ and a position and
momentum given by the relative coordinates $\vec{r}$ and $\vec{p}_{r}$. The
quantum mechanical model of our interest is based on this Hamiltonian. The
usual rules of quantization in the position representation lead to the time-
independent Schrödinger equation for our two-particles system

$\displaystyle\left[-\left(\frac{\hbar^{2}}{2\mu}\right)\vec{\nabla}^{2}+V(r)\right]\Psi(\vec{r})=E\Psi(\vec{r}).$
(2.4)

Since $V(r)$ is a rational central potential, the eigenfunctions
$\Psi(\vec{r})$ are separable into radial and angular parts, the last one
given by the spherical harmonics

$\displaystyle\Psi(\vec{r})=\frac{1}{r}u_{k,l}(r)Y_{l}^{m}(\theta,\varphi).$

The differential equation of our interest corresponds to the radial part of
(2.4) as follows

$\displaystyle\left[-\left(\frac{\hbar^{2}}{2\mu}\right)\frac{{\rm
d}^{2}}{{\rm d}r^{2}}+\frac{l(l+1)}{2\mu
r^{2}}+V(r)\right]u_{k,l}(r)=E_{k,l}u_{k,l}(r),$

where $l$ and $m$ are the usual quantum numbers for angular momentum; $k$
represents the different values of energy for fixed $l$, and it can be either
discrete or continuous. Defining the _effective radial potential_ as $V^{\rm
eff}(r)\equiv l(l+1)/\big{(}2\mu r^{2}\big{)}+V(r)$ and leaving the second
derivative in $r$ on one side, we have

$\displaystyle\left(\frac{2\mu}{\hbar^{2}}\right)\big{[}V^{\rm
eff}(r)-E_{k,l}\big{]}u_{k,l}(r)=\frac{{\rm d}^{2}}{{\rm d}r^{2}}u_{k,l}(r)$
(2.5)

at this point we define a rescaled potential $v^{\rm
eff}(r)\equiv\big{(}\frac{2\mu}{\hbar^{2}}\big{)}V^{\rm eff}(r)\ $and a
similarly rescaled energy
$\big{(}\frac{2\mu}{\hbar^{2}}\big{)}E_{k,l}\equiv\varepsilon_{k,l}$; in this
case equation (2.5) turns out to be

$\displaystyle\big{[}v^{\rm
eff}(r)-\varepsilon_{k,l}\big{]}u_{k,l}(r)=\frac{{\rm d}^{2}}{{\rm
d}r^{2}}u_{k,l}(r).$ (2.6)

In this way it is natural to define a rescaled Hamiltonian as
$\mathcal{H}\equiv\big{[}{-}\frac{{\rm d}^{2}}{{\rm d}r^{2}}+v(r)\big{]}$ in
order to recover (2.6):

$\displaystyle\mathcal{H}_{\rm eff}u_{k,l}(r)\equiv\left[-\frac{{\rm
d}^{2}}{{\rm d}r^{2}}+v^{\rm
eff}(r)\right]u_{k,l}(r)=\varepsilon_{k,l}u_{k,l}(r).$ (2.7)

The case for $l=0$ defines the non-effective potential, and is also of great
interest for our study

$\displaystyle\mathcal{H}u_{k}(r)\equiv\left[-\frac{{\rm d}^{2}}{{\rm
d}r^{2}}+v(r)\right]u_{k}(r)=\varepsilon_{k}u_{k}(r),$ (2.8)

where we have simplified $u_{k,l=0}$ and $\varepsilon_{k,l=0}$ to $u_{k}$ and
$\varepsilon_{k}$, respectively. We observe that (2.8) is a rescaled version
of

$\displaystyle\widetilde{\mathcal{H}}u_{k}(r)\equiv\left[-\left(\frac{\hbar^{2}}{2\mu}\right)\frac{{\rm
d}^{2}}{{\rm d}r^{2}}+V(r)\right]u_{k}(r)=E_{k}u_{k}(r).$ (2.9)

equations (2.7), (2.8) and (2.9) are the subject of our mathematical and
physical analysis in Section 3.

### The $\boldsymbol{12-6}$ Lennard-Jones potential and its generalizations

The 12-6 Lennard-Jones potential is usually presented in terms of two
constants $A$ and $B$

$\displaystyle V_{12-6}(r)=-{\frac{A}{r^{6}}}+{\frac{B}{r^{12}}},$ (2.10)

where the negative term $-A/r^{6}$ leads to van der Waals attractive fields
and comes from the second-order correction in perturbation theory to the
dipole-dipole interaction between two atoms [16]. The positive term $B/r^{12}$
models the short range electronic repulsion between atoms and has no
theoretical justification; it was empirically chosen because it fits
reasonably good data coming from experiments with diatomic gases [34]. An
alternative version is given by

$\displaystyle
V_{12-6}(r)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],\qquad
A=4\epsilon\sigma^{6},\qquad B=4\epsilon\sigma^{12},$ (2.11)

where $\epsilon$ is the atomic depth of the potential well, $\sigma$ is the
finite distance at which the inter-particle potential is zero and $r$ is the
distance between the particles (see Fig. 1). In mathematical terms, $\sigma>0$
and $\epsilon>0$ satisfy that $V(\sigma)=0$ and
$V\big{(}\sqrt[6]{2}\sigma\big{)}=-\epsilon$; this means that $\sigma$ is a
zero potential length and the point
$\big{(}\begin{matrix}\sqrt[6]{2}\sigma,&-\epsilon\end{matrix}\big{)}$ is the
local minimum of the potential in the interval $(0,\infty)$. It can easily be
shown that there is no other critical point in such interval. In physical
terms the well depth $\epsilon$ and the zero potential length $\sigma$ are
parameters that describe the cohesive and repulsive forces that take place in
a gas or liquid at the molecular level. $\epsilon$ measures the strength of
the attraction between pairs of molecules and $\sigma$ is the radius of the
repulsive core when two molecules collide.

Figure 1: $\frac{v(r)}{\epsilon}$ vs. $\frac{r}{\sigma}$ plot for the rescaled
$12-6$ Lennard-Jones potential given in equation (2.13) for $C=0$.

In order to explore with the differential Galois theory the integrability of
the Schrödinger equation with the Lennard-Jones potential (2.10) and other
related cases, we introduce the generalized effective version with arbitrary
powers $\nu$ and $\delta$ given in [34]

$\displaystyle
V_{\delta-\nu}(r)\equiv-{\frac{A}{r^{\nu}}}+{\frac{B}{r^{\delta}}},\qquad
V_{\delta-\nu}^{\rm eff}(r)\equiv
V_{\delta-\nu}(r)+\frac{C}{r^{2}}=-{\frac{A}{r^{\nu}}}+{\frac{B}{r^{\delta}}}+\frac{C}{r^{2}},$

where $0<\nu<\delta$, $A>0$, $B>0$, $C\geq 0$. Its rescaled version is given
by

$\displaystyle
v_{\delta-\nu}(r)\equiv\left(\frac{2\mu}{\hbar^{2}}\right)V_{\delta-\nu}(r)=-{\frac{\bar{A}}{r^{\nu}}}+{\frac{\bar{B}}{r^{\delta}}},$
(2.12) $\displaystyle v_{\delta-\nu}^{\rm
eff}(r)\equiv\left(\frac{2\mu}{\hbar^{2}}\right)V_{\delta-\nu}^{\rm
eff}(r)=-{\frac{\bar{A}}{r^{\nu}}}+{\frac{\bar{B}}{r^{\delta}}}+\frac{\bar{C}}{r^{2}},$
(2.13)
$\displaystyle\bar{A}\equiv\left(\frac{2\mu}{\hbar^{2}}\right)A,\qquad\bar{B}\equiv\left(\frac{2\mu}{\hbar^{2}}\right)B,\qquad\bar{C}\equiv\left(\frac{2\mu}{\hbar^{2}}\right)C.$

The special case for $\delta=2\nu-2$ and some of its analytic advantages has
been studied by J. Pade in [43]

$\displaystyle
V_{(2\nu-2)-\nu}(r)\equiv-{\frac{A}{r^{\nu}}}+{\frac{B}{r^{2\nu-2}}}.$

In the mentioned article, a special attention has been drawn to the $\nu=6$
case, and its ability to fit experimental data:

$\displaystyle V_{10-6}(r)\equiv-{\frac{A}{r^{6}}}+{\frac{B}{r^{10}}}.$ (2.14)

In Section 3 we will explore the interesting features of (2.14) in the realm
of SUSYQM.

### The second virial coefficient and its dependence on the potential

The virial equation of state for a gas expresses the deviation from the ideal
behavior as a power series in the density $\rho$:

$\displaystyle\frac{p}{kT}=\rho+B_{2}(T)\rho^{2}+B_{3}(T)\rho^{3}+\cdots.$

The coefficients $B_{n}(T)$ are called the virial coefficients and they are
unique real functions of the temperature. The second virial coefficient
$B_{2}(T)$ represents the most significant deviation from the ideal behavior,
since it is the prefactor in the term of order $\rho^{2}$ in the series. It is
a customary result from equilibrium statistical mechanics (see [36]) that
$B_{2}(T)$ is a radial integral of the pair-potential $v(r)$ given by

$\displaystyle B_{2}(T)=2\pi\int^{\infty}_{0}\big{(}1-{\rm
e}^{-\frac{v(r)}{kT}}\big{)}r^{2}{\rm d}r.$ (2.15)

A thorough study by Keller and Zumino of the properties of (2.15) has shown
that a unique potential function can only be obtained from $B_{2}(T)$ if the
potential behaves monotonically [28]. This is clearly not the case for the
Lennard-Jones potential and all its variants. As a result, there exists an
ambiguity in the choice of the microscopic potential, leading to the same
thermodynamic function $B_{2}(T)$. In addition to this analytic inexactness
there is also the limited range of measurements of $B_{2}(T)$ for low
temperatures. The aforementioned limitations lead to several possibilities of
choice for $v(r)$, at least from measurements of $B_{2}$, specially for the
power of the repulsive term $B/r^{\delta}$. The possibilities range from $n=9$
to $n=14$ since the early works of Lennard-Jones (see [32, 33]) and De Boer
(see [19]). We come back to this point in the next section, giving some hints
about the applicability of the $10-6$ potential for low temperatures.

### The dimensionless Schrödinger equation
and the De Boer principle of corresponding states

In 1948 J. De Boer introduced a dimensionless representation of the
Schrödinger equation employing $\sigma$ and $\epsilon$ in order to construct
dimensionless lengths and energies [18]

$\displaystyle\tilde{r}\equiv\frac{r}{\sigma},\qquad\tilde{E}\equiv\frac{E}{\epsilon},\qquad\tilde{V}\equiv\frac{V}{\epsilon}.$
(2.16)

As a result, the radial Schrödinger equation (2.9) for $l=0$ can be
transformed into the dimensionless form given by

$\displaystyle\left[-\frac{\Lambda^{2}}{2}\frac{{\rm d}^{2}}{{\rm
d}\tilde{r}^{2}}+\tilde{V}(\tilde{r})\right]u(\tilde{r})=\tilde{E}u(\tilde{r})$
(2.17)

provided that the potential $V(r)$ can be expressed in the generic form
$V(r)=\epsilon f(r/\sigma)$, where $f(r)$ is a well-defined dimensionless
interaction function and $\Lambda\equiv\hbar/(\sigma\sqrt{\mu\epsilon})$ [18].
From (2.17) we see that $\Lambda$, the so-called De Boer parameter, is the
only parameter in the equation that gives information about the particular
microscopic characteristics of the system. From this fact, De Boer was able to
formulate his principle of corresponding states, which is a “quantum”
generalization of the van der Waals law of corresponding states for classical
gases and liquids [23, 44, 50]. The De Boer principle of corresponding states
tells us that two different systems with equal value of $\Lambda$ have
identical thermodynamic properties [18]. In Section 3 we exploit this
principle in order to give an interpretation to the supersymmetric integrable
model for zero energy, we propose with the $10-6$ Lennard-Jones potential.

### Supersymmetric quantum mechanics

We implement in this work the simplest realisation of SUSYQM for one-
dimensional quantum systems [21], which includes besides the Hamiltonian
operator $H$, two _fermionic_ operators $Q^{\pm}$ or _supercharges_ such that
they commute with $H$

$\displaystyle\big{[}Q^{\pm},H\big{]}=0$ (2.18)

and satisfy the algebra

$\displaystyle\big{\\{}Q^{-},Q^{+}\big{\\}}=H,\qquad\big{(}Q^{\pm}\big{)}^{2}=0.$
(2.19)

The second relation means that $Q^{\pm}$ are nilpotent operators. A usual
representation of the algebra, given in equations (2.18) and (2.19), presents
the Hamiltonian $H$ of the system, as a diagonal two component matrix of
partner Hamiltonians $H_{\pm}$

$\displaystyle H\equiv\left(\begin{matrix}H_{+}&0\\\
0&H_{-}\end{matrix}\right),$

where $Q^{\pm}$ are $2\times 2$ diagonal matrices involving the _Ladder_
operators $A^{\pm}$

$\displaystyle Q^{-}\equiv\left(\begin{matrix}0&0\\\
A^{-}&0\end{matrix}\right),\qquad Q^{+}\equiv\left(\begin{matrix}0&A^{+}\\\
0&0\end{matrix}\right),$

such that

$\displaystyle H\equiv\big{\\{}Q^{-},Q^{+}\big{\\}}\equiv
Q^{-}Q^{+}+Q^{+}Q^{-}=\left(\begin{matrix}A^{+}A^{-}&0\\\
0&A^{-}A^{+}\end{matrix}\right)\equiv\left(\begin{matrix}H_{+}&0\\\
0&H_{-}\end{matrix}\right),$ (2.20)

and $A^{\pm}$ are defined in terms of the derivative $\frac{{\rm d}}{{\rm
d}x}$ and an arbitrary complex function $w(r)$, called the _superpotential_

$\displaystyle A^{\pm}=\mp\frac{{\rm d}}{{\rm d}r}+w(r).$ (2.21)

Since the products $A^{+}A^{-}$ and $A^{-}A^{+}$ with $A^{\pm}$ defined in
(2.21) lead to

$\displaystyle A^{+}A^{-}=-\frac{{\rm d}^{2}}{{\rm d}r^{2}}+w^{2}-\frac{{\rm
d}w}{{\rm d}r},\qquad A^{-}A^{+}=-\frac{{\rm d}^{2}}{{\rm
d}r^{2}}+w^{2}+\frac{{\rm d}w}{{\rm d}r},$

then, from (2.20) it results natural to identify

$\displaystyle H_{\pm}=-\frac{{\rm d}^{2}}{{\rm d}r^{2}}+w^{2}\pm\frac{{\rm
d}w}{{\rm d}r},$

which leads directly to a definition of the so-called partner potentials
$v_{\pm}$ given by

$\displaystyle v_{\pm}\equiv w^{2}\pm\frac{{\rm d}w}{{\rm d}r},$ (2.22)

such that

$\displaystyle H_{\pm}=-\frac{{\rm d}^{2}}{{\rm d}r^{2}}+v_{\pm}.$

Each of the two equations in (2.22) define a Riccati differential equation for
the superpotential $w$, if $v_{\pm}$ are known. Let’s recall that the
superpotential can also be found from the zero-energy base state $\psi_{0}$,
by computing $w=-\psi^{\prime}_{0}/\psi_{0}$, where $\psi_{0}$ is a solution
of the Schrödinger equation with the $v_{-}$ potential

$\displaystyle\psi^{\prime\prime}_{0}=v_{-}\psi_{0},\qquad\psi_{0}={\rm
e}^{-\int w{\rm d}r}$ (2.23)

(see Witten [51]). Riccati equations play a fundamental role in the study of
integrability in SUSYQM. For a systematic study of this subject see references
[1, 3].

### Differential Galois theory

Exact solutions of differential equation is a hard but important task in
different disciplines. Sometimes numerical methods cannot be implemented in
general, if the equation has free generic parameters. _Differential Galois
theory_ , also known as _Picard–Vessiot theory_ , is a powerful theory to
solve explicitly, in the case when it is possible, linear differential
equations.

Analogous to the concept of field in classical Galois theory, there exists the
concept _differential field_ in differential Galois theory, which is a field
satisfying the differential Leibniz rules. Similarly, a differential extension
$L$ of the differential field $K$ means that $K$ is a subfield of $L$
preserving the differential Leibniz rules. In particular for a given linear
differential with coefficients in $K$, if $C_{L}=C_{K}$ (the field of
constants of $L$ is the same field of constants of $K$) and $L$ is generated
over $K$ by a fundamental set of solutions of such differential equation, then
$L$ is called the _Picard–Vessiot extension of $K$_. Recall that the field of
constants of $K$ is defined as $C_{K}:=\\{k\in K\colon k^{\prime}=0\\}$, where
${}^{\prime}:={\rm d}/{\rm d}x$.

In the same way as we are interested in finding the roots of the polynomials
over a base field, usually $\mathbb{Q}$, using arithmetical and algebraic
conditions, we would like to have explicit solutions of differential equations
over a differential base field $K=\mathbb{C}(x)$, with field of constants
$C_{K}=\mathbb{C}$, using elementary functions and quadratures. The
differential Galois theory considers more general differential fields, but for
our purpose is enough to consider $\mathbb{C}(x)$. Thus, the differential
Galois group ($\mathrm{DGal}(L/K)$), as analogically as in the polynomial
case, is the group of all _differential automorphisms_ that restricted to the
base field coincide with the identity. Moreover if $\langle
y_{1},y_{2},\ldots,y_{n}\rangle$ is a basis of solutions of

$\displaystyle\frac{{\rm d}^{n}y}{{\rm d}x^{n}}+a_{n-1}\frac{{\rm
d}^{n-1}y}{{\rm d}x^{n-1}}+\cdots+a_{1}\frac{{\rm d}y}{{\rm
d}x}+a_{0}y=0,\qquad a_{i}\in\mathbb{C}(x),$

then for each differential automorphism $\sigma\in\mathrm{DGal}(L/K)$ there
exists a matrix $A_{\sigma}\in\mathrm{GL}(n,\mathbb{C})$ (i.e.,
$a_{ij}\in\mathbb{C}$, $1\leq i,j\leq n$ and $\det(A_{\sigma})\neq 0$) such
that

$\displaystyle\sigma(\mathbf{Y})=A_{\sigma}\mathbf{Y},\qquad\mathbf{Y}=\begin{pmatrix}z_{1}\\\
z_{2}\\\ \vdots\\\ z_{n}\end{pmatrix},$ $\displaystyle
A_{\sigma}=\begin{pmatrix}\alpha_{11}&\alpha_{12}&\ldots&\alpha_{1n}\\\
\alpha_{21}&\alpha_{22}&\ldots&\alpha_{2n}\\\ \vdots&\vdots&\vdots&\vdots\\\
\alpha_{n1}&\alpha_{n2}&\ldots&\alpha_{nn}\end{pmatrix},\qquad\mathrm{DGal}(L/K)\cong
G\subset\mathrm{GL}(n,\mathbb{C}).$

In particular,
$\mathrm{SL}(n,\mathbb{C})=\\{A\in\mathrm{GL}(n,\mathbb{C})\colon\det(A)=1\\}$.
Due to
$G=\\{A_{\sigma}\colon\sigma\in\mathrm{DGal}(L/K)\\}\subset\mathrm{GL}(n,\mathbb{C})$,
we see that $\mathrm{DGal}(L/K)$ has a faithful representation as an algebraic
group of matrices in where $G^{0}$ denotes the connected identity component of
$G$ (the biggest algebraic connected subgroup of $G$). In this terminology, we
say that a linear differential equation is _integrable in the sense of
differential Galois theory_ whether the connected identity component of its
differential Galois is a solvable group. Moreover, this definition of
integrability leads to the obtaining of solutions in closed form if and only
if $G^{0}$ is solvable, see [49] for full explanation and details. From now
on, integrable in this paper means integrable in terms of differential Galois
theory, see [47].

To accomplish our purposes, we are interested in second-order differential
equations of the form

$\displaystyle z^{\prime\prime}+az^{\prime}+bz=0,\qquad a,b\in\mathbb{C}(x).$
(2.24)

Equation (2.24) can be transformed into equations in the form

$\displaystyle y^{\prime\prime}=ry,\qquad r={a^{2}\over 4}+{a^{\prime}\over
2}-b,\qquad\text{and}\qquad z={\rm e}^{{-1\over 2}\int a{\rm d}x}y,$ (2.25)

see [10]. Jerald Kovacic developed in 1986 an algorithm to solve explicitly
second-order differential equations with rational coefficients given in the
form of equation (2.25), see [30]. In [20] another version of Kovacic’s
algorithm is presented, and it is applied to solve several second-order
differential equations with special functions as solutions. The version of
Kovacic’s algorithm presented here corresponds to [6], see also [1, 10, 11,
14].

As mentioned, Kovacic’s algorithm cannot be applied when the coefficients of
the second-order differential equations are not rational functions. Therefore
we need to transform such differential equations to apply Kovacic’s algorithm.
A possible solution to this problem was developed in [1, 3, 11], the so-called
Hamiltonian algebrization. However, we are interested in transformations that
preserve the differential Galois group (at least their connected identity
component), in other words, the transformation must be either _isogaloisian_ ,
_virtually isogaloisian_ or _strongly isogaloisian_ , see [1, 11].

One important differential equation in this work is the Whittaker’s
differential equation, which is given by

$\displaystyle\partial_{x}^{2}y=\left(\frac{1}{4}-{\kappa\over
x}+{4\mu^{2}-1\over 4x^{2}}\right)y.$ (2.26)

The Galoisian structure of this equation has been deeply studied in [45], see
also [20]. The following theorem provides the conditions of the integrability
in the sense of differential Galois theory of equation (2.26).

######  Theorem 2.1 ([45]).

The Whittaker’s differential equation (2.26) is integrable $($in the sense of
differential Galois theory$)$ if and only if either,
$\kappa+\mu\in\frac{1}{2}+\mathbb{N}$, or
$\kappa-\mu\in\frac{1}{2}+\mathbb{N}$, or
$-\kappa+\mu\in\frac{1}{2}+\mathbb{N}$, or
$-\kappa-\mu\in\frac{1}{2}+\mathbb{N}$.

The Bessel’s equation is a particular case of the confluent hypergeometric
equation and is given by

$\displaystyle\partial_{x}^{2}y+{1\over x}\partial_{x}y+{x^{2}-n^{2}\over
x^{2}}y=0.$ (2.27)

Under a suitable transformation, the reduced form of the Bessel’s equation is
a particular case of the Whittaker’s equation. Thus, we can obtain the
following well known result, see [29, p. 417] and see also [30, 40]:

######  Corollary 2.2.

The Bessel’s differential equation (2.27) is integrable $($solvable by
quadratures$)$ if and only if $n\in\frac{1}{2}+\mathbb{Z}$.

######  Definition 2.3 (Hamiltonian change of variable, [6]).

A change of variable $z=z(x)$ is called Hamiltonian if
$(z(x),\partial_{x}z(x))$ is a solution curve of the autonomous classical one
degree of freedom Hamiltonian system

$\displaystyle\partial_{x}z=\partial_{w}H,\qquad\partial_{x}w=-\partial_{z}H\qquad\textrm{with}\quad
H=H(z,w)={w^{2}\over 2}+V(z),$

for some $V\in K$.

######  Proposition 2.4 (Hamiltonian algebrization, [6]).

The differential equation

$\displaystyle\partial_{x}^{2}{r}=q(x)r$

is algebrizable through a Hamiltonian change of variable $z=z(x)$ if and only
if there exist $f$, $\alpha$ such that

$\displaystyle{\partial_{z}\alpha\over\alpha},\qquad{f\over\alpha}\in\mathbb{C}(z),\qquad\text{where}\qquad
f(z(x))=q(x),\qquad\alpha(z)=2(H-V(z))=(\partial_{x}z)^{2}.$

Furthermore, the algebraic form of the equation $\partial_{x}^{2}{r}=q(x)r$ is

$\displaystyle\partial_{z}^{2}y+{1\over
2}{\partial_{z}\alpha\over\alpha}\partial_{z}y-{f\over\alpha}y=0,\qquad
r(x)=y(z(x)).$

Next, we follow the references [1, 6, 11] to describe Kovacic’s algorithm.
Thus, to solve second-order differential equations with rational coefficients
we use should Kovacic’s algorithm, which is presented in Appendix A.

##  3 Main results

In this section we present the main contributions of this paper. First we will
show that for the usual $\nu=6$, $\delta=12$ Lennard-Jones potential, the
Schrödinger equation is non-integrable in the sense of differential Galois
theory for any value of energy. In contrast for $\delta=10$ and $\nu=6$ we
show the integrability, in the sense of differential Galois theory, as a
special case of a general theorem for $\delta=2\nu-2$ with
$\delta,\nu\in\mathbb{N}$ (see Theorem 3.3 and its subsequent remark). From
the physical point of view, the $10-6$ case is of the most remarkable
importance. Since we preserve the physically grounded $-{1/r}^{6}$ term coming
from dipole-dipole interactions and responsible of the van der Waals forces;
but we replace never the less, the rather arbitrary ${1/r}^{12}$ term
responsible for the repulsion of the particles in the many body system, and
leading to a non-integrable differential equation; with an equally arbitrary
${1/r}^{10}$ term, but leading to an integrable one. We will dedicate the
subsequent sections to show the advantages and physical interest of this
special choice (see Fig. 2 for a graphic comparison of both potentials).

Figure 2: $\frac{v(r)}{\epsilon}$ vs. $\frac{r}{\sigma}$ plot for the rescaled
$12-6$ (black) and $10-6$ (grey) Lennard-Jones potentials given in equation
(2.13) for the same molecular parameters $\sigma$, $\epsilon$ and $C=0$.

We start by considering the radial Schrödinger equation (2.7) with the
generalized effective Lennard-Jones potential (2.13)

$\displaystyle\mathcal{H}u_{k,l}(r)=\varepsilon_{k,l}u_{k,l}(r),\qquad\mathcal{H}\equiv-\frac{{\rm
d}^{2}}{{\rm d}r^{2}}+v_{\delta-\nu}^{\rm eff}(r)=-\frac{{\rm d}^{2}}{{\rm
d}r^{2}}-{\frac{\bar{A}}{r^{\nu}}}+{\frac{\bar{B}}{r^{\delta}}}+\frac{\bar{C}}{r^{2}},$
$\displaystyle 0<\nu<\delta\in\mathbb{N\subset
Z},\qquad\bar{A}>0,\qquad\bar{B}>0,\qquad\bar{C}\geq 0.$ (3.1)

Setting $\mathbb{C}(r)$ as the differential field of equation (3.1) with the
derivative $\frac{{\rm d}}{{\rm d}r}$, we set also
$\bar{A},\bar{B},\bar{C}\in\mathbb{C}$.

######  Theorem 3.1.

Schrödinger equation with original $12-6$ Lennard-Jones effective potential is
not integrable in the sense of differential Galois theory for any value of the
energy and for all $A,B\in\mathbb{C}^{*}$, $C\in\mathbb{C}$.

###### Proof.

Considering $\nu=6$ and $\delta=12$ in equation (3.1) we arrive to the
Schrödinger equation with effective original $12-6$ Lennard-Jones potential.
Now, applying the Hamiltonian change of variable $z=r^{2}$ over such
Schrödinger equation we arrive to the differential equation

$\displaystyle
u^{\prime\prime}_{k,l}+\frac{1}{2z}u^{\prime}_{k,l}+\left(\frac{A}{4z^{4}}-\frac{B}{4z^{7}}-\frac{C}{4z^{2}}+\frac{\varepsilon_{k,l}}{4z}\right)u_{k,l}=0.$

Now, the change of dependent variable

$\displaystyle u_{k,l}={\Phi_{k,l}\over\sqrt[4]{z}}$

leads to the differential equation

$\displaystyle\Phi_{k,l}^{\prime\prime}=\left({-4\varepsilon_{k,l}z^{6}+(4C-3)z^{5}-4Az^{3}+4B\over
16z^{7}}\right)\Phi_{k,l}.$ (3.2)

After applying Kovacic’s algorithm, see Appendix A, we observe that equation
(3.2) falls in case 4 for any $\varepsilon_{k,l}\in\mathbb{C}$ because there
are not suitable conditions in step 1 for case 1 and case 3. The second step
is not satisfied in case 2 because $D=\varnothing$ due to $E_{0}=\\{7\\}$,
$E_{\infty}=\\{1,2\\}$ and there are not integers satisfying the condition for
$D\neq\varnothing$. Thus we conclude that Schrödinger equation with original
$12-6$ Lennard-Jones effective potential is not integrable in the sense of
differential Galois theory for any value of the energy. ∎

### Supersymmetric quantum mechanics and the Lennard-Jones superpotential

The implementation of Hamiltonian algebrization and Kovacic’s algorithm
reaches a considerable power in the realm of supersymmetric quantum mechanics.
In fact the integrability of second-order linear equations like the radial
Schrödinger equation (3.1) subject of our study, via the Kovacic’s algorithm
is deeply related with the properties of the solutions of the associated
Riccati equation in the supersymmetric extension of the theory [1]. Taking
this as a motivation, we go further in this section and propose a
superpotential leading to the non-effective part (2.12) of (2.13) ($C=0$) as
one of two partner potentials. If we denote the superpotential in one
dimension as $w(r)$ the corresponding partner potentials are given by equation
(2.22) (see [1, 21, 51], among others)

$\displaystyle v_{\pm}(r)\equiv w^{2}(r)\pm\frac{{\rm d}w}{{\rm d}r}$ (3.3)

corresponding for each case to a Riccati equation for $w$. Knowing that
$0<\nu<\delta$ we identify terms in (3.3) with terms in (2.12) as follows

$\displaystyle w^{2}(r)={\frac{\bar{B}}{r^{\delta}},\qquad\frac{{\rm d}w}{{\rm
d}r}=\frac{\bar{A}}{r^{\nu}}}$

as a consequence we have

$\displaystyle w(r)=\pm{\frac{\sqrt{\bar{B}}}{r^{\frac{\delta}{2}}},\qquad
w(r)=-\frac{\bar{A}}{(\nu-1)r^{\nu-1}}+C}.$

A simple choice for $w(r)$ is given by

$\displaystyle w(r)\equiv-{\frac{\sqrt{\bar{B}}}{r^{\frac{\delta}{2}}}},$
(3.4)

where the following identities should hold

$\displaystyle\sqrt{\bar{B}}=\bar{A}/(\nu-1),\qquad\delta=2(\nu-1),\qquad
C=0.$ (3.5)

As a result we identify $v_{\delta-\nu}(r)$ in (2.12) with $v_{-}$ and we have
from (3.3), the following expressions for the partner potentials

$\displaystyle
v_{-}={\frac{\bar{B}}{r^{\delta}}-\frac{\bar{A}}{r^{\nu}}=v_{\delta-\nu}(r)},\qquad
v_{+}={\frac{\bar{B}}{r^{\delta}}+\frac{\bar{A}}{r^{\nu}}}.$ (3.6)

According to equation (2.23), the corresponding wave function for the zero
energy level using $v_{-}$ is given by

$\displaystyle\psi_{0}(r)={\rm
e}^{-\frac{2\sqrt{\bar{B}}}{(\delta-2)r^{\frac{\delta-2}{2}}}}={\rm
e}^{-\frac{A}{(\nu-1)(\nu-2)r^{\nu-2}}+Cr}.$ (3.7)

Figure 3: Wave function for $v_{10-6}$ with zero energy, $\bar{A}=5$,
$\bar{B}=1$ and $C=0$.

An example of wave function for this potential is given in Fig. 3. Summarizing
we conclude that expressions in (3.5) set conditions for the existence of a
superpotential in the form of (3.4), and a supersymmetric extension to
equation (3.1) with partner potentials in (3.6). The case for
$\delta=2(\nu-1)$ thus appears in a natural way, as a simple condition for
defining a supersymmetric model. The property of integrability for zero energy
of this case to be proven in Theorem 3.3, makes it an appealing model to
further explore the relation between supersymmetry and integrability already
studied in [1].

### The $\boldsymbol{10-6}$ Lennard-Jones superpotential and the De Boer
parameter

As mentioned before the case for $\nu=6$, $\delta=10$ is of particular
relevance from physical grounds. One of the aims of this work is to explore
the analytical advantages of $v_{10-6}$ in contrast to $v_{12-6}$. The $10-6$
potential in terms of the molecular parameters $\sigma$ and $\epsilon$ is
given by

$\displaystyle
V_{10-6}(r)=\alpha\epsilon\left[\left(\frac{\sigma}{r}\right)^{10}-\left(\frac{\sigma}{r}\right)^{6}\right],$
(3.8)

where $\alpha$ is chosen so that $\epsilon$ is the minimum energy (the well
depth) and $\sigma$, as mentioned before, is the value where $V_{10-6}$
vanishes. As a result we have for this case111Notice in (2.11) that
$\alpha\equiv 4$ for the $12-6$ potential. $\alpha\equiv(25/6)\sqrt{5/3}$.
Thus the rescaled $10-6$ Lennard-Jones potential reads

$\displaystyle
v_{10-6}(r)\equiv\left(\frac{2\mu}{\hbar^{2}}\right)V_{10-6}(r)=\left(\frac{2\mu\alpha\epsilon}{\hbar^{2}}\right)\left[\left(\frac{\sigma}{r}\right)^{10}-\left(\frac{\sigma}{r}\right)^{6}\right]$
(3.9)

or in the equivalent $A-B$ form, we have the following

$\displaystyle
v_{10-6}(r)=\frac{\bar{B}}{r^{10}}-\frac{\bar{A}}{r^{6}}\qquad\text{with}\qquad\bar{A}\equiv\frac{2\mu\alpha\epsilon\left(\sigma\right)^{6}}{\hbar^{2}},\qquad\text{and}\qquad\bar{B}\equiv\frac{2\mu\alpha\epsilon(\sigma)^{10}}{\hbar^{2}}.$
(3.10)

Clearly $v_{10-6}$ fulfills the condition in (3.5) for $\nu=6$; as a result
the rescaled superpotential for $\delta=10$ in (3.4) takes the form

$\displaystyle w_{10-6}(r)=-{\frac{\sqrt{\bar{B}}}{r^{5}}},$

where $\bar{A}$ $=5\sqrt{\bar{B}}$ as we easily check from (3.5). Equivalently

$\displaystyle
w_{10-6}(r)=-\frac{\bar{A}}{5r^{5}}=-\frac{2\mu\alpha\epsilon(\sigma)^{6}}{5\hbar^{2}r^{5}}=-\frac{5\mu\sqrt{5/3}\epsilon(\sigma)^{6}}{3\hbar^{2}r^{5}}.$
(3.11)

The condition $\bar{A}$ $=5\sqrt{\bar{B}}$ can be written in the following
suggestive dimensionless form

$\displaystyle\frac{\hbar^{2}}{\mu\epsilon(\sigma)^{2}}=\frac{1}{3}\sqrt{\frac{5}{3}}\approx
0.4303,$ (3.12)

where we have applied definitions in (3.10) and
$\alpha\equiv(25/6)\sqrt{5/3}$. We will call it from now on the supersymmetric
condition (for short SUSY condition) for the $10-6$ Lennard-Jones potential.
In terms of the so-called De Boer parameter
$\Lambda\equiv\hbar/(\sigma\sqrt{\mu\epsilon})$, which gives a degree of the
quantum character of the system [18], we have $\Lambda^{2}\approx\allowbreak
0.4303$ or similarly $\Lambda\approx 0.6559$. An important remark at this
point is that the SUSY condition given in the form $\bar{A}$
$=5\sqrt{\bar{B}}$ will appear again in Theorems 3.2 and 3.3, in the context
of the Martinet–Ramis theorem, that is, Theorem 2.1.

As a summary of this section, we conclude that the fulfillment of condition
(3.12), guarantees not only the solvability of the Schrödinger equation (2.8)
with the potential $v_{10-6}$ in (3.9) (through the Martinet–Ramis theorem, as
we will see) but also the existence of a superpotential given by expression
(3.11), which correspondingly leads to $v_{10-6}$ as one of the partner
potentials $v_{\pm}(r)$ defined through the Riccati equations in (3.3). The
supersymmetric model thus formulated considers a specific version of the
radial Schrödinger equation (2.9) or equivalently the rescaled form (2.8),
where we set in both equations $l=0$ for the angular momentum, and $V_{10-6}$
and $v_{10-6}$ are given in (3.8) and (3.9). We will start in the next
section, a physical analysis of the model, in the light of the De Boer
principle of corresponding states.

### The low temperature behavior of the $\boldsymbol{10-6}$ Lennard-Jones gas

Recalling the discussion in the previous section, about the dimensionless
representation (2.17) of the Schrödinger equation we start by noticing that
the $10-6$ potential (3.8) fulfills the condition $V(r)=\epsilon f(r/\sigma)$
if we identify $f(r/\sigma)$ with
$\alpha\big{[}\big{(}\frac{\sigma}{r}\big{)}^{10}-\big{(}\frac{\sigma}{r}\big{)}^{6}\big{]}$;
as a result we have

$\displaystyle\tilde{V}_{10-6}(r)\equiv\frac{V_{10-6}(r)}{\epsilon}=(25/6)\sqrt{5/3}\left[\left(\frac{\sigma}{r}\right)^{10}-\left(\frac{\sigma}{r}\right)^{6}\right]$
$\displaystyle\hphantom{\tilde{V}_{10-6}(r)}{}=(25/6)\sqrt{5/3}\left[\left(\frac{1}{\tilde{r}}\right)^{10}-\left(\frac{1}{\tilde{r}}\right)^{6}\right],$
(3.13)

where we have used the definitions in (2.16) $\tilde{r}\equiv\frac{r}{\sigma}$
and $\tilde{V}\equiv\frac{V}{\epsilon}$. The Schrödinger equation takes thus
the form in (2.17):

$\displaystyle\left[-\frac{\Lambda^{2}}{2}\frac{{\rm d}^{2}}{{\rm
d}\tilde{r}^{2}}+\tilde{V}_{10-6}(r)\right]u(r)=\tilde{E}u(r).$ (3.14)

As mentioned in Section 2, the De Boer principle of corresponding states tells
us that two different systems with equal value of $\Lambda$ have identical
thermodynamical properties [18]. In this sense the SUSY condition
$\Lambda^{2}=\hbar^{2}/\big{[}(\sigma)^{2}\mu\epsilon\big{]}=(1/3)\sqrt{5/3}\approx
0.4303$ given in (3.12); and representing a definite set of combinations of
values of the parameters $\sigma$, $\epsilon$, and $\mu$; that accounts for
$\Lambda^{2}=(1/3)\sqrt{5/3}$; is defining through the principle of
corresponding states, a specific set of physical systems with equivalent
thermodynamical properties. These systems have the special feature of being
described by a Supersymmetric potential of the form (3.11) leading to (3.14)
with the potential (3.13) as the partner potential
$V_{-}(r)\equiv\big{(}\frac{\hbar^{2}}{2\mu}\big{)}v_{-}(r)$ in (3.6) with
$\nu=6$.

We have found after a brief review of the literature, a significant
coincidence between the specific value for $\Lambda^{2}\approx 0.4303$ and the
value of $\Lambda^{2}=0.456$ reported by Miller, Nosanow and Parish [39] for a
second-order liquid to gas phase transition of a Bose–Einstein condensate at
zero temperature. Since their calculation is an approximate one, made in the
framework of the variational method; it is a worthy task (to be done
elsewhere) to investigate the advantages of our exact approach to the
calculation of properties of such many-body systems at low temperatures in the
context of the quantum extension to the principle of corresponding states.

In Fig. 4 we see a plot of the second virial coefficient calculated
numerically for both the $12-6$ and $10-6$ potentials from the integral
definition in (2.15). Relying on $B_{2}(r)$ as a quantity that gives
information of the microscopic pair-potential (with the previously mentioned
limitations) we see an asymptotic closeness of both functions for low
temperatures, that hints for the reliability of our supersymetric model with
the $10-6$ Lennard-Jones potential, near the absolute zero.

Figure 4: $B_{2}(T)$ vs. $T/\epsilon$ plot for the $12-6$ (black) and $10-6$
(grey) Lennard-Jones potentials for the same molecular parameters $\sigma$ and
$\epsilon$. The closeness of both functions for low temperatures near to
absolute zero is a hint of the reliability of the $10-6$ potential in that
region.

### Integrability of the $\boldsymbol{10-6}$ Lennard-Jones potential and its
generalization

The following result is valid for any potential $v(r)$ belonging to a
differential field.

######  Theorem 3.2.

Consider $v(r)$ belonging to a differential field $K$, then the following
statements hold

  * •

The only one change of dependent variable that allows to transform the radial
equation of the Schrödinger equation into the Schrödinger equation with
effective potential is $\varphi\colon u\mapsto\varphi(u)=ru$, where $u$ is the
solution for the radial equation and $\varphi(u)$ is the solution for the
Schrödinger equation with effective potential.

  * •

The differential Galois groups of the radial equation and Schrödinger equation
with effective potential are subgroups of $\mathrm{SL}(2,\mathbb{C})$.

  * •

The transformation $\varphi$ is strongly isogaloisian.

###### Proof.

We proceed according to each item.

  * •

Applying the transformation given in equation (2.24) and equation (2.25) we
obtain it because $2/r$ is the coefficient of the first derivative of the
radial equation after separation of variables. Thus, applying the change of
variable $\varphi\colon u\mapsto\varphi(u)=ru$ we arrive to the Schrödinger
equation with effective potential.

  * •

The Wronskian of two independent solutions of the Schrödinger equation with
effective potential is constant and constants are in the base field.
Similarly, the Wronskian of two independent solutions of the radial equation
belongs to the base field. Therefore, automorphisms over such solutions acts
by multiplication of matrices belonging to $\mathrm{SL}(2,\mathbb{C})$, that
is, $\sigma(U)=A_{\sigma}U$, $\sigma(\varphi(U))=A_{\sigma}\varphi(U)$ and
$\det(A_{\sigma})=1$. Thus, $A_{\sigma}\in\mathrm{SL}(2,\mathbb{C})$.

  * •

Applying the differential automorphism $\sigma$ over $\varphi(u)$ we observe
that $\sigma(\varphi(u))=\sigma(r)\sigma(u)=r\sigma(u)$, which implies that
differential Galois group only depends on the solutions $u$ because $r$ is in
the base field and the differential Galois group will be the same with the
same base field. Thus, the transformation $\varphi$ is strongly isogaloisian.

Thus we conclude the proof. ∎

The following result corresponds to the integrability conditions for the
$10-6$ L-J potential and its generalization.

######  Theorem 3.3.

The Schrödinger equation with $(2\nu-2)-\nu$ L-J potential, given in equation
(3.1), is integrable for zero energy in the sense of differential Galois
theory if and only if

$\displaystyle A=\pm\sqrt{B}\big{(}\pm\sqrt{1+4C}+\nu-2+2mk-4m\big{)},\qquad
m\in\mathbb{Z}.$

###### Proof.

The Schrödinger equation given in equation (3.1), with zero energy, is
transformed into the Whittaker’s differential equation (2.26) through the
change of variables

$\displaystyle u_{k,l}=\sqrt{r^{\nu-1}}\phi_{k,l},\qquad
r=\sqrt[\nu-2]{2\sqrt{B}\over(\nu-2)z}$

with parameters

$\displaystyle\kappa={A\over\sqrt{B}(2\nu-4)}\qquad\text{and}\qquad\mu={\sqrt{1+4C}\over
2\nu-4}.$

Applying Martinet–Ramis theorem we have

$\displaystyle\pm{A\over\sqrt{B}(2\nu-4)}\pm{\sqrt{1+4C}\over
2\nu-4}\in\mathbb{Z}+\frac{1}{2}.$

Assuming $m\in\mathbb{Z}$ we obtain

$\displaystyle A=\pm\sqrt{B}\big{(}\pm\sqrt{1+4C}+\nu-2+2m\nu-4m\big{)},\qquad
m\in\mathbb{Z},$

which is the integrability condition for the Schrödinger equation with
$(2\nu-2)-\nu$ L-J potential and its wave function corresponds to equation
(3.7) with $\delta=2\nu-2$. ∎

######  Remark 3.4.

We observe that Theorem 3.3 refers to the integrability in the sense of
differential Galois theory, which is not related with square integrable wave
functions. Another key point is that we are not considering energies different
than zero and excited states, this is an open problem for this generalized
potential. In particular, the theorem includes the $10-6$ L-J potential, for
$C=0$ and $\nu=6$. Therefore the Schrödinger equation with $L-J$ 10-6 is
integrable for zero energy when $A=\pm\sqrt{B}(8m+4\pm 1)$, while energies
different than zero and excited states were not considered in this paper.
Moreover, for zero energy and $m=-1$ we recover the integrability condition
obtained through SUSYQM for this potential, i.e., the Schrödinger equation
with $10-6$ L-J potential is also integrable in the sense of differential
Galois theory for $A\in\big{\\{}\pm 3\sqrt{B},\pm 5\sqrt{B}\big{\\}}$.

##  4 Final remarks and open questions

In this paper we have shown that there exist no explicit solutions of the
radial Schrödinger equation with the usual $12-6$ Lennard-Jones potential for
any value of the energy. We have proposed an alternative supersymmetric model
with a $10-6$, $v_{-}$ partner potential, that preserves the $-1/r^{6}$ van
der Waals attraction. We have found through the De Boer principle of
corresponding states, initial hints that this model could represent a low
temperature system determined by a $\Lambda^{2}\approx 0.4303$ value of the
2nd power of the De Boer parameter. We have studied possible generalizations
of the Lennard-Jones potential, where the Schrödinger equation is integrable
in the sense of differential Galois theory.

Further work can be developed looking for similar theorems of integrability in
the sense of differential Galois theory for $E\neq 0$ and excited states, for
the $10-6$ potential and other generalizations. Relations between square
integrable wave functions and solutions in closed form of SE for
generalizations of L-J potentials should be explored in further works too.

We hope that this paper can be the starting point of further works involving
SUSYQM, DGT and statistical mechanics, which are not easy topics. Although we
tried to write a readable preliminaries about these topics, we know that it
was not enough and the reader should complement with references suggested by
us, otherwise this paper could be a large paper, which was not the target.

##  Appendix A Kovacic algorithm

The version of Kovacic’s algorithm presented in this appendix is based in the
improved version given in [6]. There are four cases in Kovacic’s algorithm.
Only for cases 1, 2 and 3 we can solve the differential equation, but for the
case 4 the differential equation is not integrable. It is possible that
Kovacic’s algorithm can provide us only one solution ($\zeta_{1}$), so that we
can obtain the second solution ($\zeta_{2}$) through

$\displaystyle\zeta_{2}=\zeta_{1}\int\frac{{\rm d}x}{\zeta_{1}^{2}}.$

Notations. For the differential equation given by

$\displaystyle\partial_{x}^{2}\zeta=r\zeta,\qquad r={s\over t},\qquad
s,t\in\mathbb{C}[x],$

we use the following notations:

  1. 1)

denote by $\Gamma^{\prime}$ be the set of (finite) poles of $r$,
$\Gamma^{\prime}=\\{c\in\mathbb{C}\colon t(c)=0\\}$,

  2. 2)

denote by $\Gamma=\Gamma^{\prime}\cup\\{\infty\\}$,

  3. 3)

by the order of $r$ at $c\in\Gamma^{\prime}$, $\circ(r_{c})$, we mean the
multiplicity of $c$ as a pole of $r$,

  4. 4)

by the order of $r$ at $\infty$, $\circ(r_{\infty})$, we mean the order of
$\infty$ as a zero of $r$. That is $\circ(r_{\infty})=\deg(t)-\deg(s)$.

### The four cases

Case 1. In this case $[\sqrt{r}]_{c}$ and $[\sqrt{r}]_{\infty}$ means the
Laurent series of $\sqrt{r}$ at $c$ and the Laurent series of $\sqrt{r}$ at
$\infty$ respectively. Furthermore, we define $\varepsilon(p)$ as follows: if
$p\in\Gamma$, then $\varepsilon(p)\in\\{+,-\\}$. Finally, the complex numbers
$\alpha_{c}^{+}$, $\alpha_{c}^{-}$, $\alpha_{\infty}^{+}$,
$\alpha_{\infty}^{-}$ will be defined in the first step. If the differential
equation has no poles it only can fall in this case.

Step 1. Search for each $c\in\Gamma^{\prime}$ and for $\infty$ the
corresponding situation as follows:

  1. $(c_{0})$

If $\circ(r_{c})=0$, then

$\displaystyle[\sqrt{r}]_{c}=0,\qquad\alpha_{c}^{\pm}=0.$

  2. $(c_{1})$

If $\circ(r_{c})=1$, then

$\displaystyle[\sqrt{r}]_{c}=0,\qquad\alpha_{c}^{\pm}=1.$

  3. $(c_{2})$

If $\circ(r_{c})=2$, and

$\displaystyle
r=\cdots+b(x-c)^{-2}+\cdots,\qquad\textrm{then}\qquad[\sqrt{r}]_{c}=0,\qquad\alpha_{c}^{\pm}=\frac{1\pm\sqrt{1+4b}}{2}.$

  4. $(c_{3})$

If $\circ(r_{c})=2v\geq 4$, and

$\displaystyle
r=\big{(}a(x-c)^{-v}+\dots+d(x-c)^{-2}\big{)}^{2}+b(x-c)^{-(v+1)}+\cdots,\qquad\textrm{then}$
$\displaystyle[\sqrt{r}]_{c}=a(x-c)^{-v}+\dots+d(x-c)^{-2},\qquad\alpha_{c}^{\pm}=\frac{1}{2}\left(\pm\frac{b}{a}+v\right).$

  5. $(\infty_{1})$

If $\circ(r_{\infty})>2$, then

$\displaystyle[\sqrt{r}]_{\infty}=0,\qquad\alpha_{\infty}^{+}=0,\qquad\alpha_{\infty}^{-}=1.$

  6. $(\infty_{2})$

If $\circ(r_{\infty})=2$, and $r=\cdots+bx^{2}+\cdots$, then

$\displaystyle[\sqrt{r}]_{\infty}=0,\qquad\alpha_{\infty}^{\pm}=\frac{1\pm\sqrt{1+4b}}{2}.$

  7. $(\infty_{3})$

If $\circ(r_{\infty})=-2v\leq 0$, and

$\displaystyle
r=\big{(}ax^{v}+\dots+d\big{)}^{2}+bx^{v-1}+\cdots,\qquad\textrm{then}$
$\displaystyle[\sqrt{r}]_{\infty}=ax^{v}+\dots+d,\qquad\text{and}\qquad\alpha_{\infty}^{\pm}=\frac{1}{2}\left(\pm\frac{b}{a}-v\right).$

Step 2. Find $D\neq\varnothing$ defined by

$\displaystyle D=\left\\{n\in\mathbb{Z}_{+}\colon
n=\alpha_{\infty}^{\varepsilon(\infty)}-\sum\limits_{c\in\Gamma^{\prime}}\alpha_{c}^{\varepsilon(c)},\,\forall\,(\varepsilon(p))_{p\in\Gamma}\right\\}.$

If $D=\varnothing$, then we should start with the case 2. Now, if
$\mathrm{Card}(D)>0$, then for each $n\in D$ we search $\omega$
$\in\mathbb{C}(x)$ such that

$\displaystyle\omega=\varepsilon(\infty)[\sqrt{r}]_{\infty}+\sum\limits_{c\in\Gamma^{\prime}}\big{(}\varepsilon(c)[\sqrt{r}]_{c}+{\alpha_{c}^{\varepsilon(c)}}{(x-c)^{-1}}\big{)}.$

Step 3. For each $n\in D$, search for a monic polynomial $P_{n}$ of degree $n$
with

$\displaystyle\partial_{x}^{2}P_{n}+2\omega\partial_{x}P_{n}+\big{(}\partial_{x}\omega+\omega^{2}-r\big{)}P_{n}=0.$

If success is achieved then $\zeta_{1}=P_{n}{\rm e}^{\int\omega}$ is a
solution of the differential equation. Else, case 1 cannot hold.

Case 2. Search for each $c\in\Gamma^{\prime}$ and for $\infty$ the
corresponding situation as follows:

Step 1. Search for each $c\in\Gamma^{\prime}$ and $\infty$ the sets
$E_{c}\neq\varnothing$ and $E_{\infty}\neq\varnothing$. For each
$c\in\Gamma^{\prime}$ and for $\infty$ we define $E_{c}\subset\mathbb{Z}$ and
$E_{\infty}\subset\mathbb{Z}$ as follows:

  1. ($c_{1}$)

If $\circ(r_{c})=1$, then $E_{c}=\\{4\\}$.

  2. ($c_{2}$)

If $\circ(r_{c})=2$, and $r=\cdots+b(x-c)^{-2}+\cdots$, then
$E_{c}=\big{\\{}2+k\sqrt{1+4b}\colon k=0,\pm 2\big{\\}}$.

  3. ($c_{3}$)

If $\circ(r_{c})=v>2$, then $E_{c}=\\{v\\}$.

  4. $(\infty_{1})$

If $\circ(r_{\infty})>2$, then $E_{\infty}=\\{0,2,4\\}$.

  5. $(\infty_{2})$

If $\circ(r_{\infty})=2$, and $r=\cdots+bx^{2}+\cdots$, then
$E_{\infty}=\big{\\{}2+k\sqrt{1+4b}\colon k=0,\pm 2\big{\\}}$.

  6. $(\infty_{3})$

If $\circ(r_{\infty})=v<2$, then $E_{\infty}=\\{v\\}$.

Step 2. Find $D\neq\varnothing$ defined by

$\displaystyle D=\left\\{n\in\mathbb{Z}_{+}\colon
n=\frac{1}{2}\left(e_{\infty}-\sum\limits_{c\in\Gamma^{\prime}}e_{c}\right),\,\forall\,e_{p}\in
E_{p},\,p\in\Gamma\right\\}.$

If $D=\varnothing$, then we should start the case 3. Now, if
$\mathrm{Card}(D)>0$, then for each $n\in D$ we search a rational function
$\theta$ defined by

$\displaystyle\theta=\frac{1}{2}\sum\limits_{c\in\Gamma^{\prime}}\frac{e_{c}}{x-c}.$

Step 3. For each $n\in D$, search a monic polynomial $P_{n}$ of degree $n$,
such that

$\displaystyle\partial_{x}^{3}P_{n}+3\theta\partial_{x}^{2}P_{n}+\big{(}3\partial_{x}\theta+3\theta^{2}-4r\big{)}\partial_{x}P_{n}+\big{(}\partial_{x}^{2}\theta+3\theta\partial_{x}\theta+\theta^{3}-4r\theta-2\partial_{x}r\big{)}P_{n}=0.$

If $P_{n}$ does not exist, then case 2 cannot hold. If such a polynomial is
found, set $\phi=\theta+\partial_{x}P_{n}/P_{n}$ and let $\omega$ be a
solution of

$\displaystyle\omega^{2}+\phi\omega+{1\over
2}\big{(}\partial_{x}\phi+\phi^{2}-2r\big{)}=0.$

Then $\zeta_{1}={\rm e}^{\int\omega}$ is a solution of the differential
equation.

Case 3. Search for each $c\in\Gamma^{\prime}$ and for $\infty$ the
corresponding situation as follows:

Step 1. Search for each $c\in\Gamma^{\prime}$ and $\infty$ the sets
$E_{c}\neq\varnothing$ and $E_{\infty}\neq\varnothing$. For each
$c\in\Gamma^{\prime}$ and for $\infty$ we define $E_{c}\subset\mathbb{Z}$ and
$E_{\infty}\subset\mathbb{Z}$ as follows:

  1. $(c_{1})$

If $\circ(r_{c})=1$, then $E_{c}=\\{12\\}$.

  2. $(c_{2})$

If $\circ(r_{c})=2$, and $r=\cdots+b(x-c)^{-2}+\cdots$, then

$\displaystyle E_{c}=\big{\\{}6+k\sqrt{1+4b}\colon k=0,\pm 1,\pm 2,\pm 3,\pm
4,\pm 5,\pm 6\big{\\}}.$

  3. $(\infty)$

If $\circ(r_{\infty})=v\geq 2$, and $r=\cdots+bx^{2}+\cdots$, then

$\displaystyle E_{\infty}=\left\\{6+{12k\over m}\sqrt{1+4b}\colon k=0,\pm
1,\pm 2,\pm 3,\pm 4,\pm 5,\pm 6\right\\},\qquad m\in\\{4,6,12\\}.$

Step 2. Find $D\neq\varnothing$ defined by

$\displaystyle D=\left\\{n\in\mathbb{Z}_{+}\colon
n=\frac{m}{12}\left(e_{\infty}-\sum\limits_{c\in\Gamma^{\prime}}e_{c}\right),\,\forall\,e_{p}\in
E_{p},\,p\in\Gamma\right\\}.$

In this case we start with $m=4$ to obtain the solution, afterwards $m=6$ and
finally $m=12$. If $D=\varnothing$, then the differential equation is not
integrable because it falls in the case 4. Now, if $\mathrm{Card}(D)>0$, then
for each $n\in D$ with its respective $m$, search a rational function

$\displaystyle\theta={m\over
12}\sum\limits_{c\in\Gamma^{\prime}}\frac{e_{c}}{x-c}$

and a polynomial $S$ defined as

$\displaystyle S=\prod\limits_{c\in\Gamma^{\prime}}(x-c).$

Step 3. Search for each $n\in D$, with its respective $m$, a monic polynomial
$P_{n}=P$ of degree $n$, such that its coefficients can be determined
recursively by

$\displaystyle P_{-1}=0,\qquad P_{m}=-P,$ $\displaystyle
P_{i-1}=-S\partial_{x}P_{i}-((m-i)\partial_{x}S-S\theta)P_{i}-(m-i)(i+1)S^{2}rP_{i+1},$

where $i\in\\{0,1,\dots,m-1,m\\}$. If $P$ does not exist, then the
differential equation is not integrable because it falls in case 4. Now, if
$P$ exists search $\omega$ such that

$\displaystyle\sum\limits_{i=0}^{m}\frac{S^{i}P}{(m-i)!}\omega^{i}=0,$

then a solution of the differential equation is given by

$\displaystyle\zeta={\rm e}^{\int\omega},$

where $\omega$ is solution of the previous polynomial of degree $m$.

### Acknowledgements

P.A.-H. thanks to Universidad Simón Bolívar, Research Project _Métodos
Algebraicos y Combinatorios en Sistemas Dinámicos y Física Matemática_. He
also acknowledges and thanks the support of COLCIENCIAS through grant numbers
FP44842-013-2018 of the Fondo Nacional de Financiamiento para la Ciencia, la
Tecnología y la Innovación. E.T. wishes to thank the German Service of
Academic Exchange (DAAD) for financial support, and Professor M. Reuter at the
Institute of Physics in Uni-Mainz for stimulating discussions about this work.
Finally, the authors thank to the anonymous referees for their valuable
comments and suggestions.

## References

  * [1] Acosta-Humánez P.B., Galoisian approach to supersymmetric quantum mechanics, Ph.D. Thesis, Universitat Politècnica de Catalunya, Barcelona, 2009, arXiv:0906.3532.
  * [2] Acosta-Humánez P.B., Nonautonomous Hamiltonian systems and Morales–Ramis theory. I. The case $\ddot{x}=f(x,t)$, SIAM J. Appl. Dyn. Syst. 8 (2009), 279–297, arXiv:0808.3028.
  * [3] Acosta-Humánez P.B., Galoisian approach to supersymmetric quantum mechanics. The integrability analysis of the Schrödinger equation by means of differential Galois theory, VDM Verlag, Dr. Müller, Berlin, 2010.
  * [4] Acosta-Humánez P.B., Alvarez-Ramírez M., Blázquez-Sanz D., Delgado J., Non-integrability criterium for normal variational equations around an integrable subsystem and an example: the Wilberforce spring-pendulum, Discrete Contin. Dyn. Syst. Ser. A 33 (2013), 965–986, arXiv:1104.0312.
  * [5] Acosta-Humánez P.B., Álvarez Ramírez M., Delgado J., Non-integrability of some few body problems in two degrees of freedom, Qual. Theory Dyn. Syst. 8 (2009), 209–239, arXiv:0811.2638.
  * [6] Acosta-Humánez P.B., Blázquez-Sanz D., Hamiltonian system and variational equations with polynomial coefficients, in Dynamic Systems and Applications, Vol. 5, Dynamic, Atlanta, GA, 2008, 6–10.
  * [7] Acosta-Humánez P.B., Blázquez-Sanz D., Non-integrability of some Hamiltonians with rational potentials, Discrete Contin. Dyn. Syst. Ser. B 10 (2008), 265–293, math-ph/0610010.
  * [8] Acosta-Humánez P.B., Blazquez-Sanz D., Vargas-Contreras C.A., On Hamiltonian potentials with quartic polynomial normal variational equations, Nonlinear Stud. 16 (2009), 299–313, arXiv:0809.0135.
  * [9] Acosta-Humánez P.B., Kryuchkov S.I., Suazo E., Suslov S.K., Degenerate parametric amplification of squeezed photons: explicit solutions, statistics, means and variance, J. Nonlinear Opt. Phys. Mater. 24 (2015), 1550021, 27 pages, arXiv:1311.2479.
  * [10] Acosta-Humánez P.B., Lázaro J.T., Morales-Ruiz J.J., Pantazi C., On the integrability of polynomial vector fields in the plane by means of Picard–Vessiot theory, Discrete Contin. Dyn. Syst. Ser. A 35 (2015), 1767–1800.
  * [11] Acosta-Humánez P.B., Morales Ruiz J.J., Weil J.A., Galoisian approach to integrability of Schrödinger equation, Rep. Math. Phys. 67 (2011), 305–374, arXiv:1008.3445.
  * [12] Acosta-Humánez P.B., Pantazi C., Darboux integrals for Schrödinger planar vector fields via Darboux transformations, SIGMA 8 (2012), 043, 26 pages, arXiv:1111.0120.
  * [13] Acosta-Humánez P.B., Suazo E., Liouvillian propagators, Riccati equation and differential Galois theory, J. Phys. A: Math. Theor. 46 (2013), 455203, 17 pages, arXiv:1304.5698.
  * [14] Acosta-Humánez P.B., Suazo E., Liouvillian propagators and degenerate parametric amplification with time-dependent pump amplitude and phase, in Analysis, Modelling, Optimization, and Numerical Techniques, Springer Proc. Math. Stat., Vol. 121, Springer, Cham, 2015, 295–307.
  * [15] Braverman A., Etingof P., Gaitsgory D., Quantum integrable systems and differential Galois theory, Transform. Groups 2 (1997), 31–56, alg-geom/9607012.
  * [16] Cohen-Tannoudji C., Diu B., Lalöe F., Quantum mechanics, Vol. 1, John Wiley & Sons, New York, 1977.
  * [17] Cooper F., Khare A., Sukhatme U., Supersymmetry and quantum mechanics, Phys. Rep. 251 (1995), 267–385, hep-th/9405029.
  * [18] De Boer J., Quantum theory of condensed permanent gases. I. The law of corresponding states, Physica 14 (1948), 139–148.
  * [19] De Boer J., Michels A., Contribution to the quantum-mechanical theory of the equation of state and the law of corresponding states. Determination of the law of force of helium, Physica 5 (1938), 945–957.
  * [20] Duval A., Loday-Richaud M., Kovačič’s algorithm and its application to some families of special functions, Appl. Algebra Engrg. Comm. Comput. 3 (1992), 211–246.
  * [21] Gangopadhyaya A., Mallow J.V., Rasinariu C., Supersymmetric quantum mechanics. An introduction, World Scientific, Singapore, 2011.
  * [22] Gozzi E., Reuter M., Thacker W.D., Symmetries of the classical path integral on a generalized phase-space manifold, Phys. Rev. D 46 (1992), 757–765.
  * [23] Guggenheim E.A., The principle of corresponding states, J. Chem. Phys. 13 (1945), 253–261.
  * [24] Hansen J.-P., Verlet L., Phase transition on the Lennard-Jones system, Phys. Rev. 184 (1969), 151–161.
  * [25] Horváth G., Kawazoe K., Method for the calculation of effective pore size distribution in molecular sieve carbon, J. Chem. Eng. Japan 16 (1983), 470–475.
  * [26] Hurley A.C., Lennard-Jones J.E., Pople J.A., The molecular orbital theory of chemical valency. XVI. A theory of paired-electrons in polyatomic molecules, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 220 (1953), 446–455.
  * [27] Jorgensen W.L., Transferable intermolecular potential functions for water, alcohols and ethers. Application to liquid water, J. Amer. Chem. Soc. 103 (1981), 335–340.
  * [28] Keller J.B., Zumino B., Determination of intermolecular potentials from thermodynamic data and the law of corresponding states, J. Chem. Phys. 30 (1959), 1351–1353.
  * [29] Kolchin E.R., Differential algebra and algebraic groups, Pure and Applied Mathematics, Vol. 54, Academic Press, New York – London, 1973.
  * [30] Kovacic J.J., An algorithm for solving second order linear homogeneous differential equations, J. Symbolic Comput. 2 (1986), 3–43.
  * [31] Landau D.P., Binder K., A guide to Monte Carlo simulations in statistical physics, 3rd ed., Cambridge University Press, Cambridge, 2009.
  * [32] Lennard-Jones J.E., On the determination of molecular fields. I. From the variation of the viscosity of a gas with temperature, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 106 (1924), 441–462.
  * [33] Lennard-Jones J.E., On the determination of molecular fields. II. From the equation of state of a gas, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 106 (1924), 463–477.
  * [34] Lennard-Jones J.E., Cohesion, Proc. Phys. Soc. 43 (1931), 461–482.
  * [35] Lennard-Jones J.E., Devonshire A.F., Critical phenomena in gases - I, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 163 (1937), 53–70.
  * [36] McQuarrie D.A., Statistical mechanics, University Science Books, Sausalito, 2000.
  * [37] Mecke M., Winkelmann J., Fischer J., Molecular dynamics simulation of the liquid-vapor interface: the Lennard-Jones fluid, J. Chem. Phys. 107 (1997), 9264–9270.
  * [38] Mie G., Zur kinetischen Theorie der einatomigen Körper, Ann. Phys. 11 (1903), 657–697.
  * [39] Miller M.D., Nosanow L.H., Parish L.J., Zero-temperature properties of matter and the quantum theorem of corresponding states. II. The liquid-to-gas phase transition for Fermi and Bose systems, Phys. Rev. B 15 (1977), 214–229.
  * [40] Morales Ruiz J.J., Differential Galois theory and non-integrability of Hamiltonian systems, Progress in Mathematics, Vol. 179, Birkhäuser Verlag, Basel, 1999.
  * [41] Mulero A., Cuadros F., Isosteric heat of adsorption for monolayers of Lennard-Jones fluids onto flat surfaces, Chem. Phys. 205 (1996), 379–388.
  * [42] Olivier J.P., Modeling physical adsorption on porous and nonporous solids using density functional theory, J. Porous Mater. 2 (1995), 9–17.
  * [43] Pade J., Exact scattering length for a potential of Lennard-Jones type, Eur. Phys. J. D 44 (2007), 345–350.
  * [44] Pitzer K.S., Corresponding states for perfect liquids, J. Chem. Phys. 7 (1939), 583–590.
  * [45] Ramis J.-P., Martinet J., Théorie de Galois différentielle et resommation, in Computer Algebra and Differential Equations, Comput. Math. Appl., Academic Press, London, 1990, 117–214.
  * [46] Semenov-Tian-Shansky M.A., Lax operators, Poisson groups, and the differential Galois theory, Theoret. and Math. Phys. 181 (2014), 1279–1301.
  * [47] Singer M.F., Liouvillian solutions of $n$th order homogeneous linear differential equations, Amer. J. Math. 103 (1981), 661–682.
  * [48] Storck S., Bretinger H., Maier W.F., Characterization of micro- and mesoporous solids by physisorption methods and pore-size analysis, Appl. Catalysis A 174 (1998), 137–146.
  * [49] van der Put M., Singer M.F., Galois theory of linear differential equations, Grundlehren der Mathematischen Wissenschaften, Vol. 328, Springer-Verlag, Berlin, 2003.
  * [50] van der Waals J.D., The law of corresponding states for different substances, Proc. Kon. Nederl. Akad. Wetenschappen 15 (1913), 971–981.
  * [51] Witten E., Dynamical breaking of supersymmetry, Nuclear Phys. B 188 (1981), 513–554.