Jan Helm
Department of Electrical Engineering, Technical University Berlin, Berlin, Germany.
DOI: 10.4236/jmp.2025.1610073   PDF    HTML   XML   56 Downloads   464 Views   Citations

Abstract

This paper presents the mathematical models and calculation of states of matter dynamics (gasses, fluids, solids) and phase transitions from statistical viewpoint, with new calculation methods and results. We introduce here a new ansatz for the partition function, based on the generalized Landau theory, and apply it to two intermolecular potentials: Lennard-Jones and dipole-dipole. For these two potentials, we derive the basic thermodynamic variables, and compare the results with experimental data for Lennard-Jones substance (fluid argon) and for dipole substance (ethanol). This approach allows to derive the thermodynamic properties solely from its intermolecular potential. The paper introduces two novel methods. (1) Landau theory applied to solid-fluid-gas (fgs) systems; (2) Analytic expression for the partition function. The paper is structured as follows. Chapters 2 and 3 introduce the basics of phase transitions, and partition function. Chapter 4 describes the generalized Landau theory and its classical formulation for magnetic systems. Chapter 5 gives a detailed description of the van-der-Waals theory. Chapters 6 and 7 present the two representative substances with model and phase transitions: the Lennard-Jones substance and the dipole substance. Chapters 8 and 9 introduce the ansatz and the formulation of the partition function for fgs systems. Chapters 10 and 11 are the calculation of the radial distribution function and the partition function for the two chosen substances. Chapters 12 and 13 present the results: calculated pressure profiles, equation-of-state and characteristic points for the two substances.

Keywords

Phase Transition, Partition Function, Equation-of-State, Radial Distribution Function, Landau Theory, Intermolecular Potential, Van-der-Waals Equation, Lennard-Jones Substance, Dipole Substance

Share and Cite:

1. Introduction

In statistical mechanics, the key entity is the partition function Z, from which all important entities of a thermal solid-fluid-gas (fgs) system can be derived: free energy F, Gibbs energy G, entropy S, and pressure p in dependence of the thermal basic variables N, V, T, where the pressure determines the equation-of-state (eos).

In magnetic systems, the magnetization M plays the role of pressure, and the basic variables are N, B, T, where B is the magnetic field flux.

The phase transitions in fgs systems are described by saturation curves (e.g. fluid-gas saturation curve) characterized by derivatives of pressure ($\frac{\partial p}{\partial V}=0$ ) significant points are critical point ($\frac{{\partial }^{2}p}{\partial V^2}=0$ ) and multiple points (e.g. triple point solid-fluid-gas) i.e. branching points. Here, phase transitions are first-order transitions, i.e. first derivative of Gibbs energy, volume $V={\left(\frac{\partial G}{\partial p}\right)}_T$ , is discontinuous at phase transitions.

The phase transitions in magnetic systems are characterized by derivatives of magnetization ($\frac{\partial M}{\partial B}=0$ ).

Here, phase transitions are second-order transitions, i.e. second derivative of free (or Gibbs) energy, susceptibility ${\chi }_T={\left(\frac{{\partial }^{2}F}{\partial B^2}\right)}_T$ , is discontinuous.

