The Ornstein-Zernike Equation: three distinct approaches

Abstract

Liquid state theories are grounded in statistical mechanics, making key concepts challenging for newcomers. Important results and discussions are often scattered across various papers and textbooks, many of which address broader topics and include only a single chapter on liquid state theory. As a result, gathering the necessary foundational knowledge can be time-consuming and requires some insight in finding the right resources. This paper aims to provide a clear and accessible introduction to liquid state theory by discussing the derivation of the Ornstein-Zernike equation through three different approaches: 1) intuitive arguments based on particle interactions within the system; 2) a statistical analysis of particle fluctuations in different subvolumes, accounting for interactions between those subvolumes, as originally presented by Ornstein and Zernike; and 3) the use of correlation functions derived from statistical mechanics and their variations due to changes in an external potential. While the third approach is more complex, it offers essential insights and mathematical intuition for deeper exploration in the field.

Keywords:
Liquid state theory; Statistical mechanics; Ornstein-Zernike equation; Particle interactions; Particle fluctuations; Correlation functions

1. Introduction

The approach given by J. D. van der Waals in 1873 in his doctoral thesis On the continuity of the gaseous and liquid state , translated into english by J. S. Rowlinson, [¹] marked the first step towards developing a theory of the liquid state that accounted for not only repulsive forces, as was previously believed, but also attractive ones. The correct interpretation of an additional term proposed by Hirn in 1863, [²] as being a consequence of attractive forces and the innumerable applications of the proposed equation to various systems, awarded J. D. van der Waals the Nobel Prize in Physics in 1910 [³]. The title of the thesis is explained in his Nobel Lecture: Thus, I conceived the idea that there is no essential difference between the gaseous and the liquid state of matter…And so the idea of continuity occurred to me.

Although this approach, in its original form, does not start with a Hamiltonian, from which the results can be derived, it was later proven to be the consequence of a potential energy characterized by a hard core and a dispersion term, see references [⁴, ⁵, ⁶,]. This opened the way for investigating a rigorous treatment of the liquid state [⁷, ⁸,].

Based on the work of J. W. Gibbs, in his book on Statistical Mechanics, [⁹] several integrals involving the potential energy were developed. While studying accidental deviations of density of single substances at the critical point, by taking into account the influences of deviations in one subvolume of the system into another, Ornstein (then a student of H. A. Lorentz) and Zernike proposed in 1914, an integral equation involving the potential energy and two correlation functions. The Ornstein-Zernike equation, as to be discussed in this work, is a mathematical framework behind the original idea of van der Waals on the continuity of matter. This was only possible after the introduction of attractive and repulsive forces, described by a potential energy function.

It is important to note that Ornstein-Zernike (OZ) equation is not the only method in statistical physics to study fluids. A series of integro-differential equations firstly proposed by Yvon and later rederived by others, knows as Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy, [¹⁵] can also be employed to study dynamical properties of fluids, such as self-diffusion in dense liquids [¹⁶]. The method for solving this problem consists in truncating the BBGKY series, enabling its use to derive quantum [¹⁷] and classical kinetic theories, such as the Boltzmann equation [¹⁸]. Although both OZ and BBGKY theories can be applied to classical fluids, their applications differ: the BBGKY hierarchy is used for dynamical systems, while OZ is specific to fluids in equilibrium. In this paper, focus will be given to OZ equation.

The OZ equation expresses the relationship between two types of correlation function: the direct correlation function, which is short-ranged, and the total correlation function, which accounts for the long-ranged interactions. These two functions are extremely important, as they are related to experimental properties such as: the partial molar volume [¹⁰, ¹¹, ¹²], excess Helmholtz free energy and chemical potential [¹², ¹³,], isothermal compressibility, pressure [¹³] and the structure factor [¹⁴].

These correlation functions and the Ornstein-Zernike equation will be discussed in this paper. Initially, the OZ equation will be introduced intuitively by discussing the interactions between particles in a system. Then, a statistical reasoning follows, as in the original work, [¹⁹] for obtaining this equation. The third method presented, using functional derivatives, is more elaborate but necessary for the reader to follow the literature on the theory of liquids and solutions. In particular, the closure relations, which gives a relationship between the two types of correlation functions, can be elaborated using functional derivatives [²⁰]. Also, the classical density functional theory of fluids will use this formalism to derive important results. To make the material more understandable, different notations are going to be used for each approach

2. The N-Body Density Function

Rigorous derivations of the Ornstein-Zernike equation will be given. However, before proceeding it is important to understand the concept of n-body density functions.

In the canonical ensemble, i.e. with constant volume, temperature and number of particles, the classical probability of a system containing N identical particles to be found with momenta and coordinates given by p₁, p₂, …, p_N and r₁, r₂, …, r_N in a differential volume of the phase-space, d^pNd^rN, is given by [¹³]

(1)$$d{p}^N\left(p^N,r^N\right)=\frac{e^{-\beta H\left(p^N,r^N\right)}d{p}^{N}d{r}^N}{\int e^{-\beta H\left(p^N,r^N\right)}d{p}^{N}d{r}^N}$$

with $\beta =\frac{1}{k_{B}T}$, k_B the Boltzmann constant and T the temperature, H(p^N,r^N) being the system’s Hamiltonian. Integrating over all coordinates (r^N) and momenta (p^N), i.e. summing the probability contribution of all states, the result is 1. While studying systems in equilibrium, the primary interest is in the spacial configuration of the particles regardless of its momenta. Therefore, by integrating over the all possible values of p^N and noting that $H\left(p^N,r^N\right)={\displaystyle {\sum }_i^n\frac{p_i^2}{2m}}+E_p(r^N),$, one is left with

(2)$$d{p}^N\left(r^N\right)=\frac{e^{-\beta E_p\left(r^N\right)}d{r}^N}{\int e^{-\beta E_p\left(r^N\right)}d{r}^N}$$

with E_p(r) the potential energy function that describes the interaction between particles.

It is also possible to write the probability of some particles subset to be found at a given configuration regardless of the remaining ones. Thus, the probability of finding the particle with coordinate r₁ in dr₁, r₂ in dr₂, …, r_n in dr_n and the particles with coordinates r_{n + 1}, r_{n + 2}, …, r_N in any possible position will be,

