Global dynamics of a delayed diffusive virus infection model with cell-mediated immunity and cell-to-cell transmission

1.   Introduction

A great number of experimental and theoretical studies of soft materials have suggested that diverse microscopic structures may exist in colloidal particle suspensions. With sizes between $10$ and $1000$ nanometers, the colloidal suspensions form complex fluids more sensitive to mechanical stress than the purely molecular systems. Generally, the colloidal particles may have well defined electric charges which must be compensated by those of solvent molecules to ensure the electrical neutrality of the system. Thus the complexity of these fluids is increased by the presence of dissolved ions, surfactant or polymer molecules and other small solutes in the solvent ^[1]. Aqueous solutions of micellar dispersions, globular proteins, microemulsions and diblock copolymer systems giving rise to systems separated by microphases with periodic density modulation are particularly interesting. Synthesized polymers such as polyethylene (PE), polystyrene (PS), polymethylmethacrylate (PMMA)... are present in many chemical realizations of biopolymers ^[2,3] like cellulose, proteins, actin filaments and deoxyribonucleic acid (DNA). They are used in various forms ^[4] such as gels, rubbers, synthetic fabrics, molded pieces, etc.

Soft materials are also scientifically important for testing and improving our understanding of self-assembly hypotheses more sophisticated than is needed to describe the molecules, see e.g., ^[5]. Self-assembly is a process in which components spontaneously form ordered agammaegates without human intervention. Although much of the self-assembly processes have been performed on molecular components, interesting applications of self-assembly can be found with particles of larger sizes ^[6].

Among the particles suspended in a solvent, one can single out two distinct types. The lyophobic particles abhor the solvent and undergo a weak dissolution. These particles carry charges of same sign causing a mutual repulsion that prevents them from approaching each other. The suspensions of arsenic trisulfide (As$_{2}$S$_{3}$) or clay are typical examples, with a high negative charge that precludes the particles from coagulating. Conversely, the lyophilic particles, such as starch suspensions, feel a strong attraction for the solvent. They are characterized by a significant solvation, i.e., a great ability to agammaegate with the solvent molecules, until formation of gels when the concentrations are high enough ^[3]. The difference in behavior of these two types of colloidal suspensions is related to the process of charge compensation and to the local changes undergone by the electrolyte solution in presence of charged particles.

As regards the conditions of stability in colloidal solutions, they strongly depend on mutual interactions between particles and ions. Unlike molecular systems where the range of interaction is comparable to the particle size, in colloidal suspensions the range of interaction is smaller than the particle size. Under certain circumstances, the sticky hard sphere model ^[7], albeit crude, has been used successfully in studies on colloidal suspensions ^[8]. As a result, the nanometer-sized particles form clusters and ordered colloidal suspensions widely depending on concentration, temperature, solvent composition, chemical constituents on the colloids, etc. Between colloidal particles, the main source of attractive forces is of the van der Waals type. But the colloidal suspensions must often be stabilized by repulsive forces to prevent the agammaegation of particles of the dispersed phase ^[3]. The two main mechanisms that produce the repulsive forces are electrostatic stabilization and steric stabilization ^[9]. Electrostatic stabilization is based on the repulsion of colloids carrying similar electrical charges. Steric stabilization consists in fixing the polymers having a high affinity with the solvent, on the colloidal particles, so that the latter, is getting close, would prefer to interact with the solvent and repel each other. Contrarily, in some cases, non-adsorbed polymer solution is expelled from the thin excluded volume surrounding each colloidal particle, causing a weak attractive depletion force between the colloidal particles, which is inherently entropic ^[10]. Thus, a colloidal suspension is said to be stable when a sufficiently strong repulsive force counteracts the van der Waals attraction. In other words, the stability of a colloidal suspension results from the balance between attractive interactions and repulsive interactions which are exerted on the particles.

On the side of experiment, the confocal microscopy ^[11,12,13] and the small angle neutron scattering ^[14] reveal the formation of clusters growing in size and becoming strongly anisotropic, eventually coalescing into elongated structures. For the advance in theory of these intriguing complex fluids, the cornerstone is the definition of the interparticle potential, which is then rigorously included in simulations done using the method of the statistical mechanics ^[15].

The known interaction models for colloids encompass an attractive well at short distances just behind the collision diameter, followed by a long-range repulsion. As already mentioned, the long-range attraction favoring the agammaegation is countered by the long-range repulsion, which prohibits further growth. As an example, the double Yukawa potential ^[16] with appropriate values of the parameters has the good features for studying the colloidal suspensions ^{[17,18,19,20]}. The long-range repulsion changes drastically the phase diagram of purely attractive potentials at long distances ^[21]. Specifically, in case of long-range forces between colloids being only attractive, the colloidal suspensions exhibit a stable liquid-vapor transition depending on the potential well. The inclusion of long-range repulsion makes possible metastable fluid–fluid transition and so complicates the thermodynamic phase transitions considerably ^[17,22,23]. The phase diagram of the colloid-polymer systems depend on the ratio of the colloid particle size $\sigma$ to the gyration radius $r_{G}$ of the polymer that corresponds roughly to the position of attractive well of the potential ^[1,10]. It has been shown ^[24] that the critical point disappears when $r_{G}/\sigma < 1/3$. A similar absence of liquid–gas transition is found for proteins and large molecules as the fullerene C$_{60}$. An effective way to determine the strength and range of interparticle interactions is to study the temperature-induced phase transition, where the inverse of attractive strength of the potential is used as a measure of the temperature.

This work is devoted to the features of a certain class of soft materials such as colloids, globular proteins, microemulsions, copolymers, etc. and their modelling in terms of effective pair potential. Our study is limited to spherical particles and interactions between them which depend on the distance only, not on mutual configurations. The only exception will be the inclusion of the interactions potential between two plates. For details on these questions, we refer the reader to the literature ^[25,26,27] and specially to the book of Israelachvili ^[28] that contains many problems and discussion topics. As the statistical-mechanical simulations proceeding from microscopic definitions do eventually lead to predicting macroscopic properties, some results of static structure and thermodynamics of the cluster fluid phases will be presented in the last part of the manuscript.

2.   Colloids, globular proteins, microemulsions, copolymers

Colloidal suspensions are involved in many phenomena (precipitation, flocculation, sedimentation...) playing an important role in the formation of clusters and in the transport of certain substances ^[4,14,29,30]. It is of great interest for environment to regulate sedimentation in estuaries silted up with sand, clay and silt or to control the transport of toxic metals in water or sol. As another example, various constituents of the cytoplasm form a colloidal suspension in the living cell, the dispersion medium of which consists of water containing mineral salts and monosaccharides in solution, and the dispersed phase is formed of protein macromolecules.

Certain colloidal suspensions are rich in spheroidal agammaegates of nanometric size formed of molecules having a hydrophilic polar head directed towards the solvent and a hydrophobic tail directed towards the core. In some cases, the molecules can self-assemble in aqueous medium to form micelles, liposomes or bilayer sheets (Figure 1). Liposomes (micelles consisting of concentric lipid bilayers enclosing aqueous compartments) are increasingly used in pharmaceutical research (to carry drugs and monitor their accumulation) and in medical imaging by introducing fluorescent-labeled substances.

As already mentioned, attractive interactions promote ordered phases, whereas long-range repulsions inhibit this tendency. At low density, long-range repulsions may lead to glasses and, combined with long-range attractions, a variety of others states may be found ^[31,32] as gels, glassy states of clusters, reminiscent states of Wigner glasses. When the density is increasing, clusters can become more and more anisotropic and lead, at the end, to gelation by forming a network with chainlike local structure ^[22,33]. For small ranges of interaction compared to the particle size, simulation computations have shown the presence of a Wigner glass of clusters observed in charged colloidal dispersions for very low concentration ^[34], in analogy with the Wigner crystals produced by long-range Coulombic interactions in the electron gas ^[35].

Globular proteins in aqueous solutions have been widely studied. A protein is a biological macromolecule formed of one or more polypeptide chains, each of them consisting of small subunits (amino acids) linked together, that constitutes their elementary structure. Different chemical properties determine the spatial arrangement of the amino acids along the polypeptide chain. A special case is represented by the spherical globular proteins, which are highly soluble in the solvent (Figure 2). They play an important role in metabolism. Representative of globular proteins are albumins and globulins, abundant in animal cells and blood serum. Hemoglobin is a protein whose primary function is the transport of oxygen in humans and other vertebrates, and which is responsible for the color of red blood cells. Lysozyme is encountered in a number of secretions (tears, saliva, breast milk, mucus...) and egg white in many species of animals. Enzymes are most often proteins, consisting of amino acid chains and acting in lowering the activation energy of a chemical reaction, thereby increasing reaction rate.