It is generally accepted that the behavior of a thermal fgs system is solely determined by its inter-molecular potential (e.g. Lennard-Jones potential in fluid argon, dipole-dipole OH-binding in ethanol), $Z={\displaystyle \int {\phi }_{rdf}dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$.

The behavior of a thermal magnetic system is determined by its Landau function $\phi =M$ magnetization, and its Landau energy

$E\left(\phi ,\beta \right)=\frac{\mu \left(\beta \right)}{2}{\phi }^2+\frac{\lambda \left(\beta \right)}{4!}{\phi }^4-\beta \phi B$ .

In generalized form, the partition function Z is described by the functional integral $Z={\displaystyle \int D\phi dV\exp \left(-H\left(\overrightarrow{r},\phi ,\beta \right)\right)}$, where $H\left(\overrightarrow{r},\phi ,\beta \right)$ is the thermodynamic Hamiltonian $H\left(\overrightarrow{r},\phi ,\beta \right)=E(\phi )+\beta u\left(\overrightarrow{r},\phi \right)$.

$D\phi ={\displaystyle \sum _i\frac{\partial \phi }{\partial c_i}}d{c}_i$ for Landau function with parameters $\phi \left(c_i\right)={\displaystyle \sum _i\frac{\partial \phi }{\partial c_i}}d{c}_i$ (fgs systems).

$Z={\displaystyle \int D\phi dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$, approximately $Z\approx {\displaystyle \int {\phi }_{rdf}dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$.

$D\phi =d\phi$ for simple Landau function without parameters (magnetic systems)

$Z={\displaystyle \int d\phi dV\exp \left(-E\left(\phi \right)\right)}$

The goal in statistical mechanics is to calculate the phase transitions and their significant points analytically from the partition function, and to express thermal functions behavior at critical points $T_C$ in asymptotic form $f\left(T-T_C\right)={\left(T-T_C\right)}^{\kappa }$ , where κ is a critical exponent, e.g.

$M~{\left(T_C-T\right)}^{\beta },T\to T_C-,B=0$

$p~{\left(T_C-T\right)}^{\beta },T\to T_C-$

This is at present possible only for

• magnetic systems

with Landau energy $E\left(\phi ,\beta \right)=\frac{\mu \left(\beta \right)}{2}{\phi }^2+\frac{\lambda \left(\beta \right)}{4!}{\phi }^4-\beta \phi B$

• vdWaals eos of a real gas

with eos $p\left(v-v_0\right)\beta =\left(1-\frac{\beta }{2\left(v-v_0\right)}{\displaystyle \underset{r=\sigma }{\overset{\infty }{\int }}\left|u\left(r\right)\right|d^{3}r}\right)$ , where $v_0=\frac{2\pi }{3}{\sigma }^3$ is the

hard-core volume= vdWaals b-parameter, $\sigma$ is the hard-core radius.

In this paper, we formulate a general Landau theory for thermal systems, where the original Landau ansatz for magnetic systems and an ansatz with rdf function for fgs systems are special cases. The fgs ansatz is applied to two intermolecular potentials: Lennard-Jones and dipole-dipole.

For these two potentials, we derive the basic thermodynamic variables: partition function, free energy, pressure, then equation-of-state, saturation curve, characteristic points, and compare the results with experimental data.

This approach allows to derive the thermodynamic properties solely from its intermolecular potential, and it is carried out here in good agreement with material data for Lennard-Jones substance (fluid argon) and for dipole substance (ethanol).

We introduce here two novel methods.

• Landau theory applied to solid-fluid-gas (fgs) systems

The Landau theory of magnetic systems is transferred to fgs systems, using the representation of partition function Z as an integral over intermolecular potential and radial-distribution function (rdf). This Landau ansatz yields the correct ideal-gas equation-of-state (eos) and the vdWaals eos. Furthermore, the calculated partition function gives the correct thermodynamic data (critical point, triple point, saturation curve, melting curve, eos) for the two chosen substances.

In particular, we show that the vdWaals model follows in weak-binding approximation from the generalized Landau theory for fgs systems.

We also show that the generalized Landau theory for fgs systems yields the correct eos, saturation curve and critical/triple point for Lennard-Jones potential (e.g. fluid argon) and for polar covalent binding (= dipole-dipole potential) (e.g. ethanol).

The calculation results are in good agreement with measured parameters for fluid argon resp. for ethanol.

• Analytic expression for the partition function

We introduce here a symbolic-numeric calculation method for the integral in the partition function Z.

In general, it is impossible to solve the integral in Z in closed form with symbolic parameters.

Numerical calculation on lattice in parameter space with a fit yields a closed expression, but is prohibitively costly in computing time for more than 2 parameters.

An alternative half-analytic approach, used here, is to calculate the integral symbolically as a weighted sum on lattice, using Clenshaw-Curtis quadrature.

This yields an analytic function as a closed expression with parameters, which is large (up to 200,000 terms) and memory-consuming, but can be evaluated fast point-wise. The function converges uniformly against the true integral with O(h²), i.e. quadratically in lattice unit size h.

2. Phase Transition Basics

2.1. Phase Transitions

In general, we distinguish two types of phase transitions: first-order and second-order (Ehrenfest classification), first-order transitions are discontinuous in the first derivative of free energy F, second-order are continuous in the first derivative and discontinuous in the second derivative ([1] chap. 1.1).

The derivatives are in an order variable of the system, in solid-fluid-gas PT’s it is usually volume V or density $\rho =\frac{N}{V}$ , in magnetic PT’s it is the magnetization M.

• First-order phase transitions

These apply mostly to the solid-fluid-gas transitions, where the discontinuous

variable is the volume $V=\frac{\partial G}{\partial p}$ , or equivalently the average distance λ.

Here the average distance λ, the density ρ, specific volume $v=\frac{V}{N}$ and heat capacity C_v are discontinuous.

Transitions happen at the melting temperature T_m, and at the boiling temperature T_b, there is a latent heat $\Delta C$ at the transition temperature.

• Second-order (continuous) phase transitions

These apply mostly to the phase transitions in magnetic materials (ferromagnetic, antiferromagnetic) and in superconducting materials.

In magnetic materials, the magnetization $M={\left(\frac{\partial F}{\partial B}\right)}_T$ is continuous at Curie temperature T_C: in ferromagnetic phase transition at T_C spontaneous magnetization $M={\left(\frac{\partial F}{\partial B}\right)}_T$ becomes 0, and there is no latent heat.

The magnetic susceptibility ${\chi }_T={\left(\frac{{\partial }^{2}F}{\partial B^2}\right)}_T$ is discontinuous at T_C.

In superconducting materials, the density n_s of superconducting electrons is continuous at critical temperature T_c, where it becomes zero. The coherence length at small perturbations $\xi$ is the analogue of the susceptibility: it is discontinuous at the critical temperature T_C.

Maxwell’s equal-area rule in first-order transitions

When $F\left(\lambda ,\beta \right)$ has a maximum in λ, then pressure

$p\left(\lambda ,\beta \right)=-\frac{1}{4\pi {\lambda }^2}\frac{\partial F\left(\lambda ,\beta \right)}{\partial \lambda }$ has negative-zero-positive transition, and a region of negative pressure is unstable, therefore there is a “jump” in λ over this area, i.e. a phase transition. The transition is determined by Maxwell’s equal-area rule $p\left(v_1,T\right)=p\left(v_2,T\right)$ and ${\displaystyle \underset{v_1}{\overset{v_2}{\int }}p\left(v,T\right)dv}=0$ (equivalently $F\left(v_1,T\right)=F\left(v_2,T\right)$ )

So we have the following behavior for phase transitions:

fluid-->gas at saturation curve:

${\lambda }_1={\lambda }_f\left(T\right)\to {\lambda }_2={\lambda }_g\left(T\right)$ determined by Maxwell’s equal-area rule

$F\left(v_1,T\right)=F\left(v_2,T\right)$ , $p\left(v_1,T\right)=p\left(v_2,T\right)$ where $p\left(v,T\right)=-\frac{\partial F\left(v,T\right)}{\partial v}$

fluid-->solid at fluid-solid line,

${\lambda }_1={\lambda }_f\left(T\right)\to {\lambda }_2={\lambda }_s\left(T\right)$ determined by Maxwell’s equal-area rule

$F\left(v_1,T\right)=F\left(v_2,T\right)$ , $p\left(v_1,T\right)=p\left(v_2,T\right)$ where $p\left(v,T\right)=-\frac{\partial F\left(v,T\right)}{\partial v}$ .

Heat capacity behavior in continuous magnetic phase transitions

Heat capacity C_v often diverges in the neighborhood of T_c as

$C~\left|T-T_c\right|{}^{-\alpha }$ , with $0<\alpha <1$ mostly $\alpha \ll 1$

sometimes $C~\log \left(\frac{T_c}{\left|T-T_c\right|}\right)$

in general $C\left(T\right)~\frac{1}{\alpha }\left({\left|\frac{T-T_c}{T_c}\right|}^{-\alpha }-1\right)$ , C finite, dC/dT need not be finite.

2.2. Binding Potential

• Ionic bonding: Coulomb potential

$U\left(r\right)=\frac{k{q}_1{q}_2}{r}$ (1)

k= 2.31 × 10⁻²⁸ J·m, q_i in e₀-units.

• Ionic bonding with free electron: screened Coulomb potential

For screened Coulomb potential

$U\left(r\right)=A\exp \left(-r/l_D\right)/r$ (2)

where l_D=Debye length.

• Covalent bonding

London potential between spheres with radii R₁, R₂

$U\left(z;R_1,R_2\right)=-\frac{A}{6}\left(\frac{2R_1{R}_2}{z^2-{\left(R_1+R_2\right)}^2}+\frac{2R_1{R}_2}{z^2-{\left(R_1-R_2\right)}^2}+\log \left(\frac{z^2-{\left(R_1+R_2\right)}^2}{z^2-{\left(R_1-R_2\right)}^2}\right)\right)$ (3)

$z=R_1+R_2+r$ ,$A={10}^{-20}{\mathrm{...10}}^{-19}\text{J}$

Lennard-Jones potential

$U_{LJ}\left(r,\sigma ,\epsilon \right)=4\epsilon \left({\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right)$ (4)

for argon ${\sigma }_{LJ}\left(Ar\right)=3.4\dot{\text{A}}$, $\epsilon \left(Ar\right)=0.0123\text{eV}$ , $p_0\left(Ar\right)=\frac{\epsilon }{{\sigma }^3}=689\text{MPa}$

Morse potential

$U_M\left(r,\sigma ,\epsilon \right)=\epsilon \left({\left(1-\exp \left(-\frac{r-\sigma }{\sigma }\right)\right)}^2-1\right)$ (5)

• Dipole potential

For arbitrarily positioned dipoles, with angle θ₁ resp. θ₂ to center connection line,

$U\left(r\right)=-\frac{e^2{r}_1{r}_2{q}_1{q}_2}{r^3}\left(-\cos \left({\theta }_1-{\theta }_2\right)+3\cos \left({\theta }_1\right)\cos \left({\theta }_2\right)\right)$ (6)

when both dipoles perpendicular to connection line, inverse direction:

${\theta }_1=\pi /2$ , ${\theta }_2=-\pi /2$ , $V\left(r\right)=-\frac{e^2{r}_1{r}_2{q}_1{q}_2}{r^3}$

with thermal screening

$U\left(r\right)=-\frac{2{\left(e^2{r}_1{r}_2{q}_1{q}_2\right)}^2}{3r^6}\frac{1}{k_{B}T}$ in cgs.

2.3. Phase Transition Solid-fluid-gas

These are the typical first-order phase transitions found in nature ([2]-[5]).

• p-T-phase diagram

The corresponding pressure-temperature phase diagram p(T) consists of three regions: solid, liquid, gas (Figure 1). The saturation curve between liquid and gas region starts with the branching triple point, and ends at the critical point, where the separation of fluid and gas ceases to exist.

Figure 1. A typical p-T-phase diagram. The solid green line applies to most substances; the dashed green line gives the anomalous behavior of water [6].

• p-V-phase diagram

The corresponding pressure-volume (isotherm) phase diagram p(T) describes the two-phase transition liquid-gas along the saturation curve (Figure 2).

In the two-phase region the volume “jumps” from the liquid to the gas phase at constant pressure. The isotherm part between the minimum and the following maximum is forbidden, since the compressibility there is negative.

The transition pressure is defined by the Maxwell-area-rule, as depicted schematically below.

In the Maxwell-area-rule, the area above and below the pressure is equal, i.e. $p\left(v_1,T\right)=p\left(v_2,T\right)$ , $F\left(v_1,T\right)=F\left(v_2,T\right)$ : pressure and free energy at transition points are equal.

Figure 2. Isotherms of a gas. The red line is the critical isotherm, with critical point K. The dashed lines represent parts of isotherms which are forbidden since the gradient would be positive, giving the gas in this region a negative compressibility [6]. Maxwell’s equal-area rule: the area above and below the transition pressure is equal.

2.4. Triple Point, Critical Point

2.4.1. Triple Point

At the triple point [6], the saturation curve ends, fluid-->gas and fluid-->solid curves meet,

i.e. $p\left(v_g,T\right)=p\left(v_f,T\right)=p\left(v_s,T\right)$

The triple point can be found as the branching point of saturation curve in the $p\left(\lambda ,\beta \right)$ diagram.

2.4.2. Critical Point

At the critical point [6], only one phase exists, the isotherm has a turning point there.

$\frac{\partial p}{\partial v}=0$ , $\frac{{\partial }^{2}p}{\partial v^2}=0$

At liquid-gas critical point:

compressibility diverges ${\kappa }_T={-\frac{1}{V}\frac{\partial V}{\partial p}|}_T~{\left|T-T_c\right|}^{-\gamma }$

specific volume jump goes to zero $v=\frac{V}{N}$ , $v_g-v_l~{\left|T-T_c\right|}^{\beta }$

Heat capacity diverges $C_v={-T\frac{{\partial }^{2}F}{\partial T^2}|}_V~{\left|T-T_c\right|}^{-\alpha }$

pressure derivative is discontinuous $\left|p-p_c\right|~{\left|V-V_c\right|}^{\delta }$

The critical exponents $\alpha ,\beta ,\gamma ,\delta$ obey the interdependence laws like in the magnetic phase transitions.

2.5. Definitions of Critical Exponents in Fluid, Magnetic Systems

Critical exponents describe the thermal behavior at characteristic points in form of power law with universal (=substance-independent) exponent.

The power law has the summary form $\left|\Delta f_d\right|~{\left|\Delta f_f\right|}^{\epsilon },T\to T_c$ , where $f_f$ is a fundamental variable temperature T (resp. $\beta$ ) or volume V (resp. $\lambda$ ), and $f_d$ is a dependent variable like pressure p, density ρ, compressibility κ_T, heat capacity C_V, considered at a characteristic point with temperature $T_c$ .

In case of second-ordered derivatives like heat capacity $C_v\equiv {\left(\frac{\partial U}{\partial T}\right)}_V=k_B{\beta }^2{\left(\frac{{\partial }^2\log Z}{\partial {\beta }^2}\right)}_V$ , the power law can be divergent i.e. ε < 0.

The universal character of critical exponents results from the fact, that the partition function can always be formulated in substance-independent and dimensionless way, using as units the fundamental energy ${\epsilon }_0$ and the molecular (hard-core) radius ${\sigma }_0$ from the intermolecular potential. Examples are the Lennard-Jones fluid and the vdWaals fluid.

We distinguish 6 critical exponents in fluid/magnetic systems, which characterize the behavior of thermodynamic variables near critical temperature T_c, correlation length ξ, volume V resp. magnetic field flux B, pressure p resp. magnetization M.

They are defined as follows in Table 1 ([1] chap.5).

Table 1. Definition of critical exponents.

exponent	definition			value
α, heat capacity CB, CV	$C_B~{\alpha }^{-1}\left({\left(\left|T-T_C\right|/T_C\right)}^{-\alpha }-1\right),T\to T_C,B=0$ $C_V~{\alpha }^{-1}\left({\left(\left|T-T_C\right|/T_C\right)}^{-\alpha }-1\right),T\to T_C$			$2-\nu d$
β, magnetization M(T), density ρ(T)	$M~{\left(T_C-T\right)}^{\beta },T\to T_C-,B=0$ $\left|\frac{\rho -{\rho }_c}{{\rho }_c}\right|~{\left|\frac{T-T_C}{T_C}\right|}^{\beta },T\to T_C$			$\left(d+\eta -2\right)\nu /2$
γ, susceptibility χT, compressibility κT	${\chi }_T~{\left|T-T_C\right|}^{-\gamma },T\to T_C,B=0$ ${\kappa }_T~{\left|T-T_C\right|}^{-\gamma },T\to T_C$			$\left(2-\eta \right)\nu$
δ, magnetization M(B), pressure p(T), density ρ(T)		$M~B^{1/\delta },T\to T_C,B=0$ $\left|\frac{p-p_c}{p_c}\right|~{\left|\frac{\rho -{\rho }_C}{{\rho }_C}\right|}^{\delta },T\to T_C$	$\frac{d+2-\eta }{d-2+\eta }$
η, correlation function G(2)		$G^{\left(2\right)}\left(r\right)~1/r^{d-2+\eta },T=T_C$
ν, correlation length ξ		$\xi ~{\left|T-T_C\right|}^{-\nu }$

Values of 6 critical exponents are shown in the following Table 2.

Table 2. Values of critical exponents.

	vdWaals	LenJones	Xe	binary fluid	β-brass	4He	Fe	Ni
d dim.	3	3	1	1	1	2	3	3
α		0.11	<0.2	0.113	0.05	-0.014	-0.03	0.04
β	0.5 exp 0.325	0.328	0.35	0.322	0.305	0.34	0.37	0.358
γ	1. 1.24	1.24	1.3	1.239	1.25	1.33	1.33	1.33
δ	3 exp 4.8	4.8	4.2 ± 0.6	4.58		3.95	4.3	4.29
η		0.034	0.1 ± 0.1	0.017 ± 0.015	0.08 ± 0.07	0.021 ± 0.05	0.07 ± 0.04	0.041 ± 0.01
ν		0.63	≈0.57	0.625	0.65	0.672	0.69	0.64

3. Partition Function, Radial Distribution Function

3.1. Radial Distribution Function

Correlation function or radial distribution function (rdf) g(r) is a measure of the probability that a particle will be located a distance r from another particle ([7]-[9]), it obeys the normalization condition

$4\pi {\displaystyle \underset{0}{\overset{\infty }{\int }}dr}{r}^{2}g\left(r\right)=N-1\approx N$ , ${\displaystyle \underset{0}{\overset{\infty }{\int }}dr}{r}^{2}g\left(r\right)=\frac{V}{4\pi }$ (7)

Given a potential $U\left({\overrightarrow{r}}_1,\cdots ,{\overrightarrow{r}}_N\right)$ , we obtain the Hamiltonian

$H={\displaystyle \sum _{i=1}^{3N}\frac{p_i^2}{2m}}+U\left({\overrightarrow{r}}_1,\cdots ,{\overrightarrow{r}}_N\right)$ (8)

The corresponding partition function reads

$\begin{array}{c}Z\left(N,V,T\right)=\frac{1}{N!h^{3N}}{\displaystyle \int d^{3N}p{d}^{3N}r\exp \left(-\beta {\displaystyle \sum _{i=1}^{3N}\frac{p_i^2}{2m}}\right)\exp \left(-\beta U\left({\overrightarrow{r}}_1,\cdots ,{\overrightarrow{r}}_N\right)\right)}\\ =\frac{1}{N!h^{3N}}{\displaystyle \int d{r}_1\cdots d{r}_N\exp \left(-\beta U\left({\overrightarrow{r}}_1,\cdots ,{\overrightarrow{r}}_N\right)\right)}=\frac{Z_N}{N!h^{3N}}\end{array}$ (9)

with thermal wavelength $\zeta =\sqrt{2\pi }\hslash c\sqrt{\frac{\beta }{m{c}^2}}=\sqrt{\frac{\beta h^2}{2\pi m}}$ ,

example:

$m\left({\text{O}}_{\text{2}}\right)=16\text{GeV}$ , $\beta \left(T=300\text{K}\right)=1/0.026\text{eV}$ , $\zeta =0.14\dot{\text{A}}$

For the ideal gas $U=0$

$Z\left(N,V,T\right)=\frac{1}{N!}{\left(\frac{V}{h^3}{\left(\frac{2\pi m}{\beta }\right)}^{3/2}\right)}^N=\frac{1}{N!}{\left(\frac{V}{{\zeta }^3}\right)}^N$ (10)

The pressure becomes $p=\frac{1}{\beta }\frac{\partial \log Z}{\partial V}$ ,

$\begin{array}{c}\log Z=\left(N\left(\log V-3\log \zeta \right)-\log N!\right)\\ \approx \left(N\left(\log V-3\log \zeta \right)-N\left(\log N-1\right)\right)\end{array}$

$p\beta =\frac{\partial \log Z}{\partial V}=\frac{\partial }{\partial V}\left(N\left(\log V-3\log \zeta \right)-\log N!\right)$

because of the Stirling formula $N!\approx \sqrt{2\pi N}{\left(\frac{N}{e}\right)}^N$

$\log N!=N\left(\log N-1\right)+\frac{1}{2}\log N+\frac{1}{2}\log \left(2\pi \right)\approx N\left(\log N-1\right)$ ,

$p\beta =\frac{\partial \log Z}{\partial V}=\frac{N}{V}=\frac{1}{v}=\rho$ , or $p\beta v=1$ , $p\beta =\rho$ (11)

where $v=\frac{V}{N}={\lambda }^3=\frac{1}{\rho }$ is the specific volume, and λ is the average distance.

We obtain for the energy of ideal gas

$E=-\frac{\partial \log Z}{\partial \beta }=\frac{3}{2}\frac{N}{\beta }$ (12b)

so the specific energy (= per particle) becomes

${\rho }_E=\frac{E}{N}=\frac{3}{2\beta }$ (12a)

Average energy with potential

With the configurational partition function

$Z_N={\displaystyle \int d^3{r}_1\cdots d^3{r}_N\exp \left(-\beta U\left({\overrightarrow{r}}_1,\cdots ,{\overrightarrow{r}}_N\right)\right)}$ (13a)

we obtain

$\log Z=\log Z_N-3\log \zeta -\log N!$ (13b)

So the energy becomes

$\begin{array}{c}E=-\frac{\partial \log Z}{\partial \beta }\\ =\frac{3N}{\zeta }\frac{\partial \zeta }{\partial \beta }-\frac{1}{Z_N}{\displaystyle \int d{r}_1\cdots d{r}_{N}U\left({\overrightarrow{r}}_1,\cdots ,{\overrightarrow{r}}_N\right)\exp \left(-\beta U\left({\overrightarrow{r}}_1,\cdots ,{\overrightarrow{r}}_N\right)\right)}\\ =\frac{3}{2}\frac{N}{\beta }+\langle U\rangle \end{array}$ (14)

where $\langle U\rangle =\frac{N^2}{2V}{\displaystyle \underset{0}{\overset{\infty }{\int }}dr4\pi r^{2}u\left(r\right)g\left(r\right)}$ (grand canonical ensemble)

$E=\frac{3}{2}\frac{N}{\beta }+2\pi \frac{N^2}{V}{\displaystyle {\int }_0^{\infty }dr{r}^{2}u\left(r\right)g\left(r\right)}$

so ${\rho }_E=\frac{E}{N}=\frac{3}{2\beta }+2\pi \rho {\displaystyle {\int }_0^{\infty }dr{r}^{2}u\left(r\right)}g\left(r\right)$ , $\frac{E}{\epsilon N}=\frac{3}{2}\frac{1}{\beta \epsilon }+2\pi \rho {\displaystyle {\int }_0^{\infty }dr{r}^2{u}_0\left(r\right)g\left(r,\beta \epsilon ,\rho \right)}$ (15)

where ε is the characteristic energy, and the dimensionless potential is

$u_0\left(r\right)=\frac{u\left(r\right)}{\epsilon }$

Partition function Z in terms of g(r)

$\frac{\log Z\left(N,V,\beta \right)}{N}=-\frac{1}{N}{\displaystyle \int d{\beta }_0\frac{E}{\epsilon }}=-{\displaystyle \int d{\beta }_0\left(3/\left(2{\beta }_0\right)+2\pi \frac{1}{V}{\displaystyle {\int }_0^{\infty }dr{r}^2{u}_0\left(r\right)g\left(r,\rho ,{\beta }_0\right)}\right)}$ (16)

with the dimensionless variable ${\beta }_0=\beta \epsilon$ , reformulated it becomes

$\frac{\log Z\left(\rho ,{\beta }_0\right)}{N}=-\frac{3}{2}\log {\beta }_0-2\pi \frac{1}{V}{\displaystyle \int d{\beta }_0{\displaystyle {\int }_0^{\infty }dr{r}^2{u}_0\left(r\right)g\left(r,\rho ,{\beta }_0\right)}}$ .

In general, the temperature dependence of the rdf is weak

$g\left(r,\rho ,{\beta }_0\right)\approx g\left(r,\rho \right)$ and the average density is constant ${\rho }_0=\frac{1}{v}=\frac{N}{V}$ , so we

obtain with

$\frac{\log Z\left(\rho ,\beta \epsilon \right)}{N}\approx -\frac{3}{2}\log \beta \epsilon -\frac{\beta \epsilon }{2V}{\displaystyle \int dV{u}_0\left(r\right)g\left(r\right){\rho }_0}=-\frac{3}{2}\log \beta \epsilon -\frac{\beta \epsilon }{2}{\displaystyle \int dv{u}_0\left(r\right)g\left(r\right)}$

The pressure becomes $p=\frac{1}{\beta }\frac{\partial \log Z}{\partial V}$ ([8] chap. 9)

$\frac{1}{Z_N}\frac{\partial \log Z_N}{\partial V}=\frac{N}{V}+\frac{\beta }{3V}\langle {\displaystyle \sum _i{\overrightarrow{r}}_i\cdot {\overrightarrow{F}}_i}\rangle$

with pair potential force ${\overrightarrow{F}}_{12}=-\frac{\partial U}{\partial {\overrightarrow{r}}_{12}}=-u^{\prime }\left(\left|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2\right|\right)\left(\frac{{\overrightarrow{r}}_1-{\overrightarrow{r}}_2}{\left|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2\right|}\right)=-u^{\prime }\left(r_{12}\right)\frac{{\overrightarrow{r}}_{12}}{r_{12}}$ follows

$\frac{\beta }{3V}\langle {\displaystyle \sum _i{\overrightarrow{r}}_i\cdot {\overrightarrow{F}}_i}\rangle =\frac{\beta N^2}{6V^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}dr4\pi r^3{u}^{\prime }\left(r\right)g\left(r\right)}$

The pressure becomes

$\frac{p}{kT}=\rho -\frac{{\rho }^2}{6kT}{\displaystyle \underset{0}{\overset{\infty }{\int }}dr4\pi r^3{u}^{\prime }\left(r\right)g\left(r;kT\right)}$ , $g=g\left(r,\rho ,kT\right)$ , $g\left(r;\rho ,kT\right)={\displaystyle \sum _{j=0}^{\infty }{\rho }^j{g}_j\left(r;kT\right)}$

reformulated

$p\beta =\rho -\frac{{\rho }^2\beta \epsilon }{6}{\displaystyle \underset{0}{\overset{\infty }{\int }}dr4\pi r^3{u^{\prime }}_0\left(r\right)g\left(r,\beta \epsilon ,\rho \right)}$ (17)

Compare: vdWaals eos reads $\frac{p}{kT}=\frac{\rho }{1-b\rho }-\frac{a{\rho }^2}{kT}$ , $v=\frac{V}{N}=\frac{1}{\rho }$ , $b=2\pi {\sigma }^3/3$ , resp. $p\beta =\frac{\rho }{1-b\rho }-a{\rho }^2\beta$ ,

from which follows for a-parameter $a\left(\rho ,kT\right)=\frac{1}{6}{\displaystyle \underset{0}{\overset{\infty }{\int }}dr4\pi r^3{u}^{\prime }\left(r\right)g\left(r,\rho ,kT\right)}$

$p\beta =\rho +{\displaystyle \sum _{j=0}^{\infty }{B}_{j+2}\left(T\right){\rho }^{j+2}}$ , with $B_{j+2}\left(T\right)=-\frac{1}{6kT}{\displaystyle \underset{0}{\overset{\infty }{\int }}dr4\pi r^3{u}^{\prime }\left(r\right)g_j\left(r;kT\right)}$

3.2. Ornstein-Zernike Equation

The Ornstein-Zernike equation calculates the direct correlation function c(r), which describes the pure correlation of a molecule with a neighboring molecule at distance ${\overrightarrow{r}}_{12}$ , whereas the total correlation function g(r) (= radial distribution function) takes into account also the correlation of the neighboring molecule with a third molecule at ${\overrightarrow{r}}_{23}$ , and all following molecules in the correlation chain.

The Ornstein-Zernike equation is an integral equation of the form [10]

$h\left(r_{12}\right)=c\left(r_{12}\right)+\rho {\displaystyle \int d^3{r}_3}c\left(r_{13}\right)h\left(r_{32}\right)$ , where $h\left(r_{12}\right)=g\left(r_{12}\right)-1$ (18)

with direct correlation function c(r) and total correlation function (=radial distribution function) g(r) and with Fourier-transforms

$\widehat{h}\left(\overrightarrow{k}\right)=\frac{\widehat{c}\left(\overrightarrow{k}\right)}{1-\rho \widehat{c}\left(\overrightarrow{k}\right)}$, $\widehat{c}\left(\overrightarrow{k}\right)=\frac{\widehat{h}\left(\overrightarrow{k}\right)}{1+\rho \widehat{h}\left(\overrightarrow{k}\right)}$

Hypernetted-chain equation

Hypernetted-chain equation (HC equation) is a closure relation to solve the Ornstein-Zernike equation (HCOZ equation) which relates the direct correlation function to the total correlation function.

$\log g\left({\overrightarrow{r}}_{12}\right)+\beta u\left({\overrightarrow{r}}_{12}\right)=\rho {\displaystyle \int \left(h\left({\overrightarrow{r}}_{13}\right)-\log g\left({\overrightarrow{r}}_{13}\right)-\beta u\left({\overrightarrow{r}}_{13}\right)\right)h\left({\overrightarrow{r}}_{23}\right)}{d}^3{\overrightarrow{r}}_{23}$ (19a)

where ρ = N/V is the number density of molecules, $h\left(r\right)=g\left(r\right)-1$ , g(r) is the radial distribution function, u(r) is the direct interaction potential between pairs, β=1/k_BT, and under the integral ${\overrightarrow{r}}_{13}={\overrightarrow{r}}_{12}+{\overrightarrow{r}}_{23}$ .

HC equation yields the correlation function $g\left(r,\beta \right)$ , from the interaction potential u(r).

with $r_0=\left|{\overrightarrow{r}}_{12}\right|$ , $r_1=\left|{\overrightarrow{r}}_{23}\right|$ , $R=\left|{\overrightarrow{r}}_{13}\right|$ , and from cosine theorem $R_1^2=r_0^2+r_1^2{\cos }^2\theta +2r_0{r}_1\cos \theta \cos \varphi$ , $R^2=r_1^2{\sin }^2\theta +R_1^2=r_0^2+r_1^2+2r_0{r}_1\cos \theta \cos \varphi$ (Figure 3).

Figure 3. Distance in the HC equation.

We reformulate HC equation

$\begin{array}{l}\log g\left(r_0\right)+\beta u\left(r_0\right)\\ =\rho \left(1-g\left(r_0\right)\right){\displaystyle \int \left(\left(1-g\left(R\right)\right)-\log g\left(R\right)-\beta u\left(R\right)\right)}{r}_1{}^{2}d{r}_1\sin \theta d\theta d\varphi \end{array}$

concisely $\log g\left(r_0\right)+\beta u\left(r_0\right)-\rho \left(1-g\left(r_0\right)\right)I\left(g,\beta u\right)=0$ (19b)

$I\left(g,\beta u\right)={\displaystyle \int \left(\left(1-g\left(R\right)\right)-\log g\left(R\right)-\beta u\left(R\right)\right)}{r}_1^{2}d{r}_1\sin \theta d\theta d\varphi$

where $\rho =\frac{N}{V}$ , for ideal gas ${\rho }_0=\rho \left(T=273\text{K},p=1\text{at}\right)=2.69\times {10}^{25}{\text{m}}^{-3}$ , or in

angstrom $\dot{\text{A}}={10}^{-10}\text{m}$, ${\rho }_0=2.69\times {10}^{-5}{}^{\dot{\text{A}}}\ll 1$ and mean distance ${\lambda }_0=\lambda \left(T=273\text{K},p=1\text{at}\right)=33\dot{\text{A}}$ in the limit of small $\rho \ll 1$ $g\left(r_0\right)=\exp \left(-\beta u\left(r_0\right)\right)\approx 1-\beta u\left(r_0\right)$ .

3.3. Radial Distribution Function Lennard-Jones Potential

Lennard-Jones potential reads [11]

$u_{LJ}\left(r,\sigma ,\epsilon \right)=4\epsilon \left({\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right)$

σ is the van der Waals radius = distance at which u = 0, $\sigma =r_0$ , where r₀ = half molecule diameter, r is the distance between particles.

The correlation function (=radial distribution function rdf) is here a decaying harmonics (Figure 4), with a sharp first maximum at $r=\sigma$ , and an amplitude of ~4x the second maximum amplitude.

Figure 4. The plot of a typical radial distribution function for the monatomic Lennard-Jones liquid. (here with σ = 3.73Å and ϵ = 0.294 kcal/mol at a temperature of 111.06 K).

3.4. Radial Distribution Function and the vdWaals Equation

3.4.1. Virial Equation

The virial equation ([8] [12] chap. 1.2) is a generalization of the vdWaals equation, where the ideal gas law is expanded in a power series in specific volume v, or in

particle density $\rho =\frac{1}{v}$ .

The ideal gas law $p\beta v=\frac{p\beta }{\rho }=1$ can be formulated as an approximate invariant $E_c\left(v,\beta \right)=p\beta v$ . Now we expand E_c in a power series in specific volume v, or in particle density $\rho =\frac{1}{v}$ .

$E_c\left(v,\beta \right)=1+\frac{B_1\left(\beta \right)}{v}+\frac{B_2\left(\beta \right)}{v^2}+\cdots$

resp. $E_c\left(\rho ,\beta \right)=1+B_1\left(\beta \right)\rho +B_2\left(\beta \right){\rho }^2+\cdots$

The virial equation in integral reads

$p\beta =\rho \left(1-\frac{\rho \beta }{6}{\displaystyle \int d^{3}rg\left(\overrightarrow{r}\right)\overrightarrow{r}\cdot \nabla u\left(\overrightarrow{r}\right)}\right)$ virial equation (20a)

In comparison with vdWaals equation $p\beta =\frac{\rho }{1-\rho b}-a{\rho }^2\beta$ we see that it is

the virial expansion truncated at first order, and with a volume correction.

We can reformulate the virial equation

$p\beta =\frac{\rho }{1-\rho b}-\frac{{\rho }^2\beta }{6}{\displaystyle \int d^{3}rg\left(\overrightarrow{r}\right)\overrightarrow{r}\cdot \nabla u\left(\overrightarrow{r}\right)}$ virial with vdWaals-correction $b=2\pi {\sigma }^3/3$ (20b)

$p\beta =\frac{\rho }{1-\rho \left(2\pi {\sigma }^3/3\right)}-\frac{{\rho }^2\beta }{6}{\displaystyle \underset{\sigma }{\overset{\infty }{\int }}dr4\pi r^3{u}^{\prime }\left(r\right)g\left(r\right)}$ spherical symmetry, with vdWaals-correction (20c)

3.4.2. Van-der-Waals Equation

The vdWaals equation can be derived from the hard-sphere potential

$u_0\left(r\right)=\{\begin{array}{l}0r>\sigma \\ \infty r\le \sigma \end{array}$ , $g_0\left(r\right)=\{\begin{array}{l}1r>\sigma \\ 0r\le \sigma \end{array}$ , with excluded volume per particle $b=\frac{2}{3}\pi {\sigma }^3$ .

The partition function is

$Z\left(N,V,T\right)=\frac{1}{N!h^{3N}}{\displaystyle \int d{r}_1\cdots d{r}_N\exp \left(-\beta U\left({\overrightarrow{r}}_1,\cdots ,{\overrightarrow{r}}_N\right)\right)}=\frac{Z_N}{N!h^{3N}}$

with 0-order including own-volume correction $Z_N^{\left(0\right)}={\left(V-Nb\right)}^N$ ,

$\begin{array}{c}\log Z=\left(N\left(\log V-3\log \varsigma \right)-\log N!\right)\\ \approx \left(N\left(\log V-3\log \varsigma \right)-N\left(\log N-1\right)\right)\end{array}$

The free energy becomes $F\left(N,V,T\right)\approx -\frac{1}{\beta }\log \left(\frac{{\left(V-Nb\right)}^N}{N!{\varsigma }^{3N}}\right)-a\frac{N^2}{V}$

and specific free energy (per particle) with specific volume v and density ρ

$\begin{array}{c}f\left(v,T\right)=\frac{F\left(N,v,T\right)}{N}\approx -\frac{1}{N\beta }\log \left(\frac{N{\left(v-b\right)}^N}{N!{\varsigma }^{3N}}\right)-\frac{a}{v}\\ =-\frac{1}{\beta }\left(\log \left(\frac{v-b}{{\varsigma }^3}\right)-\left(\log N-1\right)\right)-\frac{a}{v}\end{array}$

The pressure results as

$p=-\frac{\partial F}{\partial V}=-\frac{\partial f}{\partial v}=\frac{1}{\beta \left(v-b\right)}-\frac{a\beta }{v^2}=\frac{\rho }{\beta \left(1-\rho b\right)}-a\beta {\rho }^2$ (21)

which is the vdWaals equation.

3.4.3. Derivation of the Rdf Function g(r)

The modified virial eq. reads $p\beta =\frac{\rho }{1-\rho b}-\frac{{\rho }^2\beta }{6}{\displaystyle \underset{\sigma }{\overset{\infty }{\int }}dr4\pi r^3{u}^{\prime }\left(r\right)g\left(r\right)}$ , vdWaals eq. is $p\beta =\frac{\rho }{1-\rho b}-a{\rho }^2\beta$ thus ${\displaystyle \underset{\sigma }{\overset{\infty }{\int }}dr{r}^3{u}^{\prime }\left(r\right)g\left(r\right)}=\frac{3a}{2\pi }$ , ${\lim }_{r\to \infty }g\left(r\right)=1$ , ${\lim }_{r\to \sigma +}g\left(r\right)=0$

A possible ansatz for g(r) is

$g\left(r\right)={\Theta }_H\left(r-\sigma \right){\displaystyle \sum _k\left({\alpha }_k\cos \left(\frac{2k}{\pi }\frac{r-\sigma }{\sigma }\right)\right)}$ (22)

or explicitly adapted for Lennard-Jones potential

$\begin{array}{l}{g}_{LJ}\left(r,\sigma ,dr,{\alpha }_1,k,ex\right)\\ =\left(1-\frac{1}{1+\exp \left(\frac{r-\sigma }{dr\sigma }\right)}\right)\left(1+\exp \left(-\frac{exr}{\sigma }\right)\left(1-\frac{1}{1+\exp \left(\frac{r-\sigma }{dr\sigma }\right)}\right)\left({\alpha }_1\cos \left(\frac{2k}{\pi }\frac{r-\sigma }{\sigma }\right)\right)\right)\end{array}$

Where step-up-function ${\Theta }_H\left(r,\sigma ,dr\right)=\left(1-\frac{1}{1+\exp \left(\frac{r-\sigma }{dr\sigma }\right)}\right)$ with the step radius σ, relative exponential width dr.

4. Landau Theory

4.1. Basics of Landau Theory

In the Landau theory of phase transitions, the free energy is expressed in terms of a complex order parameter field $\phi$ , where the quantity ${\left|\phi \left(r\right)\right|}^2$ is a measure of the local particle density in analogy to a quantum mechanical wave function.

The generalized Landau ansatz for the partition function is a functional integral [1]

$Z={\displaystyle \int D\phi dV\exp \left(-H\left(\overrightarrow{r},\phi ,\beta \right)\right)}$ (23)

where $H\left(\overrightarrow{r},\phi ,\beta \right)$ is the thermodynamic Hamiltonian.,

$H\left(\overrightarrow{r},\phi ,\beta \right)=E(\phi )+\beta u\left(\overrightarrow{r}\right)$, where $E\left(\phi \right)$ is the φ-induced Landau energy, $\beta u\left(\overrightarrow{r}\right)$ is the thermal intermolecular potential.

The functional integration $D\phi$ has here a precise mathematical meaning:

$D\phi ={\displaystyle \sum _i\frac{\partial \phi }{\partial c_i}}d{c}_i$ and is perfectly well-defined.

The function variable φ is a function of location $\overrightarrow{r}$ and molecular distribution parameters c_i (like molecular diameter σ, average distance λ, correlation length l_c, local distribution periods l_ak) $\phi =\phi \left(r;c_i\right)$ .

The volume-integration runs over the location differential $r^2\sin \theta drd\theta d\varphi$ .

The actual function φ₀, which yields the valid partition function of a

thermodynamic system, is found by minimization of the free energy

$F=-\frac{1}{\beta }\log Z$ in the parameters c_i in a value range Ω, which yield optimal

parameters c_i0 and the corresponding function ${\phi }_0=\phi \left(r;c_{i0}\right)$

$F\left(\phi \left(c_{i0}\right)\right)=\min \left(F\left(\phi \left(c_i\right)\right);c_i\in \Omega \right)$ .

Landau function and the rdf function: a heuristic derivation

We start with the configurational partition function $Z_N={\displaystyle \int d^3{r}_1\cdots d^3{r}_N\exp \left(-\beta U\left({\overrightarrow{r}}_1,\cdots ,{\overrightarrow{r}}_N\right)\right)}$,

and reformulate it for one particle with spherical symmetry ${\left(Z_N\right)}^{1/N}=z\left(\beta \right)={\displaystyle \int d^{3}rg\left(r\right)\exp \left(-\beta u\left(r\right)\right)}$

where $z\left(\beta \right)$ is the specific (per particle) partition function, $u\left(r\right)=\epsilon u_0\left(r\right)$ is the relative intermolecular potential with characteristic energy ε, $g\left(r\right)$ is the rdf function.

We reformulate it $z\left(\beta \epsilon \right)={\displaystyle \int dVg\left(r\right)\exp \left(-\beta \epsilon u_0\left(r\right)\right)}$ .

Now we introduce the finite-range factor $f_{\lambda }=\exp \left(-r/\lambda \right)$ with the average

distance λ, where the specific volume is $v=\frac{V}{N}={\lambda }^3$ . The finite-range factor is

approximately $f_{\lambda }=1$ for small distances and large λ in weakly bound systems like gases. $f_{\lambda }$ represents the exponentially decaying local influence within the range λ of an individual particle.

Furthermore, we must take into account the restricted volume with radius σ around the particle, where there are no particles, so we introduce the step function ${\Theta }_H\left(r,\sigma \right)$ into the integral.

Taking all together, we obtain the general Landau form for one particle

$Z={\displaystyle \int D\phi dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}\approx {\displaystyle \int \phi \left(r,\sigma ,\lambda ;c_k\right)dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$ (24a)

where the Landau function has the form

$\phi \left(r,\sigma ,\lambda ;c_k\right)={\Theta }_H\left(r,\sigma \right)\exp \left(-\frac{r}{\lambda }\right)g\left(r;c_k\right)$ with parameters $c_k$ (24b)

We know from measurements that the rdf function can be described as a exponentially damped harmonic function, so we make the general ansatz

$g\left(r;c_k\right)=\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)$ (24c)

Forms of Landau function

For first-order transitions, we make the generalized Landau ansatz for the partition function as a functional integral,

The Landau function depends on parameters $\left(c_k\right)=\left(l_c,{\alpha }_k,l_{a,k},{\alpha }_{0k}\right)$ , the

differential is $D\phi ={\displaystyle \sum \frac{\partial \phi }{\partial c_k}}d{c}_k$ .

Also, the Boltzmann-exponential contains only the intermolecular potential u, the partition function is

$Z={\displaystyle \int D\phi dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$

The ansatz for the Landau function for solid-fluid-gas transition is the λ-damped rdf function

$\begin{array}{l}\phi \left(r,\sigma ,\lambda ;c_k\right)=\exp \left(-\frac{r}{\lambda }\right){\phi }_{rdf}\left(r,\sigma ;c_k\right)\\ {\phi }_{rdf}\left(r,\sigma ;c_k\right)={\Theta }_H\left(r,\sigma \right)\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)\end{array}$ (24d)

with the soft-step-up function

${\Theta }_H\left(r,\sigma ,\Delta r\right)=\left(1-\frac{1}{1+\exp \left(\frac{r-\sigma }{\Delta r\sigma }\right)}\right)$

Calculation with this ansatz yields results, which are successfully verified for Lennard-Jones-substances and dipole-substances in chap.10.

Here, we have the typical behavior of first-order phase transitions: there is latent heat, the first derivative of free energy is discontinuous, different phases (e.g. liquid-gas) coexist at transition point.

For second-order transitions, like in magnetic systems, the Landau function has no parameters $\left(c_k\right)=\left\{\right\}$ , so the differential is $D\phi =d\phi$ .

Also, the partition function contains a φ-generated kinetic energy in the exponential.

$Z={\displaystyle \int D\phi dV\exp \left(-E\left(\phi \right)\right)}$ (25a)

$E\left(\phi \right)$ is the energy generated by $\phi$ , and is a generalization of the kinetic energy $\frac{\mu }{2}{\phi }^2$ , in Landau’s original ansatz

$E\left(\phi ,\beta \right)=\frac{\mu \left(\beta \right)}{2}{\phi }^2+\frac{\lambda \left(\beta \right)}{4!}{\phi }^4$ (25b)

For magnetic systems, $\phi =M_s=\frac{M}{N}$ specific magnetization, and the thermal potential is $-\beta B\phi$ , with interaction $E\left(\phi ,\beta \right)=\frac{\mu \left(\beta \right)}{2}{\phi }^2+\frac{\lambda \left(\beta \right)}{4!}{\phi }^4-\beta \phi B$

Without interaction, the free energy is

$F\left(T,V,N\right)=-\frac{1}{\beta }\log Z\approx \frac{1}{\beta }E\left(\phi ,\beta \right)=\frac{1}{\beta }\left(\frac{\mu \left(\beta \right)}{2}{\phi }^2+\frac{\lambda \left(\beta \right)}{4!}{\phi }^4\right)$ ,

In this way, we obtain the original Landau ansatz

$F\left(T,\phi \right)=F_0-a\left(T\right){\phi }^2+\frac{b\left(T\right)}{2}{\phi }^4$ , $a\left(T\right)=\frac{1}{\beta }\frac{\mu \left(\beta \right)}{2}$ , $b\left(T\right)=-\frac{\lambda \left(\beta \right)}{12\beta }$

At critical temperature T_c: $\mu \left(\beta ={\beta }_c=\frac{1}{k_B{T}_c}\right)=0$ , $\mu \left(\beta \right)={\mu }_0\left(1-\frac{{\beta }_c}{\beta }\right)$ , ${\beta }_c=\frac{1}{k_B{T}_c}$

$a\left(T\right)\approx a_0\left(T-T_c\right),T\to T_c$ , $a_0=-{\mu }^{\prime }\left({\beta }_c\right){\beta }_c^2{k}_B\left(T-T_c\right)$ , $b\left(T_c\right)=-\frac{\lambda \left({\beta }_c\right)}{12{\beta }_c}$

With the assumptions $b\left(T\right)=b_0,T\to T_c$ and $a\left(T\right)\approx a_0\left(T-T_c\right),T\to T_c$ near the critical temperature T_c, minimizing the free energy with respect to $\phi$ requires

$\frac{\partial F}{\partial \phi }=-2a\left(T\right)\phi +2b\left(T\right){\phi }^3=0$

The solutions either $\phi =0$ , or

${\phi }^2=\frac{a}{b}=\frac{a_0}{b_0}\left(T-T_c\right)$ , $\phi \propto {\left|T-T_c\right|}^{1/2}$

This yields for free energy the typical second-order transition:

$F-F_0=\{\begin{array}{l}\frac{a_0^2}{2b_0}{\left(T-T_c\right)}^{2}T<T_c\\ 0T>T_c\end{array}$ free energy, and its first derivative are continuous, $c_v=-T\frac{{\partial }^{2}F}{\partial T^2}=\{\begin{array}{l}\frac{a_0^2}{b_0}TT<T_c\\ 0T>T_c\end{array}$ specific heat (second derivative) is discontinuous.

4.2. The Landau-Ising Model

The dimensionless Landau-Ginzburg Hamiltonian density in d dimensions reads ([1] chap.7)

$H_{LG}=\frac{1}{2}{\alpha }^2{\left|\nabla \cdot \overrightarrow{\phi }\right|}^2+\frac{1}{2}{\mu }^2{\left|\overrightarrow{\phi }\right|}^2+\frac{1}{4!}{\lambda }^2{\left|\overrightarrow{\phi }\right|}^4$ (26a)

where ϕ = magnetization in the Ising model, with external field B it becomes

$H_{LG}=\frac{1}{2}{\alpha }^2{\left|\nabla \cdot \overrightarrow{\phi }\right|}^2+\frac{1}{2}{\mu }^2{\left|\overrightarrow{\phi }\right|}^2+\frac{1}{4!}{\lambda }^2{\left|\overrightarrow{\phi }\right|}^4-\beta \overrightarrow{B}\cdot \overrightarrow{\phi }$ (26b)

with parameters:

characteristic length α, dimensionless ${\mu }^2\propto \frac{T-T_c}{T_c}$ whose temperature variation drives phase transition, (positive) d-volume λ > 0.

The partition function is $Z_{LG}={\displaystyle \int D\phi \exp \left(-{\displaystyle \int d^{d}rH\left(\phi \right)}\right)}$ (functional integral), the simplest approximation is $Z_{LG}=A\exp \left(-{\displaystyle \int d^{d}rH\left({\phi }_0\right)}\right)$ , with A constant, ${\phi }_0$ minimizing the integral ${\displaystyle \int d^{d}rH\left({\phi }_0\right)}$ the simplified integral becomes ${\displaystyle \int d^{d}rH({\phi }_0)}=V\left(const+\frac{1}{2}{\mu }^2{\left|{\overrightarrow{\phi }}_0\right|}^2+\frac{1}{4!}{\lambda }^2{\left|{\overrightarrow{\phi }}_0\right|}^4-\beta \overrightarrow{B}\cdot {\overrightarrow{\phi }}_0\right)$, where V = d-volume of the system.

The free energy density becomes

$f=-\frac{1}{\beta V}\log Z=\frac{const}{\beta }+\frac{1}{2\beta }{\mu }^2{\phi }_0^2+\frac{1}{4!\beta }{\lambda }^2{\phi }_0^4-B{\phi }_0$ , free energy

$F=V\left(\frac{const}{\beta }+\frac{1}{2\beta }{\mu }^2{\phi }_0^2+\frac{1}{4!\beta }{\lambda }^2{\phi }_0^4-B{\phi }_0\right)$ (26c)

Near T = T_c critical temperature, ${\mu }^2\propto \left(T-T_c\right)$ , for B = 0 minimization solution $\frac{\partial H_{LG}}{\partial \phi }=0$ gives ${\mu }^2{\phi }_0+\frac{1}{3!}{\lambda }^2{\phi }_0^3=0$ , from this follows ${\phi }_0^2=\{\begin{array}{l}0{\mu }^2\ge 0\\ -\frac{6{\mu }^2}{\lambda }{\mu }^2<0\end{array}$ , ${\phi }_0\propto {\left|T-T_c\right|}^{1/2}$ .

The free energy density becomes $f=\frac{1}{\beta }\left(const+\{\begin{array}{l}0T>T_c\\ -\frac{3{\mu }^4}{2\lambda }T<T_c\end{array}\right)$ and the internal energy density becomes

$u=\frac{d\left(\beta f\right)}{d\beta }=\{\begin{array}{l}0T>T_c\\ -\frac{3{\mu }^2}{\lambda }\frac{d{\mu }^2}{d\beta }+\frac{3{\mu }^4}{2{\lambda }^2}\frac{d\lambda }{d\beta }T<T_c\end{array}$ ,

is continuous at T_c (${\mu }^2=0$ ), so there is no latent heat.

The specific heat is discontinuous at T_c

${k_B{c}_v|}_{{\mu }^2=0}=-{\beta }^2{\frac{du}{d\beta }|}_{{\mu }^2=0}=\{\begin{array}{l}0T=T_c+\\ \frac{3{\beta }^2}{\lambda }{\left(\frac{d{\mu }^2}{d\beta }\right)}^{2}T=T_c-\end{array}$

With field B, minimization condition reads ${\mu }^2{\phi }_0+\frac{1}{3!}{\lambda }^2{\phi }_0^3-\beta B=0$ , which is a reduced cubic ${\phi }_0^3+\frac{6{\mu }^2}{\lambda }{\phi }_0-\frac{6\beta }{\lambda }B=0$ , Cardano form $x^3+3px+2q=0$ ,

$u=\sqrt[3]{-q+\sqrt{q^2+p^3}}$ , $v=\sqrt[3]{-q-\sqrt{q^2+p^3}}$

$u=\sqrt[3]{\frac{3B}{\lambda }+\sqrt{\frac{9B^2}{{\lambda }^2}+\frac{8{\mu }^6}{{\lambda }^3}}}$ , $v=\sqrt[3]{\frac{3B}{\lambda }-\sqrt{\frac{9B^2}{{\lambda }^2}+\frac{8{\mu }^6}{{\lambda }^3}}}$

with solution ${\phi }_{0,1}=u+v$ .

For $T=T_c$ , ${\mu }^2=0$ , ${\phi }_{0,1}^3=\frac{6\beta }{\lambda }B$ , so critical exponent δ = 1/3, and for susceptibility ${\chi }_T=\frac{\partial {\phi }_0}{\partial B}$ we obtain $\left({\mu }^2+\lambda {\phi }_0^2/2\right){\chi }_T=\beta$ , so

${{\chi }_T|}_{B=0}=\{\begin{array}{l}\frac{\beta }{{\mu }^2}T>T_c\\ \frac{\beta }{2{\left|\mu \right|}^2}T<T_c\end{array}$

At the critical temperature T_c, the susceptibility has a singularity.

5. Van-der-Waals Theory

5.1. Basic vdWaals Theory

The van der Waals equation for real gasses reads ([7]-[9])

$p=\frac{RT}{v_m-b}-\frac{a}{v_m^2}$ , $p=\frac{RT{\rho }_m}{1-b{\rho }_m}-a{\rho }_m^2$ , $v=1/{\rho }_m$ , ${\rho }_m=\rho \frac{1}{N_A}=\frac{N/N_A}{V}=\frac{n_m}{V}$ (27a)

where ${\rho }_m$ is the molar particle density, p is pressure, T is temperature, and v_m = VN_A/N is molar volume, N_A is the Avogadro constant, V is the volume, and N is the number of molecules, R = N_A/k is the universal gas constant, k is the Boltzmann constant, and a and b are substance-specific constants.

Another form of the vdWaals equation is

$p=\frac{N{k}_{B}T}{V-bN}-a\frac{N^2}{V^2}$

In specific form (per particle) the vdWaals eos becomes

$p=\frac{1}{\beta \left(v-b\right)}-\frac{a}{v^2}$ (27b)

The constant a expresses the strength of the molecular interactions, it has dimension [p V²], or [E V]. The constant b denotes the excluded self-volume of a particle.

Ideal gas law is the limit a = 0, b = 0

$p=\frac{RT}{v_m}$ (28a)

Another form ideal gas

$pV=N{k}_{B}T$ (28b)

specific per particle $pv=k_{B}T=1/\beta$ , or $pv\beta =1$ in terms of mean distance $\lambda =v^{1/3}$ and inverse thermal energy β we obtain

$p=\frac{1}{\beta {\lambda }^3}=\frac{1}{\beta v}$ (28c)

The isothermal compressibility becomes

${\kappa }_T=\frac{v_m^2{\left(v_m-b\right)}^2}{RT{v}_m^3-2a{\left(v_m-b\right)}^2}$

and coefficient of thermal expansion,

$\alpha =\frac{R{v}_m^2\left(v_m-b\right)}{RT{v}_m^3-2a{\left(v_m-b\right)}^2}$

in the limit $v_m\to \infty$ : $\alpha =\frac{1}{T}$ , ${\kappa }_T=\frac{v_m}{RT}$ .

The vdWaals parameters can be calculated from molecular parameters $\sigma$ = molecule diameter, $\epsilon$ = characteristic energy of the inter-molecular potential.

• b-parameter $b=\frac{2\pi }{3}{\sigma }^3$

• a-parameter generally $a=Ib\epsilon$ , I = dimensionless factor, depends on the form of $u/\epsilon$ -cut-off Lennard-Jones potential: $u_{LJ}\left(r,\sigma ,\epsilon \right)=4\epsilon \left({\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right)$

$a=\frac{4e}{3\left(3e-8\right)}\epsilon {\sigma }^3=23.4\epsilon {\sigma }^3$ [13]

$\min \left(u_{cLJ}\left(r,1,1\right),r\right)=u_{cLJ}\left(r_{\min }=1.22\right)=-0.58$ , min. energy $u_0=0.58\epsilon$ , Morse potential $u_M\left(r,\sigma ,\epsilon \right)=\epsilon \left({\left(1-\exp \left(-\frac{r-\sigma }{\sigma }\right)\right)}^2-1\right)$

$a=\frac{35e}{8\left(3e-8\right)}\epsilon {\sigma }^3=76.8\epsilon {\sigma }^3$ [13]

Derivation of van der Waals equation

From the first-order partition function [14]

$Z\left(N,V,T\right)=\frac{V^N}{N!{\zeta }^{3N}}\left(1+\frac{\beta N^2}{2V}{\displaystyle \int u\left(r\right)d^{3}r}+\cdots \right)$ , where $u\left(r\right)$ is the inter-molecular potential, $\zeta$ thermal wavelength, we obtain for pressure $p=-\frac{\partial F}{\partial V}=\frac{1}{\beta }\frac{\partial \log Z}{\partial V}\approx \frac{N{k}_{B}T}{V}\left(1+\frac{N\beta }{2V}{\displaystyle \int u\left(r\right)d^{3}r}\right)$ , in variables $\lambda ,\beta$

$\begin{array}{c}p\left(\lambda ,\beta \right)=-\frac{\partial F}{\partial V}=\frac{1}{\beta }\frac{\partial \log Z}{\partial V}=\frac{1}{\beta }\frac{1}{4\pi {\lambda }^2}\frac{1}{Z}\frac{\partial Z\left(\lambda ,\beta \right)}{\partial \lambda }\\ \approx \frac{1}{{\lambda }^3\beta }\left(1+\frac{1}{2{\lambda }^3}{\displaystyle \int u\left(r\right)d^{3}r}\right)\end{array}$ (29a)

The equation-of-state (eos) reads

$p\left(v-v_0\right)\beta =\left(1-\frac{\beta }{2\left(v-v_0\right)}{\displaystyle \underset{r=\sigma }{\overset{\infty }{\int }}\left|u\left(r\right)\right|d^{3}r}\right)$ (29b)

where $v_0=\frac{2}{3}\pi {\sigma }^3$ is the own volume of the molecule, and $u\left(r\right)$ is the attractive (negative) energy density of the intermolecular potential, ${\displaystyle \underset{r=\sigma }{\overset{\infty }{\int }}\left|u\left(r\right)\right|d^{3}r}=v_0{u}_0$ , with u₀ = mean inter-molecular energy.

We obtain the vdWaals eos in the form $p=\frac{1}{\beta \left(v-b\right)}-\frac{a}{{\left(v-b\right)}^2}$ (29c)

where $b=v_0=\frac{2}{3}\pi {\sigma }^3$ , $a=\frac{1}{2}{\displaystyle \underset{r=\sigma }{\overset{\infty }{\int }}\left|u\left(r\right)\right|d^{3}r}=\frac{v_0{u}_0}{2}$ .

5.2. Mean Distance in Ideal Gas

The probability to find a particle at the distance from the origin between r and r + dr is

$P_N\left(r\right)dr=4\pi r^{2}dr{\left(1-\frac{4\pi r^3}{3V}\right)}^{N-1}=\frac{3}{a}{\left(\frac{r}{a}\right)}^{2}dr{\left(1-{\left(\frac{r}{a}\right)}^3\frac{1}{N}\right)}^{N-1}$

where we substituted $\frac{1}{V}=\frac{3}{4\pi N{a}^3}$ , and a is the mean distance.

Finally, taking the $N\to \infty$ limit we obtain

$P\left(r\right)=\frac{3}{a}{\left(\frac{r}{a}\right)}^2\exp \left(-{\left(\frac{r}{a}\right)}^3\right)$

The distribution peaks at $r_{\max }=\sqrt[3]{\frac{2}{3}}a=0.874a$ .

5.3. Liquid-gas Transition

The vdWaals eos is $p=\frac{kT}{v-b}-\frac{a}{v^2}$ , $v=\frac{V}{N}$

The critical temperature T_c in the liquid-gas transition results from

$\frac{dp}{dv}=0$ , $\frac{d^{2}p}{d{v}^2}=0$ , follows $p_c{\left(v-v_c\right)}^3=0$

$k{T}_c=\frac{8a}{27b}$ , $v_c=3b$ , $p_c=\frac{a}{27b^2}$

With reduced variables $T_r=\frac{T}{T_c}$ , $v_r=\frac{v}{v_c}$ , $p_r=\frac{p}{p_c}$

vdWaals equation becomes universal

$p_r=\frac{8}{3}\frac{T_r}{v_r-1/3}-\frac{3}{v_r^2}$ (29d)

and also and the compressibility ratio is universal $\frac{p_c{v}_c}{k_B{T}_c}=\frac{3}{8}=0.375$ (Figure 5).

Figure 5. [14] Below are shown four isotherms of the universal vdWaals equation in relative coordinates with the spinodal curve $\frac{\partial p}{\partial v}=0$ (black dash-dot curve) and the saturation curve (red dash-dot curve). The critical point lies at the turning point $\frac{{\partial }^{2}p}{\partial v^2}=0$ on the orange isotherm. The saturation curve (left wing=fluid, right wing=gas) left (low volume) wing ends at the triple point, its points are determined by Maxwell’s equal-area rule $F\left(v_1,T\right)=F\left(v_2,T\right)$ , $p\left(v_1,T\right)=p\left(v_2,T\right)$ , where $p\left(v,T\right)=-\frac{\partial F\left(v,T\right)}{\partial V}$ .

Saturation curve

Extended principle of corresponding states has been suggested in which $p_r=p\left(v_r,T_r,\varphi \right)$ ,

where ϕ is a substance-dependent dimensionless parameter, $\varphi =\frac{p_c{v}_c}{k_B{T}_c}$

better candidate is $\omega =-{\log }_{10}\left(p_r\left(T_r=0.7\right)\right)$

The approximate saturation curve is [15] (Figure 6)

$\log \left(p_{rs}\right)=5.37\left(1-\frac{1}{T_r}\right)+\omega \left(\text{7}\text{.49}-11.18T_r^3+3.69T_r^6+17.93\log T_r\right)$

with $\omega =-{\log }_{10}\left(p_r\left(T_r=0.7\right)\right)$

Figure 6. The family of saturation curves, showing the vdW curve as a member (blue curve). The blue dots are calculated from Lekner’s solution [16].

6. Lennard-Jones Substance

6.1. Lennard-Jones Fluid

The Lennard-Jones Potential is given by the following equation [17] [18]:

$u_{LJ}\left(r,\sigma ,\epsilon \right)=4\epsilon \left({\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right)$

σ is the van der Waals radius = distance at which u = 0, $\sigma =r_0$ , where r₀ = molecule diameter, r is the distance between particles.

Cut-off Lennard-Jones function

Infinite potential at r = 0 is unrealistic, much better is the corresponding cut-off potential.

$u_{cLJ}\left(r,\sigma ,\epsilon \right)=4\epsilon \left(\frac{1}{\frac{1}{2}+{\left(\frac{r}{\sigma }\right)}^{12}}-\frac{1}{1+{\left(\frac{r}{\sigma }\right)}^6}\right)$ , normalized $u_{cLJ}\left(r,\sigma =1,\epsilon =1\right)$ has the form [13] $\min \left(u_{cLJ}\left(r,1,1\right),r\right)=u_{cLJ}\left(r_{\min }=1.22\right)=-0.58$ (Figures 7-9).

Figure 7. Normalized Lennard-Jones potential [13].

Figure 8. Lennard-Jones fluid isobars p(v,T): (constant-pressure curves) [18].

Figure 9. Phase diagram of the Lennard-Jones substance [19], star = critical point, circle indicates the vapor-liquid-solid triple point, triangle indicates the vapor-solid (fcc)-solid (hcp) triple point, solid lines indicate coexistence lines of two phases, dashed lines indicate the vapor-liquid spinodal.

The most important characteristic points of the Lennard-Jones potential are the critical point and the vapor-liquid-solid triple point.

The critical point parameters are

$\epsilon {\beta }_{cr}=1/1.321=0.755$ , $v_{cr}=1/0.316{\sigma }^3=3.16{\sigma }^3$ , ${\lambda }_{cr}={\left(v_{cr}\right)}^{1/3}=1.47\sigma$ ,

$p_{cr}=0.129p_0=6.47\text{MPa}$ , $p_0\left(Ar\right)=\frac{\epsilon }{{\sigma }_{LJ}^3}=50.1\text{MPa}$ ,

The triple point parameters are

$T_{tr}=83.8\text{K}$ , density fluid, solid ${\rho }_{tr,f}=0.845{\sigma }^{-3}$ , ${\rho }_{tr,s}=0.961{\sigma }^{-3}$

$\epsilon {\beta }_{tr}=1/0.69=1.45$ , $v_{tr,f}={\sigma }^3/0.845=1.18{\sigma }^3$ , ${\lambda }_{tr}={\left(v_{tr}\right)}^{1/3}=1.06\sigma$ ,

$p_{tr}=0.0012p_0=0.060\text{MPa}$ , $p_0\left(Ar\right)=\frac{\epsilon }{{\sigma }_{LJ}^3}=50.1\text{MPa}$ ,

gas (300K, 1bar): $\rho =2.7\times {10}^{19}{\text{cm}}^{-3}$ , $v=1/\rho =3.7\times {10}^{-20}{\text{cm}}^3$ , $\lambda =v^{1/3}=3.3\times {10}^{-7}\text{cm}=33\dot{\text{A}}$,

mean free path $l_{free}=52\dot{\text{A}}$

Comparison with measured values in LJ-Argon (see below)

triple point $T_t=83.7\text{K}$ , $p_{tr}=0.68\text{bar}=0.068\text{MPa}=1.36\times {10}^{-3}{p}_0$ $\beta \epsilon \left(T_{tr}\right)=1.67$

critical point $T_c=150.8\text{K}$ , $p_c=48.3\text{bar}=4.83\text{MPa}=0.0964p_0$ $\beta \epsilon \left(T_c\right)=0.938$

Measured values for argon are as follows [20]-[22].

${\sigma }_R\left(Ar\right)=3.72\dot{\text{A}}$, ${\sigma }_{LJ}\left(Ar\right)=3.4\dot{\text{A}}$, $\epsilon \left(Ar\right)=0.0123\text{eV}$ , $p_0\left(Ar\right)=\frac{\epsilon }{{\sigma }^3}=689\text{MPa}$

$10\text{bar}=1\text{MPa}={10}^6\text{J}/{\text{m}}^{\text{3}}={10}^6\times 6.242\times {10}^{18}\text{eV}/{\text{m}}^{\text{3}}=6.242\times {10}^{-6}\text{eV}/{}^{\dot{\text{A}}}$, $p_0\left(Ar\right)=\frac{\epsilon }{{\sigma }_{LJ}^3}=0.312\times {10}^{-3}\text{eV}/{}^{\dot{\text{A}}}=0.501\times {10}^2\text{MPa}=50.1\text{MPa}$

$k_B\times 1\text{K}=0.026\text{eV}/300=0.000086\text{eV}$ , $\epsilon \left(Ar\right)=142.095\times 0.000086=0.0123\text{eV}$ , ${\sigma }_{LJ}\left(Ar\right)=0.336\text{nm}=3.36\dot{\text{A}}$

$T_m\left(p=1\text{bar}\right)=83.8\text{K}$ , $\beta =1/\left(0.026\times 83.8/300\right)=1/0.00726=137.7{\text{eV}}^{-1}$ , $\beta \epsilon \left(T_m\right)=1.69$

$T_b\left(p=1\text{bar}\right)=87.3\text{K}$ , $\beta =1/\left(0.026\times 87.3/300\right)=1/0.00757=132.2{\text{eV}}^{-1}$ , $\beta \epsilon \left(T_b\right)=1.62$

triple point $T_t=83.7\text{K}$ , $p_t=0.68\text{bar}=0.068\text{MPa}=1.36\times {10}^{-3}{p}_0$ $\beta \epsilon \left(T_t\right)=1.67$

critical point $T_c=150.8\text{K}$ , $p_c=48.3\text{bar}=4.83\text{MPa}=0.097p_0$ $\beta \epsilon \left(T_t\right)=0.978$

6.2. Lennard-Jones Radial Distribution Function

The measured and calculated rdf functions of argon LJ-substance is shown in Figure 10, Figure 11.

Figure 10. [23] The first comparison between the rdf of a real liquid from neutron scattering and the rdfs obtained independently by Scott and Bernal from RCPs of about 1000 equal spheres; period $l_a=0.75$ , distribution weight $g\left(r=1\right)=2.35$ , where $g\left(r=\infty \right)=1$ , correlation length $l_c=0.68$ .

Figure 11. [24] The radial distribution functions of solid (T = 50 K), liquid (T = 80 K), and gaseous argon (T = 300 K). The radii are given in reduced units of the molecular diameter (σ = 3.822Å).

Approximate parameters of radial distribution function

$f_{rdf}\left(r,\theta ,\sigma ;c_k\right)=\Theta \left(r,\sigma \right)\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)f_{\theta }\left(\theta ,{\alpha }_{0k}\right)$

gas: $\alpha =0.5$ , extinction factor extfactor (r = 0.75/8) = 1.37, $l_c=0.75/8/0.318=0.29$ , $l_a=0.75$

fluid: $\alpha =2$ , period $l_a=0.75$ , extfactor(r = l_a/2) = 4.54, correlation length $l_c=l_a/1.49=0.50$

solid fcc: $r=1:{\alpha }_1+{\alpha }_2=5.3$ , $r=1.37:-{\alpha }_1-0.83{\alpha }_2=0.25$ , period $l_{a1}=0.75$ , $l_{a2}=0.63$ , extfactor(r = l_a) = 2.08, $l_c=l_{a1}/\log \left(2.08\right)=1.36l_{a1}=1.02$ , ${\alpha }_2=1.47$ , ${\alpha }_1=3.83=2.6{\alpha }_2$

solid fcc theory: $a=5.47\dot{\text{A}}=1.61{\sigma }_{LJ}$, $l_{a1}=a/2=0.80{\sigma }_{LJ}$ , $l_{a2}=a/\left(2\sqrt{2}\right)=0.57{\sigma }_{LJ}$ , ${\alpha }_1=n\left(r_1\right)/2=3$ , ${\alpha }_2=n\left(r_2\right)/2=2$

solid bcc theory: $a=5.47\dot{\text{A}}=1.61{\sigma }_{LJ}$, $l_{a1}=a\sqrt{3}/2=0.86a=1.39{\sigma }_{LJ}$ , ${\alpha }_1=n\left(r_1\right)/2=4$ .

6.3. Lennard-Jones Phase Transition Points and Parameters Argon

The used parameters are:

$\epsilon =0.0123\text{eV}$ , $\sigma =3.4\dot{\text{A}}$, $p_0\left(Ar\right)=\frac{\epsilon }{{\sigma }_{LJ}^3}=50.1\text{MPa}$ , $b=\frac{2\pi }{3}{\sigma }^3$ , $a=\frac{e}{\pi \left(3e-8\right)}\epsilon {\sigma }^3=5.59\epsilon {\sigma }^3$

Ar tr vapor-liquid-solid triple point

$T=83.7\text{K}$ , $\beta \epsilon =1.67$ ,$\lambda =1.05\sigma$ , $p=1.36\times {10}^{-3}{p}_0$ ,

$\rho =\frac{1}{v}=\frac{1}{{\lambda }^3}=0.87{\sigma }^{-3}$

fluid 80K: $\alpha =2$ , $l_a=0.75$ , $l_c=0.50$ .

Ar cr critical point

$T=150.8\text{K}$ , $\beta \epsilon =0.978$ ,$\lambda =1.47\sigma$ , $p=0.0964p_0$ , $\rho =\frac{1}{v}=\frac{1}{{\lambda }^3}=0.317{\sigma }^{-3}$

$\alpha =0.5$ , $l_c=0.29$ , $l_a=0.75$

Ar gas 300 K, 1 bar

$v=1/\rho =3.7\times {10}^{-20}{\text{cm}}^3$ , $\lambda =v^{1/3}=26\times {10}^{-7}\text{cm}=26\dot{\text{A}}$, $l_{free}=52\dot{\text{A}}$

Ar phase transitions

fluid-->gas at saturation curve:

$\lambda ={\lambda }_f\left(v_1,T\right)=\left\{\mathrm{0.65...0.9}\right\}\sigma \to \lambda ={\lambda }_g\left(v_2,T\right)$ determined by Maxwell’s equal-area rule

$F\left(v_1,T\right)=F\left(v_2,T\right)$ , $p\left(v_1,T\right)=p\left(v_2,T\right)$ where $p\left(v,T\right)=-\frac{\partial F\left(v,T\right)}{\partial V}$

fluid-->solid at fluid-solid line:

$\lambda ={\lambda }_f\left(v_1,T\right)=\left\{\mathrm{0.65...0.9}\right\}\sigma \to \lambda ={\lambda }_s\left(v_2,T\right)$ determined by Maxwell’s equal-area rule

$F\left(v_1,T\right)=F\left(v_2,T\right)$ , $p\left(v_1,T\right)=p\left(v_2,T\right)$ where $p\left(v,T\right)=-\frac{\partial F\left(v,T\right)}{\partial V}$ .

7. Dipole Substance

7.1. Dipole-dipole Interaction

• Parallel dipoles [25]

$V\text{}\left(r\right)=\frac{2{\mu }_1{\mu }_2}{r^3}$ in cgs, where $e^2=1.43\times {10}^{-9}\text{eV}\cdot \text{m}=14.3\text{eV}\cdot \dot{\text{A}}$ dipole moment ${\mu }_l=e{r}_l{q}_l$ , $l=1,2$

• Arbitrarily positioned dipoles [25] with angle θ₁ resp. θ₂ to center connection line,

$V\left(r\right)=-\frac{e^2{r}_1{r}_2{q}_1{q}_2}{r^3}\left(-\cos \left({\theta }_1-{\theta }_2\right)+3\cos \left({\theta }_1\right)\cos \left({\theta }_2\right)\right)$ ,

when both dipoles perpendicular to connection line, inverse direction:

${\theta }_1=\pi /2$ , ${\theta }_2=-\pi /2$ , $V\left(r\right)=-\frac{e^2{r}_1{r}_2{q}_1{q}_2}{r^3}$

with thermal screening

$V\left(r\right)=-\frac{2{\left(e^2{r}_1{r}_2{q}_1{q}_2\right)}^2}{3r^6}\frac{1}{k_{B}T}$ in cgs

The singularity-free cutoff-dipole potential is

$V_{cDD}\left(r,\sigma ,\Delta r,\theta \right)=V_0\left({\Theta }_L\left(r,\sigma ,\Delta r\right)-\frac{{\Theta }_H\left(r,\sigma ,\Delta r\right)}{{\left(r/\sigma \right)}^3+\Delta r}\cos \left(\theta \right)\right)$

where ${\Theta }_H\left(r,\sigma ,\Delta r\right)=\left(1-\frac{1}{1+\exp \left(\frac{r-\sigma }{\Delta r\sigma }\right)}\right)$ is the soft-step-up-function,

and ${\Theta }_L\left(r,\sigma ,\Delta r\right)=\frac{1}{1+\exp \left(\frac{r-\sigma }{\Delta r\sigma }\right)}$ is the soft-step-down-function.

7.2. Dipole Substance Ethanol

The following Table 3 contains the relevant data for ethanol.

Table 3. Phase data ethanol [26].

Critical pressure	6.25	MPa = MN/m2
Critical temperature	513.9	K
Critical volume	169	cm3/mol
Triple point pressure	4.3 × 10−10	MPa = MN/m2
Triple point temperature	150	K
Spec. volume (liquid)	58.7	cm3/mol

Solid ethanol

Ethanol’s crystal structure can involve different conformations of the ethanol molecules, including gauche and trans. The molecules are held together by hydrogen bonds.

The crystal structure of ethanol at 87 K. C₂H₅OH, monoclinic mP, a = 5.377 (4), b = 6.882 (5), c = 8.255 (8) Ȧ,

β = 102.2 (1)˚, V = 298.6 Ȧ³ at 87 K, Z = 4, Dx = 1.025 g cm⁻³.

$\sigma \left(Eth\right)=2.6\dot{\text{A}}$, $l_{a1}=5.37\dot{\text{A}}\approx 2\sigma \left(Eth\right)=5.2\dot{\text{A}}$, $l_{a2}=6.88\dot{A}=1.28l_{a1}$, $l_{a3}=8.25\dot{\text{A}}=1.54l_{a1}$

${\alpha }_1={\alpha }_2={\alpha }_3=6$

The phase diagram of ethanol is shown in Figure 12.

Figure 12. [27] A schematic phase diagram of ethanol. TP, BP, and CP are the triple point, the boiling point, and the critical point, respectively. The gray region above CP represents supercritical ethanol SCE. The broken line within SCE is the ridge of density fluctuations.

The rdf function of ethanol is shown in the following Figure 13.

Figure 13. rdf function for ethanol = ethanol and water-ethanol interactions at three ethanol concentrations x_et = 0.2, 0.22, 0.5 at T = 300K [28].

7.3. Dipole Substance Ethanol Transition Points

Maximum energy OH hydrogen bonding potential is V₀ = 0.216 eV ([26] [29]),

hydrogen bonding dipole interaction potential is $V\left(r,\sigma \right)=-\frac{V_0}{{\left(\frac{r}{\sigma }\right)}^3}$ with the

non-singular cutoff-dipole potential

$V_{cDD}\left(r,\sigma ,\Delta r\right)=V_0\left({\Theta }_L\left(r,\sigma ,\Delta r\right)-\frac{{\Theta }_H\left(r,\sigma ,\Delta r\right)}{{\left(r/\sigma \right)}^3+\Delta r}\right)$

where ${\Theta }_H\left(r,\sigma ,\Delta r\right)=\left(1-\frac{1}{1+\exp \left(\frac{r-\sigma }{\Delta r\sigma }\right)}\right)$ is the soft-step-up-function, and ${\Theta }_L\left(r,\sigma ,\Delta r\right)=\frac{1}{1+\exp \left(\frac{r-\sigma }{\Delta r\sigma }\right)}$ is the soft-step-down-function and the ε-parameter, in analogy to LJ-potential $\epsilon =V_0/4$

${\epsilon }_h\left(Eth\right)=V_0/4=0.0542\text{eV}=4.4\epsilon \left(Ar\right)$

The used parameters are:

$\sigma \left(Eth\right)=v^{1/3}=2.6\dot{\text{A}}$, $p_0=\frac{{\epsilon }_h\left(Eth\right)}{\sigma {\left(Eth\right)}^3}=\frac{0.0542\text{eV}}{5.2{}^{\dot{\text{A}}}}=3.85\times {10}^{-4}\text{eV}\cdot {}^{\dot{\text{A}}}=61.75\text{MPa}$, ${\epsilon }_h\left(Eth\right)=V_0/4=0.0542\text{eV}$

Critical point

$T=513\text{K}$ , $\beta =1/0.0445{\text{eV}}^{-1}$ $\beta {\epsilon }_h\left(Eth\right)=0.0542/0.0445=1.218$

$1\text{MPa}=6.242\times {10}^{-6}\text{eV}\cdot {}^{\dot{\text{A}}}$, $\text{eV}\cdot {}^{\dot{\text{A}}}=1.602\times {10}^5\text{MPa}$

$p_{cr}=6.25\text{MPa}=\frac{6.25}{61.75}{p}_0=0.101p_0$

$v_{cr}=169{\text{cm}}^{\text{3}}/\text{mol}=2.80\times {10}^{-22}{\text{cm}}^3=\frac{169\times {10}^{24}}{6.022\times {10}^{23}}{}^{\dot{\text{A}}}=280.6{}^{\dot{\text{A}}}$,

${\lambda }_{cr}={\left(v_{cr}\right)}^{1/3}=6.55\dot{\text{A}}=2.5\sigma$, ${\rho }_{cr}=1/{2.5}^3=0.064$

Triple point

$T=150\text{K}$ , $\beta =1/0.013{\text{eV}}^{-1}$ $\beta {\epsilon }_h\left(Eth\right)=0.0542/0.013=4.17$

$p_{tr}=4.3\times {10}^{-10}\text{MPa}=\frac{4.3\times {10}^{-10}}{61.75}{p}_0=0.696\times {10}^{-11}{p}_0$

$v_{fl}=58.7{\text{cm}}^{\text{3}}/\text{mol}=\frac{58.7\times {10}^{24}}{6.022\times {10}^{23}}{}^{\dot{\text{A}}}=9.75{}^{\dot{\text{A}}}={\left(2.13\dot{\text{A}}\right)}^3$, ${\lambda }_{fl}=2.13\dot{\text{A}}=0.48\sigma$,

$v_{fl}=587{\text{cm}}^{\text{3}}/\text{mol}=\frac{587\times {10}^{24}}{6.022\times {10}^{23}}{}^{\dot{\text{A}}}=97.5{}^{\dot{\text{A}}}={\left(4.6\dot{\text{A}}\right)}^3$, ${\lambda }_{fl}=4.6\dot{\text{A}}=1.05\sigma$

${\rho }_{fl}=1/0.11=9.$

melting (p=1bar) $T_m=159\text{K}$ , $\beta {\epsilon }_h\left(Eth\right)=3.93$ .

8. Landau Theory for Solid-Fluid-Gas: Ansatz and Basic Properties

8.1. Ansatz

We start with a generalized Landau ansatz for the partition function as a functional integral [1]

$Z={\displaystyle \int D\phi dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$

where $u\left(\overrightarrow{r}\right)$ is the intermolecular potential, and $\phi$ is a test function from a parameterized function family $\phi =\phi \left(\overrightarrow{r},\lambda ;c_k\right)$ with parameters $c_k$ , the average distance $\lambda =v^{1/3}$ , and the inverse thermal energy $\beta =\frac{1}{kT}$ .

The functional integral is expanded with differentials: volume $dV=2\pi r^2\sin \theta d\theta dr$ (assuming cylindrical symmetry) and functional

$D\phi ={\displaystyle \sum \frac{\partial \phi }{\partial c_k}}d{c}_k$ , and integrated accordingly.

The test function family must encompass radial distribution functions (rdf) $f_{rdf}\left(\overrightarrow{r};c_k\right)$ , since they describe the distribution of available states, we have the

ansatz $\phi =\phi \left(\overrightarrow{r};c_k\right)=f_1\left(\lambda ,\beta ;c_k\right)f_{rdf}\left(\overrightarrow{r};c_k\right)$

where $f_1\left(\overrightarrow{r};\lambda ,\beta ;c_k\right)$ is a function of the average distance $\lambda =v^{1/3}$ (resp. average volume per particle v), the inverse thermal energy $\beta$ , and the parameters $c_k$ .

The potential has the form $u\left(\overrightarrow{r},\sigma ,u_0\right)$ , where $\sigma$ is the particle radius, resp. the repulsive hardcore-radius, and $u_0=4\epsilon$ is the maximum potential energy, $\epsilon$ is the characteristic energy in the Lennard-Jones potential.

The simplest form for $f_1\left(\overrightarrow{r};\lambda ,\beta ;c_k\right)$ is exponential damping with $\lambda$ (states outside the average distance are of no importance) : $f_1\left(\overrightarrow{r};\lambda ,\beta ;c_k\right)=\exp \left(-r/\lambda \right)$ , it is shown to be correct, since it yields the ideal gas law in the limit of zero interaction u = 0.

The form of the rdf function is based on the known form for the three aggregates (gas, fluid, solid)

$f_{rdf}\left(r,\theta ,\sigma ;c_k\right)={\Theta }_H\left(r,\sigma \right)\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)f_{\theta }\left(\theta ,{\alpha }_{0k}\right)$ , $k=1,\cdots ,n_k$ (30a)

where $l_c$ is the correlation length, $l_{a,k}$ are the lattice constants of the underlying solid ($n_k=1,2,3$ or more, depending on the molecule symmetry), ${\alpha }_k$ are the corresponding amplitudes, ${\Theta }_H\left(r,\sigma \right)$ is step-up function,

$\phi \left(r,\theta ,\lambda ,\sigma ;c_k\right)=\exp \left(-\frac{r}{\lambda }\right)f_{rdf}\left(r,\theta ,\sigma ;c_k\right)$

The angular factor $f_{\theta }\left(\theta ,{\alpha }_{0k}\right)$ is the axial part, which is introduced for a non-radial (axially symmetric) potential, e.g. for the dipole potential, where θ is the polar angle, and the angle factor becomes simply$f_{\theta }\left(\theta ,{\alpha }_{0k}\right)=\left({\alpha }_{00}+{\alpha }_{01}\cos \theta \right)$ .

The idealized form of the rdf function for the three aggregates is then (radial part)

• gas $l_c=0$ ,

$f_{rdf}\left(r,\sigma ;c_k\right)=\Theta \left(r,\sigma \right)$ (30b)

pure unitstep-up function Θ, no correlation

• fluid $0<l_c~\sigma$ ,

$f_{rdf}\left(r,\sigma ;c_k\right)=\Theta \left(r,\sigma \right)\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)$ (30c)

• crystalline solid $l_c=\infty$ ,

$f_{rdf}\left(r,\sigma ;c_k\right)=\Theta \left(r,\sigma \right)\left(1+\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)$ (30d)

pure harmonic (periodic with lattice constants $l_{a,k}$ ).

With this ansatz we obtain the general Landau function

$\begin{array}{l}{\phi }_s\left(r,\theta ,\sigma ,\lambda ,l_c,{\alpha }_k,l_{a,k},{\alpha }_{0k}\right)\\ ={\Theta }_H\left(r,\sigma \right)\exp \left(-r/\lambda \right)\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)f_{\theta }\left(\theta ,{\alpha }_{0k}\right)\end{array}$ (30e)

with parameters $c_k=\left(l_c,{\alpha }_{k_1},l_{a,k_1},{\alpha }_{0k_2}\right)$ .

We can skip the step-function in the general Landau function, if we limit the range of r to $r\ge \sigma$ , so the general Landau function becomes

$\begin{array}{l}\phi \left(r,\theta ,\sigma ,\lambda ,l_c,{\alpha }_k,l_{a,k},{\alpha }_{0k}\right)\\ =\exp \left(-r/\lambda \right)\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)f_{\theta }\left(\theta ,{\alpha }_{0k}\right)\end{array}$ (30f)

The partition function becomes then

$\begin{array}{l}Z\left(\lambda ,\beta ,\sigma ,\epsilon ;c_k\right)\\ ={\displaystyle \int \left({\displaystyle \underset{r=\sigma }{\overset{\lambda }{\int }}{\displaystyle \underset{\theta =0}{\overset{\pi }{\int }}2\pi r^2\sin \theta d\theta dr{\displaystyle \sum _k\frac{\partial \phi \left(r,\theta ,\sigma ,\lambda ;c_k\right)}{\partial c_k}\exp \left(-\beta u\left(r,\theta ,\sigma ,\epsilon \right)\right)}}}\right)d{c}_{k}d\lambda }\end{array}$ (31)

We can set $\sigma =1,\epsilon =1$ , and measure radius in $unit\left(r\right)=\sigma$ , and energy in $unit\left(u\right)=\epsilon$ , and skip $\sigma ,\epsilon$ in the following.

We obtain the minimal parameters $c_k=\left(l_c,{\alpha }_{k_1},l_{a,k_1},{\alpha }_{0k_2}\right)$ by minimization of free energy $F=-\frac{1}{\beta }\log Z$ or equivalently by maximization of Z.

The minimal parameters $c_k$ determine the rdf $f_{rdf}\left(r,\theta ,\sigma ;c_k\right)$ in the partition function.

The other method of fixing the parameters is to solve the Hypernetted-Chain-Ornstein-Zernicke equation in order to find the rdf (see B.3).

The rdf function is uniquely determined by the intermolecular potential $u\left(r,\theta ,\sigma ,\epsilon \right)$ , and the same is valid for the rdf function $f_{rdf}\left(r,\theta ,\sigma ;c_k\right)$ in Z resulting from the maximization of Z, therefore the two resulting parameter sets $c_k$ should be identical. This is indeed the case for the two intermolecular potentials considered here: Lennard-Jones potential and dipole potential, as shown in chap. 12 and chap. 13.

8.2. Calculation

Calculation methods

In thermodynamics, the main task is the evaluation of multiple indefinite integrals with parameters.

Such integrals are in most cases not solvable in closed form, i.e. as expressions in parameters and variables using elementary functions (e.g. sin, exp) and special functions (e.g. Bessel functions, Gamma function).

The usual ansatz with a cut-off power series expansion (for exp) does not work, because the Boltzmann exponent $\beta u\left(r,\theta \right)$ is not necessarily small.

Basically, two calculation methods are applicable here.

• Numerical evaluation on parameter lattice

Here, the integral is evaluated pointwise on a multidimensional lattice for parameters c_k and variables ($r,\theta ,\lambda ,\beta$ ) and then interpolated or fitted (in polynomials or trigonometric functions) to obtain a closed expression. The calculation time grows as $N^{n_v+n_p}$ , where N is the number of lattice points in one dimension, n_v and n_p is the number of variables (here n_v = 4) and the number of parameters (here n_p = 5, ..., 8). In practice, one needs lattices with N = 400 - 1000 points in order to obtain reliable results. This is viable only with massive parallel processing on supercomputers.

• Numeric-symbolic evaluation on parameter lattice

Here, one replaces the parametric indefinite integral by a symbolic quadrature sum on a linear equidistant lattice, using some quadrature formula, e.g. Clenshaw-Curtis quadrature with weights w_k. An indefinite integral is represented by multiplying each summand by a stepup-function $\Theta \left(X-x\right)$ , where x is the integration variable and X is the upper limit variable.

In this way, the summation is variable, the lattice points $x_k>X$ do not contribute.

The transformation scheme is for the simplest quadrature rule:

${\displaystyle \underset{x=x_0}{\overset{X}{\int }}f\left(x;p_k\right)}dx\to {\displaystyle \sum _{j=1}^{j=N-1}{w}_j\Theta \left(X-x_j\right)f\left(\frac{x_{j+1}+x_j}{2};p_k\right)\frac{x_{j+1}-x_j}{N}}$

We obtain a closed symbolic expression in parameters p_k and integration limit X, which is easy to handle, but needs a large memory.

The expression value converges with the simplest Clenshaw-Curtis quadrature

against the true value with an error $\Delta I\approx \frac{h^2}{2}m\left(\left|f^{\prime }\right|\right)$ , where h is the lattice stepsize

and $m\left(\left|f^{\prime }\right|\right)$ is the average value of the integrand derivative.

On a 3.5 GHz workstation with 12 parallel processes, the calculation time for $p_{cDD}\left(\lambda ,\beta \right)$ with a DD-potential is 100s, the calculation of $p_{cDD}\left(\lambda ,\beta \right)$ on a N = 80 lattice takes about 10 h, and needs about 30GB memory per process, total memory needed is about 400 GB.

Numerical discretization error

The numerical discretization error in λ and β is approximately equal to the stepsize of the lattice.

An estimated error is $\Delta \lambda \approx 0.05$ , $\Delta \beta \approx 0.1$ for the value range $1\le \lambda \le 2$ $1\le \beta \le 5$ , where the characteristic points lie, which yields average relative errors $\frac{\Delta \lambda }{\lambda }\approx 0.03$ , $\frac{\Delta \beta }{\beta }\approx 0.03$ .

The relative error of p can be estimated from the ideal gas eos $p=\frac{1}{\beta {\lambda }^3}$ , namely $\frac{\Delta p}{p}=3\frac{\Delta \lambda }{\lambda }+\frac{\Delta \beta }{\beta }=0.12$ , which is roughly in agreement with the observed numerical error.

9. Partition Function for Ideal Gas and vdWaals-Gas

For the ideal gas, the potential is zero, and the rdf is the unitstep function.

In the following, we set $\sigma =1,\epsilon =1$ , and measure radius in $unit\left(r\right)=\sigma$ , and energy in $unit\left(u\right)=\epsilon$ .

The Landau function for the ideal gas becomes simply

${\phi }_{idg}\left(r,\lambda \right)=\exp \left(-\frac{r}{\lambda }\right)$ , where we set for the volume (per particle) $v={\lambda }^3$ and the partition function is

$Z_{idg}\left({\lambda }_1\right)={\displaystyle \underset{\lambda =1}{\overset{{\lambda }_1}{\int }}\left({\displaystyle \underset{r=1}{\overset{r=\lambda }{\int }}\frac{\partial {\phi }_{idg}\left(r,\lambda \right)}{\partial \lambda }4\pi r^2}dr\right)}d\lambda$ (32a)

It can be evaluated analytically with the result

$Z_{idg}\left(\lambda \right)={\left(-\frac{64\pi }{3e}{\lambda }^3+4\exp \left(-\frac{1}{\lambda }\right)\pi \lambda \left(1+2\lambda +2{\lambda }^2\right)\right)|}_1^{\lambda }$ (32b)

For the pressure we obtain

$\begin{array}{c}{p}_{idg}\left(\lambda ,\beta \right)=\frac{1}{3\beta {\lambda }^2}\frac{\partial \log Z_{idg}\left(\lambda \right)}{\partial \lambda }\\ =\frac{1}{\beta {\lambda }^3}\frac{\left(6e-16\exp \left(1/\lambda \right)\right)+{\lambda }^{-2}e\left(1/\lambda +3+6\lambda \right)}{\left(6e-16\exp \left(1/\lambda \right)\right)+{\lambda }^{-2}e\left(3+6\lambda \right)}\end{array}$ (33a)

and the limit $\lambda \to \infty$ ,

$p_{1c}\left(\lambda ,\beta \right)=\frac{1}{\beta {\lambda }^3}\left(1+\frac{1}{{\lambda }^3}\frac{e}{6e-16}\right)\approx \frac{1}{\beta \left({\lambda }^3-\frac{e}{6e-16}\right)}=\frac{1}{\beta \left(v-v_0\right)}$ (33b)

which is the ideal gas law with a van-der-Waals self-volume $v_0=\frac{e}{6e-16}=8.78$ .

The naive expectation value for the self-volume is a cube with side length $2\sigma =2$ , and self-volume $v_{00}=2^3=8$ , so the naive assessment is quite good.

9.1. vdWaals Gas with Lennard-Jones Potential

A vdWaals gas is ideal gas with a (weak) potential correction.

The pressure obeys the vdWaals eos

$p=\frac{1}{\beta \left(v-b\right)}-\frac{a}{v^2}$

In case of the Lennard-Jones potential we have the partition function

$\phi \left(r,\lambda \right)=\exp \left(-\frac{r}{\lambda }\right)$

$Z_{LJ}\left({\lambda }_1\right)={\displaystyle \underset{\lambda =1}{\overset{{\lambda }_1}{\int }}\left({\displaystyle \underset{r=1}{\overset{r=\lambda }{\int }}\frac{\partial \phi \left(r,\lambda \right)}{\partial \lambda }4\pi r^2\exp \left(-\beta u_{cLJ}\left(r,1,1\right)\right)dr}\right)}d\lambda$ (34)

and the pressure

$p_{LJ}\left(\lambda ,\beta \right)=\frac{1}{3\beta {\lambda }^2{Z}_{LJ}\left(\lambda ,\beta \right)}\frac{\partial Z_{LJ}\left(\lambda \right)}{\partial \lambda }$ (35)

which can calculated numerically on a λ-β-lattice [13].

The ideal gas invariant $E_{ig}\left(p,\lambda ,\beta \right)=p\beta {\lambda }^3$ is equal 1, for the vdWaals gas it reads

$E_{vdW}\left(p,\lambda ,\beta \right)=\frac{1}{1-b/{\lambda }^3}-\frac{a\beta }{{\lambda }^3}$ (36)

For the vdWaals-LJ gas-fluid it has the form (Figure 14).

Figure 14. State equation of vdWaals-Lennard-Jones gas-fluid.

We fit $E$ with the vdWaals eos $E_{vdW}\left(p,\lambda ,\beta \right)$ and obtain the result for the vdWaals parameters

$a=15.6$ , $b=29$ .

9.2. vdWaals Gas with Dipole Potential

In case of the dipole potential we have the partition function

$Z_{DD}\left({\lambda }_1\right)={\displaystyle \underset{\lambda =1}{\overset{{\lambda }_1}{\int }}\left({\displaystyle \underset{\theta =0}{\overset{\theta =\pi }{\int }}{\displaystyle \underset{r=1}{\overset{r=\lambda }{\int }}\frac{\partial \phi \left(r,\lambda \right)}{\partial \lambda }2\pi r^2\sin \theta \exp \left(-\beta u_{cDD}\left(r,1,1,\theta \right)\right)}drd\theta }\right)}d\lambda$ (37)

and the pressure

$p_{DD}\left(\lambda ,\beta \right)=\frac{1}{3\beta {\lambda }^2{Z}_{DD}\left(\lambda ,\beta \right)}\frac{\partial Z_{DD}\left(\lambda \right)}{\partial \lambda }$ (38)

which can calculated numerically on a λ-β-lattice [13].

For the vdWaals-dipole gas-fluid the invariant E has the form [13] (Figure 15).

Figure 15. State equation of vdWaals-dipole gas-fluid.

Fitting E with $E_{vdW}\left(p,\lambda ,\beta \right)$ yields the result for the vdWaals parameters

$a=48.6$ , $b=92.7$ .

10. Calculation of the Radial Distribution Function from Hypernetted-Chain-Ornstein-Zernicke Equation

10.1. Potentials

The general potential function has the form $u\left(r,\sigma ,\epsilon \right)$ , where σ is the hard-core radius, and ε is the characteristic energy, and $u_0=4\epsilon$ is the maximal energy.

The hard-core part of the potential is made non-singular and step-like by a cut-off.

We have the following typical intermolecular potentials

• cut-off Lennard-Jones potential (LJ)

$u_{cLJ}\left(r,\sigma ,\epsilon \right)=4\epsilon \left(\frac{1}{\frac{1}{2}+{\left(\frac{r}{\sigma }\right)}^{12}}-\frac{1}{1+{\left(\frac{r}{\sigma }\right)}^6}\right)$ (Figure 16),

Figure 16. Plot $u_{cLJ}\left(r,1,1\right)$ .

• dipole potential (DD)

$u_{cDD}\left(r,\sigma ,\epsilon ,\theta \right)=4\epsilon \left({\Theta }_L\left(r,\sigma ,0.1\right)-\left(1-{\Theta }_L\left(r,\sigma ,0.1\right)\right)u_{DD}\left(r,\sigma ,\theta ,0.1\right)\right)$

with the soft-stepdown function ${\Theta }_L\left(r,\sigma ,\Delta r\right)=\frac{1}{1+\exp \left(\frac{r-\sigma }{\Delta r\sigma }\right)}$ , and the pure cut-off dipole potential $u_{DD}\left(r,\sigma ,\theta ,\Delta r\right)=\frac{\cos \theta }{\Delta r+{\left(\frac{r}{\sigma }\right)}^3}$ (Figure 17)

Figure 17. Plot $u_{cDD}\left(r,1,1,\pi /3\right)$ .

• Coulomb potential

$u_{cCB}\left(r,\sigma ,\epsilon \right)=4\epsilon \left({\Theta }_L\left(r,\sigma ,0.1\right)-\left(1-{\Theta }_L\left(r,\sigma ,0.1\right)\right)u_{CB}\left(r,\sigma ,0.1\right)\right)$

with the pure cut-off Coulomb potential $u_{CB}\left(r,\sigma ,\Delta r\right)=\frac{1}{\Delta r+\frac{r}{\sigma }}$ (Figure 18)

Figure 18. Plot $u_{cCB}\left(r,1,1\right)$ .

In analogy to the Lennard-Jones potential, we distinguish te maximum energy V₀, and the characteristic energy ε = V₀/4.

The intermolecular potential determines the physical properties (like phase transitions) and the structure (crystal form) of a substance.

We deal here mainly with three potentials and corresponding substances: Lennard-Jones (LJ), dipole (D), and ionic (C).

• LJ-substance with LJ-potential

A typical LJ-substance is argon, or more generally noble gases. Here we use argon as the standard representative of a LJ-substance.

• D-substance with dipole potential

Many organic fluids are dipole substances, a typical example are alcohols, the best known being methanol and ethanol with a hydrogen bond OH-dipole potential (OH-bond energy V₀ ≈ 0.26eV) [27] [29].

Water is a special case, since it has two potentials: OH-dipole potential and covalent LJ-potential; accordingly, water consists of two liquids: high-density covalent-bound phase and low-density OH-dipole-bound phase with tetraedric structure (this phase is dominant in water ice) [30].

Here we use ethanol and methanol as the standard representatives of a dipole-substance.

For solid methanol we have the crystal structure:

-solid orthorombic: different xyz-lattice constants $l_{a1},l_{a2},l_{a3}$ , all angles ${\beta }_k=\pi /2$

For solid ethanol, the z-axis is slightly skewed, we have the crystal structure

-solid monoclinic: different xyz-lattice constants $l_{a1},l_{a2},l_{a3}$ , ${\beta }_3=102˚$

Since the solid ethanol structure is approximately orthorombic, we assume for the solid dipole-substance the orthorombic crystalline structure.

• ionic substance with Coulomb potential

Typical ionic substances are salts, e.g. sodium chloride NaCl. They have much higher melting points, because the Coulomb potentiall is much stronger than the dipole potential.

10.2. Phase Transitions

For the three potentials represented by three substances

• Lennard-Jones substance argon

• dipole substance ethanol

• ionic substance sodium chloride NaCl

we obtain the following phase transition data.

• Lennard-Jones substance Argon

$\epsilon =0.0123\text{eV}$ , $\sigma =3.4\dot{\text{A}}$, $p_0\left(Ar\right)=\frac{\epsilon }{{\sigma }_{LJ}^3}=50.1\text{MPa}$

Triple point

$T=83.7\text{K}$ , $\beta \epsilon =1.67$ ,$\lambda =1.05\sigma$ , $p=1.36\times {10}^{-3}{p}_0$ ,

$\rho =\frac{1}{v}=\frac{1}{{\lambda }^3}=0.87{\sigma }^{-3}$

Critical point

$T=150.8\text{K}$ , $\beta \epsilon =0.978$ , $\lambda =1.47\sigma$ ,$p=0.0964p_0$ , $\rho =\frac{1}{v}=\frac{1}{{\lambda }^3}=0.317{\sigma }^{-3}$

• Dipole substance ethanol

$\sigma \left(Eth\right)=v^{1/3}=2.6\dot{\text{A}}$, $p_0=\frac{{\epsilon }_h\left(Eth\right)}{\sigma {\left(Eth\right)}^3}=\frac{0.065\text{eV}}{17.6{}^{\dot{\text{A}}}}=3.7\times {10}^{-3}\text{eV}\cdot {}^{\dot{\text{A}}}=591\text{MPa}$,

${\epsilon }_h\left(Eth\right)=V_0/4=0.065\text{eV}$ , $V_0=0.26\text{eV}=V\left(OH\right)$

Critical point

$T=513\text{K}$ , $\beta {\epsilon }_h\left(Eth\right)=0.0542/0.0445=1.46$

$p_{cr}=6.25\text{MPa}=\frac{6.25}{591}{p}_0=0.0105p_0$

$v_{cr}=280.6{}^{\dot{\text{A}}}$, ${\lambda }_{cr}={\left(v_{cr}\right)}^{1/3}=6.55\dot{\text{A}}=2.5\sigma$

Triple point

$T=150\text{K}$ , $\beta {\epsilon }_h\left(Eth\right)=0.065/0.013=5.0$

$p_{tr}=4.3\times {10}^{-10}\text{MPa}=\frac{4.3\times {10}^{-10}}{591}{p}_0=0.73\times {10}^{-14}{p}_0$

$v_{fl}=9.75{}^{\dot{\text{A}}}$, ${\lambda }_{fl}=2.13\dot{\text{A}}=0.85\sigma$, ${\lambda }_{tr}\approx 2.5\dot{\text{A}}\approx 1\sigma$

• Ionic substance sodium chloride

$V_0\left(\text{NaCl}\right)=4.27\text{eV}$ , $\epsilon \left(\text{NaCl}\right)=V_0/4=1.067\text{eV}$ , $\sigma \left(\text{NaCl}\right)=2.8\dot{\text{A}}$

Triple point

$T_{tr}\approx T_m=1074\text{K}$ , $\beta =1/0.00721\text{eV}$ , $\beta \epsilon =1.067/0.00721=148$

Critical point

$T_{cr}\approx T_b=1738\text{K}$ , $\beta =1/0.1726\text{eV}$ , $\beta \epsilon =1.067/0.1726=6.18$

10.3. Calculation of the Radial Distribution Function for Lennard-Jones Fluid

The Hypernetted-Chain-Ornstein-Zernicke equation (HCOZ) reads

$\begin{array}{l}\log g\left(r_0\right)+\beta u\left(r_0\right)\\ =\rho \left(1-g\left(r_0\right)\right){\displaystyle \int \left(\left(1-g\left(R\right)\right)-\log g\left(R\right)-\beta u\left(R\right)\right)}{r}_1^{2}d{r}_1\sin \theta d\theta d\varphi \end{array}$

where $R^2=r_1^2{\sin }^2\theta +R_1^2=r_0^2+r_1^2+2r_0{r}_1\cos \theta \cos \varphi$

concisely $\log g\left(r_0\right)+\beta u\left(r_0\right)-\rho \left(1-g\left(r_0\right)\right)I\left(g,\beta u\right)=0$

$I\left(g,\beta u\right)={\displaystyle \int \left(\left(1-g\left(R\right)\right)-\log g\left(R\right)-\beta u\left(R\right)\right)}{r}_1^{2}d{r}_1\sin \theta d\theta d\varphi$

We make the general ansatz for the radial distribution function (rdf)

$g\left(r,\sigma ,lc,{\alpha }_k,l{a}_k\right)={\Theta }_H\left(r,\sigma ,\Delta r\right)\left(1+\Delta g\left(r,\sigma ,lc,{\alpha }_k,l{a}_k\right)\right)$

$\Delta g\left(r,\sigma ,lc,{\alpha }_k,l{a}_k\right)=\exp \left(-\frac{r}{lc}\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l{a}_k}\right)}\right)$

where ${\Theta }_H\left(r,\sigma ,\Delta r\right)=\left(1-\frac{1}{1+\exp \left(\frac{r-\sigma }{\Delta r\sigma }\right)}\right)$ is the soft-step-up-function with $\Delta r\ll \sigma$ .

Here lc is the damping length and ${\alpha }_k,l_{a_k}$ describe the amplitudes and the periods of the harmonics involved in the rdf.

The harmonics are determined by the crystalline geometry of the solid.

Theoretical values are (lattice unit a, normally $a=\sigma$ )

• solid fcc: $l_{a1}=a/2$ , $l_{a2}=a/\left(2\sqrt{2}\right)$ , ${\alpha }_1=n\left(r_1\right)/2=3$ , ${\alpha }_2=n\left(r_2\right)/2=2$

• solid bcc: $l_{a1}=a\sqrt{3}/2=0.86a$ , ${\alpha }_1=n\left(r_1\right)/2=4$

• solid hdp (hexagonal dense pack): $l_{a1}=a$ , neighbors $n\left(r_1\right)=6$ , ${\alpha }_1=n\left(r_1\right)/2=3$

• solid pc (primitive cubic):): $l_{a1}=a$ , $l_{a2}=a\sqrt{2}$ , neighbors $n\left(r_1\right)=6$ $n\left(r_2\right)=4$ , ${\alpha }_1=n\left(r_1\right)/2=3$ , ${\alpha }_2=n\left(r_2\right)/2=2$

For Lennard-Jones, we choose the fcc-ansatz, i.e. n = 2 harmonics, ${\alpha }_k,l_{a_k}$ $k=1,2$

The integral becomes

$I\left(r,\sigma ,lc,{\alpha }_k,l{a}_k\right)={\displaystyle \int \left(\left(1-g\left(R\right)\right)-\log g\left(R\right)-\beta u\left(R\right)\right)}{r}_1^{2}d{r}_1\sin \theta d\theta d\phi$

we can set σ = 1, ε = 1, and unit(r) = σ, unit(energy) = ε, then we reformulate HC-equation

$\begin{array}{l}0=E\left(r_0,\beta ,\rho ,lc,{\alpha }_k,l{a}_k\right)\\ \equiv \log g\left(r_0,lc,{\alpha }_k,l{a}_k\right)+\beta u_0\left(r_0\right)-\rho \left(1-g\left(r_0,lc,{\alpha }_k,l{a}_k\right)\right)I\left(r_0,\beta ,lc,{\alpha }_k,l{a}_k\right)\end{array}$

as a minimization problem on lattice L(r)

${\displaystyle \sum _{r_{0i}\in L\left(r\right)}\left|E\left(r_{0i},\beta ,\rho ,lc,{\alpha }_k,l{a}_k\right)\right|}=\min \left(lc,{\alpha }_k,l{a}_k\right)$

where $u_0\left(r_0\right)=u\left(r_0,\sigma =1,\epsilon =1\right)$ with a function solution $g_L\left(r,\rho ,\beta ;c_k\right)$ with coefficients, which are functions of the parameters $\left(lc,{\alpha }_k,l{a}_k\right)$

$c_k\left(\beta ,\rho \right)=\left(lc\left(\beta ,\rho \right),{\alpha }_k\left(\beta ,\rho \right),l{a}_k\left(\beta ,\rho \right)\right)$

In the limit $r_0\to \infty$ we have

$E\left(\right)=\Delta g\left(r_0\right)+\beta u_0\left(r_0\right)-\rho \Delta g\left(r_0\right)\left(2\Delta g\left(r_0\right)+\beta u_0\left(r_0\right)\right)I_0\left(r_0\right)$

$0=\Delta g\left(r_0\right)\left(1-\rho \beta u_0\left(r_0\right)I_0\left(r_0\right)\right)+\beta u_0\left(r_0\right)-2\rho \Delta g{\left(r_0\right)}^2{I}_0\left(r_0\right)$

follows $\Delta g\left(r_0\right)\approx -\frac{\beta u_0\left(r_0\right)}{1-\rho \beta u_0\left(r_0\right)I_0\left(r_0\right)}$ , which is an extended limit for small $\rho \ll 1$ .

With the above ansatz for the Lennard-Jones radial distribution function with $\sigma =1,\epsilon =1$ .

$g\left(r,\sigma ,lc,{\alpha }_k,l{a}_k\right)$ , and error minimization of the HCOZ equation we obtain the numerical solution

$g_L\left(r,\beta ,\rho \right)$ , equivalently ${\tilde{g}}_L\left(r,\beta ,\lambda \right)$ with $\rho =1/{\lambda }^3$ .

The results for rdf function are the following [13].

The dependence on β is weak, the dependence on ρ is negligible (the sporadic bumps are numerical artefacts) (Figure 19):

Figure 19. rdf function for LJ-substance in dependence of $\beta$ and $\rho$ for r = 1.4, HCOZ equation.

A typical r-profile is (here fluid near the critical point $\beta =1$ , $\rho =0.4$ , $\lambda =1.35$ ) (Figure 20).

Figure 20. rdf function r-profile for $\beta =1$ , $\rho =0.4$ ,$\lambda =1.35$ .

As in other calculations, there is a sharp minimum at $r=1.5\sigma$ .

10.4. Calculation of the Radial Distribution Function for Dipole Fluid

The Hypernetted-Chain-Ornstein-Zernicke equation (HCOZ) reads

$\begin{array}{l}\log g\left(r_0,{\theta }_0\right)+\beta u\left(r_0,{\theta }_0\right)\\ =\rho \left(1-g\left(r_0,{\theta }_0\right)\right){\displaystyle \int \left(\left(1-g\left(R,{\theta }_0\right)\right)-\log g\left(R,{\theta }_0\right)-\beta u\left(R,{\theta }_0\right)\right)}{r}_1^{2}d{r}_1\sin {\theta }_{1}d{\theta }_{1}d\varphi \end{array}$

where $R^2=r_1^2{\sin }^2{\theta }_1+R_1^2=r_0^2+r_1^2+2r_0{r}_1\cos {\theta }_1\cos \varphi$

concisely $\log g\left(r_0,{\theta }_0\right)+\beta u\left(r_0,{\theta }_0\right)-\rho \left(1-g\left(r_0,{\theta }_0\right)\right)I\left(g,\beta u\right)=0$

$I\left(g,\beta u\right)={\displaystyle \int \left(\left(1-g\left(R,\theta \right)\right)-\log g\left(R,\theta \right)-\beta u\left(R,\theta \right)\right)}{r}_1^{2}d{r}_1\sin {\theta }_{1}d{\theta }_{1}d\varphi$

We make the general ansatz for the radial distribution function (rdf)

$g\left(r,\sigma ,lc,{\alpha }_k,l{a}_k,\theta ,{\alpha }_0\right)={\Theta }_H\left(r,\sigma ,\Delta r\right)\left(1+\Delta g\left(r,\sigma ,lc,{\alpha }_k,l{a}_k,\theta ,{\alpha }_0\right)\right)$

$\Delta g\left(r,\sigma ,lc,{\alpha }_k,l{a}_k,\theta ,{\alpha }_0\right)=\exp \left(-\frac{r}{lc}\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l{a}_k}\right)}\right)\left(\cos \theta +{\alpha }_0\right)$

where ${\Theta }_H\left(r,\sigma ,\Delta r\right)=\left(1-\frac{1}{1+\exp \left(\frac{r-\sigma }{\Delta r\sigma }\right)}\right)$ is the soft-step-up-function with $\Delta r\ll \sigma$ .

Here lc is the damping length and ${\alpha }_k,l_{a_k}$ describe the amplitudes and the periods of the harmonics involved in the rdf, θ is the polar angle variable, α₀ is the amplitude of the spherically symmetric part of rdf.

We have here two-dimensional positional space with two variables $r,\theta$ .

The harmonics are determined by the crystalline geometry of the solid.

For the dipole potential, we make the rdf ansatz for orthorombic (see above), i.e. n = 3 harmonics, ${\alpha }_k,l_{a_k}$ $k=1,2,3$ .

With the above ansatz for the dipole radial distribution function with $\sigma =1,\epsilon =1$ $g\left(r,\sigma ,lc,{\alpha }_k,l{a}_k,\theta ,{\alpha }_0\right)$ , and error minimization of the HCOZ equation we obtain the numerical solution $g_D\left(r,\theta ,\rho ,\beta ;c_k\right)$ , equivalently ${\tilde{g}}_D\left(r,\theta ,\beta ,\lambda ;c_k\right)$ with $\rho =1/{\lambda }^3$ [13], with coefficients $c_k\left(\beta ,\rho \right)=\left(lc\left(\beta ,\rho \right),{\alpha }_k\left(\beta ,\rho \right),l{a}_k\left(\beta ,\rho \right),\theta \left(\beta ,\rho \right),{\alpha }_0\left(\beta ,\rho \right)\right)$ , which are functions of the parameters $\left(lc,{\alpha }_k,l{a}_k,\theta ,{\alpha }_0\right)$ .

The results for rdf function are as follows [13] (Figures 21-23).

The dependence on $\beta$ is significant for small values only, there is practically no dependence on $\rho$ , resp. $\lambda$ .

Figure 21. rdf function for D-substance in dependence of $\beta$ and $\rho$ for r = 1.2, θ = 0.78, HCOZ equation.

We show calculated specific rdf profiles for D-substance (average distance $\lambda$ in σ, inverse thermal energy $\beta$ in 1/ε).

• at critical point in $r,\theta$

Figure 22. rdf function r-θ-profile for $\beta =1.5$ , $\lambda =2.5$ .

• a typical rdf profile near critical point is ($\beta =1.1$ , $\rho =0.6$ , $\lambda =1.18$ )

Figure 23. rdf function r-profile for $\beta =1.1$ ,$\lambda =1.18$ , $\theta =0.78$ .

10.5. Minimization of Free Energy for Lennard-Jones Fluid

We obtain the minimal parameters $c_k=\left(l_c,{\alpha }_{k_1},l_{a,k_1}\right)$ by minimization of free energy $F=-\frac{1}{\beta }\log Z$ pointwise on $\left(\lambda ,\beta \right)$ lattice and fitting $c_k\left(\lambda ,\beta \right)$ in the form

$c_k\left(\lambda ,\beta \right)=\left(lc\left(\lambda ,\beta \right),{\alpha }_k\left(\lambda ,\beta \right),l{a}_k\left(\lambda ,\beta \right)\right)$

Here $\lambda$ is the average distance and $\beta =\frac{1}{k_{B}T}$ is the Boltzmann factor.

The minimization on lattice is very time-consuming, so it was made with low relative precision (only 0.01) in order to reduce calculation time. The calculation time for a 81 × 81 point lattice with 12 parallel processes was t = 8 h, with memory 12 × 7 + 30 GB.

The results for the parameters $c_k\left(\lambda ,\beta \right)$ are in the same range as from the HCOZ calculation, but differ somewhat in the $\left(\lambda ,\beta \right)$ -dependence [13] (Figure 24, Figure 25).

Figure 24. rdf function for LJ-substance in dependence of $\beta$ and $\lambda$ for r = 1., free energy minimization.

However, the rdf function near the saturation curve and the critical point are very similar, as well as the results for the free energy F.

A typical r-profile is (here fluid near the critical point $\beta =1$ , $\rho =0.4$ , $\lambda =1.35$ ), compare the HCOZ profile in chap. 10.3.

Figure 25. rdf function r-profile for $\beta =1$ , $\lambda =1.35$ .

10.6. Minimization of Free Energy for Dipole Fluid

[13] We obtain the minimal parameters $c_k=\left(l_c,{\alpha }_{k_1},l_{a,k_1},{\alpha }_0\right)$ by minimization of

free energy $F=-\frac{1}{\beta }\log Z$ pointwise on $\left(\lambda ,\beta \right)$ lattice and fitting $c_k\left(\lambda ,\beta \right)$ in

the form

$c_k\left(\lambda ,\beta \right)=\left(lc\left(\lambda ,\beta \right),{\alpha }_k\left(\lambda ,\beta \right),l{a}_k\left(\lambda ,\beta \right),{\alpha }_0\left(\lambda ,\beta \right)\right)$

Here $\lambda$ is the average distance and $\beta =\frac{1}{k_{B}T}$ is the Boltzmann factor.

The results for the parameters $c_k\left(\lambda ,\beta \right)$ are in the same range as from the HCOZ calculation, but differ somewhat in the $\left(\lambda ,\beta \right)$ behavior. However, the rdf function near the saturation curve and the critical point are very similar, as well as the results for the free energy F.

The minimization was made with low relative precision (only 0.01) in order to reduce calculation time. Still, the calculation time for a 81x81 point lattice with 12 parallel processes was t = 83h, with memory 12x8 + 31 GB.

The results for the parameters $c_k\left(\lambda ,\beta \right)$ are in the same range as from the HCOZ calculation, but differ somewhat in the $\left(\lambda ,\beta \right)$ behavior (Figure 26, Figure 27).

Figure 26. rdf function for D-substance in dependence of $\beta$ and $\lambda$ for r = 1.2, θ = 0.78, free energy minimization.

However, the rdf function near the saturation curve and the critical point are very similar, as well as the results for the free energy F.

A typical rdf profile near critical point is ($\beta =1.1$ , $\rho =0.6$ , $\lambda =1.18$ ), compare the HCOZ profile in chap. 10.4.

Figure 27. rdf function r-profile for $\beta =1.1$ , $\lambda =1.18$ , $\theta =0.78$ .

11. Partition Function and Free Energy with Calculated Radial Distribution Function

We start with the general Landau function ($\sigma =1,\epsilon =1$ )

$\begin{array}{l}\phi \left(r,\theta ,\lambda ,l_c,{\alpha }_k,l_{a,k},{\alpha }_{0k}\right)\\ =\exp \left(-r/\lambda \right)\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-1}{l_{a,k}}\right)}\right)\right)\left({\alpha }_{00}+{\alpha }_{01}\cos \theta \right)\end{array}$ (39a)