(3)$$d{p}^{(n)}(r^n)=\frac{d{r}_{1}d{r}_2\text{…}d{r}_n}{Q_N}\times \int e^{-\beta E_p(r_1,r_2,\text{…},r_N)}d{r}_{n+1}d{r}_{n+2}\text{…}d{r}_N$$

in which Q_N = ∫e^−E_p(rN)dr_N is called configuration integral. Since the particles are identical, one can have any of those N particles in dr₁, N − 1 in dr₂, …, N−n + 1 in dr_n. Therefore, this multiplicity of possibilities must be taken into account by multiplying equation (3) by

(4)$$N(N-1)\text{…}(N-n+1)=\frac{N!}{(N-n)!}$$

Then, the probability of finding a subset of particles at a given configuration is

(5)$$d{p}_g^{(n)}(r^n)=\frac{N!}{(N-n)!}\frac{d{r}_{1}d{r}_2\text{…}d{r}_n}{Q_N}\times \int e^{-\beta E_p(r_1,r_2,\text{…},r_N)}d{r}_{n+1}d{r}_{n+2}\text{…}d{r}_N$$

This quantity is known as the generic differential probability.

The generic differential probability, if divided by dr₁dr₂…dr_n, is the probability density, ${\rho }^{(n)}(r_1,r_2,\text{…},r_n)=\frac{d{p}_g^{(n)}(r^n)}{d{r}_{1}d{r}_2\text{…}d{r}_n}$ and will have the dimension of number of particles per volume, justifying the variable name. Therefore,

(6)$${\rho }^{(n)}(r_1,r_2,\text{…},r_n)=\frac{N!}{(N-n)!}\frac{1}{Q_N}\int e^{-\beta E_p(r_1,r_2,\text{…},r_N)}d{r}_{n+1}d{r}_{n+2}\text{…}d{r}_N$$

Some particular cases of probability density are

1. One-body probability density:

   (7)$${\rho }^{(1)}(r_1)=\frac{N}{Q_N}\int e^{-\beta E_p(r_1,r_2,\text{…},r_N)}d{r}_{2}d{r}_3\text{…}d{r}_N$$

2. Two-body probability density:

   (8)$${\rho }^{(2)}(r_1,r_2)=\frac{N(N-1)}{Q_N}\times \int e^{-\beta E_p(r_1,r_2,\text{…},r_N)}d{r}_{3}d{r}_4\text{…}d{r}_N$$

3. Three-body probability density:

   (9)$${\rho }^{(3)}(r_1,r_2,r_3)=\frac{N(N-1)(N-2)}{Q_N}\times \int e^{-\beta E_p(r_1,r_2,\text{…},r_N)}d{r}_{4}d{r}_5\text{…}d{r}_N$$

These definitions will be useful in the discussion on liquids. In particular, ρ⁽²⁾(r₁,r₂) will be related to the radial distribution function as

(10)$$g(r_1,r_2)=\frac{{\rho }^{(2)}(r_1,r_2)}{{\rho }^{(1)}(r_1){\rho }^{(1)}(r_2)}$$

Also, if the fluid is isotropic, i.e. the particle configurations are spherically symmetric, one can translate the origin to one of the particles position, r_j. Since the origin is located at particle j position, the system potential energy will only depends on the remaining particle positions. Thus, the one-body probability density is rewritten as

(11)$$\begin{array}{lll}{\rho }^{(1)}(r_1)\hfill & =\hfill & \frac{N{\displaystyle \int e^{-\beta E_p(r^{\prime }2,\cdots ,r^{\prime }N)}}d{r}^{\prime }{}_2\dots d{r}^{\prime }N}{N{\displaystyle \int e^{-\beta E_p(r^{\prime }2,\cdots ,r^{\prime }N)}}d{r}^{\prime }{}_2\dots d{r}^{\prime }N}\hfill \\ \hfill & =\hfill & \frac{N{\displaystyle \int e^{-\beta E_p(r^{\prime }2,\cdots ,r^{\prime }N)}}d{r}^{\prime }{}_2\dots d{r}^{\prime }N}{V{\displaystyle \int e^{-\beta E_p(r^{\prime }2,\cdots ,r^{\prime }N)}}d{r}^{\prime }{}_2\dots d{r}^{\prime }N}\hfill \end{array}$$

with $\rho =\frac{N}{V}$ being the system density.

3. Total and Direct Correlation Functions

While studying systems containing N interacting particles, it is intuitive to infer that correlations between these particles will exist. There will be two correlation functions to be discussed in the following sections:

1. Direct correlation function, c(r₁,r₂).

2. Total correlation function, h(r₁,r₂).

These were defined considering any two particles in the system, but there are analogous definitions for several variables. The direct correlation function acts between near neighbors and it is short-ranged. Its range is approximately the same as that of the potential energy between particles. In contrast, the total correlation function is long-ranged and will take into account influence of both near and distant neighbors.

A third correlation functions, related to h(r₁, r₂), is defined as,

(12)$$g(r_1,r_2)=h(r_1,r_2)+1$$

and is known as radial distribution function. These correlation functions are of central importance in liquid state theory, since from these it is possible to calculate thermodynamic quantities [¹², ¹³] and the structure factor: properties that can be measured experimentally [²¹, ²², ²³,]. For clarity reasons, these properties will be discussed below considering a spherically symmetric systems, e.g. a homogeneous monoatomic fluid.

From the total correlation function, it is possible to calculate the structure factor for the fluid as

(13)$$S(q)=1+4\pi \rho \int r^{2}h(r)\frac{\sin (qr)}{qr}dr$$

with q the reciprocal of r. A detailed derivation of equation (13) is given in reference [¹⁴]. Also, thermodynamic properties, as for example, the virial pressure can be calculated as [¹³]

(14)$$\frac{\beta p^v}{\rho }=1-\frac{\beta \rho }{6}\int \left[rg(r)\frac{d{E}_p(r)}{dr}\right]dr$$

The isothermal compressibility is calculated as follows [¹³]