A number of techniques have been designed to produce proteins crystals. Systematic investigation of the process is difficult due to the slow time scale: the growth of a regular crystal may take weeks or months. Aside from the crystallization of dense colloidal suspensions, spherical clusters of proteins may also be formed at lower concentrations. Recent experiments indicate that clusters formed in dilute solutions exhibit no tendency of bulk phase transition ^[36,37]. In contrast, a liquid–liquid phase separation occurs when a protein solution is cooled below its cloud point temperature, i.e., the temperature above which the aqueous solution remains practically homogeneous and transparent, whereas below that temperature it separates into two coexisting liquids, one is rich in proteins and the other is poor in proteins. The onset of the formation of droplets via liquid–liquid phase separation (sometimes called coacervation) takes only few minutes. Note that liquid–liquid transition may be metastable because of the long time of equilibration; the demixion into two distinct liquid phases separated by a sharp meniscus may take many days. The clusters formation is different from the critical nucleus occurring during the vapour–liquid phase separation.

Of particular interest is the spontaneous microphase formation in systems of spherically symmetric particles lying at the interface of two simple liquids. These systems, known as microemulsions, are of industrial importance because they produce very low interfacial tensions (by cosolubilizing certain elements in the medium) and of scientific interest because the formation of their microstructures is not yet well understood. A microemulsion is a thermodynamically stable dispersion of two immiscible liquids, optically isotropic and stabilized by one surfactant (amphiphilic molecules). Since the typical sizes of mesophases (of few nanometers) are smaller than the light wavelength, the dispersions are generally transparent. In contrast, the macroemulsions made of large droplets of diameter greater than 1 $\mu$m, in the gravitational field, are thermodynamically unstable in the sense that the droplets sediment or float, depending on the densities of the dispersed phases.

A microemulsion contains a polar solvent (water), a nonpolar solvent (oil) and a surfactant, composed of long amphiphilic molecules having at one end a hydrophilic polar segment and at the other end a hydrophobic segment. Usually, a second amphyphilic component (cosurfactant) is introduced in order to ensure the flexibility of the interface and favour the microemulsion formation. By its tendency to locate between water and oil regions, the surfactant reduces the interfacial tension by several orders of magnitude and imposes a topological order of the microemulsions, but microemulsions remain geometrically disordered in the interstitial domains ^[38]. The role of the surfactant is to leave the dispersion macroscopically homogeneous and to arrange the mesophases on a microscopic scale. The microstructure has been examined by various techniques, e.g., light scattering, conductivity, diffusion by NMR, etc. as a function of the temperature ^[39,40,41]. It seems now admitted that, at low volume fraction of one component (water or oil) in the other, the picture of microemulsion structure is that of discrete swollen micelles of which the long-range order remains unpredictable. At intermediate volume fraction of water and oil, microemulsions may organize themselves into mesophases spatially periodic, separated by equilibrium bicontinuous structures like sheets between water-rich and oil-rich domains ^[38,42,43].

As in microemulsions, many exotic geometries have been observed in diblock copolymers consisting of two types of monomers. The monomers A and B are arranged such that there is a chain of each monomer grafted together to form a single copolymer chain (Figure 3). Ternary mixtures of A, B homopolymers and AB diblock copolymers are similar in some respects to mixtures of water, oil and surfactant. The copolymer acts as a surfactant for the A and B homopolymer mixture, decreasing the surface tension between the A-rich and B-rich domains, as it accumulates on the interface, and maintaining a topological order of the domains. Above a certain temperature, the large collection of diblock copolymers is equally distributed in a disordered state, like in liquids. Below that temperature, when entropy is dominated by enthalpy, the diblock copolymers form regular structures common to the microemulsions ^[44,45,46]. Despite this similarity, the diblock copolymer is generally not a solubilizer of homopolymers as efficient as good amphiphile in the oil and water mixture.

3.   The Gouy–Chapman model for charged particles

In the colloidal suspensions, the repulsive part of the potential generated by the electrostatic forces plays a more significant role than for other simple liquids. The interface between a charged particle and the solvent is the seat of a potential difference that, in the simplest model, gives rise to an electrical double layer: on one side of the interface the charge is clearly positive while on the other side it is clearly negative. Therefore, the electrical double layer is equivalent to a plane capacitor. The first quantitative study of the electrical double layer has been achieved by Helmholtz (1879) who considered a solid immersed in a solution. The study showed that the charged particles do not form a rigid layer, because of the thermal agitation. Later on, Gouy ^[47] and Chapman ^[48] independently developed a rational model showing that the charged particles farthest from the interface may come and go more freely than those near the interface. This is the diffuse layer model. Note in passing that the Gouy–Chapman model is based on a similar model developed by Debye and Hückel ^[49] for strong electrolytes.

The Gouy–Chapman model deals with the case of colloidal particles surrounded by a diffuse layer. Like the Debye–Hückel model, the Gouy–Chapman model is based on the solution of the simultaneous Boltzmann and Poisson equations that, together, make possible the calculation of the electric potential $\Phi({\bf{r}})$ in terms of the coordinates ${\bf{r}}$ and ionic concentration $c$ in the solution.

The Boltzmann equation results from statistical considerations but it could be otherwise supported by a thermodynamic argument, as a direct consequence of the equilibrium of a system of charged particles ^[28]. In such a system not subject to an interparticle interaction, it is judicious to consider the electrochemical potential ($\mu _{i}+z_{i}e\Phi _{i}),$ where $\mu _{i}\left( = \mu ^{0}+k_{B}T\ln c\right)$ is the chemical potential of the component $i$ and $z_{i}e\Phi _{i}$ the electrical work. If $c_{A}$ and $c_{B}$ are supposed to be the concentrations of the same constituent at two distinct points of the solution, and $\Phi _{A}$ and $\Phi _{B}$ the electric potentials at these respective positions, the thermodynamic equilibrium is given by the equality of electrochemical potentials:

Moreover, assuming that $\mu _{A}^{0} = \mu _{B}^{0}$ and that the electric potential in the unperturbed region, say at the point $B$, is $\Phi _{B} = 0,$ Eq 1 simplifies to a form that is precisely the Boltzmann distribution:

Thus, the Boltzmann distribution provides the variation law of the concentration of a constituent into the diffuse layer, as a function of the concentration $c_{B}$ in the unperturbed region and electric potential $\Phi _{A}$ at the point where the concentration $c_{A}$ is calculated.

Concerning the Poisson equation, which is a consequence of Coulomb's law, it connects the electric charge density $\rho({\bf{r}})$ with the electric potential $\Phi({\bf{r}})$ for each point of the diffuse layer. The local form of the Poisson equation is $\nabla^{2}\Phi({\bf{r}}) = -\rho({\bf{r}})/(\varepsilon_0\varepsilon)$, where $\varepsilon_0\varepsilon$ is the dielectric constant of the suspension and $\nabla^2$ the Laplacian operator in Cartesian coordinates.

In the Gouy–Chapman model, the colloidal particles are treated as parallelepipeds of thin thickness (Figure 4). For ensuring that the equipotentials in the diffuse layer are parallel to the larger surface of the parallelepiped, the Poisson equation at one dimension is used:

If the electric charge is assumed to be uniformly distributed on the surface of the colloidal particle, the charge density $\rho (x)$ to be introduced in Eq 3 is:

where $n_{i}(x)$ is the number of ions of species $i$, per unit volume, given by the Boltzmann distribution (Eq 2), and $z_{i}$ the charge of ions of species $i$ assigned of the corresponding sign. The parameter $n_{0i}$ is the number of ions, per unit volume, in the region not perturbed electrically far away from the diffuse layer, and $\Phi (x)$ the electric potential at the point where $n_{i}(x)$ is calculated. By combining Eqs 4 and 3, the differential equation giving the variation of electric potential, as a function of the abscissa $x$, reads:

Without changing the generality of the calculation, we assume that the numbers of ions of two species are equal ($n_{0+} = n_{0-} = n_{0}$), so that the previous differential equation becomes:

On multiplying both sides of the above equation by $2\frac{d\Phi(x)}{dx}$, the first integration is easily performed. The result is:

We then carry out the changes of variables $\alpha = {e}/(2k_{B}T)$ and $\kappa = \alpha\left[{2n_{0}k_{B}T}/({\varepsilon_0\varepsilon}) \right]^{ \frac{1}{2}}$, as well as $\left[1+\exp(\alpha\Phi_0)\right] = A$ and $\left[ 1-\exp(\alpha\Phi_0)\right] = B$, where $\Phi_0$ is the electric potential on the surface of the colloidal particle. With the boundary conditions at infinity: $\Phi = 0$ and ${d\Phi(x)}/{dx} = 0$, the integration of Eq 7 reduces^* to:

^* Noting that the integration constant is $C = -2$, the term in braces of Eq 7 is the perfect square:

and the integral can be easily done by separation of variables.

As expected, in the case of a negatively charged colloidal particles (clay suspension), the electrical potential is always negative, with the value $\Phi_{0}$ at $x = 0$. Equation 8 can be then used to find the asymptotic form of the electric potential $\Phi \left(x\right)$ :