where we inserted the angular factor $f_{\theta }\left(\theta ,{\alpha }_{0k}\right)=\left({\alpha }_{00}+{\alpha }_{01}\cos \theta \right)$ of the dipole potential, and the parameters are $c_k=\left(l_c,{\alpha }_{k_1},l_{a,k_1},{\alpha }_{0k_2}\right)$ .

The general partition function is then

$\begin{array}{l}Z\left({\lambda }_1,\beta ;c_k\right)\\ ={\displaystyle \underset{c_k}{\int }{\displaystyle \underset{\lambda =1}{\overset{\lambda ={\lambda }_1}{\int }}\left({\displaystyle \underset{r=1}{\overset{\lambda }{\int }}{\displaystyle \underset{\theta =0}{\overset{\pi }{\int }}2\pi r^2\sin \theta d\theta dr{\displaystyle \sum _k\frac{\partial \phi \left(r,\theta ,\lambda ;c_k\right)}{\partial c_k}\exp \left(-\beta u\left(r,\theta \right)\right)}}}\right)d\lambda }d{c}_k}\end{array}$ (39b)

integrated over parameters $c_k$ in within a suitable range.

and the general free energy is $F\left(\lambda ,\beta ;c_k\right)=-\frac{1}{\beta }\log Z\left(\lambda ,\beta ;c_k\right)$ .

