Generalized Lennard-Jones Potentials, SUSYQM and Differential Galois Theory

In this paper we start with proving that the Schr\"odinger 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.


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ν−2)−ν, 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ν − 2) − ν. 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.

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 (| r 1 − r 2 |) is given by 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 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 where M is the total mass of the system, µ is called the reduced mass. The Hamiltonian in the new coordinates takes the form where p r and p R are the canonical momenta conjugated respectively to the coordinates r = | r 1 − r 2 | and R = | 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 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 µ and a position and momentum given by the relative coordinates r and 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 Since V (r) is a rational central potential, the eigenfunctions Ψ( r) are separable into radial and angular parts, the last one given by the spherical harmonics The differential equation of our interest corresponds to the radial part of (2.4) as follows 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 eff (r) ≡ l(l + 1)/ 2µr 2 + V (r) and leaving the second derivative in r on one side, we have 2µ 2 V eff (r) − E k,l u k,l (r) = d 2 dr 2 u k,l (r) (2.5) at this point we define a rescaled potential v eff (r) ≡ 2µ 2 V eff (r) and a similarly rescaled energy 2µ 2 E k,l ≡ ε k,l ; in this case equation (2.5) turns out to be v eff (r) − ε k,l u k,l (r) = d 2 dr 2 u k,l (r). (2.6) In this way it is natural to define a rescaled Hamiltonian as H ≡ − d 2 dr 2 + v(r) in order to recover (2.6): The case for l = 0 defines the non-effective potential, and is also of great interest for our study where we have simplified u k,l=0 and ε k,l=0 to u k and ε k , respectively. We observe that (2.8) is a rescaled version of equations (2.7), (2.8) and (2.9) are the subject of our mathematical and physical analysis in Section 3.

The 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 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 where is the atomic depth of the potential well, σ 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, σ > 0 and > 0 satisfy that V (σ) = 0 and V 6 √ 2σ = − ; this means that σ is a zero potential length and the point 6 √ 2σ, − is the local minimum of the potential in the interval (0, ∞). It can easily be shown that there is no other critical point in such interval. In physical terms  the well depth and the zero potential length σ are parameters that describe the cohesive and repulsive forces that take place in a gas or liquid at the molecular level. measures the strength of the attraction between pairs of molecules and σ is the radius of the repulsive core when two molecules collide.
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 ν and δ given in [34] where 0 < ν < δ, A > 0, B > 0, C ≥ 0. Its rescaled version is given by The special case for δ = 2ν − 2 and some of its analytic advantages has been studied by J. Pade in [43] In the mentioned article, a special attention has been drawn to the ν = 6 case, and its ability to fit experimental data: (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 ρ: 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 ρ 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 pairpotential v(r) given by 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 δ . 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 σ and in order to construct dimensionless lengths and energies [18] r As a result, the radial Schrödinger equation (2.9) for l = 0 can be transformed into the dimensionless form given by provided that the potential V (r) can be expressed in the generic form V (r) = f (r/σ), where f (r) is a well-defined dimensionless interaction function and Λ ≡ /(σ √ µ ) [18]. From (2.17) we see that Λ, 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 Λ 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 ± or supercharges such that they commute with H and satisfy the algebra The second relation means that Q ± are nilpotent operators. A usual representation of the algebra, given in equations (2.18) 20) and A ± are defined in terms of the derivative d dx and an arbitrary complex function w(r), called the superpotential then, from (2.20) it results natural to identify which leads directly to a definition of the so-called partner potentials v ± given by such that Each of the two equations in (2.22) define a Riccati differential equation for the superpotential w, if v ± are known. Let's recall that the superpotential can also be found from the zero-energy base state ψ 0 , by computing w = −ψ 0 /ψ 0 , where ψ 0 is a solution of the Schrödinger equation with the v − potential (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 In the same way as we are interested in finding the roots of the polynomials over a base field, usually Q, using arithmetical and algebraic conditions, we would like to have explicit solutions of differential equations over a differential base field K = C(x), with field of constants C K = C, using elementary functions and quadratures. The differential Galois theory considers more general differential fields, but for our purpose is enough to consider C(x). Thus, the differential Galois group (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 y 1 , y 2 , . . . , y n is a basis of solutions of d n y dx n + a n−1 then for each differential automorphism σ ∈ DGal(L/K) there exists a matrix A σ ∈ GL(n, C) (i.e., a ij ∈ C, 1 ≤ i, j ≤ n and det(A σ ) = 0) such that 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 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 secondorder 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 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).
The Bessel's equation is a particular case of the confluent hypergeometric equation and is given by for some V ∈ K.
Furthermore, the algebraic form of the equation 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.