showing that $\Phi(x)$ decreases exponentially, tending to $0$ by negative values when $x$ goes to infinity (Figure 5). It follows from this function that the parameters $\alpha$ and $\kappa$ have the dimensions of the inverse of an electrical potential and the inverse of a length, respectively: recall that $e$ is expressed in C, $k_{B}T$ in J, $n_{0}$ in m$^{-3}$ and $\varepsilon_0\varepsilon$ in C$^{2}{\cdot}$N$^{-1}{\cdot}$m$^{-2}$. As a consequence of Eq 9, we define the pair potential $u(x)$ between two particles as the product of the electrical potential $\Phi(x)$ by the charge of the particle $-Q$, namely:

The cation and anion distributions in the diffuse layer are easy to determine by eliminating $\Phi(x)$ between Eqs 8 and 4. The distributions calculated here correspond to the case of a clay suspension for which the surface is negatively charged, for given values of $n_0$ and $\phi_0$. If we note that

the distribution of anions (Eq 4, with $z_{i} = -1)$ can be written with the help of Eq 8 as:

and the distribution of cations (with $z_{i} = +1)$ is:

It should be noted that the curves representative of the distributions $n_{+}(x)$ and $n_{-}(x)$ tend both toward the asymptotic limit $n_{0}$ when $x$ goes to infinity, each curve approaching the asymptote from a different side because $n_{+} > n_{-},$ $dn_{+}/dx < 0$ and $dn_{-}/dx > 0$ (with $\Phi < 0$ and $d\Phi /dx > 0,$ see Eq 9). Schematic representations of the ionic concentrations as a function of the distance are shown in Figure 6. At low distances, the two curves are strongly asymmetric with respect to the asymptote $n_{0}$. For a colloidal particle negatively charged, this means that the concentration of cations in the diffuse layer is greater than the concentration of anions. In this specific case, the negative charge of the clay suspension is balanced by an excess of cations and a deficit of anions, with respect to the concentration $n_{0}$ of the unperturbed solution. For ions in the solution of same sign as the colloidal particle, it is said that the adsorption is negative, while for ions of opposite sign to the colloidal particle the adsorption is positive.

The previous expressions of the electric potential and ionic concentrations correspond to the case where the colloidal particles are far from each other, so that their mutual interactions are neglected. If the flat and parallel faces of two identical particles approach one another, the diffuse layers interfere with each other and the electrical potential is modified. The electric potential between two particles is portrayed by the full line in Figure 5.

4.   Repulsive potential of Verwey and Overbeek

The method used by Verwey and Overbeek ^[50] to determine the expression of the repulsive part of the potential is the same as that of Gouy ^[47] and Chapman ^[48], with the difference that the colloidal particles are spherical. Over the last seventy years, the stability of colloidal suspensions has been extensively studied with the repulsive potential developed by Derjaguin and Landau ^[51] as well as by Verwey and Overbeek ^[50]. This potential was originally constructed to describe the effective interactions between spherical colloidal particles in an aqueous monovalent electrolyte solution $(n_{0+} = n_{0-} = n_{0})$, which has given rise to the concept of electrostatic screening. In this section, we present a neat demonstration of this repulsive potential due to Denton ^[52], which is based on an analogy between the metallic liquids and the colloidal suspensions ^[53,54]. Indeed, if the macroions^† are assumed to be tantamount to the metallic ions and the microions to the electrons, it is logical to continue the analogy by using the formalism of the screening theory (linear response theory) to derive the effective potential between the macroions.

^† The terms frequently used in literature are macroions for the colloidal particles and microions for the ions in solution. These are the ones that will be used here.

The system studied by Denton consists of $N_{M}$ macroions and $N_{m}$ microions occupying a volume $V$, at temperature $T$. The macroions and microions are immersed in a solvent constituted by a uniform continuous medium of dielectric constant $\varepsilon_{0}\varepsilon$. The macroions are charged hard spheres of diameter $\sigma$, carrying a charge $-Ze$ uniformly distributed on the surface. Note that the charge of macroions is negative, contrary to the metallic ions. All the point microions carry the same charge $ze$. Therefore, the overall electrical neutrality condition reduces to $ZN_{M} = zN_{m}.$

By analogy with the liquid metals, the Hamiltonian of the system may be evaluated with the energies of interaction of the different components ^[15]. The potential for the direct interaction between two macroions is expressed as:

and the potential between two microions is:

Regarding the potential between a macroion and a microion, it is chosen under the following form proposed by van Roij and Hansen ^[55]:

The potential is Coulombic outside the macroion and constant inside. It may be noted that this potential differs from the Ashcroft pseudopotential used in liquid metals ^[56], where the value zero in the ionic core is justified by the presence of core electrons. As expected, the parameter $\alpha$ must be determined afterwards, within the condition that the microionic density profile in the heart of the macroion is zero.

After defining the potentials between the different entities, the response function theory may now be used in order to investigate how the system reacts to the perturbation caused by the presence of macroions. This involves calculating the linear response function $\chi^{(1)}(q)$ of the classical one component plasma (OCP), composed of a mixing of macroions and microions. Unlike the method employed to determine the response function $\chi^{(1)}(q)$ for liquids metals, for the classical OCP the general formalism ^[57,58] must be used, which connects the response function $\chi^{(1)}(q)$ with the structure factor $S(q)$ by the formula:

where $n_{m}$ is the effective density of microions, and $c(q)$ the Fourier transform of the direct correlation function $c(r)$ that will be defined in the following.

The calculation of the effective density $n_{m}$ of microions is performed in considering only the volume $V_{m}$ available to the microions, i.e., the volume $V$ of the system excluding the macroion volume ($\frac{\pi \sigma ^{3}}{6}N_{M})$ not occupied by the microions. So, the expression of the effective density is:

As regards the function $c(q)$, it can be estimated with the mean spherical approximation (MSA), proposed by Lebowitz and Percus ^[59]. The MSA of statistical mechanics amounts to identifying the direct correlation function $c(r)$ as the scaled potential $u_{mm}(r)$ by means of the formula $c(r)\simeq -\beta u_{mm}(r)$. Thus the function $c(q)$ that is the Fourier transform of $c(r)$ reads:

The substitution of the Eq 19 into Eq 17 leads to the linear response function $\chi^{(1)}(q)$:

where:

With the response function $\chi^{(1)}(q)$, it becomes possible to determine the contribution of the potential induced by the microions, $u_{ind}(r)$. If the notations of the liquid metal are transposed to the colloidal suspensions, $u_{ind}(r)$ reads:

where $u_{Mm}(q)$ is the Fourier transform of the potential $u_{Mm}(r)$ defined by Eq 14. The calculation of $u_{Mm}(q)$ is conducted as that of the Ashcroft form factor. Its outcome is:

Note the presence of the unknown parameter $\alpha$ in the previous relation. As it has been indicated, its determination is made possible by requiring that the microionic density profile cancels out in the core of macroions. To ensure this condition, one can calculate first the density profile of microions in reciprocal space with the standard relation:

where ${\bf{R}}_{l}$ is the position of center of mass of the $l$-th macroion. The substitution of $\chi^{(1)}(q)$ and $u_{Mm}(q)$ in Eq 24 gives $n_{m}(q),$ namely:

In order to determine the density profile of the microions in real space, the inverse Fourier transform of $n_{m}(q)$ is carried out. However, to simplify the calculations, Denton^[52] suggested to single out the macroion at ${\bf{R}}\! = \!0$ and neglect the other terms in the sum over $l$. The result of the inverse Fourier transform is^‡:

^‡ The Fourier transform and inverse Fourier transform are defined by the expressions:

Taking into account that the density profile is zero when $r < \sigma/2$, the expression of the unknown parameter $\alpha$ reduces to:

Consequently, if $\alpha$ is inserted in Eq 23, the expression of $u_{Mm}(q)$ becomes:

Then, the potential $u_{ind}(r)$ induced by the microions is obtained by inserting the expressions of $\chi ^{(1)}(q)$ (Eq 20) and $u_{Mm}(q)$ (Eq 28) into (Eq 22), namely:

The calculation of the inverse Fourier transform is simple though tedious. The result is:

Finally, by adding the induced contribution $u_{ind}(r)$ to the direct contribution $u_{MM}(r),$ given by Eq 14, we obtain the repulsive potential between macroions beyond $\sigma$ in the explicit form:

This expression is the same as that found by Verwey and Overbeek with the Boltzmann and Poisson equations. While the Gouy–Chapman repulsive potential (Eq 10) between the parallel faces of two charged particles, $u(x)$, is a purely Coulomb one, the Verwey-Overbeck repulsive potential (Eq 31) between two spherical colloidal particles, $u_{rep}(r)$, is screened Coulomb potential decreasing more rapidly than $u(x).$ In addition, both potentials are very sensitive to the value of the salt concentration in the solution^§, decreasing with $r$ much more rapidly when $\kappa$ is large, i.e., when the effective density and the valence of the microions are more important.

^§ The salt concentration $c$ is expressed in moles per liter (M) and the average number of ions $n_{0}$ in ion per m$^{3}$. These two quantities are related to each other by the formula:

where $\mathbb{N}$ is the Avogadro number.

5.   Attractive potential of Hamaker