The fundamental thermodynamic are the average distance $\lambda =v^{1/3}$ (resp. average volume per particle v), and the inverse thermal energy $\beta =\frac{1}{kT}$ .

Now, we obtain the actual free energy $F\left(\lambda ,\beta \right)$ by minimization of $F\left(\lambda ,\beta ;c{}_k\right)$ in the parameters $c_k$ within a suitable parameter range.

Alternatively, we can solve the Hypernetted-Chain-Ornstein-Zernicke equation in order to find the correct rdf with the corresponding parameters (see B.3), and this is the procedure, which we choose here.

11.1. LJ Potential

The Lennard-Jones potential as the form [13] $u_{cLJ}\left(r\right)=4\left(\frac{1}{\frac{1}{2}+r^{12}}-\frac{1}{1+r^6}\right)$

We start with the general partition function above,

$Z\left(\lambda ,\beta ;c_k\right)={\displaystyle \int \left({\displaystyle \underset{r=1}{\overset{\lambda }{\int }}4\pi r^{2}dr{\displaystyle \sum _k\frac{\partial \phi \left(r,\theta ,\lambda ;c_k\right)}{\partial c_k}\exp \left(-\beta u_{cLJ}\left(r\right)\right)}}\right)d{c}_{k}d\lambda }$ (40a)