(15)$$\frac{\rho {\chi }_T}{\beta }=1+\rho \int h(r)dr$$

Another information that can be obtained from the radial distribution function is the number of particles contained in a spherical volume around the reference particle:

(16)$$N={\int }_{R_0}^{R}4\pi r^2\rho g(r)dr$$

in which R₀ and R are the smaller and larger radii of a spherical shell, respectively. From equation (16) it is possible to see that the radial distribution function gives the deviation of the local density with respect to the bulk fluid density, ρ.

4. Consecutive Influence of Direct Correlation

Consider a system containing three particles, as represented in Figure 1. The total correlation between particles 1 and 2 can be written as the direct correlation between 1 and 2 and an indirect correlation in which particle 1 interacts with 3 and 3 with 2.

Figure 1
A system of three (left) and four (right) interacting particles with the direct and indirect correlations indicated by the full and dashed arrows, respectively.

Now, one needs to sum over all the possible positions for particle 3. However, since there is a probability of finding 3 at position r₃ regardless of the remaining particle positions, the indirect correlation also needs to take this into account. Thus, one accounts the indirect correlation through the following integral ∫c(r₁, r₃)ρ⁽¹⁾(r₃)c(r₃, r₂)dr₃. For an isotropic fluid ρ⁽¹⁾(r₃) = ρ and the total correlation is written as

(17)$$h12=c12+\rho \int c13c32d3$$

with fij = f(r_i, r_j) and dj = dr_j. A similar reasoning can be carried out for a system containing four interacting particles, as represented in Figure 1. It is important to note at this point that although the interactions 1 → 4 → 2 and 1 → 4 → 3 → 2 were omitted in the figure, those are still considered. As for example, in the indirect correlation considering three particles, r₃ in ∫c(r₁, r₃)ρ⁽¹⁾(r₃)c(r₃, r₂)dr₃ will take into account the positions of both particles labeled as 3 and 4 in the figure.

To make the influence of the particles more clear, the notation 1 → 3 → 2 is modified to $1\to \underset{\_}{3}\to 2$, indicating that the integration is to be performed under the third particle. Using this notation one can continue the influence of a fourth particles as,

(18)$$\begin{array}{lll}h12\hfill & =\hfill & c12+\rho \stackrel{1\to \underset{\_}{3}\to 2}{\overbrace{{\displaystyle \int }c13c32d3}}+{\rho }^2\stackrel{1\to \underset{\_}{3}\to \underset{\_}{4}\to 2}{\overbrace{{\displaystyle \int }{\displaystyle \int }c13c34c42d3d4}}+\dots \hfill \\ \hfill & =\hfill & c12+\rho \int c13\underset{h32}{\underbrace{\left(c32+\rho {\displaystyle \int }c34c42d4\right)}}d3\hfill \\ \hfill & =\hfill & c12+\rho \int c13h32d3\hfill \end{array}$$

in which equation (17) was used to introduce the correlation h32. More particles and interactions can be added,

(19)$$\begin{array}{lll}h12\hfill & =\hfill & c12+\rho \stackrel{1\to \underset{¯}{3}\to 2}{\overbrace{{\displaystyle \int c13c32d3}}}+{\rho }^2\stackrel{1\to \underset{¯}{3}\to \underset{¯}{4}\to 2}{\overbrace{{\displaystyle \int {\displaystyle \int c13c34c42d3d4}}}}+{\rho }^3\stackrel{1\to \underset{¯}{3}\to \underset{¯}{4}\to \underset{¯}{5}\to 2}{\overbrace{{\displaystyle \int {\displaystyle \int {\displaystyle \int c13c34c45c52d3d4d5+\dots }}}}}\hfill \\ \hfill & =\hfill & c12+\rho \int c13\times \underset{h32}{\underbrace{\left(c32+\rho \int c34c42d4{\rho }^2\int \int c34c45c52d4d5+\dots \right)d3}}\hfill \\ \hfill & =\hfill & c12+\rho \int c13h32d3\hfill \end{array}$$

resulting in the same equation. This result is known as the Ornstein-Zernike equation, generally written as

(20)$$h(r_{\text{1}},r_2)=c(r_{\text{1}},r_2)+\rho {\displaystyle \int h(r_{\text{1}},r)c(r_{\text{2}},r)}dr$$

5. The Ornstein-Zernike Equation from Statistical Fluctuations

Consider a volume, V, subdivided into several volumes dv_k. It is assumed that the average number of particles in the system is N, that is ⟨N⟩. Thus, $N={\displaystyle {\sum }_k{N}_k}$ and $〈N〉=〈{\displaystyle {\sum }_k{N}_k}〉={\displaystyle {\sum }_k〈N_k〉}$. The variance in the number of particle for this system is given by

(21)$$\begin{array}{lll}var(N)\hfill & =\hfill & 〈{\left(N-〈N〉\right)}^2〉\hfill \\ \hfill & =\hfill & 〈{\left({\displaystyle {\sum }_k{N}_k}-{\displaystyle {\sum }_k〈N_k〉}\right)}^2〉\hfill \\ \hfill & =\hfill & 〈{\displaystyle {\sum }_k\left(N_k-〈N_k〉\right){\displaystyle {\sum }_j\left(N_j-〈N_j〉\right)}}〉\hfill \\ \hfill & =\hfill & {\displaystyle {\sum }_k{\displaystyle {\sum }_j〈\mathrm{\Delta }{N}_k\mathrm{\Delta }{N}_j〉}}\hfill \\ \hfill & =\hfill & {\displaystyle {\sum }_k〈\mathrm{\Delta }{N}_k^2〉}+{\displaystyle {\sum }_k{\displaystyle {\sum }_{j\ne k}〈\mathrm{\Delta }{N}_k{N}_j〉}}\hfill \end{array}$$

with ΔN_k = N_k − N_k being the fluctuations in the number of particles for the subvolume dv_k. Therefore, the variance in the system total number of particles is equals to the sum of variance in the number of particles in each subvolume and the covariance between the number of particles fluctuations for every pair of subvolumes.