Beside the repulsive potential of Verwey and Overbeek ^[50], which depends on the salt concentration in solution, there is an attractive potential attributed to microscopic forces acting between the atoms of the colloidal particles. Remember that polar molecules form permanent dipoles tending to align themselves under the action of the van der Waals attractive forces. As regards the nonpolar molecules (and atoms), they form instantaneous dipoles capable of inducing dipoles in neighboring atoms, which are the cause of dispersion forces of London ^[60]. For the molecules and atoms, the van der Waals and London forces extend to distances that do not exceed one nanometer.

When two colloidal particles approach one another, each atom belonging to the first particle interacts with all the atoms of the second particle in a manner more or less additive. As we shall see, the sum of all these microscopic interactions is likely to generate an attractive force between the colloidal particles, called the Hamaker force ^[61]. Whereas the London forces between the atoms have ranges of one nanometer, the Hamaker forces between the colloidal particles extend over a range of hundred nanometers while decreasing more slowly. This range is comparable to that of the repulsive Verwey and Overbeek forces.

In order to derive the expression for the Hamaker force, we first calculate the potential between an atom of the colloidal particle 1 and all the atoms of the nearest colloidal particle 2 (Figure 7). Taking account of the potential between two neighboring atoms resulting from the formation of two instantaneous dipoles ^[60], $u(r) = -C_6/r^6$, we may compute the potential between the atom $1$ of the particle $1$ and the particle $2$ by performing the sum over all the atoms of the particle $2$:

To obtain the result, the discrete sum over the atoms is replaced by an integral over the volume of particle $2$ by considering the number of atoms per unit volume $\rho$ and the volume element $d\tau$:

Then, the integration element $d\tau$ in the above equation is taken to be an infinitesimal volume between two concentric spherical caps intercepted by the angle $2\theta$, of which the volume element is:

If the atom of interest (atom $1$) is close to particle $2$, it can be taken that it sees a flat surface. Therefore the potential between atom $1$ and particle $2$ is reduced to the following expression:

It should be noted that this potential decreases as $1/r^{3}$ whereas the potential between two dipoles decreases as $1/r^{6}$. But to determine the potential between two colloidal particles, it remains to carry out the sum over all the atoms of particle $1$, that is tantamount to integrating over the volume of particle $1$. If $S$ denotes the area of the two opposite faces, the integration element is $dV = Sdr$ and the potential between two colloidal particles, per unit area, is equal to the integral of Eq 35 taken between $r$ and infinity:

where $A\ ( = \pi ^{2}\rho C_{6})$ is the Hamaker constant. Equation 36 yields the attractive potential, per unit area, between two colloidal particles, the surfaces of which are flat and face to face. The variation law $1/r^{2}$ should be noted. Note also that the Hamaker constant $A$ has the dimensions of energy. It varies little with the nature of particles and solvent, and its numerical value is of the order of the thermal energy $k_{B}T$. For example, it has the value $A = 6.3 \times 10^{-21}$ J for quartz–water–quartz interfaces and $A = 10^{-18}$ J for copper–vacuum–copper interfaces. The values of the Hamaker constant may be found in literature for all kinds of surfaces separated by a vacuum or a solvent ^[62].

There are alternative Hamaker formulas for other shapes of particles but they are more complicated ^[28,63]. We only mention the expression of Hamaker's potential between two spheres of diameter $\sigma$:

When $r > \sigma$, Eq 37 is written in the approximate form:

so that, for $r\gg \sigma$, the potential behaves as:

and for $r\simeq \sigma ^{+},$ as:

When two spherical particles with diameters $\sigma _{1}$ and $\sigma _{2}$ are close to each other, the following formula may be used:

In contrast to the repulsive potential of Verwey and Overbeck, the attractive potential of Hamaker between two spherical colloidal particles decreases as $r^{-2}$ and depends very little on the nature of the electrolyte, while they extend over comparable ranges.

6.   Casimir potential

It has been mentioned that the dispersion forces have their origin in the fluctuations of atomic dipole moments. The process is as follows. The electric field created by the dipole moment of an atom, propagating at the speed of light, induces a dipole moment on close neighboring atoms, the direction of which is the same as the first dipole moment. The attraction between them is maximum if the two dipole moments are aligned, as in two polar molecules. But in the case of two atoms far away from each other, there is a misalignment of the instantaneous dipole moments. This results in a reduction of the dispersion forces because the electric field emitted by the second atom arrives at the first one with a certain phaseshift. This process gives rise to a delayed effect of dispersion forces ^[64] hardly affecting the interaction between the atoms. By contrast, it can significantly alter the interactions between colloidal particles when the separation distances are greater.

The first application of the delayed effect to the calculation of the interactions between colloidal suspensions has been studied with Eq 36 required for two parallelepipedical colloidal particles located face to face, with the delay function $f(p)$ as ^[65]:

where $f(p)$ is an empirical function varying between $1$ and $0$ when the parameter $p$ increases from $1$ to infinity. The parameter $p$ yields a relative measure of the distance $r$ separating the plates with respect to the wavelength $\lambda$ of the radiation of the electromagnetic field. With the parameter $p$ defined by the relation $p = 2{\pi }r/{\lambda }$, the empirical function $f(p)$ is expressed as:

Incidentally, it can be noted that the dispersion forces cause the appearance of a particular type of forces acting in the vacuum. By investigating the fluctuations of the electromagnetic field, Casimir ^[66] showed that a macroscopic force must be exerted between two parallel plates (conductive and uncharged) in vacuum, in equilibrium with the photons of the vacuum. Qualitatively the phenomenon might be explained as follows. The conductive plates form a resonant cavity, imposing the boundary conditions $L = p\lambda$ on the photons, where $L$ is the distance between the two surfaces, $\lambda$ the wavelength of a photon and $p$ an integer. Among the photons contained in the enclosure, only those with a wavelength less than $L$ are confined between the plates. Therefore, the deficit of photons between the plates gives rise to a reduction in pressure and an attraction of the plates.

In quantum mechanics, the energy levels of a harmonic oscillator are given by the relation:

where $\hbar \omega$ $= h\nu$, and $n$ is the number of photons present in the mode of pulsation $\omega$. However, at zero temperature, although there is no photon in the enclosure, the energy of the electromagnetic quantum vacuum at zero-point ($n = 0)$ is $\frac{1 }{2}\hbar\omega$, generating a pressure at the origin of a macroscopic force between the plates.

One of the laws of black body radiation indicates that the number of photons in an enclosure is proportional to the cube of the temperature. Thus, at zero temperature no real photon is present.

It should be mentioned that the functional form of the pressure created between the plates, at zero temperature, can be obtained by a simple dimensional analysis. Noting that the effect is of quantum and relativistic origin, it may be expected that the pressure $F/S$ (in N$\cdot$m$^{-2}$) depends on the universal constants $\hbar$ (in J$\cdot$s) and $c$ (in m$\cdot$s$^{-1}$), as well as the spacing $L$ (in m) between the plates ^[67,68]. The equation of dimensions binding these four variables is:

The exact formula ^[69] of the pressure indicates that $F/S$ follows a universal law independent of the nature of the conductors, the expression of which is:

The sign of the pressure shows that the Casimir force is attractive and that the potential $u_{att}^{C}$, per unit area, is of the form:

Expressions for the Casimir force have been obtained for alternative configurations, and the temperature dependence has been also studied. As an example, it has been shown that the expression for the attractive force exerted between a plate and a sphere of radius $R$ is $F = (\pi ^{3}/360)R(\hbar c/L^{3})$. Between a plate and a sphere of radius $R = 98\, { \mu }$m, the force is $F = 33 \times 10^{-12}$ N when $L = 200$ nm. Although this force is very small, its magnitude is comparable to that of the forces involved in biological systems, which are perfectly measurable by atomic force microscopy. Inspired by the classical Casimir effect, an additional mechanism for the confinement of density fluctuations in an electrolyte solution by the walls has been recently investigated ^[70].

7.   Interparticle interaction potentials

The most common expression used for the potential between two spherical colloidal particles of diameter $\sigma$ consists of (ⅰ) the hard-sphere potential $u_{HS}(r)$ of diameter $\sigma$, (ⅱ) the purely repulsive potential $u_{rep}(r)$ resulting of the diffuse layer (Eq 31) and (ⅲ) the attractive potential $u_{att}(r)$ caused by the dispersion forces (Eq 37). Together, these three contributions make up the DLVO (Derjaguin-Landau-Verwey-Overbeek) potential, the form of which is:

It is clear that the stability of colloidal suspensions is accounting for the competition between the attractive and repulsive contributions to the DLVO potential. Figure 8 displays a schematic representation of the DLVO potential for two different electrolytes. If the electrolyte concentration is important, $\kappa$ is large and the decay of the repulsive part of the DLVO potential is rapid. For large concentrations in salt ($c\simeq 1$ M), the attractive part of the DLVO potential prevails, which helps to explain the irreversible agammaegation of macroions in such solutions. Therefore, the repulsive part of the DLVO potential should be overshadowed by the attractive part, thus reducing the potential barrier. For lower salt concentrations, the repulsive part is increasing and gives rise to a higher potential barrier ($\sim k_{B}T$) that prevents the contact of the macroions. This explains the reversible liquid–gas transition observed in electrolyte solutions of salt concentration less than $10^{-3}$ M. For salt concentrations even lower, the attractive Hamaker potential is partially masked by the repulsive contribution. In fact, all experiments with electrolyte solutions where the salt concentration is about $10^{-5}$ to $10^{-6}$ M have shown that the DLVO potential has a long repulsive barrier ^[71,72]. Nevertheless, it should be emphasized that the tail of the DLVO potential for long distances is always negative because the attractive Hamaker potential decreases less rapidly than the repulsive Verwey and Overbeck potential.