where $c_k=\left(l_c,{\alpha }_{k_1},l_{a,k_1}\right)$ and insert the parameters from the HCOZ solution for the Lennard-Jones potential from chap.10: $c_k=c_k\left(u_{cLJ}\right)$ .

The result is the partition function of the LJ-potential $Z_{cLJ}\left(\lambda ,\beta \right)$ and the corresponding free energy $F_{cLJ}\left(\lambda ,\beta \right)=-\frac{1}{\beta }\log Z_{cLJ}\left(\lambda ,\beta \right)$

A plot of $F_{cLJ}\left(\lambda ,\beta \right)$ is shown below [13], Figure 27.

At the left edge we see the “bulge” of the LJ-saturation curve fluid-gas, running between the triple point and the critical point (see chap.2) $1.\le \beta \le 1.6$ . It is much better recognizable in the pressure and equation-of-state plots in the next chap.12, chap.13.

Figure 27. A plot of free energy of Lennard-Jones substance $F_{cLJ}\left(\lambda ,\beta \right)$ in dependence of $\lambda$ and $\beta$ .

11.2. Dipole Potential

The dipole potential as the form [13]

$u_{cDD}\left(r,\sigma ,\epsilon ,\theta \right)=4\epsilon \left(step10\left(r,\sigma ,0.1\right)-\left(1-step10\left(r,\sigma ,0.1\right)\right)u_{DD}\left(r,\sigma ,\theta ,0.1\right)\right)$