The correlation between fluctuations on different volume elements can be expressed in two distinct ways. The first one relates the fluctuations in volume k with those occurring in i, taking into account the direct and indirect correlations through a total correlation function, h(r_ik):

(22)$$\text{Δ}{N}_k=h(r_{ik})\rho d{v}_k\text{Δ}{N}_i$$

The term ρdv_k will scale the result taking into account the system density and the volume of k. The second way relates the fluctuations in volume k with those occurring in all other subvolumes directly:

(23)$$\begin{array}{lll}\mathrm{\Delta }{N}_k\hfill & =\hfill & {\displaystyle \sum _{j\ne k}c\left(r_{j}k\right)\rho d{v}_k\mathrm{\Delta }{N}_j}\hfill \\ \hfill & =\hfill & c(r_{ik})\rho d{v}_k\mathrm{\Delta }{N}_i+{\displaystyle \sum _{j\ne k,i}c\left(r_{j}k\right)\rho d{v}_k\mathrm{\Delta }{N}_j}\hfill \end{array}$$

in which the terms for subvolume i were made explicit. Since the fluctuations observed in k must be equal regardless of the which expression is used, equations (22) and (23) must be equal. Thus,

(24)$$h(r_{ik})\rho d{v}_k\text{Δ}{N}_i=c(r_{ik})\rho d{v}_k\text{Δ}{N}_i+{\displaystyle \sum _{j\ne k,i}c(r_{jk})}\rho d{v}_k\text{Δ}{N}_i$$

It is possible to use equation (22) on the righthand side of (24) to express the dependence of ΔN_j on the fluctuation in i:

(25)$$h(r_{ik})\rho d{v}_k\mathrm{\Delta }{N}_i=c\left(r_{ik}\right)\rho d{v}_k\mathrm{\Delta }{N}_i+{\displaystyle \sum _{j\ne k,i}c(r_{jk})}\rho d{v}_{k}h(r_{ij})\rho d{v}_j\mathrm{\Delta }{N}_i$$

dividing both sides by ρdvkΔN_i, one arrives at

(26)$$h(r_{ik})=c\left(r_{ik}\right)+\rho {\displaystyle \sum _{j\ne k,i}c(r_{jk})}h(r_{ij})d{v}_j$$

which is the discrete form of the Ornstein-Zenike equation. Taking the limit of dv_k → 0, this result can be expressed as an integral and equation (20) is recovered.

Statistical considerations of the fluctuations, as presented in the present derivation, were used in the original proof [¹⁹, ²⁴].

6. Fluids Under an External Potential

Although fluids are homogeneous in the bulk, interactions with any interface, such as confining walls or a different phase, will perturb the liquid and inhomogeneities will appear [¹⁵]. Therefore, analyzing fluids under the influence of an external potential is particularly insightful, as it allows one to recover the properties of homogeneous fluids by taking the limit of a vanishing external potential.

6.1. First variation of ln⁡(Z_N)

Considering the total external potential as the sum of its individual contributions acting on each particle in the system, one has

(27)$$w(r^N)=\sum _{i=1}^{N}w(r_i)$$

with r_i being the coordinates of the i-th atom and ^rN denoting the set of coordinates for the N particles in the system. Thus, the partition function can be written as

(28)$$Z_N=\int e^{-\beta {\sum }_i^N\frac{{\text{p}}_i^2}{2m}}d{p}^N\int e^{-\beta \left(E_p(r^N)+{\sum }_i\mathrm{}w(r_i)\right)}d{r}^N$$

As before, integrating over the momenta

(29)$$Z_N=\frac{1}{N!}{V}^N{\text{Λ}}^{3N}\int e^{-\beta \left(E_p(r^N)+{\sum }_i\mathrm{}w(r_i)\right)}d{r}^N$$

with $\text{Λ}=\sqrt{\frac{2\pi m{k}_{B}T}{h^2}}$, in which m is the particle mass and h the Planck’s constant. The Helmholtz free energy is calculated from the partition function as

(30)$$F=-k_{B}T\ln \left(Z_N\right)$$

Therefore, now it is possible to evaluate changes in the free energy due to variations on the external potential. Considering that the volume is constant, one has the variation for the partition function’s logarithm given by

(31)$$\delta \ln \left(Z_N\right)=-\frac{\beta }{Q_N}\int \sum _j\left(\delta w(r_j)\right)e^{-\beta \left(E_p(r^N)+{\sum }_i\mathrm{}w(r_i)\right)}d{r}^N$$

Noting that the external potential is the same for all particles, thus ${\sum }_{j=1}^N\left(\delta w(r_j)\right)=N\delta w(r_1)$, and comparing with equation (7), this result can be rewritten as

(32)$$\begin{array}{lll}\delta \ln \left(Z_N\right)\hfill & =\hfill & -\beta \frac{N}{Q_N}{\displaystyle {\int }_{r_1}\delta }w(r_1)\times \left[{\displaystyle \int e^{-\beta \left(E_p(r^N)+{\displaystyle {\sum }_{i}w}(r_i)\right)}}d{r}_2,\dots ,d{r}_N\right]d{r}_1\hfill \\ \hfill & =\hfill & -\beta {\displaystyle \int {\rho }^{(1)}}(r_1;w)\delta w(r_1)d{r}_1\hfill \end{array}$$

introducing the ρ⁽¹⁾(r₁;w) as the one-body probability density function when an external potential is acting on the particle. Thus, the Helmholtz free energy variation at a given temperature is

(33)$$\delta F=\int {\rho }^{(1)}(r_1;w)\delta w(r_1)d{r}_1$$

The formal definition of the functional derivative can be expressed as [¹⁵]

(34)$$\delta F=\int A[u;x]\delta u(x)dx$$

with x the coordinate, u(x) a function with dependence on x and $A[u;x]=\frac{\delta F}{\delta u(x)}$. Thus, by comparing equations (34) and (33) one has, taking the limit of a vanishing external potential,

(35)$${\rho }^{(1)}(r_1)={\left[\frac{\delta F}{\delta w(r_1)}\right]}_{w\to 0}$$

Therefore, the one-body probability density also measures the response of the free energy to a small variation in the system’s Hamiltonian.

6.2. Second variation of ln⁡(Z_N)