Among the mesoscopic systems, polymers are also objects studied with great attention. In solution, the polymers can be found under the form of chains of macromolecules ($i$) unfolded and entangled with each other or ($ii$) folded in on themselves, forming slightly inter-penetrable spherical objects. The classic model of a spherical polymer, developed by Daoud and Cotton ^[73], consists of a sphere of radius $R_{c}$ (called the crown radius), inside which the density profile of the monomers follows a power law established by the authors. In turn, the sphere is surrounded by a diffuse layer made of $f$ monomers, such as hair stood on head, where the density profile is characterless.

As regards the interactions acting between spherical polymers, the model of Asakura and Oosawa ^[74] is usually employed. Initially, these authors studied the interaction between two large spheres immersed in a fluid of small particles not interacting between themselves. They provided the explanation for the following mechanism: when the distance separating two large spheres is less than the diameter of small particles, the latter are expelled from the intermediate region between the large spheres, producing an anisotropy of the local pressure and an attraction force between the large spheres that gives rise to the depletion potential.

Later, Pincus ^[75] developed a potential describing the interaction between a polymer ball and a wall. This potential has the universal logarithmic form, $u(r)\sim -k_{B}T\ln(r/R_{c})$, valid at any distances. More recently, a purely repulsive potential has been proposed to represent the interaction between two spherical polymers ^[76]. Its expression is of the form:

where $f$ is the number of monomers at the periphery of the sphere of radius $R_{c}$. Outside the sphere, the potential is made of a screened Coulomb potential, whereas inside the sphere the logarithmic form of interaction with a wall is preserved.

Regarding the molecular fluids, an early study has been done to investigate the behavior of carbon-60 and the possible existence of its liquid phase ^[77]. The pair potential between two molecules of C$_{60}$ is composed of the sum of all interatomic potentials, the interatomic interaction being modelled by the Lennard–Jones potential. This potential is characterized by a strong repulsion and an attractive part substantially shorter than that of the Lennard–Jones potential. It clearly demonstrates that the existence of the liquid phase of C$_{60}$ is strictly related to the extension of the attractive part of the potential.

Other potentials with an attraction at intermediate range and a repulsion at long range can also produce peculiar phase diagrams. Unlike to the DLVO potential, the tail of these potentials for large distances is slightly positive. Such potentials are particularly useful for describing the behavior of colloidal suspensions, spherical polymers, protein solutions, etc. In this regard, one may infer that the competition between the attractive and repulsive parts, on different length scales, is responsible for the formation of inhomogeneous fluid phases, such as gels in non-equilibrium systems ^[78], microphase separation ^[17], clusters in equilibrium systems ^[19], etc. In solutions of charge-stabilized colloids and non-adsorbing polymers, the long-range repulsion is due to the charged colloidal particles, while the long-range attraction is caused by the solvent containing non-adsorbing polymers responsible for the depletion effect. It should be noted that the issue concerning the competition between long-range attractive and long-range repulsive parts of the potential had been already raised by Lebowitz and Penrose ^[21] fifty years ago. From an exact theory, the authors noticed that the presence of a long-range repulsive potential could bring the system to break down into droplets or froth, without causing a first order liquid–vapor transition. While this potential makes sense for complex liquids (charged colloidal particles, spherical polymers, globular proteins...), it is irrelevant for simple liquids.

To end this discussion, we would like to mention the influence of the general form of interparticle potential on the phase diagram of a fluid of identical particles. The respective roles of the repulsive and attractive contributions of the Lennard–Jones potential on the stability limits of the phase diagram of fluids are known since long time. They deal essentially with the liquid–vapor and liquid–solid transitions. However, in recent years, it has been shown that some empirical potentials, hardly more complicated than the Lennard–Jones potential, can produce phase diagrams radically different from those of simple liquids, giving rise to new phase transitions. Despite their spherical symmetry, such potentials can also generate inhomogeneous fluid phases ^[10,79].

Much effort is currently devoted to understanding the agammaegation properties and to predict a variety of locally inhomogeneous structural states at equilibrium by looking at the features of potential models. In addition to the square-well-linear model ^{[80,81,82,83]} and other piecewise-continuous models ^[84,85,86], here are three models often used: double Kac potential ^{[19,87,88,89]}, double Yukawa potential ^{[17,19,90,91,92]} and Lennard–Jones-like plus screened Coulomb potential ^{[78,93,94,95,96,97]}.

The double Kac model is formed by the hard-sphere potential followed by two Kac potentials, of which the overall analytical expression is ^[87,89]:

where $\beta\ ( = 1/k_{B}T)$ is the inverse temperature, $\sigma$ the hard-sphere diameter, $\gamma _{a}$ and $\gamma _{r}$ characterize the extent of the interactions and $\varepsilon _{a}\gamma _{a}^{3}$ and $\varepsilon _{r}\gamma _{r}^{3}$ their strengths. The choice of these parameters describing short-range attraction and long-range repulsion should respect the relations $\varepsilon _{r}\gamma _{r}^{3} <$ $\varepsilon _{a}\gamma _{a}^{3}\ll 1$ and $\gamma _{a} >$ $\gamma _{r} > 0.$ Another potential frequently used is the double Yukawa model ^[16] of the following form ^[18,93]:

Its strictly positive parameters $z_{i}$ and $K_{i}$ characterize, respectively, the range and the magnitude of the attractive and repulsive contributions. The smaller $z_{i}$, the greater the extension of the interaction. By a judicious choice of these parameters it is easy to generate a deep short-range attraction and a small repulsive barrier at long range. The formation and growth of equilibrium clusters in suspensions of weakly charged colloidal particles and/or small non-adsorbing polymers may be also studied with the potential composed of the Lennard–Jones and the repulsive Yukawa potentials (LJY):

with $\alpha > 6$ to ensure a very long-range attraction. Different choices of the parameters $\varepsilon$ and $\alpha$ as well as $\xi$ and $A$ lead to different scenarios such as crystal made of finite-sized clusters, lamellar phases ^[98] and elongated structures forming a connected network (gel) ^[94].

8.   Phase transitions

It is fascinating that a small variation of a parameter, like temperature, causes in some cases a spectacular modification of the system and the appearance of new phases of matter whose collective behavior bears little resemblance to that made of few particles. The liquid state theory aims to identify both range and strength of the interactions among particles in terms of experimentally accessible parameters, such as molecular length, density and temperature. But it is rare to find an interacting many-particle system likely to be described easily by a microscopic approach near the critical points. More useful is the macroscopic approach having its origin in the thermodynamic theory of critical phenomena, which ignores the interactions between particles. Since all macroscopic properties can be deduced from the free energy, their dramatic changes involved in phase transitions correspond to singularities in the free energy. As with any thermodynamic system, there is a competition between entropy and enthalpy, whose balance corresponds to the minimization of free energy. Hence, the study of phase transitions is related to finding the origin of these singularities and to the identification of new phases.

Among the numerous phase transitions, the liquid–gas transition is well known with the line in the ($p, T$) plane (Figure 9a) that ends at the critical point ($T_{c}$, $p_{c}$), and the coexistence interval containing a mixture of gas and liquid of densities $\rho _{G} = 1/v_{G}$ and $\rho _{L} = 1/v_{L}$, at $T < T_{c},$ in the ($p$, $\rho$) plane (Figure 9b). Since the liquid–gas coexistence curve $p(T)$ terminates at a critical point, there is a possibility to convert liquid to gas continuously without crossing the phase transition line, that is not possible for the liquid–solid transition line. Another interesting phase transition occurs in the Ising ferromagnetic system. The system is in a paramagnetic (disordered) phase for a temperature greater than the critical temperature (Curie temperature $T_{c})$, while it remains in a ferromagnetic (ordered) phase below $T_{c}$ (Figure 9c). If a magnetic field $H$ is applied to a ferromagnetic material (Figure 9d), a discontinuity occurs in magnetization density $m = M/V$ as the magnetic field $H$ goes to zero, for $T < T_{c}$. From these figures, it can be seen that the isotherms $p(\rho)$ and $H(m)$ have much in common due to the similarity between $H$ and $p\ $ and between $m$ and $\rho$. Note that jumps in the magnetization $m$ and the density $\rho$ decrease when the temperature increases towards $T_{c}$.