Main results
In this section we present the main contributions of this paper. First we will show that for the usual ν = 6, δ = 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 δ = 10 and ν = 6 we show the integrability, in the sense of differential Galois theory, as a special case of a general theorem for δ = 2ν − 2 with δ, ν ∈ 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).
We start by considering the radial Schrödinger equation (2.7) with the generalized effective Lennard-Jones potential (2.13) Setting C(r) as the differential field of equation (3.1) with the derivative d dr , we set alsō A,B,C ∈ 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 ∈ C * , C ∈ C.
Proof . Considering ν = 6 and δ = 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 Now, the change of dependent variable After applying Kovacic's algorithm, see Appendix A, we observe that equation (

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 secondorder 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 (  as a consequence we have A simple choice for w(r) is given by where the following identities should hold As a result we identify v δ−ν (r) in (2.12) with v − and we have from (3.3), the following expressions for the partner potentials According to equation (2.23), the corresponding wave function for the zero energy level using v − is given by 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 δ = 2(ν − 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 10 − 6 Lennard-Jones superpotential and the De Boer parameter As mentioned before the case for ν = 6, δ = 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 σ and is given by where α is chosen so that is the minimum energy (the well depth) and σ, as mentioned before, is the value where V 10−6 vanishes. As a result we have for this case 1 α ≡ (25/6) 5/3. Thus the rescaled 10 − 6 Lennard-Jones potential reads where we have applied definitions in (3.10) and α ≡ (25/6) 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 Λ ≡ /(σ √ µ ), which gives a degree of the quantum character of the system [18], we have Λ 2 ≈ 0.4303 or similarly Λ ≈ 0.6559. An important remark at this point is that the SUSY condition given in the formĀ = 5 √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 ± (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 10 − 6 Lennard-Jones gas As mentioned in Section 2, the De Boer principle of corresponding states tells us that two different systems with equal value of Λ have identical thermodynamical properties [18]. In this sense the SUSY condition Λ 2 = 2 / (σ) 2 µ = (1/3) 5/3 ≈ 0.4303 given in (3.12); and representing a definite set of combinations of values of the parameters σ, , and µ; that accounts for Λ 2 = (1/3) 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) ≡ 2 2µ v − (r) in (3.6) with ν = 6. We have found after a brief review of the literature, a significant coincidence between the specific value for Λ 2 ≈ 0.4303 and the value of Λ 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.

Integrability of the 10 − 6 Lennard-Jones potential and its generalization
The following result is valid for any potential v(r) belonging to a differential field. • The differential Galois groups of the radial equation and Schrödinger equation with effective potential are subgroups of SL(2, C).
• The transformation ϕ is strongly isogaloisian.  Figure 4. B 2 (T ) vs. T / plot for the 12 − 6 (black) and 10 − 6 (grey) Lennard-Jones potentials for the same molecular parameters σ and . 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.
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 ϕ : u → ϕ(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 SL(2, C), that is, σ(U ) = A σ U , σ(ϕ(U )) = A σ ϕ(U ) and det(A σ ) = 1. Thus, A σ ∈ SL(2, C).
• Applying the differential automorphism σ over ϕ(u) we observe that σ(ϕ(u)) = σ(r)σ(u) = rσ(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 ϕ 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.
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 Applying Martinet-Ramis theorem we have Assuming m ∈ Z we obtain which is the integrability condition for the Schrödinger equation with (2ν − 2) − ν L-J potential and its wave function corresponds to equation (3.7) with δ = 2ν − 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 ν = 6. Therefore the Schrödinger equation with L − J 10-6 is integrable for zero energy when A = ± √ B(8m + 4 ± 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 ∈ ±3 √ B, ±5 √ B .

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 Λ 2 ≈ 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 = 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.

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 (ζ 1 ), so that we can obtain the second solution (ζ 2 ) through Step 2. Find D = ∅ defined by D = n ∈ Z + : n = α ε(∞) ∞ − c∈Γ α ε(c) c , ∀ (ε(p)) p∈Γ .
If D = ∅, then we should start with the case 2. Now, if Card(D) > 0, then for each n ∈ D we search ω ∈ C(x) such that Step 3. For each n ∈ D, search for a monic polynomial P n of degree n with ∂ 2 x P n + 2ω∂ x P n + ∂ x ω + ω 2 − r P n = 0.
If success is achieved then ζ 1 = P n e ω is a solution of the differential equation. Else, case 1 cannot hold. Case 2. Search for each c ∈ Γ and for ∞ the corresponding situation as follows: Step 1. Search for each c ∈ Γ and ∞ the sets E c = ∅ and E ∞ = ∅. For each c ∈ Γ and for ∞ we define E c ⊂ Z and E ∞ ⊂ Z as follows: (c 1 ) If •(r c ) = 1, then E c = {4}.
If D = ∅, then we should start the case 3. Now, if Card(D) > 0, then for each n ∈ D we search a rational function θ defined by θ = 1 2 c∈Γ e c x − c .
If P n does not exist, then case 2 cannot hold. If such a polynomial is found, set φ = θ + ∂ x P n /P n and let ω be a solution of Then ζ 1 = e ω is a solution of the differential equation. Case 3. Search for each c ∈ Γ and for ∞ the corresponding situation as follows: Step 1. Search for each c ∈ Γ and ∞ the sets E c = ∅ and E ∞ = ∅. For each c ∈ Γ and for ∞ we define E c ⊂ Z and E ∞ ⊂ Z as follows: (c 1 ) If •(r c ) = 1, then E c = {12}.
Step 2. Step 3. Search for each n ∈ D, with its respective m, a monic polynomial P n = P of degree n, such that its coefficients can be determined recursively by P −1 = 0, P m = −P, where i ∈ {0, 1, . . . , 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 ω such that m i=0 S i P (m − i)! ω i = 0, then a solution of the differential equation is given by where ω is solution of the previous polynomial of degree m.