The analysis of variations on ln (Z_N) due to changes in the external potential will be extended to the second order, δ² ln Z. Noting from equation (32) that $\frac{\delta \ln \left(Z_N\right)}{\delta w_1}=-\beta {\rho }^{(1)}(r_1;w)$, the desired results can be acquired by carrying out δρ⁽¹⁾(r₁;w). Thus, taking the functional derivative of ρ⁽¹⁾(r₁;w):

(36)$$\begin{array}{lll}\delta {\rho }^{(1)}(r_1;w)\hfill & =\hfill & N\delta \text{}\left(\text{}\frac{1}{Q_N}\int e^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i)}d{r}_{2}d{r}_3\dots d{r}_N\text{}\text{}\right)\hfill \\ \hfill & =\hfill & \frac{N{Q}_N\delta \left({\displaystyle \int e^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}}d{r}_{2}d{r}_3\dots d{r}_N\text{}\right)-N\left(\delta Q_N\right){\displaystyle \int e^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}}d{r}_{2}d{r}_3\dots d{r}_N}{Q_N^2}\hfill \\ \hfill & =\hfill & \frac{N\delta \left({\displaystyle \int e^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}}d{r}_{2}d{r}_3\dots d{r}_N\right)}{Q_N}-\frac{N(\delta Q_N)\text{}{\displaystyle \int \text{}}{e}^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_i\text{}}w(r_i))}\text{}d{r}_{2}d{r}_3\dots d{r}_N}{Q_N^2}\hfill \\ \hfill & =\hfill & P-S\hfill \end{array}$$

Developing the first term in the last equality, P,

(37a)$$P=\frac{N}{Q_N}\delta \left({\displaystyle \int e^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}}d{r}_{2}d{r}_3\dots d{r}_N\right)$$

(37b)$$=-\beta \frac{N}{Q_N}{\displaystyle \int \left({\displaystyle \sum _{i=1}\delta }{w}_i\right)}{e}^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}d{r}_{2}d{r}_3\dots d{r}_N$$

(37c)$$=-\beta \frac{N}{Q_N}{\displaystyle \int \delta }{w}_1{e}^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}d{r}_{2}d{r}_3\dots d{r}_N-\beta \frac{N}{Q_N}{\displaystyle \int \left({\displaystyle \sum _{i=2}\delta }{w}_i\right)}{e}^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}d{r}_{2}d{r}_3\dots d{r}_N$$

(37d)$$=-\beta \frac{N}{Q_N}{\displaystyle \int \delta }{w}_1{e}^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}d{r}_{2}d{r}_3\dots d{r}_N-\beta \frac{N(N-1)}{Q_N}{\displaystyle \int \delta }{w}_2\left({\displaystyle \int e^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}}d{r}_3\dots d{r}_N\right)d{r}_2$$

(37e)$$=-\beta \delta w_1{\rho }^{(1)}(r_1;w)-\beta {\displaystyle \int \delta }{w}_2{\rho }^{(2)}(r_1,r_2)d{r}_2$$

(37f)$$=-\beta {\displaystyle \int \delta }{w}_2{\rho }^{(1)}(r_1;w)\delta (r_1-r_2)d{r}_2-\beta {\displaystyle \int \delta }{w}_2{\rho }^{(2)}(r_1,r_2)d{r}_2$$

in which δw(r_i) has been shortened to δw_i for brevity in the notation, the first term in the summation of equation (37b) was made explicit in equation (37c) and the delta function property f(r₁) = ∫δ(r₁ − r₂)f(r₂)dr₂ was used on the first term of equation (37f). To develop S, the variation of Q_N with respect to w will be acquired first:

(38)$$\delta Q_N=-\beta N{\displaystyle \int \delta }{w}_2{e}^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}d{r}_{1}d{r}_2\dots d{r}_N$$

Using this result on the second term:

(39)$$\begin{array}{lll}S\hfill & =\hfill & \frac{-N(-\beta N{\displaystyle \int \delta }{w}_2{e}^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}d{r}_{1}d{r}_2\dots d{r}_N)\times {\displaystyle \int e^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}}d{r}_{2}d{r}_3\dots d{r}_N}{Q_N^2}\hfill \\ \hfill & =\hfill & \beta N\frac{{\displaystyle \int \delta }{w}_2{e}^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}d{r}_{1}d{r}_2\dots d{r}_N}{Q_N}\times N\frac{{\displaystyle \int e^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}}d{r}_{2}d{r}_3\dots d{r}_N}{Q_N}\hfill \\ \hfill & =\hfill & \beta N\frac{{\displaystyle \int \delta }{w}_2{e}^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}d{r}_{1}d{r}_2\dots d{r}_N}{Q_N}\times {\rho }^{(1)}(r_1;w)\hfill \\ \hfill & =\hfill & \beta \frac{{\displaystyle {\int }_{r_2}\delta }{w}_2\left(N{\displaystyle \int e^{-\beta (E_p(r_1,r_2,\dots ,r_N)+{\displaystyle {\sum }_{i}w}(r_i))}}d{r}_{1}d{r}_3\dots d{r}_N\right)d{r}_2}{Q_N}\times {\rho }^{(1)}(r_1;w)\hfill \\ \hfill & =\hfill & \beta {\displaystyle \int \delta }{w}_2{\rho }^{(1)}(r_2;w){\rho }^{(1)}(r_1;w)d{r}_2\hfill \end{array}$$

Combining these results in equation (36)

(40)$$\begin{array}{lll}{\delta }^2\text{ln}Z\hfill & =\hfill & -\beta \delta {\rho }^{(1)}(r_1;w)\hfill \\ \hfill & =\hfill & -\beta \left[-\beta \int \delta w(r_2){\rho }^{(1)}(r_1;w)\delta (r_1-r_2)d{r}_2-\beta \int \delta w_2(r_2){\rho }^{(2)}(r_1,r_2)d{r}_2+\beta \int \delta w_2(r_2){\rho }^{(1)}(r_2;w){\rho }^{(1)}(r_1;w)d{r}_2\right]\hfill \\ \hfill & =\hfill & {\beta }^2{\displaystyle \int \delta w(r_2)}\left[{\rho }^{(1)}(r_1;w)\delta (r_1-r_2)+{\rho }^{(2)}(r_1,r_2)-{\rho }^{(1)}(r_2;w){\rho }^{(1)}(r_1;w)\right]d{r}_2\hfill \end{array}$$