In spite of apparent similarities between these two phase transitions, they differ in their specific features. Generally, to distinguish two phases of a material one uses an order parameter, which is equal to zero in the disordered phase and non-zero in the ordered phase (usually the low temperature one). For the liquid–gas transition, the order parameter ($\rho -\rho _{c})$ is defined as the difference between the density $\rho$ and the critical density $\rho _{c},$ and for the ferromagnetic transition, the role of the order parameter is played by the magnetization density $m$. In the case of the liquid–gas transition, the order parameter has the dimensionality $n = 1,$ since the density is a scalar, while for the ferromagnetic transition, the order parameter has the dimension $n = 3$ because of the vectorial character of the magnetization $M$. Figure 10 shows the phase diagram of the Ising ferromagnetic system representing the magnetization density $m$ as a function of the temperature $T$ and the magnetic field $H.$ There are two ways in which the transition may occur as $H$ varies. For temperatures $T < T_{c},$ the system undergoes a first-order phase transition with magnetization passing abruptly from a positive value to a negative value. In contrast, near the critical point, the system undergoes a continuous or second-order phase transition with magnetization varying continuously. It is said that a second-order phase transition is accompanied by a spontaneous symmetry breaking in which the system chooses to be in either an up or down-spin phase.

When a system is divided in a small number of domains, keeping the thermodynamic variables constant, the macroscopic properties of each domain remain usually the same as the whole. However, as the number of domains becomes important in such a way that the size of each domain becomes less than a certain length $\xi$, the properties of each domain may start to considerably differ over the system. The typical length scale $\xi$ is known as the correlation length, i.e., the length over which the fluctuations of microscopic variables are correlated. At the critical point, the system is divided into locally ordered domains scattered in the bulk, whose mean size corresponds to the correlation length $\xi$, and mean lifetime to the correlation time $\tau.$ The experiments indicate that the correlation length and the correlation time increase more and more as the critical point is approached from the low temperatures side. In particular, it is observed that the correlation length $\xi$ remains finite for a first-order phase transition, but becomes infinite for a second-order one. When $T$ tends to $T_{c}$, the correlation length diverges according to the simple form $\xi \sim (T_{c}-T)^{-\nu }$, where $\nu$ is known as a critical exponent. As a consequence of large fluctuations, which become correlated over all distances near the critical point, the system is in a single phase, since the disordered phase has a stability similar to the ordered phase.

It is a matter of fact that the singular behavior of various thermodynamic functions, in the vicinity of the critical point, is governed by a set of critical exponents for transitions as different as those of liquid–gas and ferromagnetic. In addition to the exponent $\nu,$ the most commonly encountered critical exponents are: ($i$) $\beta$ which describes the order parameter $m$ below $T_{c}$ as $m\sim (T_{c}-T)^{\beta }$, ($ii$) $\gamma$ describing the susceptibility (i.e., the response of the order parameter $m$ to its conjugate field $H)$ $\chi = \left. \frac{\partial m}{ \partial H}\right\vert _{H = 0}\sim \left\vert T-T_{c}\right\vert ^{-\gamma }$, ($iii$) $\alpha$ describing the heat capacity (i.e., the thermal response function) $C = \left. \frac{\partial E}{\partial T}\right\vert \sim \left\vert T-T_{c}\right\vert ^{-\alpha }$, etc. It is remarkable that the critical phenomena of very different physical systems show the same temperature dependence, with roughly similar values of the critical exponents. For the liquid–gas transition $\alpha$ ${\simeq }$ $0.11$, $\beta$ ${\simeq }$ $0.32$, $\gamma$ ${\simeq }$ $1.24$, $\nu$ ${\simeq }$ $0.63$ and for the ferromagnetic transition $\alpha$ ${\simeq }$ $-0.1$, $\beta$ ${\simeq }$ $0.34$, $\gamma$ ${\simeq }$ $1.4$, $\nu$ ${\simeq }$ $0.7$, related by $\alpha +2\beta +\gamma \simeq 2.$ Detailed studies of this universality applied to critical phenomena (see, e.g., ^[99]) go beyond the scope of this review.

The divergence of the correlation length near second-order phase transition indicates that the thermodynamical properties are insensitive to microscopic details of the system. This is a motivation to search for a phenomenological description of critical phenomena capable of describing a wide variety of model systems. In the next section, we introduce such phenomenological approach to second-order phase transitions, known as the Landau-Ginzburg theory, in which the order parameter goes to zero continuously. It should be noted that this theory is quite general in a way that it may be used for a great number of systems with various dimensionality: three-dimensional systems (fluids, ferromagnetic systems...), two-dimensional systems (films, emulsions...), one-dimensional systems (polymers, chains of atoms weakly coupled...).

9.   Phenomenological approach: Landau-Ginzburg theory

The Landau-Ginzburg theory ^[100] has been first developed to investigate the charged superfluids; its way of operation is the minimization of the free-energy functional. If applied to complex fluids near the critical point, the expansion of the free-energy functional can be written as an integral over all space of an appropriate function of the order parameter and its spatial deviations as:

where the order parameter, noted $m({\bf{r}}) = \rho ({\bf{r}})-\rho _{c}$, is the difference between the local density $\rho ({\bf{r}})$ and the critical density $\rho _{c}$. It may be noted that the odd powers of $m({\bf{r}})$ are absent in the integrand (free-energy density) by virtue of the conservation law of particles. Besides, the third and fourth terms are the square gradient and square Laplacian of the order parameter, respectively. Depending on the system under study, other terms may be added or dropped in the expansion. For instance, in homogeneous fluids, all the spatial derivatives of the order parameter are neglected even if $m({\bf{r }})$ is subject to fluctuations at the microscopic level^||. On the contrary, in inhomogeneous fluids the spatial derivatives are expected. As an example, the Ornstein–Zernike (OZ) theory may be rederived ^[101] with Eq 53, when only the coefficients $a$ and $c$ are kept. By stability consideration, both parameters $a$ and $c$ must be positive, however, $a$ vanishes at the critical point.

^|| This local inhomogeniety assumption has been discussed by Landau in his theory of the spatial correlation functions, see Landau LD, Lifshitz EM (1968) Statistical Physics, London: Pergamon.

In other circumstances, as in microemulsions ^[42,102] and homopolymer-diblock copolymer blends ^[103,104], the terms pertinent to a study of the critical behavior are those involving the coefficients $a,$ $c$ and $d$, with a positive value of $a$ in order that the free energy has a minimum. In that case, it has been argued that the stability condition imposes $c^{2}-4ad < 0$, so that $d$ is positive and $c$ may be positive or negative. It should be mentioned that at the change of sign of $c,$ by varying an appropriate parameter, the system undergoes a phase separation into a bulk phase or a spatially modulated phase. The condition $c = 0$, referred to as the Lifshitz point, has been first introduced in the context of ferromagnetic systems ^[105], but it is also known to exist in other systems like liquid crystals ^[106], microemulsions ^[42,102], block copolymers ^[103,104].

In the Landau-Ginzburg theory, it is convenient to write the free-energy functional $F\left[m({\bf{r}})\right]$ in an alternative form by expanding the order parameter and its spatial deviations in Fourier series as:

After substitution of these relations in Eq 53 and integration, all the terms of the form $m_{{\bf{q}}}m_{{\bf{q}}^{\prime }}\exp \left[i({\bf{q+q}}^{\prime }){\bf{\cdot r}})\right]$ vanish when ${\bf{q}} ^{\prime }\neq -{\bf{q}}$ and are equal to $\left| m_{{\bf{q}} }\right|^{2}$ when ${\bf{q}}^{\prime } = -{\bf{q}}$, so that $F\left[m_{ {\bf{q}}}\right]$ reads:

It should be noted that each term of the sum in Eq 57 depends only on single density modes $m_{{\bf{q}}}$ contributing additively to the free energy. Thus, the different $m_{{\bf{q}}}$ are statistically independent, i.e., noninteracting. Because of the truncation of the Fourier series, Eq 57 is only valid for small $q,$ where the wavelength $1/q$ is great compared to the interparticle distance.

10.   Density fluctuations near the critical point

As an application of the Landau-Ginzburg phenomenological approach, we perform a quantitative description of the density fluctuations in self-assembly structures such as microemulsions and diblock copolymers. These systems are homogeneous and disordered over macroscopic length scales, but they are inhomogeneous over mesoscopic length scales ^[38]. The presence of fluctuations in the systems is provided by the emergence of a peak in the static structure factor located between the long-wavelength limit and the principal peak of diffraction. These fluctuations give rise to an instability of the homogeneous phase leading to microphase separation.

It should be pointed out that the free energy (Eq 57) may be written as a function of the density fluctuations under the form ^[101] :

where $S(q)$ is the static structure factor, i.e., the linear response of the system with respect to the particle density. This relation results from the fluctuation-dissipation theorem meaning that the structure factor $S(q)$ of the system is directly proportional to the susceptibility $\chi (q)$ according to the formula $\chi ^{-1}(q) = -(\delta^2F\left[m_{{\bf{q}}} \right] /\delta m_{{\bf{q}}}^{2})$. The comparison of Eq 57 with Eq 58 allows to deduce the following expression for the structure factor at small $q$:

