Statistical theory of fluids with a complex electric structure: Application to solutions of soft-core dipolar particles
SStatistical theory of fluids with a complex electric structure: Application tosolutions of soft-core dipolar particles
Yury A. Budkov
1, 2, a)1)
Tikhonov Moscow Institute of Electronics and Mathematics,National Research University Higher School of Economics, Tallinskaya st. 34,123458 Moscow, Russia G.A. Krestov Institute of Solution Chemistry of the Russian Academy of Sciences,Akademicheskaya st. 1, 153045 Ivanovo, Russia
Based on the thermodynamic perturbation theory (TPT) and the Random phaseapproximation (RPA), we present a statistical theory of solutions of electrically neu-tral soft molecules, every of which is modelled as a set of sites that interact witheach other through the potentials, presented as the sum of the Coulomb potentialand arbitrary soft-core potential. As an application of our formalism, we formu-late a general statistical theory of solution of the soft-core dipolar particles. Forthe latter, we obtain a new analytical relation for the screening function. As a spe-cial case, we apply this theory to describing the phase behavior of a solution of thedipolar particles interacting with each other in addition to the electrostatic potentialthrough the repulsive Gaussian potential – Gaussian core dipolar model (GCDM).Using the obtained analytic expression for the total free energy of the GCDM, weobtain the liquid-liquid phase separation with an upper critical point. The developedformalism could be used as a general framework for the coarse-grained description ofthermodynamic properties of solutions of macromolecules, such as proteins, betaines,polypeptides, etc. a) [email protected] a r X i v : . [ c ond - m a t . s o f t ] M a r . INTRODUCTION As is well known, in contrast to simple fluids whose molecules interact with each otherthrough the potentials, always containing a short-range repulsive component (hard-core) ,the intermolecular interaction in complex macromolecular fluids can be described by soft-corepotentials . The simplest soft-core potential, which can describe the effective interactionbetween the centers of mass of polymer coils in dilute and semi-dilute polymer solutions isthe Gaussian core potential . So far, such a Gaussian core model (GCM) of fluid with theGaussian intermolecular potential of interaction has been extensively studied for the bulkphase, as well as at the interfaces . The most surprising result obtained by Louis et al. is that thermodynamic properties of the GCM at sufficiently small energetic parameter ofthe repulsive interaction in a liquid state region can be described with good accuracy in awide range of state parameters within the Random phase approximation (RPA). Compar-ing the analytical calculation with the Monte-Carlo simulation results, the authors showedthe applicability of the RPA for the whole liquid state region of the system phase diagram.Thus, one can conclude that the availability of an accurate analytical equation of state forGCM allows us to use this model as a reference system for describing within the thermo-dynamic perturbation theory (TPT) more complex soft-core macromolecular fluids whosemolecules carry additional microscopic quantities, such as multipole moments, polarizability,quadrupolarizability, etc .One can especially stress, among others, a number of macromolecules (such as proteins,polypeptides and betaines) that have some important applications and under certain con-ditions acquire a large dipole moment ( (cid:39) − D ) . However, in order to describetheoretically the electrostatic interactions between such big dipolar macromolecules, it isincorrect to model them as the point-like dipoles (or even dipolar hard spheres) similar topolar molecules in the theory of simple dipolar fluids due to a very big dipole length( ∼ − nm ). Instead, it is necessary to consider a dipolar molecule as a set of chargedcenters (sites), separated by a fixed or fluctuating distance . So far, only several analyt-ical models of dipolar fluids taking into account the internal electric structure of dipolarmolecules and only three of them discuss possible applications to thermodynami-cal description of soft matter.In paper , the authors proposed a statistical theory of conformational behavior of a2ipolar polymer chain. In contrast to the previous theoretical models of solutions of dipolarpolymer chains , modeling the electrostatic interaction of monomers as an interactionof point-like dipoles, the authors considered the monomeric units of dipolar chain as hardcharged dumbbells . Accounting properly for the conformational entropy, the excluded vol-ume interactions of hard dumbbells and their electrostatic interactions, the authors achievedgood agreement between the theoretical and MD simulation results for the conformationalbehavior of dipolar chain. The authors pointed out that their model might be used forqualitative describing the conformational behavior of polyzwitterions, such as betaines inpolar solvents, folding of the proteins with the polar residues, weak polyelectrolytes andpolyelectrolyte-based gels. In paper , the authors formulated a field-theoretic methodologyfor describing the thermodynamic behavior of polarizable soft matter. The model is basedon fluid elements, referred to as beads, which can carry a net monopole charge at their cen-ter of mass and a permanent or induced dipole, described through a Drude-type distributedcharge approach . Beads of different types can be mixed or linked into polymer chains withan arbitrary mechanism of chain flexibility . The authors claim that their general fieldtheoretic model can be relevant for the description of a vast range of soft matter systems,such as polyelectrolytes, ionic liquids, polymerized ionic liquids, ionomers, block copolymers,among others. A nonlocal field-theoretic model has been recently proposed in our paper ,which is a generalization of formerly formulated a field-theoretical model of dipolar fluidwith a fixed distance between the charged groups of dipolar particles. We have developed amodel of dipolar particles as dimers of oppositely charged point-like polar groups, separatedby a fluctuating distance, distributed in accordance with a certain probability distributionfunction g ( r ) . Using the proposed nonlocal field-theoretical model, we have derived a non-linear integro-differential equation with respect to the mean-field electrostatic potential,generalizing the Poisson-Boltzmann-Langevin equation for the point-like dipoles obtainedin papers and generalized for the case of a polarizable solvent in . We have used theobtained equation in its linearized form to derive expressions for the mean-field electrostaticpotential of the point-like test ion and its solvation free energy in a salt-free solution, as wellas in a solution in the presence of salt ions. As the main result of that paper, we obtained ageneral expression for the bulk electrostatic free energy of the ion-dipole mixture within theRandom phase approximation (RPA). In addition, we pointed out that the developed theorycould be applied to thermodynamic description of a solution of dipolar biomacromolecules3nd showed that nonlocal effects, related to the nonzero length of dipolar particles mustbe important for the aqueous solutions of proteins with the dipole length ≈ nm (dipolemoment ≈ D ) that is quite accessible for the real protein aqueous solutions (see, forinstance, ). We would like to note, however, that formulated theory could be applied todescription of thermodynamic behavior of solutions of the dipolar macromolecules only atrather small concentrations, since it does not take into account the excluded volume interac-tions of the dipolar particles. In order to describe the behavior of thermodynamic functionsof macromolecular solutions in a wide range of concentration, it is necessary to take intoaccount not only the electrostatic interactions, but also the soft-core repulsive interactions.In papers the authors formulated a field-theoretical model of electrolyte solution whoseions interact with each other via the Coulomb and soft-core (Yukawa) repulsive potentials.The authors applied this theory to description of ionic liquid in the bulk, as well as at theneutral dielectric interface with account for the finite size of ions (modelled by the soft-core repulsive potential) within the variational method. However, theoretical description ofthe thermodynamic properties of soft matter systems with more complex internal electricstructure than that of the ions (dipolar and quadrupolar particles, for instance) requiresformulating a model of soft-core molecules that are modelled as a set of spatially correlatedsites interacting with each other through the Coulomb and soft-core potentials. In presentpaper, we will formulate such a field theory, based on TPT and RPA for description of ther-modynamic properties of solutions of soft-core molecules with a complex electric structure.As an illustration, we will show how the developed formalism can be applied to descriptionof the phase behavior of the solution of soft-core dipolar particles.The rest of the article is organized as follows. In sec. II, we propose a general field-theoretic formalism based on TPT and RPA for a theoretical description of thermodynamicproperties of solutions of soft-core particles with an arbitrary internal electric structure. InSec. III, we apply the developed formalism to the formulation of the general theory of soft-core dipolar particles, generalizing the previously formulated theory . Sec. IV considers thesimplified Gaussian-core dipolar model (GCDM), which is a particular case of the generalmodel, discussed in Sec. III. Sec. V presents some numerical results regarding the phasebehavior of solution of Gaussian core dipolar particles and Sec. VI contains some concludingremarks and further prospects. 4 I. GENERAL FORMALISM
In this section, based on the thermodynamic pertubation theory (TPT) and the Randomphase approximation (RPA), we will formulate a general field theoretical approach to thedescription of solutions of electrically neutral soft particles with a complex internal electricstructure.
A. Derivation of excess free energy of fluids with complex electric structure
Let us consider a solution of N electrically neutral molecules which we will model asa set of l charged sites with charges q α ( α = 1 , .., l ), so that l (cid:80) α =1 q α = 0 . Let us assumethat in each molecule, the positions r ( α ) of sites are random variables, described by theprobability distribution function g ( r (1) , .., r ( l ) ) , which in general case is determined by thequantum nature of the molecule and have to satisfy the following normalization condition (cid:90) .. (cid:90) d r ( α ) ..d r ( α l − ) g ( r (1) , .., r ( l ) ) = 1 , (1)where α i (cid:54) = α k for all i (cid:54) = k . As in our previous paper , we consider the case of animplicit solvent, modelling it as a continuous dielectric medium with the constant dielectricpermittivity ε . Let us also assume that apart from the Coulomb interactions between sitesthere are also some "volume" non-electrostatic interactions, described by the l × l matrix ofsoft-core potentials V αγ ( r ) .Therefore, bearing in mind all the above model assumptions, we can write the configura-tion integral of the solution as follows Q = (cid:90) d Γ exp [ − βH − βH el ] , (2)where the integration measure is (cid:90) d Γ( · ) = (cid:90) .. (cid:90) N (cid:89) i =1 d r (1) i ..d r ( l ) i V g ( r (1) i , .., r ( l ) i )( · ) , (3) β = ( k B T ) − is the inverse thermal energy. The Hamiltonian of non-electrostatic interactionsis H = 12 (cid:90) d r (cid:90) d r (cid:48) (cid:88) α,γ ˆ n α ( r ) V αγ ( r − r (cid:48) )ˆ n γ ( r (cid:48) ) = 12 (ˆ nV ˆ n ) , (4)5here ˆ n α ( r ) = N (cid:80) i =1 δ (cid:16) r − r ( α ) i (cid:17) is the microscopic concentration of the sites of α th type and V αγ ( r − r (cid:48) ) is the matrix of the non-electrostatic interaction potentials between the sites ofmolecules. The Hamiltonian of the electrostatic interactions is H el = 12 (cid:90) d r (cid:90) d r (cid:48) ˆ ρ ( r ) G ( r − r (cid:48) ) ˆ ρ ( r (cid:48) ) = 12 ( ˆ ρG ˆ ρ ) , (5)where ˆ ρ ( r ) = l (cid:80) α =1 q α ˆ n α ( r ) is the microscopic charge density of molecules satisfying the elec-troneutrality condition, i.e. (cid:82) d r ˆ ρ ( r ) = 0 ; G ( r − r (cid:48) ) = 1 / (4 πεε | r − r (cid:48) | ) is the Green functionof the Poisson equation.Further, using the standard Hubbard-Stratonovich transformation exp (cid:20) − β ρG ˆ ρ ) (cid:21) = (cid:90) D ϕC exp (cid:20) − β ϕG − ϕ ) + iβ ( ˆ ρϕ ) (cid:21) , (6)we arrive at the following functional representation of the configuration integral Q = (cid:90) D ϕC exp (cid:20) − β ϕG − ϕ ) (cid:21) (cid:90) d Γ exp [ − βH ] exp [ iβ ( ˆ ρϕ )] , (7)where the following short-hand notations ( ϕG − ϕ ) = (cid:90) d r (cid:90) d r (cid:48) ϕ ( r ) G − ( r , r (cid:48) ) ϕ ( r (cid:48) ) , ( f ϕ ) = (cid:90) d r f ( r ) ϕ ( r ) and the normalization constant of the Gaussian measure C = (cid:90) D ϕ exp (cid:20) − β ϕG − ϕ ) (cid:21) have been introduced. Considering a solution of molecules without Coulomb interactions asa reference system, one can write configuration integral (7) in the standard TPT form: Q = Q (cid:90) D ϕC exp (cid:20) − β ϕG − ϕ ) (cid:21) (cid:104) exp [ iβ ( ˆ ρϕ )] (cid:105) , (8)where Q = (cid:90) d Γ exp [ − βH ] (9)is the configuration integral of the reference system, which will be calculated below and thesymbol (cid:104) ( · ) (cid:105) = 1 Q (cid:90) d Γ exp [ − βH ] ( · ) (10)6enotes the average over the statistics of the reference system. Using the Random phaseapproximation (RPA), i.e. truncating the cumulant expansion in the integrand exponentby the second term, we arrive at Q ≈ Q (cid:90) D ϕC exp (cid:20) − β ϕG − ϕ ) − β ϕS c ϕ ) (cid:21) , (11)where S c ( r − r (cid:48) ) = (cid:104) ˆ ρ ( r ) ˆ ρ ( r (cid:48) ) (cid:105) = (cid:88) α,γ q α q γ S αγ ( r − r (cid:48) ) (12)is the correlator of microscopic charge density with the structure matrix S αγ ( r − r (cid:48) ) of thereference system, whose elements can be calculated as follows S αγ ( r − r (cid:48) ) = (cid:104) δ ˆ n α ( r ) δ ˆ n γ ( r (cid:48) ) (cid:105) , (13)where the local fluctuations of microscopic concentrations of the sites δ ˆ n α ( r ) = ˆ n α ( r ) −(cid:104) ˆ n α ( r ) (cid:105) are introduced.Therefore in the RPA, the configuration integral takes the form Q ≈ Q (cid:90) D ϕC exp (cid:20) − β ϕG − ϕ ) (cid:21) = Q exp (cid:20) tr (ln G − ln G ) (cid:21) , (14)where the symbol tr ( A ) denotes the trace of integral operator, so that for the kernel A ( r , r (cid:48) ) we have the following relation tr ( A ) = (cid:90) d r A ( r , r ) (15)and G − ( r − r (cid:48) ) = G − ( r − r (cid:48) ) + βS c ( r − r (cid:48) ) (16)is the renormalized reciprocal Green function of the solution medium in RPA. We would liketo note that trace of operator A with a translation invariant kernel A ( r , r (cid:48) ) = A ( r − r (cid:48) ) canbe calculated as tr ( A ) = V (cid:90) d k (2 π ) ˜ A ( k ) , (17)where V is the total volume of system and ˜ A ( k ) = (cid:90) d r A ( r ) e − i kr (18)is the Fourier-image of function A ( r ) . 7herefore, taking into account the expression (16) and definitions of trace and subtractingthe electrostatic self-energy of molecules from the final expression, we arrive at the totalexcess free energy F ex = F + V k B T (cid:90) d k (2 π ) (cid:18) ln (cid:18) κ ( k ) k (cid:19) − κ ( k ) k (cid:19) , (19)where κ ( k ) = 1 k B T εε (cid:88) α,γ q α q γ S αγ ( k ) (20)is the screening function and S αγ ( k ) are the structure factors of the reference system.It should be noted that equation (20) has been obtained first in work within the differentapproach. B. Calculation of free energy and structure factors of the reference systemwithin RPA
In order to calculate the excess free energy (19) of solution, it is necessary to calculatethe excess free energy F and the structure factors S αγ ( k ) of the reference system. For thispurpose, we determine the following auxiliary configuration integral of the reference system Q [Ψ] = (cid:90) d Γ exp (cid:20) − β n V ˆ n ) − β ( δ ˆ n Ψ) (cid:21) , (21)where ( δ ˆ n Ψ) = (cid:82) d r (cid:80) α δ ˆ n α ( r )Ψ α ( r ) with the auxiliary external potentials Ψ α ( r ) . The struc-ture factors can be calculated as S αγ ( r − r (cid:48) ) = 1 β δ δ Ψ α ( r ) δ Ψ γ ( r (cid:48) ) Q [Ψ] Q [0] (cid:12)(cid:12)(cid:12)(cid:12) Ψ=0 . (22)Further, using the Hubbard-Stratonovich transformation e − β (ˆ n V ˆ n ) = (cid:90) D Φ C V e − β ( Φ V − Φ ) + iβ (ˆ n Φ) , (23)where (cid:0) Φ V − Φ (cid:1) = (cid:90) d r (cid:90) d r (cid:48) (cid:88) α,γ Φ α ( r ) V − αγ ( r − r (cid:48) )Φ γ ( r (cid:48) ) , (ˆ n Φ) = (cid:90) d r (cid:88) α ˆ n α ( r )Φ α ( r ) , and C V = (cid:82) D Φ e − β ( Φ V − Φ ) is the renormalization constant of the Gaussian measure, weobtain the following representation of (21): Q [Ψ] = (cid:90) D Φ C V e − β ( Φ V − Φ ) (cid:90) d Γ e iβ (ˆ n Φ) − β ( δ ˆ n Ψ) , (24)8here the reciprocal matrix operator V − is determined by (cid:90) d r (cid:48)(cid:48) (cid:88) λ V αλ ( r − r (cid:48)(cid:48) ) V − λγ ( r (cid:48)(cid:48) − r (cid:48) ) = δ αγ δ ( r − r (cid:48) ) . (25)Further, rewriting the fluctuating field variable as Φ( r ) = Φ + Λ( r ) , where Φ = V (cid:82) d r Φ( r ) , and performing the integration over the constant field Φ in (24), we obtain Q [Ψ] = Q ( MF )0 e β ( n Ψ) (cid:90) D Λ C V e − β ( Λ V − Λ ) (cid:90) d Γ e iβ (ˆ n [Λ+ i Ψ]) , (26)where Q ( MF )0 = exp (cid:20) − β n V n ) (cid:21) (27)is the mean-field approximation for the configuration integral of the reference system withthe average concentrations n α ( r ) = (cid:104) ˆ n α ( r ) (cid:105) of the charged sites. Now, performing thefollowing shift Λ → Λ − i Ψ , we obtain Q [Ψ] = Q ( MF )0 e β ( n Ψ)+ β ( Ψ V − Ψ ) (cid:90) D Λ C V e − β ( Λ V − Λ ) + iβ ( Λ V − Ψ ) Q N [Λ] , (28)where Q [Λ] = (cid:90) .. (cid:90) d r (1) ..d r ( l ) V g ( r (1) , .., r ( l ) ) exp (cid:34) iβ l (cid:88) α =1 Λ α ( r ( α ) ) (cid:35) (29)is the one-particle partition function. In the thermodynamic limit N → ∞ , V → ∞ , N/V → n , we obtain Q [Ψ] = Q ( MF )0 e β ( n Ψ)+ β ( Ψ V − Ψ ) (cid:90) D Λ C V e − β ( Λ V − Λ ) + iβ ( Λ V − Ψ ) + W [Λ] , (30)where we have introduced the following functional W [Λ] = n (cid:90) .. (cid:90) d r (1) ..d r ( l ) g ( r (1) , .., r ( l ) ) (cid:32) exp (cid:34) iβ l (cid:88) α =1 Λ α ( r ( α ) ) (cid:35) − (cid:33) . (31)Taking into account that (cid:82) d r Λ( r ) = 0 , expanding the functional W [Λ] in the power-lawseries and truncating it by the quadratic term, we obtain Q [Ψ] ≈ Q ( MF )0 e β ( n Ψ)+ β ( Ψ V − Ψ ) (cid:90) D Λ C V e − β ( Λ D − Λ ) + iβ ( Λ V − Ψ )= Q ( RP A )0 e β ( n Ψ)+ β ( Ψ V − Ψ ) − β ( Ψ V − DV − Ψ ) (32)where Q ( RP A )0 = Q ( MF )0 exp (cid:20) T r (ln D − ln V ) (cid:21) (33)9s the reference system configuration integral in the RPA and operator D is determined bythe relation D − αγ ( r − r (cid:48) ) = V − αγ ( r − r (cid:48) ) + βW αγ ( r − r (cid:48) ) , (34)where the operator W has the following form W αγ ( r − r (cid:48) ) = n ( δ αγ δ ( r − r (cid:48) ) + (1 − δ αγ ) g αγ ( r − r (cid:48) )) (35)with the structure functions g αγ ( r − r (cid:48) ) = (cid:90) .. (cid:90) d r (1) ..d r ( α − d r ( α +1) ..d r ( γ − d r ( γ +1) ..d r ( l ) ×× g ( r (1) , .., r ( α − , r , r ( α +1) , .., r ( γ − , r (cid:48) , r ( γ +1) , .., r ( l ) ) (36)which describe nothing more than the probability distribution functions of the distancebetween α th and γ th sites. The symbol T r ( .. ) denotes the trace of integral matrix operatorin accordance with the following definition T r ( A ) = (cid:90) d r (cid:88) α A αα ( r , r ) . (37)Using relation (22), we obtain the following result for the structure factor S = 1 β (cid:0) V − − V − DV − (cid:1) , (38)which after some simple transformations can be rewritten as follows S − αγ ( r − r (cid:48) ) = W − αγ ( r − r (cid:48) ) + βV αγ ( r − r (cid:48) ) (39)or in the Fourier-representation S − αγ ( k ) = W − αγ ( k ) + β ˜ V αγ ( k ) . (40)Therefore, using the defintion (37) of the trace of integral matrix operator and takinginto account that det ˜ D ( k )det ˜ V ( k ) = det S ( k )det W ( k ) , we obtain the free energy of the reference system in the RPA as follows F = − k B T ln Q ( RP A )0 = V n (cid:88) α,γ ˜ V αγ (0) − V k B T (cid:90) d k (2 π ) ln det S ( k )det W ( k ) . (41)The first term in the right hand side of eq. (41) determines the mean-field approximationfor the contribution of excluded volume interactions to the total free energy, whereas thesecond term – contribution of the Gaussian fluctuations of the potential energy near itsaverage value. 10 II. GENERAL THEORY OF SOFT DIPOLAR PARTICLES
Now, we turn to the application of our general theory to describing dipolar particles. Weassume that each dipolar particle is a dimer of the oppositely charged sites with charges ± q ,separated by the fluctuating distance. As in our earlier paper , we introduce the probabilitydistribution function g ( r ) of the distance between the charged sites. The relation for thestructure factor of the dipolar particles solution in the RPA takes the following form ˆ S − ( k ) = ˆ W − ( k ) + β ˜ V ( k ) = n (1 − g ( k )) + βV ( k ) − g ( k ) n (1 − g ( k )) + βV ( k ) − g ( k ) n (1 − g ( k )) + βV ( k ) n (1 − g ( k )) + βV ( k ) , (42)where ˜ V ( k ) = ˜ V ( k ) ˜ V ( k )˜ V ( k ) ˜ V ( k ) (43)is the matrix of Fourier-images of non-electrostatic interaction potentials between sites ofdipolar particles, whereas the structure matrix has the following form W ( k ) = n ng ( k ) ng ( k ) n , (44)where g ( k ) = (cid:90) d r g ( r ) e − i kr (45)is the characteristic function corresponding to the probability distribution function g ( r ) .Further, using relations (20) and (44), we obtain the following new analytical expressionfor the screening function for the solution of soft-core dipolar particles κ ( k ) = κ D (1 − g ( k )) (cid:16) βn (cid:16) ˜ V ( k ) + ˜ V ( k ) + 2 ˜ V ( k ) (cid:17) (1 + g ( k )) (cid:17) βn (cid:16) ˜ V ( k ) + ˜ V ( k ) + 2 ˜ V ( k ) g ( k ) (cid:17) + β n (1 − g ( k )) (cid:16) ˜ V ( k ) ˜ V ( k ) − ˜ V ( k ) (cid:17) , (46)where κ D = r − D = 2 q n/ ( εε k B T ) is the square of the inverse Debye radius r D , attributedto the charged sites of dipolar particles. Note that in the absence of the excluded volumeinteractions ( ˜ V αγ ( k ) = 0 ) between the sites, the expression (46) will transform into that hasbeen obtained earlier in . 11 V. THEORY OF GAUSSIAN CORE DIPOLAR PARTICLES
In order to obtain the analytical results for the excess free energy of a dipolar particlessolution, we will simplify the general model, considering only the case of Gaussian coredipolar model (GCDM). Namely, let us consider a reference system with ˜ V ( k ) = ˜ V ( k ) > and ˜ V ( k ) = ˜ V ( k ) = 0 . Essentially, such a reference system describes a set of soft particleswith a point-like ” fantom ” particle, grafted to their center of mass. It should be noted thatthis theoretical model could be applied for coarse-grained description of the macromoleculesof dipolar polypeptides or betaines, dissolved in water. Moreover, such theoretical modelcould be applied to coarse-grained description of zwitterionic liquids whose molecules areusually comprised of the covalently bounded big organic cation and small anion .Thus, in that case we obtain the following relation for the screening function: κ ( k ) = κ D (1 − g ( k )) Q ( k ) , (47)where Q ( k ) = 1 − βn ˜ V ( k )1 + βn ˜ V ( k ) (1 − g ( k )) (48)is the auxiliary function. The electrostatic free energy can be written as follows F el = V k B T (cid:90) d k (2 π ) (cid:18) ln (cid:18) κ D k (1 − g ( k )) Q ( k ) (cid:19) − κ D k (1 − g ( k )) Q ( k ) (cid:19) . (49)As in our previous paper , to describe a distribution function g ( r ) of the distance betweencharged groups of dipolar particles, we use the Yukawa-type distribution function g ( r ) = 32 πl r exp (cid:104) −√ r/l (cid:105) , (50)which determines the following characteristic function g ( k ) = 11 + k l . (51)Using the expression (51), following ’mean-field’ approximation Q ( k ) ≈ − βn ˜ V (0)1 + βn ˜ V (0) (1 − g ( k )) = 1 − α − g ( k )) , (52)and taking integral (49), we arrive at the new analytical expression for the electrostatic freeenergy of solution of the dipolar particles F el = − V k B Tl θ ( y, α ) , (53)12here θ ( y, α ) = (cid:18) − α (cid:19) σ ( y ) − √ (cid:0) ( α + 4)( y + 2 y + y √ y ) + 8(1 + √ y ) (cid:1) π (cid:0) √ y (cid:1) ++ 3 √ (cid:104) αy √ y ) (cid:105) (cid:0) y (cid:0) √ y ( α + 2) (cid:1) + 2 y + 4(1 + √ y ) (cid:1) π (1 + √ y ) . (54)Here, α = βn ˜ V (0) / (1 + βn ˜ V (0)) < , y = l / r D is the strength of dipole-dipole interactionsand the auxiliary function σ ( y ) = √ π (cid:0) y ) / − − y (cid:1) (55)has been introduced.As in the case of point-like dipolar particles , one can obtain two limiting regimes forthe electrostatic free energy F el V k B T = − √ e ln π ( εε k B T ) (cid:16) − α + α (cid:17) , y (cid:28) − πr D (cid:0) − α − (cid:0) − e α/ (cid:1)(cid:1) , y (cid:29) . (56)As in the theory without the repulsive interactions between dipolar particles, the firstregime ( r D (cid:29) l ) corresponds to the situation, when the electrostatic interactions of chargedgroups manifest themself as the Keesom interactions of dipoles, separated by rather bigdistances; the second regime ( r D (cid:28) l ) corresponds to the case when the charged groupsof the dipolar particles do not "feel" connectivity anymore and, thus, behave as the freelymoving ions in regular electrolyte solution, participating in the screening of electrostaticinteractions. As is seen from relations (56), in both limiting regimes, the occurrence ofan additional repulsive interaction between dipolar particles results in a decrease in theabsolute value of the electrostatic free energy with respect to the case of dipolar particleswithout excluded volume interactions . Such a behavior can be explained by the fact thatoccurrence of the finite volume of the dipolar particles prevents their mutual approachingto the small distances that leads to decrease in the electrostatic free energy in its absolutevalue.Using the general expression (41), we obtain the following simple relation for the excessfree energy of the reference system, which has the form of the free energy of a one-component13uid with a soft central potential of the intermolecular interaction with a positive Fourier-image F = V n ˜ V (0)2 + V k B T (cid:90) d k (2 π ) (cid:16) ln (cid:16) βn ˜ V ( k ) (cid:17) − βn ˜ V ( k ) (cid:17) , (57)from which we have subtracted the self-energy of molecules. Choosing the GCM as a refer-ence system, i.e. using the Gaussian potential for modeling the non-electrostatic interactionsof dipolar particles V ( r ) = V exp (cid:2) − r /R (cid:3) (58)and taking into account that ˜ V ( k ) = π / R V exp [ − k R / , we obtain the following ex-pression for the reference system free energy F N k B T = ξ − βV F ( ξ )) , (59)where ξ = π / βnV R , R is the effective size of the Gaussian core and V is the energeticparameter, showing the strength of the repulsive non-electrostatic interactions; the auxiliaryfunction F ( ξ ) = 43 √ π (cid:90) dt ( − ln t ) / t ξt (60)is introduced, which can be expressed via the polylogarithm function. We would like tonote that the free energy (59) describes the thermodynamic properties of the GCM in theliquid-state region with good accuracy at sufficiently small interaction parameter V . V. NUMERICAL RESULTS AND DISCUSSIONS
Before turning to numerical calculations, we introduce the following set of dimensionlessparameters: ˜ n = nR , ˜ T = 4 πεε Rk B T /q , ˜ V = 4 πεε RV /q , ˜ l = l/R . Due to the fact thatfor the real dipolar macromolecules (such as betaines or polypeptides) the effective dipolelength l is order of the effective size of core R , we assume in further that ˜ l = 1 . In whatfollows we will consider only sufficiently small values of energetic parameter V of Gaussianpotential to satisfy the applicability condition of the RPA for the GCM .Fig. 1 demonstrates dependences of the chemical potential, expressed in thermal energy k B T units of the GCDM on the reduced concentration ˜ n at different reduced temperatures ˜ T . As one can see, at a sufficiently high temperature, the chemical potential behaves in a14 igure 1. Dependences of chemical potential µ = k B T ln( nR ) + ∂f ex /∂n ( f ex is the excess freeenergy density) on the reduced concentration ˜ n at different reduced temperatures: ˜ T = 0 . (solidline), ˜ T = 0 . (dotted line) and ˜ T = 0 . (dashed line). The data are shown for ˜ V = 1 and ˜ l = 1 . standard way – monotonically increases with the concentration. However, when the tem-perature drops below some threshold value, a van der Waals loop occurs on the chemicalpotential, indicating the liquid-liquid phase separation. Therefore, the present theory pre-dicts that sufficiently strong electrostatic interactions between GCDM in a solution mediumcan provoke a liquid-liquid phase separation. Fig.2 shows a typical phase equilibrium curve(binodal) for solution of the GCDM with an upper critical point. As is seen, a stronglyasymmetric binodal of the GCDM is very similar to the binodal of the restricted primitivemodel (RPM) of electrolyte solution (see, for instance, ). We would also like to note thatthe critical parameters for the GCDM have the same order of value that are the criticalparameters for the RPM. We would like to mention that system of Gaussian coredipolar particles with sufficiently large dipolar length is quite similar to systemsof soft charged dumbbells which phase behavior also well reproduces the crit-ical behavior of the RPM . The similar critical behavior is realized in thesystems of hard charged dumbbells . It is worth noting that quite similar coexistence15urves take place in the salt-free polyelectrolyte solutions (see, for instance, ). It can besuggested that such broad and asymmetric binodals should be realized in the systems wherethe phase behavior is determined by the long-range Coulomb interactions.
Figure 2. Binodal of the liquid-liquid phase equilibria of a solution of the Gaussian core dipolarparticles. The data are shown for ˜ V = 1 and ˜ l = 1 . The critical parameters for this case are ˜ T c = 0 . and ˜ n c = 0 . . We would also like to discuss how the soft core repulsive interactions between the dipolarparticles affect the electrostatic free energy of the GCDM. Fig. 3 demonstrates dependencesof the reduced electrostatic free energy density ˜ f el = βf el R on the concentration ˜ n atdifferent parameters ˜ V . As is seen, an increase in parameter ˜ V leads to a substantialdecrease in the absolute value of the electrostatic free energy at rather big concentrations. Asit was already pointed out in sec. III, such a behavior of the electrostatic free energy can beeasily interpreted. Indeed, making repulsive interaction parameter V bigger, we increase theeffective excluded volume of dipolar particles, that, in turn, prevents their interpenetrationand, consequently, decreases the contribution of the electrostatic interactions to the totalfree energy. 16 igure 3. Reduced electrostatic free energy density ˜ f el = βf el R as a function of the reducedconcentration ˜ n at the fixed temperature ˜ T = 1 and different Gaussian interaction parameters: ˜ V = 0 (solid line), ˜ V = 0 . (dotted line), and ˜ V = 1 (dashed line). The data are shown for ˜ l = 1 . VI. CONCLUDING REMARKS AND PROSPECTS
We have formulated a field theoretical approach, based on the thermodynamic perturba-tion theory and random phase approximation for theoretical description of thermodynamicproperties of solutions of electrically neutral soft-core molecules with an arbitrary internalelectric structure, modelled as a set of charged sites, interacting with each other through theCoulomb and soft-core potentials. We have obtained the general relation for the excess freeenergy, taking into account internal structure of molecules and electrostatic and soft-coreinteractions between them. As an illustration, we have applied the developed general theoryto theoretical description of a phase behavior of solution of the dipolar particles interactingwith each other through the electrostatic and Gaussian potentials – Gaussian-core dipolarmodel, which could be considered as a primitive model of dipolar polypeptides or betaines.For the latter model, we have derived an analytical expression for the total excess free energyof solution and, basing on it, described the phase behavior of system. We have predicted17hat at sufficiently low temperature, when the rather strong electrostatic dipole-dipole at-tractive interactions take place, the liquid-liquid phase separation occurs. We have shownthat the binodal (phase equilibrium curve) has an upper critical point and strongly asym-metric shape. We have established that increase in the energetic parameter of the Gaussianpotential at rather large concentrations of the dipolar particles makes the electrostatic freeenergy less negative.In conclusion, we would like to discuss further prospects of the developed theoreticalbackground. First of all, our theory can be applied to coarse-grained modelling of the phasebehavior of real aqueous solutions of biomacromolecules in salt-free solutions. Secondly, likethe previously formulated simpler field-theoretical statistical models , present theory canbe easily extended for the case of presence of salt ions in the bulk solution. The latter willallow us to study an influence of concentration of salt on the discussed in present paperthe liquid-liquid phase separation. Thirdly, the proposed theoretical background allows usto formulate the classical density functional theory for solutions of the Gaussian core dipo-lar particles, extending the density functional theory, proposed in papers for the simpleGaussian core model. Such a theory will allow us to describe the concentration profilesof soft-core dipolar molecules at the soft interfaces (with lipid or polymeric membranes),which is relevant for biological and industrial applications. Finally, it is necessary to verifypresented here theoretical predictions by means of the Monte Carlo or MD computer sim-ulations.
As it was shown in paper by means of MD computer simulation, atsufficiently low temperatures the system of soft charged dumbbells undergoes aliquid-liquid phase separation. It was obtained that critical parameters for softcharged dumbbells are very close to those are realized for the restricted primi-tive model. Due to the fact that the system of soft charged dumbbells is quitesimilar to studied here model, one can also expect in the MD simulation theliquid-liquid phase separation for the Gaussian core dipolar model, predictedby present theory. All these further investigations are the subject of the forthcomingpublications. 18
CKNOWLEDGMENTS
The publication was prepared within the framework of the Academic Fund Program atthe National Research University Higher School of Economics (HSE) in 2019-2021 (grant No19-01-088) and by the Russian Academic Excellence Project "5-100". The author thanksI.Ya. Erukhimovich for the motivating discussions.19
EFERENCES Andersen H C, Weeks J D, and Chandler D 1971
Phys. Rev. A (4) 1597 Andersen H C, Weeks J D, and Chandler D 1972
J. Chem. Phys. Likos C N 2001
Physics Reports Molotilin T Y, Maduar S R, and Vinogradova O I 2018
Phys. Rev. E Stillinger F H and Weber T A 1978 J. Chem. Phys. (8) 3837 Louis A A, Bolhuis P G, and Hansen J P 2000
Phys. Rev. E Archer A J, Evans R 2001
Phys. Rev. E Archer A J, Evans R 2002
J. Phys.: Condens. Matter Lang A, et al. 2000
J. Phys.: Condens. Matter Nogovitsyn E A, Gorchakova E S, and Kiselev M G 2007
Russ. J. Phys. Chem. Baeurle S A, Efimov G V, and Nogovitsyn E A 2006
Europhys. Lett. (3) 378 Baeurle S A, Efimov G V, and Nogovitsyn E A 2006
J. Chem. Phys. Frydel D and Levin Y 2018
J. Chem. Phys. Canchi D R and Garcia A E 2013
Annu. Rev. Phys. Chem. Haran G 2012
Curr. Opin. Struct. Biol. Lowe A B and McCormick C L 2002
Chem. Rev. Kudaibergenov S, et al 2006
Advances in Polymer Science Heldebrant D J, et al 2010
Green Chemistry Felder C E, Prilusky J, Silman I, and Sussman J L 2007
Nucleic Acids Res. Fedotova M V, Dmitrieva O A 2016
Amino Acids Kumar R and Fredrickson G H 2009
J. Chem. Phys. Eiberwiser A, et al.
J. Phys. Chem. B Onsager L 1936
J. Am. Chem. Soc. (8) 1486 Kirkwood J 1939
J. Chem. Phys. Hoye J S and Stell G 1976
J. Chem. Phys. Morozov K I 2007
J. Chem. Phys. Nienhuis G and Deutch J M 1972
J. Chem. Phys. Levin Y 1999
PRL (6) 1159 Zhuang B, Wang Z-G 2018
J. Chem. Phys. Coalson R D, Duncan A 1996
J. Phys. Chem. Abrashkin A, Andelman D and Orland H 2007
PRL Budkov Yu A, Kolesnikov A L and Kiselev M G 2015
EPL Budkov Y A, Kolesnikov A L, Kiselev M G 2016
J. Chem. Phys. Budkov Y A, et al 2018
Electrochimica Acta McEldrew M, et al 2018
J. Phys. Chem. Lett. Budkov Yu A 2018
J. Phys.: Condens. Matter (34) 344001 Budkov Yu A 2019
Journal of Molecular Liquids Chandler D 1977
J. Chem. Phys. Buyukdagli S and Ala-Nissila T 2013
Phys. Rev. E Buyukdagli S and Ala-Nissila T 2013
J. Chem. Phys. Buyukdagli S and Blossey R 2014
J. Phys.: Condens. Matter Buyukdagli S and Blossey R 2014
J. Chem. Phys. Gordievskaya Y D, Budkov Y A, Kramarenko E Y 2018
Soft Matter Martin J M, et al. 2016
J. Chem. Phys. Budkov Yu A and Kolesnikov A L 2016
Eur. Phys. J. E Budkov Yu A, Kalikin N N and Kolesnikov A L 2017
Eur. Phys. J. E Budkov Yu A, Kolesnikov A L and Kiselev M G 2015
J. Chem. Phys. Kolesnikov A L, Budkov Yu A, Basharova E and Kiselev M G 2017
Soft Matter Schiessel H and Pincus P 1998
Macromolecules Kundu P and Dua A 2014
J. Stat. Mech. Cherstvy A G 2010
J. Phys. Chem. B Dussi S, Rovigatti L, Sciortiono F 2013
Molecular Physics Buyukdagli S, Achim C V and Ala-Nissila T 2011
Journal of Statistical Mechanics: Theoryand Experiment Buyukdagli S and Ala-Nissila T 2012
J. Chem. Phys. Kubo R 1962
J. Phys. Soc. Jap. (7) 1100 Podgornik R 1989
J. Chem. Phys. Khokhlov A R, Khachaturian K A 1982
Polymer (12) 1742 Borue V Yu and Erukhimovich I Ya 1988
Macromolecules Brilliantov N V 1993
Phys Rev E (6) 4536 Lue L 2006
Fluid Phase Equilibria Kopanichuk I V, et al. 2016
Fluid Phase Equilibria Efimov G V, Nogovitsin E A 1996
Physica A Zubarev D N 1954
Dokl. Akad. Nauk SSSR
757 (in Russian) Iukhnovskii I R, Golovko M F 1980
Statistical theory of classical equilibrium systems (Kiev,Izdatel’stvo Naukova Dumka, In Russian) Zakharov A Yu, Loktionov I K 1999
Theoretical and Mathematical Physics (1) 532 Fisher M E and Levin Y 1993
Phys. Rev. Lett. Levin Y and Fisher M E 1996
Physica A Braun H, Hetschke R 2013
Phys. Rev. E Daub C D, Patey G N, Camp P J 2003
J. Chem. Phys. (5) 7952 Orkoulas G, Kumar S K, and Panagiotopoulos A J 2003
Phys. Rev. Lett. (4) 048303-1 Budkov Yu A, et al. 2015
J. Chem. Phys.142