Therefore, one has

(41)$$\frac{{\delta }^2\ln Z}{\delta w_2\delta w_1}=-\beta \frac{\delta {\rho }^{(1)}(r_1;w)}{\delta w_2}={\beta }^2\left[{\rho }^{(1)}(r_1;w)\delta (r_1-r_2)+{\rho }^{(2)}(r_1,r_2)-{\rho }^{(1)}(r_2;w){\rho }^{(1)}(r_1;w)\right]$$

From equation (41) it is seem that the two-body correlation function is related to the second variation of the free energy.

7. The Ornstein-Zernike Equation From the Density Functional Derivatives

Here, the Ornstein-Zernike equation will be derived from equation (41) and the inverse relation for the functional derivative, [²⁵] given by

(42)$$\int \frac{\delta {\rho }^{(1)}(r_1;w)}{\delta w(r)}\frac{\delta w(r)}{\delta {\rho }^{(1)}(r_2;w)}dr=\delta (r_1-r_2)$$

The inverse $\frac{\delta w(r)}{\delta {\rho }^{(1)}(r_2;w)}$ was derived by Percus [²⁰] and it is given by

(43)$$\frac{\delta w(r_2)}{\delta {\rho }^{(1)}(r_1;w)}=-\frac{1}{\beta }\left(\frac{\delta (r_1-r_2)}{{\rho }^{(1)}(r_1;w)}-c(r_1,r_2)\right)$$

The detailed derivation of this result can be found in reference [¹⁵]. However, since it requires carrying out an analogous procedure considering the grand canonical ensemble, it will not be discussed here.

Thus, using equations (41), (42) and (43):

(44)$$\begin{array}{lll}I\hfill & =\hfill & \int \left[{\rho }^{(1)}(r_1;w)\delta (r_1-r)+{\rho }^{(2)}(r_1,r)-\left.{\rho }^{(1)}(r;w){\rho }^{(1)}(r_1;w)\right]\left[\frac{\delta (r_2-r)}{{\rho }^{(1)}(r_2;w)}-c(r_2,r)\right]dr\right.\hfill \\ \hfill & =\hfill & I_1+I_2+I_3+I_4+I_5+I_6\hfill \\ \hfill & =\hfill & \delta (r_1-r_2)\hfill \end{array}$$

Carrying out each term in equation (44) separately,

(45)$$\begin{array}{ll}{I}_1\hfill & ={\displaystyle \int {\rho }^{(1)}}(r_1;w)\delta (r_1-r)\frac{\delta (r_2-r)}{{\rho }^{(1)}(r_2;w)}dr\hfill \\ \hfill & =\frac{{\rho }^{(1)}(r_1;w)}{{\rho }^{(1)}(r_2;w)}\delta (r_1-r_2)\hfill \\ I_2\hfill & =-{\displaystyle \int {\rho }^{(1)}}(r_1;w)\delta (r_1-r)c(r_2,r)dr\hfill \\ \hfill & =-{\rho }^{(1)}(r_1;w)c(r_2,r_1)\hfill \\ I_3\hfill & ={\displaystyle \int {\rho }^{(2)}}(r_1,r)\frac{\delta (r_2-r)}{{\rho }^{(1)}(r_2;w)}dr=\frac{{\rho }^{(2)}(r_1,r_2)}{{\rho }^{(1)}(r_2;w)}\hfill \\ I_4\hfill & =-{\displaystyle \int {\rho }^{(2)}}(r_1,r)c(r_2,r)dr\hfill \\ I_5\hfill & =-{\displaystyle \int {\rho }^{(1)}}(r;w){\rho }^{(1)}(r_1;w)\frac{\delta (r_2-r)}{{\rho }^{(1)}(r_2;w)}dr\hfill \\ \hfill & =-{\rho }^{(1)}(r_1;w)\hfill \\ I_6\hfill & ={\displaystyle \int {\rho }^{(1)}}(r;w){\rho }^{(1)}(r_1;w)c(r_2,r)dr\hfill \end{array}$$

in which the relation ∫f(x,y)δ(x − w)δ(y − w)dw = f(x,y)δ(x − y) was used in the first equation, while the equation, ${\int }_{-\infty }^{\infty }\delta (x-a)f(x,y)dx=f(a,y)$, in the other equations. Gathering the terms in I₁ + I₃ + I₅ + I₂ + I₆ + I₄,

(46)$$\stackrel{\delta (r_1-r_2)}{\overbrace{\frac{{\rho }^{(1)}(r_1;w)}{{\rho }^{(1)}(r_2;w)}\delta (r_1-r_2)}}+\frac{{\rho }^{(2)}(r_1,r_2)}{{\rho }^{(1)}(r_2;w)}-{\rho }^{(1)}(r_1;w)-{\rho }^{(1)}(r_1;w)c(r_2,r_1)+{\displaystyle \int ({\rho }^{(1)}(}r;w){\rho }^{(1)}(r_1;w)-{\rho }^{(2)}(r_1,r)c(r_2,r))dr=\delta (r_1-r_2)$$

Therefore,

(47)$$\frac{{\rho }^{(2)}(r_1,r_2)}{{\rho }^{(1)}(r_2;w)}-{\rho }^{(1)}(r_1;w)-{\rho }^{(1)}(r_1;w)c(r_2,r_1)+{\displaystyle \int ({\rho }^{(1)}(}r;w){\rho }^{(1)}(r_1;w)-{\rho }^{(2)}(r_1,r)c(r_2,r))dr=0$$

Rearranging this result, one has

(48)$$\frac{1}{{\rho }^{(1)}(r_2;w)}({\rho }^{(2)}(r_1,r_2)-{\rho }^{(1)}(r_2;w){\rho }^{(1)}(r_1;w))={\rho }^{(1)}(r_1;w)c(r_2,r_1)+{\displaystyle \int (-{\rho }^{(1)}(}r;w){\rho }^{(1)}(r_1;w)+{\rho }^{(2)}(r_1,r)c(r_2,r))dr$$