A point of particular interest is the correlation of the density fluctuations in complex systems near the critical point. To gain an insight into the density fluctuations in fluids, the total correlation function $h(r)$ may be investigated by X-rays or light scattering experiments through the structure factor $S(q)$, that is related to the Fourier transform of $h(r)$ as follows:

The microscopic nature of the density fluctuations is also studied by the Ornstein–Zernike relation whose the FT reads (Eq 17):

where $c(q)$ is the FT of direct correlation function $c(r)$. This shows that the correlation functions depend on the set of the arbitrary phenomenological parameters $a$, $c$ and $d$, which are non-universal functions of microscopic interactions. The understanding of the physical origin of long-range character of $h(r)$ is a major problem of complex fluids near the critical point. The slower the correlations attenuation the greater is the characteristic size of fluctuations. In an attempt to explain the asymptotic behavior of $h(r)$ for large $r$, one might calculate the inverse FT of Eq 59.

As already mentioned, the OZ theory may be rederived by taking account only of the terms with the coefficients $a$ and $c$ in Eq 59. In this well-known theory, which predicts the large isothermal compressibility $\chi _{T}$ at the critical point, the structure factor for small values of $q$ is:

Far away from the critical point, the structure factor at $q = 0$ is $S(0) = \left(\beta \rho a\right) ^{-1} = \rho \chi _{T}/\beta$, where $\chi _{T}$ is the isothermal compressibility. But at the critical point, $a$ tends to zero as a consequence of the divergence of the isothermal compressibility. Taking the inverse FT of the Lorentzian representation of Eq 62, presumably good for small wave vectors, provides the total correlation function where the leading term at large distances is:

with the correlation length $\xi = \sqrt{c/a}$. This is in agreement with the fluctuation theorem for the isothermal compressibility:

which predicts that $h(r)$ must vanish asymptotically at large distances more rapidly than $1/r^{3}$. The thermodynamic stability in the OZ theory requires positive values of $a$ and $c$, and at the critical point, where $a\, {\rightarrow }\, 0$, the correlation length tends to infinity, so that $h(r)$ is no longer exponentially damped at the critical point. Hence the asymptotic behavior of $h(r)$ given by Eq 63, which is only correct for moderately large values of the isothermal compressibility, away from the critical point, whereas it behaves as $1/r$ at the critical point.

The Landau-Ginzburg expansion of the free-energy functional (Eq 57) is now used to investigate the structure factors and the total correlation function of the microemulsions ^[42] and diblock copolymers ^[107,108]. In this case, $a$, $c$ and $d$ are arbitrary phenomenological parameters fit accurately on experimental scattering data of these systems up to a small microscopic cutoff compared to the inverse of the interparticle distance.

The structure factor is a crucial property of systems measured by X-ray or neutron scattering and calculated from the interparticle potential in various ways. Remember that the structure factor $S(q)$ describes the arrangement of the particles, in the reciprocal space, for any system. In the case of simple liquids, the principal peak in $S(q)$, at $q = q_{p}$, gives an indication on the mean distance between particles. More relevant for the purpose is the value of $S(q)$, at $q = 0$, that increases drastically when the first-order liquid–vapor transition occurs, due to the increase of the density fluctuations in vapor state. By contrast, in more complex systems, a pre-peak may emerge in the structure factor at the small wave vector $q_{c}$ ($0 < q_{c} < q_{p})$, when the long-range repulsion in the pair potential is sufficiently large, which is the signature of the formation of clusters in the system. In certain circumstances, the pre-peak in $S(q)$ grows, indicating a disordered-pattern of clusters, with the possibility to display spontaneous modulated phases when the pre-peak is diverging. This transition (modulated phases) as observed in $S(q)$, at $q = q_{c}$, has an analogy with the liquid–vapor transition predicted by the divergence of $S(q)$ at $q = 0$. Indeed, both transitions are the consequence of instabilities of the uniform phase with respect to the density fluctuations. In simple liquids with long-range attraction of the van der Waals type, the long-wavelength fluctuations dominate, while in complex liquids, with a long-range slightly repulsive interaction, the density fluctuations of the microscopic wavelength $2\pi /q_{c}$ are predominant, leading to the microphase separation.

Additional information about the density fluctuations in disordered phases may be extracted from the structure factor, judging by whether the discriminant $\Delta$ of the polynomial in Eq 59 changes its sign. With the conditions $\Delta = c^{2}-4ad > 0$, $a > 0$ and $d > 0$, the inverse FT of Eq 59 and straightforward algebra leads to total correlation function ^[23,42]:

with two correlation lengths:

With the conditions $\Delta = c^{2}-4ad < 0$, $a > 0$ and $d > 0$, the inverse FT of Eq 59 reads:

with two correlation lengths:

Depending on the values of the parameters $a$, $c$ and $d$, the total correlation function $h(r)$ may present either a decay resulting from a combination of two Yukawa functions, with the correlation lengths $\xi _{1}$ and $\xi _{2}$, or a damped oscillatory behavior with a periodic spatial variation $\delta$, generating modulated phases in the system.

To obtain more information on how the copolymer in mixtures of homopolymers (or the surfactant in microemulsions) tends to order the mesophases and also to order itself, it is useful to introduce the disorder line and the Lifshitz line. The disorder line is the locus of points in the phase diagram at which the oscillatory behavior appears in the total correlation function $h(r)$. It occurs when ($c^{2}-4ad)$ is negative and $c$ has a specific negative value. It should be noted that the parameter $c$ depends on the copolymer density, while the question arises whether it varies with the temperature ^[109]. The total correlation function is presented in Figure 11 for a system where $c$ is variable and $a = d = 3$. The discriminant ($c^{2}-4ad)\ $is negative when $6 > c > -6.$ It can be seen that the faster attenuation of the total correlation function $h(r)$ is for $c = 10$ whereas its oscillatory behavior begins to dominate for $c < -4$.

The Lifshitz line is the locus of points at which the peak in the structure factor $S(q)$ just starts to move out of zero wave vector towards a nonzero value. It appears when $\partial S^{2}(q)/\partial q^{2}|_{q = 0} = 0$, i.e., $c = 0$. In other words, it corresponds to the point at which the structure factor (Figure 12) ceases to decay monotonously, when $c\! < \!0$, in contrast to the disorder line which indicates the point at which the oscillatory component (Figure 11) appears in the total correlation function $h(r)$. It should be noted that the position of the peak in $S(q)$ moves to large wave vectors when $c$ decreases and divergences for $c\! = \!-6$.

In the context of diblock copolymer blends, Eq 67, with $\Delta < 0$, is needed to describe the disordered phase with the correlation length $\xi$ and the wavelength $\delta$. From an experimental point of view, the oscillation in $h(r)$ reflects the tendency of the copolymer to order the A and B monomers in space. It depends on the concentration of diblock copolymers ^[109,110] and the temperature ^[108]. At low copolymer content, the parameter $c$ is positive and the total correlation function decreases monotonously with the distance. As the copolymer density increases and the temperature remains constant, one expects the wavelength $\delta$ to decrease and to enforce the order over shorter distances. When the amount of copolymer continues to increase, $c$ changes its sign at crossing the Lifshitz line. It remains negative for larger copolymer content up to the occurrence of the disorder line, at the divergence of the peak in $S(q)$ when $\Delta = (c^{2}-4ad) = 0$. The divergence of the subsidiary peak of $S(q)$ at $q = q_{c}$ indicates a strong increase of the fluctuations and an instability of the homogeneous phase leading to microphase separation. The instability of the homogeneous phase is realized at the wavelength value of $\delta = 2\pi /q_{c}$.

11.   Structure factor of complex liquids

In Section 10, the structure factor was determined from a phenomenological approach in terms of the three parameters $a$, $b$ and $c$ to be adjusted on the experiment data. In this section, we determine the structure factor from the microscopic theory of liquids, as a function of the pair potential.

We first consider the calculation of the structure factor $S(q)$, and then investigate the features of the pre-peak on $S(q)$ (magnitude and position). When the MD simulation is used, the accuracy of results strongly depends on the potential model. In contrast, when the liquid theory methods are employed, the results depend also on the approximation used. For instance, with the simplest random phase approximation (RPA), the structure factor reads ^[57]:

where $S^{HS}(q)$ is the hard-sphere structure factor whose the expression is known accurately under the analytical form ^[111,112]:

with

and

The function $\rho c_{lr}(q)$ in Eq 69 is the FT of the direct correlation function $\rho c_{lr}(r)\simeq -\rho \beta u_{lr}(r),$ ($u_{lr}(r)$ being the long-range part of the potential), which is defined as:

If $u(r)$ is taken to be the temperature-independent double Yukawa potential (Eq 51), the FT of the direct correlation function, $\rho c_{lr}(q)$, is expressed directly in terms of the model parameters under the analytical form:

For mathematical convenience, $c_{lr}(q)$ has been calculated analytically with the extension of $u(r)$ inside the core ^[87,93]. However, in the RPA, $c_{lr}(r)$ inside the core is sometimes taken to be $c_{lr}(r < \sigma) = \lambda \left(K_{a}-K_{r}\right)$, where $\lambda$ is a fitting parameter. Different choices for $c_{lr}(r < \sigma)$ are possible by varying the parameter $\lambda$, among them $\lambda = 0$ or $\lambda = 1$ are often used. It should be noted that with an appropriate polynomial for $c_{lr}(r < \sigma)$ inside the core, whose parameters are chosen by determining the free energy, the optimized random phase approximation (ORPA) brings an appreciable improvement to the structure factor. As an example, the ORPA has been proven to be an efficient tool to predict the particularities of the structure factor of the divalent liquid metals ^[113]. It should be mentioned that the structure factor determined by RPA (with $c_{lr}(r)$ rigorously equal to $-\beta u(r)$ at large distances) is not expected be as accurate as the predictions from simulation or from integral equations calculations, because $c_{lr}(r)$ in the hard core is not equal to that of the HS potential, which must be corrected by an additional contribution ^[114]. Recently, we have developed a semianalytical "reverse" approach to link structure and microscopic interactions in two-Yukawa competing fluids ^[115].

Figure 13 displays the structure factor $S(q)$ for the double Yukawa potential with $z_{a} = 1$, $z_{r} = 0.5$ and different values of $K_{a}$, $K_{r}$ and $\rho$. It can be seen that $S(q)$ develops a pre-peak at $q_{c}$ when $K_{a}^{-1}$ is small, indicating a propensity towards cluster formation in the liquid ^[92]. Simultaneously, in the phase diagram $K_{a}^{-1}(\rho)$, we observe a curve, denoted $\lambda$-line by Archer et al. ^[17], that is generated by the locus of points where $S(q_{c})$ diverges. This curve marks the boundary between the homogeneous phase (above) and the modulated phase (below). It may be compared with the binodal, for which $p(\rho _{L)} = p(\rho _{G})$ at a given temperature, and the spinodal, for which the isothermal compressibility diverges ($\chi _{T} = { S(0)/\rho k_{B}T}$). Both the binodal and spinodal lines have the vortex at the critical point $C$. On the other hand, the $\lambda$-line intersects the binodal at two points $A$ and $B,$ indicating that the liquid–vapor transition is preempted by a microphase separation in the range of densities between $A$ and $B$. This phenomenon happens for a particular choice of potential parameters giving rise to the competition between the short-range attraction and long-range repulsion. It should be noted that the location of the $\lambda$-line is very sensitive to the choice of the approximation used to calculate the structure factor.

For a pair potential of the van der Waals type, consisting of a short-range repulsion and a long-range attraction, the theory predicts a first-order transition between vapor and liquid, at the critical temperature. Different phase behaviors may exist when the ranges of the repulsive and attractive interactions are interchanged, giving a small long-range repulsion beyond the attractive well. The possibility that a delicate balance between the short-range attraction and the long-range repulsion may stabilize equilibrium clusters has been discussed in numerous theoretical works ^[116,117]. At low densities, there are clear-cut evidences for the existence of a fluid phase of clusters, whose size and shape are expected to be sensitive to the respective ranges and strengths of the attractive and repulsive contributions to the pair potential. Upon increasing the density, the competition between the repulsive and attractive contributions may stabilize new particle arrangements. In colloidal suspensions, the presence of long-range repulsive interactions strongly alter the formation of clusters giving rise to stable chainlike clusters coalescing into elongated structures, eventually forming a connected network at gelation ^[33,98]. In other systems such as aqueous surfactants and mixtures of block copolymers, competing interactions may lead to the formation of modulated phases ^{[39,41,118,119]}.

It is frequently argued that the emergence of a pre-peak in $S(q)$ indicates the presence of a generic Intermediate Range Order (IRO) within the fluid. An IRO is defined by a locally inhomogeneous region of a system characterized by the appearance of agammaegates (not necessarily at equilibrium) with sizes divided into dimers, trimers and so on, rather than clusters of preferred finite size. A cluster phase at equilibrium refers to a stable situation where the system does not evolve towards either a homogeneous fluid or a completely resolved phase separation, except in the case of a temperature or density change.

Much has been learned about the structural and thermodynamic properties in fluids with competitive SALR interactions, both for disordered cluster fluid states and modulated (ordered) cluster phases, at intermediate densities. But the origin and the formation of cluster fluid phases, at low density, are much less well studied.

Many approaches have been used to examine the formation of cluster fluid phases by various methods: molecular dynamics and Monte-Carlo simulations ^[78,91], self-consistent integral equation approach ^[19], density functional theory ^[17], molecular thermodynamic model ^[120]. Structural criteria have been used ^{[121,122,123]} to identify the boundary between different regimes starting from a homogeneous state fluid up to the IRO. Most of these calculations were carried out with the double Yukawa potential. Besides, interesting calculations ^[96,97] have been performed recently with a potential modeled with a 100-50 Lennard-Jones potential plus a screened electrostatic repulsive interaction of Yukawa's type, for systems at low densities. With this particular pair potential, the results yield promising information on the equilibrium cluster fluid phase.

12.   Concluding remarks

In this article we consider soft materials in which the interactions at intermediate and large distances are dominated by repulsive contributions. So, the long-range attraction favoring the agammaegation is frustrated by the long-range repulsion, which prohibits further growth and changes drastically the phase diagram of purely attractive potentials. Insofar as the long-range interactions between colloids are only attractive, it is well known that soft materials exhibit a conventional liquid–vapor transition; however, in presence of a long-range repulsion, a metastable fluid–fluid transition ^[124,125] (above the liquid–vapor transition) occurs that complicates the description of the thermodynamic phase transitions ^[17,22,23]. Taking account of the variety of soft materials, this article is not intended to be complete. It must be regarded as only an introduction to the study of the interactions in such materials.

In colloidal polymer mixtures, the long-range repulsion (LR) is usually attributed to the weak screened charge carried by the colloids, while the long-range attraction (SA) arises from the depletion forces generated by nonadsorbing polymers. The expression of the repulsive part of the interparticle potential has been first determined by Verwey and Overbeek ^[50] but it has been shown, in Section 4, that it can be derived ^[52] by an analogy between the colloidal suspensions and the metallic liquids. Concerning the expression of the attractive part of the interparticle potential ^[61], presented in Section 5, it is ascribed to the sum of van der Waals forces between atoms of the colloidal particles. Notwithstanding its simplicity and its relevance, the interparticle (SALR) potential is often used under an empirical form (Section 6), able to generate long-range attractive and long-range repulsive contributions, such as the double Kac potential and the double Yukawa potential. These one-component models for interactions between solute particles, in which the solvent degrees of freedom are integrated implicitly, involve only spherically symmetric pair interactions and provide the impetus to better understand their origin. The role of attractive and repulsive contributions to the interparticle potential has been carefully investigated in globular protein solutions ^[14] and star-polymer systems ^[126,127], specially to describe the spontaneous microphase formation in systems of spherical symmetric particles, where stable clusters of a specific agammaegation number develop in the solution of monomers ^[91,120]. Clustering in these systems is important not only for theoretical studies but also for developing new materials such as nanoparticles ^[128,129], protein solutions ^[14] and other in the process of biomineralisation ^[130] by which living organisms produce minerals.

The phenomenon of microphase separation observed in microemulsions and diblock copolymers, also described by SARL potentials, is common to other physical systems exhibiting diverse microstructures at two and three dimensions, known as modulated phases. Most of the observed phenomena and available experimental data (lamellar structures, diverse arrays of cylinders and spheres...) have stimulated a strong theoretical effort ^[5,98,117]. The microscopic inhomogeneities in microemulsions or diblock copolymers are related to the chemical constituents, the range of volume fractions ^[110] and the temperature ^[108]. From phenomenological (Section 10) and microscopic (Section 11) approaches, it is now understood that frustration, arising from competing interactions over different distances, plays a crucial role in microphase separation. Another point we want to address is an investigation of the features of the total correlation function $h(r)$ and the structure factor $S(q)$ of the systems under study, for which the Ornstein–Zernicke picture is of limited use. In the context of diblock copolymer blends, the disordered phase is characterized by the correlation length $\xi$ and the wavelength $\delta$ that give rise to the oscillations in the total correlation function $h(r)$. As the volume fraction increases and the temperature remains constant, one expects the wavelength $\delta$ to decrease and to enforce the order over shorter distances. For larger copolymer content, the peak in the structure factor $S(q)$ at $q = q_{c}$ diverges. This divergence associated to the correlations of density fluctuations means an instability of the homogeneous phase leading to microphase separation.

The present review provides ample evidence of the confidence of potential models with competing attractive and repulsive interactions to describe a variety of structures in complex fluids of components larger than molecules. In conjunction with appropriate potential models, real systems such as colloidal dispersions, microemulsions, diblock copolymer blends and many others formed by self-assembly provide interesting tools and potential exciting opportunities for studying intriguing phenomena in generating nanostructures.

Acknowledgments

I would like to thank Jean-Marc Bomont for valuable discussions on this work and Andrei Postnikov for its careful reading of the manuscript.

Conflict of interest

No potential conflict of interest was reported by the author.