with the soft-stepdown function ${\Theta }_L\left(r,\sigma ,\Delta r\right)=\frac{1}{1+\exp \left(\frac{r-\sigma }{\Delta r\sigma }\right)}$ , and the pure cut-off dipole potential $u_{DD}\left(r,\sigma ,\theta ,\Delta r\right)=\frac{\cos \theta }{\Delta r+{\left(\frac{r}{\sigma }\right)}^3}$ .

We start again with the general partition function above,

$\begin{array}{l}Z\left(\lambda ,\beta ;c_k\right)\\ ={\displaystyle \int \left({\displaystyle \underset{r=1}{\overset{\lambda }{\int }}{\displaystyle \underset{\theta =0}{\overset{\pi }{\int }}2\pi r^2\sin \theta d\theta dr{\displaystyle \sum _k\frac{\partial \phi \left(r,\theta ,\lambda ;c_k\right)}{\partial c_k}\exp \left(-\beta u_{cDD}\left(r,\theta \right)\right)}}}\right)d{c}_{k}d\lambda }\end{array}$ (40b)

where $c_k=\left(l_c,{\alpha }_{k_1},l_{a,k_1},{\alpha }_{0k_2}\right)$

and insert the parameters from the HCOZ solution for the dipole potential from chap.10: $c_k=c_k\left(u_{cDD}\right)$ .