From equation (10), one has ρ⁽²⁾(r₁, r₂) = ρ⁽¹⁾(r₂;w)ρ⁽¹⁾(r₁;w)g(r₁, r₂;w) and considering the limit of a vanishing external potential, in which ρ⁽¹⁾(r) = ρ, this result reduces to ρ⁽²⁾(r₁, r₂) = ρ²g(r₁, r₂).Thus,

(49)$$\frac{1}{\rho }({\rho }^{2}g(r_1,r_2)-{\rho }^2)=\rho c(r_2,r_1)+{\displaystyle \int ({\rho }^{2}g(}{r}_1,r)-{\rho }^2)c(r_2,r)dr$$

As the total correlation function is defined as h(r) = g(r) − 1, one has

(50)$$h(r_1,r_2)=c(r_1,r_2)+\rho {\displaystyle \int h}(r_1,r)c(r_2,r)dr$$

which is the Ornstein-Zernike equation. Noting that equation (50) has two unknowns, it is necessary to define an additional relation called closure. Some examples of common closures used are Percus-Yevick, [²⁶] hypernetted chain, [²⁷, ²⁸] Kovalenko-Hirata [²⁹] and the PSE-n [³⁰].

8. Example of Lennard-Jones Fluid

A Lennard-Jones fluid consists of particles that interact through the Lennard-Jones potential, which is defined as:

(51)$$E_p(r)=4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]$$

with r = |r₂ − r₁| the distance between two particles, ϵ the depth of the potential well and σ the distance for which E_p(r = σ) = 0. A simple fluid composed of argon atoms will be discussed in this section. The parameters for argon atom are [³¹] $\frac{ϵ}{k_B}=119.8$ K and σ = 3.405 Å.

The hypernetted chain (HNC) closure relations, given by

(52)$$h(r)=e^{-\beta E_p(r)+h(r)-c(r)}-1$$

will be considered in this section for solving equation (50). One approach for acquiring this solution has been described in reference [³²]. As discussed previously, it is possible to see in Figure 2 that the direct correlation function is short-ranged while the total correlation function is long-ranged.

Figure 2
Solutions of Ornstein-Zernike equation for the Lennard-Jones potential with ρ = 0.02125Å⁻³ and T = 85 K.

From the discussion in Section From the discussion in Section 3 and the results for h(r), it is evident that there is a region (for small values of r) where particles are not allowed, known as the core region. Similarly, it is possible to identify regions where the fluid is locally denser than the bulk (around the maxima positions) and others where the fluid is depleted (around the minima positions). Therefore, from g(r), one can extract information about the solvation shells in the fluid and from Equation (16), the number of particles in these regions can be determined.

The result for S(q), calculated using equation (13), is given in Figure 3. It is possible to see that the calculated first peak position is in good agreement with experiment, while the second and third are slightly shifted. Also, the first peak intensity was underestimated. This shows that the HNC closure is not very appropriated for describing a Lennard-Jones fluid. There are other choices of closure that can be tested against experiment and the development of new closures is still an open field of research [¹², ³², ³³].

Figure 3
Calculated structure factor compared to experimental results [³¹].

From these results, some thermodynamic properties will be calculated. Using the HNC closure relation, an analytical expression can be derived for the excess chemical potential, μ^HNC [²⁹, ³⁴]. Table 1 presents the results for μ^HNC, the virial pressure (p_v), the isothermal compressibility (χ_T) and the number of atoms in the first solvation shell, along with their corresponding mathematical expressions.

Table 1
Results for excess chemical potential, virial pressure and isothermal compressibility calculated from OZ/HNC solutions.

Property	Expression	Calculated value
$\beta {\mu }^{HNC}$	$\rho \int \left[\frac{h{\left(r\right)}^2}{2}-C\left(r\right)-\frac{h\left(r\right)C\left(r\right)}{2}\right]dr$	3.44
$\frac{\beta p^v}{\rho }$	$1-\frac{\beta \rho }{6}\int \left[rg\left(r\right)\frac{d{E}_p\left(r\right)}{dr}\right]dr$	4.91
$\frac{\rho {\chi }_T}{\beta }$	$1+\rho {\displaystyle \int h(r)}dr$	0.097
N first solvation shell	$4\pi \rho {\int }_0^{R_{\text{first minimum}}}{r}^{2}g\left(r\right)dr$	11.44

The calculation of these properties is essential not only for predicting experimental outcomes but also for evaluating the performance of the closure relation against experimental data or molecular simulation results. As for example, the calculated number of atoms in the first solvation shell is approximately equal to the number of nearest neighbors (N = 12) in the face-centered cubic structure of crystalline argon [¹⁵]. Although the HNC closure does not give a perfect match for the structure factor, it does approximates well the number of atoms in the first solvation shell. Thus, the same closure may be good for predicting some properties while not as good for others.

9. Conclusion

The fundamental concepts of density and correlation functions were introduced and discussed, providing the foundation for understanding the Ornstein-Zernike equation, which was derived through three distinct approaches.

The first approach is the most intuitive and accessible, introducing the concepts of total and direct correlation functions based on particle interactions. This method relies on a straightforward reasoning, making it ideal for introductory undergraduate courses where the focus is on building an intuitive understanding of correlations in particle systems. The discussion emphasizes simplicity over mathematical rigor, which is suitable for students who are new to the field.

The second approach, while still appropriate for undergraduate students, introduces statistical concepts such as variance and covariance. This method connects the correlation functions to statistical measures, providing a richer physical and mathematical interpretation. It bridges the gap between the first method and more rigorous statistical frameworks, making it suitable for more advanced undergraduate courses.

The third approach is the most elaborated, requiring familiarity with functional derivatives and prior knowledge of statistical thermodynamics. It highlights the connection between microscopic (Hamiltonian) and macroscopic (Helmholtz free energy) quantities through correlation functions, offering a more comprehensive understanding of liquid state theory. This method is best suited for advanced students, such as graduate students, who are prepared to tackle the mathematical complexity and theoretical depth required to engage with modern research in the field.

Numerical solutions can be obtained using different methods, with the simplest approach involving the direct and inverse Fourier transforms, as discussed in references [³², ³⁵]. A faster and more stable alternative, called method of direct inversion in iterative subspace (MDIIS), is discussed in reference [³⁶].

Also, the application of the Ornstein-Zernike equation to a simple Lennard-Jones fluid was discussed. This integral equation theory can be extended to investigate more complex systems, such as molecular fluids [³⁷] and the solvation thermodynamics of molecules in solution [³⁸].

Acknowledgments

We would like to thank CNPq for the financial support.

Referencias

• ^[1] J.D. van der Waals and J.S. Rowlinson, On the Continuity of the Gaseous and Liquid State (Dover, New York, 2004).
• ^[2] G.A. Hirn, Théorie mécanique de la chaleur par G.-A. Hirn: Conséquences philosophiques et métaphysiques de la thermodynamique analyse élémentaire de l’univers (Gauthier-Villars, Paris, 1868).
• ^[3] J.D. van der Waals, The Equation of State for Gases and Liquids, available in: https://www.nobelprize.org/prizes/physics/1910/waals/lecture/, accessed in 20/01/2025.
  » https://www.nobelprize.org/prizes/physics/1910/waals/lecture/
• ^[4] J.P. Braga, Físico-Química – Aspectos Moleculares e Fenomenológicos (Editora UFV, Viçosa, 2002).
• ^[5] J.P. Braga, Termodinâmica Estatística de Átomos e Moléculas (Livraria da Física, São Paulo, 2021).
• ^[6] S.R.A. Salinas, Introdução à Física Estatística (EDUSP, São Paulo, 2013).
• ^[7] K.E. Gubbins, Molecular Physics 111, 3666 (2013).
• ^[8] J.H. Hildebrand, Annu. Rev. Phys. Chem. 32, 1 (1981).
• ^[9] J.W. Gibbs, Elementary Principles in Statistical Mechanics: Developed with Especial Reference to the Rational Foundations of Thermodynamics (C. Scribner’s Sons, New York, 1902).
• ^[10] T. Imai, M. Kinoshita and F. Hirata, J. Chem. Phys. 112 , 9469 (2000).
• ^[11] T. Luchko, N. Blinov, G.C. Limon, K.P. Joyce and A. Kovalenko, J. Comput.-Aided Mol. Des. 30 , 1115 (2016).
• ^[12] T. Tsednee and T. Luchko, Phys. Rev. E, 99, 032130 (2019).
• ^[13] L.L. Lee, Molecular Thermodynamics of Nonideal Fluids (Butterworth-Heinemann, London, 2016).
• ^[14] F.S. Carvalho and J.P. Braga, J. Mol. Model. 28 , 99 (2022).
• ^[15] J.P. Hansen and I.R. McDonald, Theory of Simple Liquids (Academic Press, New York, 2006).
• ^[16] P.C. Jordan and A.D. Greenberg, J. Chem. Phys. 56, 2181 (1972).
• ^[17] M. Brokate and M.Y. Rasulova, Phys. Part. Nucl. 41, 1014 (2010).
• ^[18] Y.L. Klimontovich, Statistical Theory of Open Systems. Volume 1: A Unified Approach to Kinetic Description of Processes in Active Systems (Springer Science+Business Media, Dordrecht, 1995).
• ^[19] L.S. Ornstein and F. Zernike, Proc. Acad. Sci. Amsterdam 17, 793 (1914).
• ^[20] J.K. Percus, Phys. Rev. Lett. 8, 462 (1962).
• ^[21] F.S. Carvalho and J.P. Braga, J. Mol. Model. 26, 193 (2020).
• ^[22] F.S. Carvalho, J.P. Braga, M.O. Alves and C.E.M. Gonçalves, Theor. Chem. Acc. 139, 29 (2020).
• ^[23] F.S. Carvalho and J.P. Braga, Braz. J. Phys. 50, 489 (2020).
• ^[24] A.E. Rodriguez and R.E. Caligaris, Teoria Estadística de Fluidos Simples en Equilibrio (Organización de los Estados Americanos, Washington, 1987).
• ^[25] R. Balescu, Statistical Dynamics: Matter Out Of Equilibrium (Imperial College Press, London, 1997).
• ^[26] J.K. Percus and G.J. Yevick, Phys. Rev. 110, 1 (1958).
• ^[27] J.M.J. van Leeuwen, J. Groeneveld and J. De Boer, Physica 25, 7 (1959).
• ^[28] E. Meeron, J. Math. Phys. 1, 165 (1960).
• ^[29] A. Kovalenko and F. Hirata, J. Chem. Phys. 110, 10095 (1999).
• ^[30] S.M. Kast and T. Kloss, J. Chem. Phys. 129, 236101 (2008).
• ^[31] J.L. Yarnell, M.J. Katz, R.G. Wenzel and S.H. Koenig, Phys. Rev. A 7, 2130 (1973).
• ^[32] F.S. Carvalho and J.P. Braga, Phys. A 576, 126065 (2021).
• ^[33] R.E.A. Goodall and A.A. Lee, Soft Matter 17, 5393 (2021).
• ^[34] S.J. Singer and D. Chandler, Mol. Phys. 55, 621 (1985).
• ^[35] F.S. Carvalho and J.P. Braga, J. Mol. Liq. 367, 120504 (2022).
• ^[36] A. Kovalenko, S. Ten-no and F. Hirata, J. Comput. Chem. 20, 928 (1999).
• ^[37] J.R.P.F.R. Fischer and S.R.H. Krienke, Mol. Phys. 93, 925 (1998).
• ^[38] T. Luchko, S. Gusarov, D.R. Roe, C. Simmerling, D.A. Case, J. Tuszynski and A. Kovalenko, J. Chem. Theory Comput. 6, 607 (2010).

Publication Dates

• Publication in this collection
  28 Feb 2025
• Date of issue
  2025

History

• Received
  15 Nov 2024
• Reviewed
  22 Jan 2025
• Accepted
  01 Feb 2025