The result is the partition function of the LJ-potential $Z_{cDD}\left(\lambda ,\beta \right)$ and the corresponding free energy $F_{cDD}\left(\lambda ,\beta \right)=-\frac{1}{\beta }\log Z_{cDD}\left(\lambda ,\beta \right)$

A plot of $F_{cDD}\left(\lambda ,\beta \right)$ is shown below [13] Figure 28.

At the left edge in front we see the DD-saturation curve fluid-gas, running between the triple point and the critical point (see chap.2) $1.5\le \beta \le 5$ , at the triple point the branches fluid-solid and solid-gas are clearly visible.

Figure 28. A plot of free energy of dipole substance $F_{cDD}\left(\lambda ,\beta \right)$ in dependence of $\lambda$ and $\beta$ .

12. Pressure Profiles

The pressure is derived from the partition function using $v={\lambda }^3$

$p\left(\lambda ,\beta \right)=-\frac{\partial F}{\partial v}=\frac{1}{\beta }\frac{\partial \log Z}{\partial v}=\frac{1}{\beta Z}\frac{1}{3{\lambda }^2}\frac{\partial Z}{\partial \lambda }$ after inserting the parameters from the HCOZ solution for the potential.

Since the calculation of $p\left(\lambda ,\beta \right)$ is time-consuming (up to about 20 s on a 3.5 GHz work station), it is calculated and evaluated as an array on an equidistant $\lambda \beta$ -lattice.

The formula for pressure simplifies, when we take into account the integration over λ in Z, and keep the dependence on parameters c_k:

$\begin{array}{l}p\left(\lambda ,\beta ;c_k\right)\\ =\frac{1}{\beta Z\left(\lambda ,\beta ;c_k\right)}\frac{1}{3{\lambda }^2}\left({\displaystyle \underset{r=1}{\overset{\lambda }{\int }}{\displaystyle \underset{\theta =0}{\overset{\pi }{\int }}2\pi r^2\sin \theta d\theta dr{\displaystyle \sum _k\frac{\partial \phi \left(r,\theta ,\lambda ;c_k\right)}{\partial c_k}\exp \left(-\beta u\left(r,\theta \right)\right)}}}\right)\end{array}$ (41)

12.1. LJ Potential

• Direct calculation of rdf function from HCOZ equation

We insert the parameters from the HCOZ solution for the Lennard-Jones potential from chap.10 $c_{c,k}=c_{c,k}\left(u_{cLJ}\right)$ into the above formula for $p\left(\lambda ,\beta ;c_{c,k}\right)$ .

The result is the pressure of the LJ-potential $p_{cLJ}\left(\lambda ,\beta \right)$ [13] (Figure 29, Figure 30).

Figure 29. A plot of the pressure of Lennard-Jones substance $p_{cLJ}\left(\lambda ,\beta \right)$ from HCOZ-eq. in dependence of $\lambda$ and $\beta$ .

Here the saturation curve fluid-gas in the range $\lambda \approx 1.6$ $1\le \beta \le 1.7$ is clearly recognizable.

The temperature range of the fluid is small $\Delta \beta =0.7$ , which is typical for LJ-fluids (e.g. for fluid argon with temperature range $\Delta T_f=67\text{K}$ ) because the covalent binding LJ-potential is so weak.

A typical saturation curve profile $p\left(\lambda ,\beta =const\right)$ at constant temperature is shown below ($\beta =1.375$ ).

Figure 30. Profile of the pressure of Lennard-Jones substance $p\left(\lambda ,\beta =1.375\right)$ in dependence of $\lambda$ .

The profile has in two regions the characteristic vdWaals-form “humps” with negative minimum and positive maximum in p.

The first hump is the solid-fluid transition at $\lambda =1.40$ , $p=0.052$ , the second hump is the fluid-gas transition at $\lambda =1.62$ , $p=0.014$

In the fluid-gas phase transition the average distance λ (or equivalently specific volume v) jumps $v_1={\lambda }_1^2\to v_2={\lambda }_2^3$ across the instable region at constant p, obeying the Maxwell rule (equal area below and above), i.e. $p\left(v_1,T\right)=p\left(v_2,T\right)$ , $F\left(v_1,T\right)=F\left(v_2,T\right)$ .

• Minimization of free energy F

We obtain the pressure of the LJ-potential $p_{mLJ}\left(\lambda ,\beta \right)$ [13] (Figure 31, Figure 32).

Figure 31. A plot of the pressure of LJ-substance $p_{cLJ}\left(\lambda ,\beta \right)$ from F-minimization in dependence of $\lambda$ and $\beta$ .

The strong bulge of the saturation curve at $\lambda =1.5$ is still there, compared to the HCOZ result.

A typical saturation curve profile $p\left(\lambda ,\beta =const\right)$ at constant temperature is shown below [13] ($\beta =1.375$ ).

Figure 32. Profile of the pressure of Lennard-Jones substance $p\left(\lambda ,\beta =1.375\right)$ in dependence of $\lambda$ .

12.2. Dipole Potential

• Direct calculation of rdf function from HCOZ equation

We insert the parameters from the HCOZ solution for the dipole potential from chap.10: $c_{c,k}=c_{c,k}\left(u_{cDD}\right)$ into the above formula for $p\left(\lambda ,\beta ;c{}_{c,k}\right)$ .

The result is the pressure of the dipole potential $p_{cDD}\left(\lambda ,\beta \right)$ [13] (Figure 33, Figure 34).

Figure 33. A plot of the pressure of dipole substance $p_{cDD}\left(\lambda ,\beta \right)$ from HCOZ-eq. in dependence of $\lambda$ and $\beta$ .

The saturation curve fluid-gas in the range $1.35\le \lambda \le 1.46$ $1.46\le \beta \le 5$ is recognizable at the left edge, followed by the solid-gas evaporation curve.

The temperature range of the fluid is $\Delta \beta \approx 3.5$ , which is much larger than for LJ-fluids (e.g. for ethanol with temperature range $\Delta T_f=363\text{K}$ ) because the dipole-potential is much stronger.

A typical saturation curve profile $p\left(\lambda ,\beta =const\right)$ at constant temperature is shown below [13] ($\beta =1.715$ ).

Figure 34. Profile of the pressure of dipole substance $p\left(\lambda ,\beta =1.715\right)$ in dependence of $\lambda$ .

As above, the profile has the characteristic vdWaals-form with negative minimum and positive maximum in p.

The transition fluid-gas is at $\lambda =1.29$ , $p=0.079$ .

• Minimization of free energy F

We obtain the pressure of the dipole potential $p_{mDD}\left(\lambda ,\beta \right)$ [13] (Figure 35, Figure 36).

Figure 35. A plot of the pressure of dipole substance $p_{cDD}\left(\lambda ,\beta \right)$ from F-minimization in dependence of $\lambda$ and $\beta$ .

The strong bulge of the saturation curve at $\lambda =1.5$ is still there, compared to the HCOZ result.

A typical saturation curve profile $p\left(\lambda ,\beta =const\right)$ at constant temperature is shown below [13] ($\beta =1.715$ ).

Figure 36. Profile of the pressure of dipole substance $p\left(\lambda ,\beta =1.715\right)$ in dependence of $\lambda$ .

13. Equation-of-state, Characteristic Points

On the saturation curve we have the following behavior:

fluid-->gas at saturation curve:

${\lambda }_1={\lambda }_f\left(\beta \right)\to {\lambda }_2={\lambda }_g\left(\beta \right)$ determined by Maxwell’s equal-area rule with $v={\lambda }^3$

$F\left({\lambda }_1,\beta \right)=F\left({\lambda }_2,\beta \right)$ , $p\left({\lambda }_1,\beta \right)=p\left({\lambda }_2,\beta \right)$ where $p\left(\lambda ,\beta \right)=-\frac{\partial F\left(\lambda ,\beta \right)}{\partial v}$ fluid-->solid at fluid-solid line, ${\lambda }_1={\lambda }_f\left(\beta \right)\to {\lambda }_2={\lambda }_s\left(\beta \right)$ determined by Maxwell’s equal-area rule $F\left({\lambda }_1,\beta \right)=F\left({\lambda }_2,\beta \right)$ , $p\left({\lambda }_1,\beta \right)=p\left({\lambda }_2,\beta \right)$ .

The saturation curve runs between the triple point and the critical point.

• Triple point

At the triple point, the saturation curve ends, fluid-->gas and fluid-->solid curves meet, i.e. $p\left({\lambda }_g,\beta \right)=p\left({\lambda }_f,\beta \right)=p\left({\lambda }_s,\beta \right)$ .

The triple point can be found as the branching point of saturation curve in the $p\left(\lambda ,\beta \right)$ diagram.

• Critical point

At the critical point, only one phase exists, the isotherm has a turning point there.

$\frac{\partial p}{\partial v}=0$ , $\frac{{\partial }^{2}p}{\partial v^2}=0$

The critical point can be determined from $p\left(\lambda ,\beta \right)$ by solving numerically the

two equations, or by minimizing numerically $w_1{\left(p\left(\lambda ,\beta \right)\right)}^2+w_2{\left(\frac{\partial p\left(\lambda ,\beta \right)}{\partial \lambda }\right)}^2$

with suitable weights $w_1,w_2$ .

Since the calculation of $p\left(\lambda ,\beta \right)$ is time-consuming, a numerical solution takes too much time [14].

It is by far preferable to inspect the $p\left(\lambda ,\beta =const\right)$ -profiles and $p\left(\lambda =const,\beta \right)$ -profiles and find the end of the saturation curve in the $p\left(\lambda ,\beta \right)$ diagram.

• Ideal gas invariant

For the ideal gas, the equation-of-state can be written in the form $p\beta {\lambda }^3=1$

That means, the function $E_c\left(\lambda ,\beta \right)=p\left(\lambda ,\beta \right)\beta {\lambda }^3$ is constant, $E_c\left(\lambda ,\beta \right)=1$

In general case of a gas-fluid, this function $E_c\left(\lambda ,\beta \right)$ , called here ideal-gas-invariant, provides important insight into its thermal behavior, and is equivalent to its equation-of-state (eos).

• Saturation curve (fluid-gas)

The saturation curve (boiling curve) runs between the triple point and the critical point. At the critical point, the fluid-gas boundary vanishes. The triple point is a branching point, where the saturation curve divides into the fluid-solid and the fluid-gas curve.

All boundary curves obey the condition $\frac{\partial p}{\partial v}=0$ . The critical point obeys additionally the condition $\frac{{\partial }^{2}p}{\partial v^2}=0$ .

The pressure $p\left(\lambda ,\beta \right)$ contains exponentials and harmonic polynomials, the condition is an algebraic-exponential equation where only the real roots are admissible, so the solution contains branching points, and endpoints (where a real solution ceases to exist).

Furthermore, at boundary curves the eos “jumps” over regions with negative pressure according to the Maxwell rule, and therefore there is latent heat, i.e. a jump in free energy.

• Critical exponents

We consider here only critical exponents of pressure and density, because higher derivatives like specific heat are not precise enough on a small lattice used here.

The critical exponent δ is defined by the summary formula

$\left|\Delta p\right|~{\left|\Delta \rho \right|}^{\delta },T\to T_{C-}$

It is calculated below for the LJ-fluid.

The critical exponent β₀ (not to be confounded with the inverse thermal energy $\beta =1/k_{B}T$ ) is defined by the summary formula

$\left|\Delta \rho \right|~{\left|\Delta T\right|}^{{\beta }_0},T\to T_{C-}$

It is calculated below for the dipole fluid.

13.1. LJ Potential

• Direct calculation of rdf function from HCOZ equation [13]

The ideal gas invariant $E_{cLJ}\left(\lambda ,\beta \right)=p_{cLJ}\left(\lambda ,\beta \right)\beta {\lambda }^3$ for the LJ-potential has the form (Figure 37).

Figure 37. Equation-of-state of LJ-substance from HCOZ-eq.

Clearly visible is the saturation curve fluid-gas in the range $\lambda \approx 1.6$ $1\le \beta \le 1.6$ , followed by the solid-gas evaporation curve in the range $1.7\le \beta \le 4$ .

Inspection of the p-profiles yields the following characteristic points.

• critical point

λ = 1.62, p = 0.016, β = 1.035, visible in the β = 1.035 profile, with the approximate turning point at λ = 1.84, measured values (argon): β_cr = 0.98, λ_cr = 1.47, p_cr

= 0.097p₀, with the characteristic pressure $p_0=\frac{\epsilon }{{\sigma }^3}=50.1\text{MPa}$ for argon.

In the following two λ = const-profiles one sees the changing profile form: the turning point disappears, signaling the end of the saturation curve (Figure 38).

Figure 38. Critical point pressure profiles of LJ-substance from HCOZ-eq.

• triple point

λ = 1.58, p = 0.004, β = 1.63, measured (argon): β_tr = 1.67, λ_tr ≈ 1.4, p_tr = 0.0014p₀.

The following two β = const profiles show the disappearance of the transition hump at λ = 1.58, Figure 39.

Figure 39. Triple point pressure profiles of LJ-substance from HCOZ-eq.

• saturation curve (fluid-gas)

For the HCOZ-solution we obtain the saturation curve from the pressure profile: see Figure 40.

Figure 40. Saturation curve of LJ-substance from HCOZ-eq.

• critical exponent δ

The fit to the curve $\Delta p\left(\lambda ,\beta ={\beta }_{cr}\right)$ at the critical point yields the value δ = 2.93, for the vdWaals fluid we have δ = 3., for the LJ-fluid δ = 4.8, see Figure 41.

Figure 41. Critical exponent curve of LJ-substance from HCOZ-eq.

• Minimization of free energy F [13]

We obtain for the ideal gas invariant $E_{cLJ}\left(\lambda ,\beta \right)=p_{mLJ}\left(\lambda ,\beta \right)\beta {\lambda }^3$ for the LJ-potential the form (Figure 42).

Figure 42. Equation-of-state of LJ-substance from F-minimization.

The comb of the saturation curve fluid-gas in the range $\lambda \approx 1.6$ is still there.

Inspection of the p-profiles yields the following characteristic points.

• critical point

λ = 1.56, p = 0.027, β = 1.035, visible in the β = 1.035 profile

• triple point

λ = 1.67, p = 0.0097, β = 1.63

13.2. Dipole Potential

• Direct calculation of rdf function from HCOZ equation [13]

The ideal gas invariant $E_{cDD}\left(\lambda ,\beta \right)=p_{cDD}\left(\lambda ,\beta \right)\beta {\lambda }^3$ for the DD-potential has the form (Figure 43).

Figure 43. Equation-of-state of dipole substance from HCOZ-eq.

At the left is the saturation curve fluid-gas in the range $1.35\le \lambda \le 1.46$ , $1.46\le \beta \le 5$ , branching off is the solid-gas evaporation curve in the range $5\le \beta \le 6$

Inspection of the p-profiles yields the following characteristic points.

λ = 1.46, β = 1.545 ± 0.1, (λ) = 0.225, measured (ethanol) p_cr = 0.051, λ = 1.49, β = 1.46.

The following two β = const profiles show the disappearance of the vdWaals turning point for β < 1.545 (Figure 44).

Figure 44. Critical point pressure profiles of dipole substance from HCOZ-eq.

• triple point

λ = 1.35, β = 5.16 + −0.2, p = 0.004

measured p_tr = 0.051, λ_tr = 1.3, β_tr = 5.

The following two β = const profiles show the disappearance of the transition hump above β = 5, signaling the end of the saturation curve (Figure 45).

Figure 45. Triple point pressure profiles of dipole substance from HCOZ-eq.

• saturation curve (fluid-gas)

For the HCOZ-solution we obtain the saturation curve from the pressure profile, see Figure 46.

Figure 46. Saturation curve for dipole substance from HCOZ-eq.

• critical exponent β₀

The fit to the curve $\Delta \rho \left(\lambda ={\lambda }_{cr},\beta \right)$ at the critical point yields the value β₀ = 0.46, for the vdWaals fluid we have

β₀ = 0.5, for the LJ-fluid β₀ = 0.33, see Figure 47.

Figure 47. Critical exponent curve for dipole substance from HCOZ-eq.

• Minimization of free energy F [13]

We obtain for the ideal gas invariant $E_{mDD}\left(\lambda ,\beta \right)=p_{mDD}\left(\lambda ,\beta \right)\beta {\lambda }^3$ for the DD-potential the form (Figure 48).

Figure 48. Equation-of-state of dipole substance from F-minimization.

Inspection of the p-profiles yields the following characteristic points.

• critical point

λ = 1.505, β = 1.545, p(λ) = 0.034

• triple point

λ = 1.35, β = 5.16, p = 0.03

14. Conclusions

The Landau theory of magnetic systems can be generalized to describe solid-fluid-gas systems as follows, $Z={\displaystyle \int D\phi dV\exp \left(-H\left(\overrightarrow{r},\phi ,\beta \right)\right)}$, where $H\left(\overrightarrow{r},\phi ,\beta \right)=E\left(\phi \right)+\beta u\left(\overrightarrow{r}\right)$ is the thermodynamic Hamiltonian, and $E\left(\phi \right)$ is the φ-induced Landau energy, $\beta u\left(\overrightarrow{r}\right)$ is the thermal intermolecular potential.

(first-order) fgs-systems: $\phi \left(r,\theta ,\lambda ,\sigma ;c_k\right)=\exp \left(-\frac{r}{\lambda }\right)f_{rdf}\left(r,\theta ,\sigma ;c_k\right)$ , radial distribution function $f_{rdf}\left(r,\theta ,\sigma ;c_k\right)$

(second-order) magnetic systems: $\phi =M$ magnetization

We have the following schematics.

Solid-fluid-gas-solid (fgs), first order transition

basic variables T, V, N

$\left(\frac{\partial F}{\partial v}\right)$ discontinuous, lat. heat $\Delta C={\left(\frac{\partial F}{\partial v}\right)}_{T=T_{+}}-{\left(\frac{\partial F}{\partial v}\right)}_{T=T_{-}}$

$Z={\displaystyle \int D\phi dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$,$D\phi ={\displaystyle \sum \frac{\partial \phi }{\partial c_k}}d{c}_k$ , $H\left(\overrightarrow{r},\phi ,\beta \right)=\beta u\left(\overrightarrow{r}\right)$

approximately $Z\approx {\displaystyle \int \phi dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$

$\left(c_k\right)=\left(l_c,{\alpha }_k,l_{a,k},{\alpha }_{0k}\right)$

$\phi ~\rho =\frac{1}{v}$ , specific volume $\nu$

$\begin{array}{l}\phi \left(r,\sigma ,\lambda ;c_k\right)=\exp \left(-\frac{r}{\lambda }\right){\phi }_{rdf}\left(r,\sigma ;c_k\right)\\ {\phi }_{rdf}\left(r,\sigma ;c_k\right)={\Theta }_H\left(r,\sigma \right)\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)\end{array}$

${\Theta }_H$ step-up function.

Minimization of free energy yields the solution

$F\left(T,V,N\right)=-\frac{1}{\beta }\log Z$ , $F\left(T,V,N\right)=\min$ $\to$ $\phi ={\phi }_0$

equivalently $Z\left(T,V,N\right)=\max$ $\to$ $\phi ={\phi }_0$ .

Magnetic, second order transition

basic variables T, M, N

$\left(\frac{{\partial }^{2}F}{\partial {\phi }^2}\right)$ discontinuous, $\Delta C=0$

$Z={\displaystyle \int D\phi dV\exp \left(-E\left(\phi \right)\right)}$ ,$D\phi =d\phi$ ,$\left(c_k\right)=\left\{\right\}$ , $H\left(\overrightarrow{r},\phi ,\beta \right)=E\left(\phi \right)$

$\phi =M$ magnetization, ${\phi }_s=\frac{M}{N}$ → $v$ fgs

φ-induced Landau energy =

$E\left(\phi ,\beta ,B\right)=\frac{\mu \left(\beta \right)}{2}{\phi }^2+\frac{\lambda \left(\beta \right)}{4!}{\phi }^4-\beta \phi B$

Minimization of free energy yields the solution

$\begin{array}{c}F\left(T,V,N\right)=-\frac{1}{\beta }\log Z\approx -\frac{1}{\beta }\left(E\left(\phi ,\beta \right)-\beta B\phi \right)\\ =V\left(\frac{const}{\beta }+\frac{1}{2\beta }{\mu }^2{\phi }^2+\frac{1}{4!\beta }{\lambda }^2{\phi }^4-B\phi \right)\end{array}$

T_c = critical temperature, ${\beta }_c=\frac{1}{k_B{T}_c}$ , ${\mu }^2\left({\beta }_c\right)=2{\mu }^{\prime }\left({\beta }_c\right)\left(\beta -{\beta }_c\right)$ , $F\left(T,V,N\right)=\min$ $\to$ $\phi ={\phi }_0$

Two solution methods: minimization of F, calculation of rdf function from HCOZ equation

• Minimization of F in parameters $c_k$

Solid-fluid-gas, first order transition: minimization of F

$Z\left(\beta ,\lambda ;c_k\right)={\displaystyle \int D\phi dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$, $D\phi ={\displaystyle \sum \frac{\partial \phi }{\partial c_k}}d{c}_k$

$F\left(\beta ,\lambda ;c_k\right)=-\frac{1}{\beta }\log Z\left(\beta ,\lambda ;c_k\right)$ , ${\lambda }^3=v$

pointwise minimization on $\left(\lambda ,\beta \right)$ lattice

$F\left(\beta ,\lambda ;c_k\right)=\min \left(c_k\right)$ $\leftrightarrow$ $Z\left(\beta ,\lambda ;c_k\right)=\max \left(c_k\right)$

→ $c_k\left(\lambda ,\beta \right)$

→ minimal free energy $F\left(\lambda ,\beta ,\sigma ,\epsilon ;c_k=c_k^{\left(0\right)}\right)$

• Calculation of rdf-function rom HCOZ equation

$F\left(\beta ,\lambda ;c_k\right)=-\frac{1}{\beta }\log Z\left(\beta ,\lambda ;c_k\right)$ , ${\lambda }^3=v$

$Z={\displaystyle \int D\phi dV\exp \left(-\beta u\left(\overrightarrow{r}\right)\right)}$,$D\phi ={\displaystyle \sum \frac{\partial \phi }{\partial c_k}}d{c}_k$

$\phi \left(r,\theta ,\lambda ,\sigma ;c_k\right)=\exp \left(-\frac{r}{\lambda }\right)f_{rdf}\left(r,\theta ,\sigma ;c_k\right)$

$\phi \left(r,\theta ,\lambda ,\sigma ;c_k\right)=\exp \left(-r/\lambda \right)\left(1+\exp \left(-r/l_c\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l_{a,k}}\right)}\right)\right)f_{\theta }\left(\theta ,{\alpha }_{0k}\right)$

parameters $\left(c_k\right)=\left(l_c,{\alpha }_k,l_{a,k},{\alpha }_{0k}\right)$

$\begin{array}{l}Z\left(\lambda ,\beta ,\sigma ,\epsilon ;c_k\right)\\ ={\displaystyle \int \left({\displaystyle \underset{r=\sigma }{\overset{\lambda }{\int }}{\displaystyle \underset{\theta =0}{\overset{\pi }{\int }}2\pi r^2\sin \theta d\theta dr{\displaystyle \sum _k\frac{\partial \phi \left(r,\theta ,\sigma ,\lambda ;c_k\right)}{\partial c_k}\exp \left(-\beta u\left(r,\theta ,\sigma ,\epsilon \right)\right)}}}\right)d{c}_{k}d\lambda }\end{array}$

Solve HypernettedChain-Ornstein-Zernicke (HCOZ) equation for rdf function $g\left(R\left(r_0,r_1,\theta ,\varphi \right);\sigma ,c_k\right)=f_{rdf}\left(r,\theta ,\sigma ;c_k\right)$

concisely $\log g\left(r_0\right)+\beta u\left(r_0\right)-\rho \left(1-g\left(r_0\right)\right)I\left(g,\beta u\right)=0$

$I\left(g,\beta u\right)={\displaystyle \int \left(\left(1-g\left(R\right)\right)-\log g\left(R\right)-\beta u\left(R\right)\right)}{r}_1{}^{2}d{r}_1\sin \theta d\theta d\varphi$

$R\left(r_0,r_1,\theta ,\varphi \right)=\sqrt{r_1^2{\sin }^2\theta +R_1^2=r_0^2+r_1^2+2r_0{r}_1\cos \theta \cos \varphi }$

→solution $f_{rdf}\left(r,\theta ,\sigma ;c_k=c_k^{\left(0\right)}\right)$

→free energy $F\left(\lambda ,\beta ,\sigma ,\epsilon ;c_k=c_k^{\left(0\right)}\right)$

Results Lennard-Jones substance (fluid argon) with HCOZ eq.

Potential $u_{cLJ}\left(r,\sigma ,\epsilon \right)=4\epsilon \left(\frac{1}{\frac{1}{2}+{\left(\frac{r}{\sigma }\right)}^{12}}-\frac{1}{1+{\left(\frac{r}{\sigma }\right)}^6}\right)$

Ansatz rdf $g\left(r,\sigma ,c_k\right)={\Theta }_H\left(r,\sigma ,\Delta r\right)\left(1+\Delta g\left(r,\sigma ,c_k\right)\right)$

$\Delta g\left(r,\sigma ,c_k\right)=\exp \left(-\frac{r}{lc}\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l{a}_k}\right)}\right)$

solution $g_L\left(r,\rho ,\beta ;c_k\right)$ with coefficients, $c_k\left(\beta ,\rho \right)=\left(lc\left(\beta ,\rho \right),{\alpha }_k\left(\beta ,\rho \right),l{a}_k\left(\beta ,\rho \right)\right)$ , $\rho =1/v={\lambda }^{-3}$

Pressure $p_L\left(\lambda ,\beta ;c_k\left(\beta ,\lambda \right)\right)=-\frac{\partial F}{\partial v}=\frac{1}{\beta }\frac{\partial \log Z}{\partial v}=\frac{1}{\beta Z}\frac{1}{3{\lambda }^2}\frac{\partial Z}{\partial \lambda }$

The ideal gas invariant $E_L\left(\lambda ,\beta \right)=p_L\left(\lambda ,\beta \right)\beta {\lambda }^3$

saturation curve fluid-gas in the range $\lambda \approx 1.6$ $1\le \beta \le 1.6$ , followed by the solid-gas evaporation curve in the range $1.7\le \beta \le 4$

• critical point

λ_cr = 1.62, p = 0.016, β_cr = 1.035, visible in the β = 1.035 profile, with the approximate turning point at λ = 1.84, measured values (argon): β_cr = 0.98, λ_cr = 1.47, p_cr = 0.097p₀.

• triple point

λ_tr = 1.58, p = 0.004, β_tr = 1.63, measured (argon): β_tr = 1.67, λ_tr ≈ 1.4, p_tr = 0.0014p₀

• saturation curve

• critical exponent δ

The fit to the curve $\Delta p\left(\lambda ,\beta ={\beta }_{cr}\right)$ , $\left|\Delta p\right|~{\left|\Delta \rho \right|}^{\delta },T\to T_{C-}$ , at the critical point yields the value δ = 2.93, for the vdWaals fluid we have δ = 3, for the LJ-fluid δ = 4.8.

Results dipole substance (ethanol) with HCOZ eq.

Potential $u_{cDD}\left(r,\sigma ,\epsilon ,\theta \right)=4\epsilon \left({\Theta }_L\left(r,\sigma ,0.1\right)-\left(1-{\Theta }_L\left(r,\sigma ,0.1\right)\right)\cos \theta /\left(0.1+{\left(r/\sigma \right)}^3\right)\right)$ ,

Ansatz rdf $g\left(r,\theta ,\sigma ,c_k\right)={\Theta }_H\left(r,\sigma ,\Delta r\right)\left(1+\Delta g\left(r,\sigma ,c_k\right)\right)$

$\Delta g\left(r,\theta ,\sigma ,c_k\right)=\exp \left(-\frac{r}{lc}\right)\left({\displaystyle \sum _k{\alpha }_k\cos \left(2\pi \frac{r-\sigma }{l{a}_k}\right)}\right)\left(\cos \theta +{\alpha }_0\right)$

solution $g_D\left(r,\theta ,\rho ,\beta ;c_k\right)$ with coefficients, $c_k\left(\beta ,\rho \right)=\left(lc\left(\beta ,\rho \right),{\alpha }_k\left(\beta ,\rho \right),l{a}_k\left(\beta ,\rho \right),{\alpha }_0\left(\beta ,\rho \right)\right)$ , $\rho =1/v={\lambda }^{-3}$

Pressure $p_D\left(\lambda ,\beta ;c_k\left(\beta ,\lambda \right)\right)=-\frac{\partial F}{\partial v}=\frac{1}{\beta }\frac{\partial \log Z}{\partial v}=\frac{1}{\beta Z}\frac{1}{3{\lambda }^2}\frac{\partial Z}{\partial \lambda }$

The ideal gas invariant $E_D\left(\lambda ,\beta \right)=p_D\left(\lambda ,\beta \right)\beta {\lambda }^3$

saturation curve fluid-gas in the range $1.35\le \lambda \le 1.46$ , $1.46\le \beta \le 5$ , branching off is the solid-gas evaporation curve in the range $5\le \beta \le 6$

• critical point

λ_cr = 1.46, β_cr = 1.545 ± 0.1, (λ) = 0.225, measured (ethanol) p_cr = 0.051, λ_cr = 1.49, β_cr = 1.46

• triple point

λ_tr = 1.35, β_tr = 5.16 ± 0.2, p = 0.004,measured p_tr = 0.051, λ_tr = 1.3, β_tr = 5.

• saturation curve (fluid-gas)

• critical exponent β₀

The fit to the curve $\Delta \rho \left(\lambda ={\lambda }_{cr},\beta \right)$ , $\left|\Delta \rho \right|~{\left|\Delta T\right|}^{{\beta }_0},T\to T_{C-}$ at the critical point yields the value β₀ = 0.46, for the vdWaals fluid we have β₀ = 0.5, for the LJ-fluid β₀ = 0.33.