Lifshitz theory of wetting films at three phase coexistence: The case of ice nucleation on Silver Iodide (AgI)
AAPS/123-QED
Lifshitz theory of wetting films at three phase coexistence:The case of ice nucleation on Silver Iodide
Juan Luengo-M´arquez and Luis G. MacDowell
Departamento de Qu´ımica F´ısica, Universidad Complutense de Madrid (Dated: February 9, 2021)
Abstract
Hypothesis:
As a fluid approaches three phase coexistence, adsorption may take place by thesuccessive formation of two intervening wetting films. The equilibrium thickness of these wettinglayers is the result of a delicate balance of intermolecular forces, as dictated by an underlyingsurface potential. The van der Waals forces for the two variable adsorption layers may be formu-lated exactly from Dzyaloshinskii-Lifshitz-Pitaevskii theory, and analytical approximations may bederived that extent well beyond the validity of conventional Hamaker theory.
Calculations:
We consider the adsorption equilibrium of water vapor on Silver Iodide where bothice and a water layers can form simultaneously and compete for the vapor as the triple point isapproached. We perform numerical calculations of Lifshitz theory for this complex system and workout analytical approximations which provide quantitative agreement with the numerical results.
Findings
At the three phase contact line between AgI/water/air, surface forces promote growthof ice both on the AgI/air and the water/vapor interfaces, lending support to a contact nucleationmode of AgI in the atmosphere. Our approach provides a framework for the description of adsorp-tion at three phase coexistence, and allows for the study of ice nucleation efficiency on atmosphericaerosols.
Keywords: Adsorption, Wetting, Phase coexistence, Surface thermodynamics, Van der Waals forces, LifshitzTheory, Hamaker constant, Heterogeneous nucleation, Ice, Silver Iodide a r X i v : . [ phy s i c s . c h e m - ph ] F e b . INTRODUCTION The preferential adsorption of gases and liquids on an inert solid substrate is ubiquitous incolloidal science, and is often accompanied by the formation of thick wetting films that spanfrom a few nanometers to several microns as two phase coexistence is approached . Thebehavior of the equilibrium layer formed can be fully characterized by an interface potential, g ( h ), which measures the free energy as a function of the film thickness, h . For thick films, g ( h ) is often dominated by long range van der Waals forces that can be described rigorouslywith the celebrated Dzyaloshinskii-Lifshitz-Pitaevskii (DLP) theory .A more complex situation arises when the adsorbed wetting film can further segregateand form a new layer between the solid substrate and the mother phase as three phasecoexistence is approached . This can be quite generally the case for substrates in contactwith multicomponent mixtures. Examples include the industrially relevant formation ofclathrate hydrates from aqueous solutions of oil or carbon dioxide , biologically relevantaqueous solutions of Dextran and Bovine Serum Albumin ; or theoretically relevant modelsolutions such as ethanol/n-alkane or lutidine/water mixtures . The wetting problemnow becomes considerably more complex, as the system can exhibit two thick layers of size,say, l and d , respectively, that are bounded by the substrate and the mother phase and canfeed one from the other depending on the prevailing thermodynamic conditions.In practice, this complex scenario can be realized for a very relevant one componenttest system, namely, atmospheric supercooled water vapor in close proximity to the triplepoint . As ice nucleates on the surface of inorganic aerosols, the resulting ice / watervapor interface is exposed. This surface is actually a complex system exhibiting a thinpremelted water layer, usually referred as Quasi-Liquid Layer (QLL). Its properties largelycondition several phenomena related to ice, such as the electrification of storm clouds, frostand snowflakes formation, or ice skating . Hereby the surface Van der Waals forces playa crucial role in the stabilization of a thick QLL .On the other hand, silver iodide may be used as a nucleus for ice formation , withapplications in cloud seeding to induce rainfall over wide areas. The conventional beliefwas that the capability of the AgI to influence the growth of ice was connected to theirlattice match. Nonetheless, recent studies have found substances with similar structuresthat do not promote nucleation. Thus they relate the faculty of the AgI to serve as ice2ucleating agent to the charge distribution in the substrate surface and to the Van derWaals forces . The present work aims at elucidating the role of Van der Waals forcesin the ice nucleating activity of silver iodide, by describing precisely the Van der Waalsinteractions in the framework of the Lifshitz theory, applied to the AgI - Ice - Liquid Water- Air system (Fig. 1).FIG. 1: Geometry of layered planar systems studied in this work. An inert substrate, AgI,in contact with water vapor (air), close to three phase coexistence can exhibit interveningphases of ice, with thickness l and water, with thickness d . Since the system is subject tomaterial equilibrium each phase feeds from the other. The equilibrium layer thickness of l and d are dictated by the bulk free energies and the surface potential. For the sake ofgenerality and conciseness of notation, we denote the phases AgI, ice, water, air as ’1’, ’2’,’3’, ’4’, respectively.The conventional view is that the van der Waals interactions between two surfaces resultfrom the summation of forces between pairs of particles. The interaction coefficients are thenadded into the Hamaker coefficient, A Ham , and the interaction for a layered planar systemprovides an energy A Ham / (12 πh ), being h the distance between the interacting media . TheLifshitz theory of the Van der Waals forces goes beyond the pairwise summation, by treatingthe Van der Waals forces through continuum thermodynamic properties of the polarizableinteracting media . As a particularly significant result of the detailed electromagnetictreatment, there appears a range of distances where the leading order interactions become retarded and decay as 1 /h . Accordingly, the effective interaction coefficient itself, A Ham ,is no longer a constant, but becomes a function of the film thickness as well. Then thegeneralized expression for the Van der Waals free energy becomes : g V dW ( h ) = − A Ham ( h )12 πh (1)In the study of ice growing at the AgI / air interface, the simple description based on inter-3ctions between AgI and air across ice is not valid. Because of the existence of premelting,we also expect the formation of a fourth layer of water in between ice and air. The extensionof Lifshitz theory to four media of which one is variable has been known for a long time .However, in the atmosphere the water layer can grow at the expense of the ice layer, de-pending on the prevailing conditions. Therefore both the thickness of the ice layer ( l ) andthat of the premelted layer ( d ) must be taken as variables. The results for the free energy g ( l, d ) as derived from the Lifshitz theory in the case that two of the layers are variable wasanticipated recently without proof . Here we provide the full derivation in the supplemen-tary section, discuss some of its implications and show how the apparently awkward generalresult may be approximated accurately in analytic form.In the next section we formulate the problem of three phase adsorption equilibrium gener-ally, and discuss the conditions for wetting. In Section III we concisely present the theoreticalframework of van der Waals forces in the modern formulation and generalize the resultsto a system of two layers of variable thickness sandwiched between infinite bodies. We thenshow that the very complicated expression is amenable to an approximate but very accurateanalytical representation. In Section IV we briefly discuss the parameterization of dielectricproperties for AgI, water and ice that is required as input into DLP theory. In Section V wepresent our results. These are structured as follows: first, we test the numerical approxima-tions performed to achieve closed analytical expressions. Secondly, the theory is exploited todescribe the van der Waals forces of the complex system AgI/Ice/water/air, and their rolein the ice nucleating efficiency of AgI. The conclusions of our work are presented in SectionVI. II. ADSORPTION EQUILIBRIUM AT THREE PHASE COEXISTENCE
Consider a fluid phase, say, medium ’4’, in contact with an inert substrate, say, medium’1’. Furthermore, consider that phase, ’4’, extends up to a finite but very large distance L → ∞ . Accordingly, ’4’ serves as a heat and mass reservoir that fixes the temperatureand chemical potential of the full system (of course, if ’4’ is a multicomponent mixture,it fixes the chemical potential of each of its components). Quite generally, the density ofphase ’4’ in the vicinity of the substrate is not that found in the bulk phase well away fromthe wall. Particularly, consider that phase 4 is approaching three phase coexistence, such4hat two additional bulk phases ’2’ and ’3’ are slightly metastable. In a bulk system inthe thermodynamic limit, slightly metastable means that these phases are not observed atall. However, close to an inert substrate, surface forces can change this situation, as oneof the metastable phases could preferentially adsorb between the substrate and the motherphase ’4’. By the same token, once, say, phase ’2’ has adsorbed preferentially in betweenthe substrate and phase ’4’, the third phase ’3’, could adsorb preferentially in between ’2’and ’4’, leading to a two layer system of phases ’2’ and ’3’ in between the substrate ’1’ andthe mother phase ’4’.The question is then what sets the equilibrium film thickness of the intervening layers,’2’ and ’3’, with thickness l and d , respectively.Since we assume the system is at fixed temperature and chemical potential as dictatedby the semi-infinite phase ’4’, the equilibrium states will be such that the grand free energyis a minimum . In the mood of capillarity theory, we assume the total free energy is thatof infinitely large bulk systems, plus the cost to form each of the interfaces. FollowingDerjaguin, however, we also need to account for the effective interaction between interfacesseparated by a finite distance, via a generalized surface potential g ( l, d ). Accordingly,the total surface free energy of the system is expressed as: ω ( l, d ) = − p l − p d − p ( L − l − d ) + γ + γ + γ + g ( l, d ) (2)where p i is the bulk pressure of phase i at the fixed temperature and chemical potential ofphase 4, and γ ij is the surface tensions between phase i and phase j . The first three termsaccount for the bulk free energy of the system; the next three correspond to the free energyto form the interfaces separating infinitely large bulk phases; and finally, g ( l, d ) accountsfor the missing interactions die to the finite extent of phases 2 and 3. By this token it followsimmediately that the surface potential is defined such that in the limit l → ∞ and d → ∞ , g = 0.It is convenient to write this result as an excess over the bulk free energy of a systemfilled with phase 4 only, so that:∆ ω ( l, d ) = − ( p − p ) l − ( p − p ) d + γ + γ + γ + g ( l, d ) (3)This general expression corrects previous results by additive constants . Here, it rationalizesin a nutshell the adsorption equilibrium of a three phase system. In practice, this does not5hange the equilibrium condition for l and d , which are obtained by equating to zero thepartial derivatives of the free energy with respect to l and d . In practice, here we will beconcerned with the situation where the system is exactly at bulk three phase coexistence,such that p = p = p . The above equation then allows us to generalize the condition forwetting at three phase coexistence.If Eq.3 has an extremal point at small l = l and d → ∞ , the system has only twointerfaces with cost γ and γ , so that it must hold:lim d →∞ g ( l , d ) = γ − γ − γ (4)This provides the wetting condition for phase 2 intervening between the substrate and phase’3’. Likewise, if Eq.3 has extrema at l → ∞ and finite d = d , the system forms interfaceswith a cost γ and γ , and we find:lim l →∞ g ( l, d ) = γ − γ − γ (5)which provides the wetting condition for phase 3 intervening between phase 2 and 4. Finally,if Eq.3 has an extrema at finite and small l = l and d = d (including the case where both l = d = 0), we find: g ( l , d ) = γ − γ − γ − γ (6)which corresponds to the condition of infinite adsorption layers of 2 and 3 intervening be-tween phase 1 and phase 4.The above results serve to illustrate the crucial significance of the surface potential g ( l, d ).Not only it dictates the allowed equilibrium values of layer thickness, but it also embodies allallowed wetting conditions in the system. The surface potential consists of contributions ofdifferent nature. First, a structural contributions , which is short range, as it decays in thelengthscale of the bulk correlation length of a few molecular diameters, conveys informationon the packing correlations between the hard core of the molecules; Second, van der Waalscontributions, which are long range and result from spontaneous electromagnetic fluctuationsof the media. Additionally, charged systems will also have electrostatic contributions, witha decay that is given by the Debye screening length.In the next section we discuss the calculation of van der Waals contributions, which, inthe absence of electrolytes are expected to dominate the long range behavior of the system.6 II. ANALYTIC APPROXIMATION FOR SURFACE VAN DER WAALS FORCESA. Lifshitz theory of Van der Waals forces
The starting point of this section is the general result for the surface free energy embodiedbetween two semiinfinite planar bodies, L and R , separated by an arbitrary number of layers( m, n, · · · ) due to van der Waals forces: g Lm...R = k B T π ∞ (cid:88) n =0 (cid:48) (cid:90) ∞ ρ dρ ln( D ERm...L D MLm...R ) (7)where D E and D M are the dispersion relations of the standing waves of electric and mag-netic modes in the system. The integral is performed over transverse components of themomentum, and the sum is performed over an infinite set of discrete Matsubara frequencies. k B is Boltzmann’s constant and we assume set temperature T = 273 .
15 K to the triple pointof water.The dispersion relations depend on the specific geometry of the system. For a layeredplanar system composed of two bulk bodies (say, L=1 and R=3) separated by a dielectric(m=2) with thickness ’ h ’, the result is well known: D E,M = 1 − ∆ E,M ∆ E,M e − ρ h (8)Assuming from now on that all magnetic susceptibilities are equal to one , the ∆ ij functionshave the form ∆ Eij = ρ j − ρ i ρ j + ρ i ∆ Mij = (cid:15) i ρ j − (cid:15) j ρ i (cid:15) i ρ j + (cid:15) j ρ i (9)where ρ i = ρ + (cid:15) i ξ n c , and (cid:15) i , the dielectric function of the medium i = 1 , , ξ n = 2 πk B T (cid:126) n, n = 0 , , ... (10)with c, the velocity of light and (cid:126) , Planck’s constant in units of angular frequency.In this work, we deal with the Van der Waals free energy between two bulk media (AgIand air), separated by two intervening layers (ice and water), with variable thicknesses ’ l ’and ’ d ’, respectively (Fig.1). The corresponding dispersion relation of the four media systemwas anticipated by Esteso et al. , and is demonstrated here in the supplementary material . The result is: D E,M = 1 − ∆ E,M ∆ E,M e − ρ l − ∆ E,M ∆ E,M e − ρ d − ∆ E,M ∆ E,M e − ρ l − ρ d (11)The subscripts 1-4 indicate, respectively, AgI, ice, liquid water and air, but we keep through-out this section this notation for the sake of generality.The related result for the interaction between a plate coated with a layer of fixed thickness l and another plate at a variable distance h has been known for a long time (c.f. ).Here, Eq.7 together with Eq.11 generalize this result for the case were both l and d arevariable. Accordingly, the layer thickness l and d stand on the same footing, and the freeenergy g ( l, d ) is now symmetrical with respect to the interchange of l and d , as expected.The result also satisfies a number of desirable physical properties, and is consistent withexpectations for the limiting cases where the layers either become infinitely thick or vanishaltogether.Firstly, one notices that in the limit where both d and l tend to infinity, D E,M →
1, sothat the surface potential vanishes, as implied in the discussion of the previous section.In the limiting case where d → ∞ , one readily finds from Eq.11 that D E,M → D E,M . Asa result, it follows: lim d →∞ g ( l, d ) = g ( l ) (12)and likewise, for l → ∞ , lim l →∞ g ( l, d ) = g ( d ) (13)With some additional algebraic work, we also find for the opposite limit of vanishing thick-ness that: lim d → ( g ( l, d ) − g ( d )) = g ( l ) (14)lim l → ( g ( l, d ) − g ( l )) = g ( d ) (15)These limiting cases serve as a check of consistency for the numerical calculation of g ( l, d ).As evidenced by the above results, g ( l, d ) contains information on g ( l ) and g ( d ).This can be illustrated explicitly upon linearization of the logarithmic term in Eq.7. Theresulting expression can be readily interpreted as given by: g ( l, d ) = g ( l ) + g ( d ) + ∆ g ( l, d ) (16)8ccordingly, g ( l, d ) may be given as a sum of g ( l ) and g ( d ), plus a correction∆ g ( l, d ) which accounts for the indirect interaction of the two macroscopic bodies acrossthe layers.Despite this simplification, the expressions for the three media potential still remain verydifficult to interpret intuitively. In the next section we derive analytical formula which allowto interpret g ( l, d ) easily and serve also as an accurate means to calculate the free energyefficiently.The results extend the theory for the calculation of interface potentials for three mediareported recently . In order to make the derivation self contained, we first summarize theresults for g from Ref. and then extend these ideas to the calculation of the correctionfunction ∆ g ( l, d ). B. Analytical approximations for the three media contributions
We start from the exact form of the 3-media contributions, g ( h ), where ’ h ’ now standsfor the thickness of layer ’2’. g ( h ) = k B T π ∞ (cid:88) n =0 (cid:48) (cid:90) ∞ ρ dρ ln( D E D M ) (17)Of course, now the dispersion relation would be that given in Eq. 8. Linearization of thelogarithm as implied in Eq.16, together with the change of variables ρ m → ρ and x → hρ m ,provides ( supplementary material 2 ): g ( h ) = − k B T πh ∞ (cid:88) n =0 (cid:48) (cid:90) ∞ r n x dx R ( n, x ) e − x (18)where R ( n, x ) = ∆ E ∆ E + ∆ M ∆ M , while the lower limit of the integral is r n = 2 h √ (cid:15) ξ n /c .Notice that the ∆ E,Mij functions are redefined in terms of x i instead of ρ i , given that x i = x + ( (cid:15) i − (cid:15) )(2 hξ n /c ) .Although this expression remains rather cumbersome, it has been shown recently thatvery accurate analytical approximations may be obtained using problem adapted one-pointGaussian quadrature rules . The one-point Gaussian quadrature (detailed in supplementarymaterial 3 ) allows the transformation of an integral (cid:82) f ( x ) w ( x ) dx without known primitive,into the product of a simple integral (cid:82) w ( x ) dx and the function f ( x ) evaluated at the9uadrature point, x . Essentially, this corresponds to the application of the mean valuetheorem, with the one-point Gaussian quadrature rule exploited as a means to estimate x .
1. First Gaussian Quadrature Approximation
First, let us separate the n = 0 contribution of the summation at 18. In this summandthe value of the Matsubara frequency (equation 10) is 0, which simplifies the expressions.Hence we can solve straightforwardly the integral and state g ( h ) = g ξ n =0123 ( h ) + g ξ n > ( h ),with g ξ n =0123 ( h ) = − k B T πh (cid:18) (cid:15) − (cid:15) (cid:15) + (cid:15) (cid:19) (cid:18) (cid:15) − (cid:15) (cid:15) + (cid:15) (cid:19) (19)Afterwards, we can obtain g ξ n > ( h ) by applying the one-point Gaussian quadrature ( sup-plementary material 4 ) for the rest of the summation, employing w ( x ) = x e − x as a weightfunction and f ( x ) = R E + R M = R ( n, x ). The results materialize in terms of r n , with thequadrature point x = (2 + 2 r n + r n ) / (1 + r n ) g ξ n > ( h ) = − k B T πh ∞ (cid:88) n =1 R ( n, x )(1 + r n ) e − r n (20)We denote the result of Eq.19 and Eq.20 as the First Gaussian Quadrature Approximation(FGQA).
2. Second Gaussian Quadrature Approximation
In order to simplify further the above result, our purpose is to turn the summation of Eq.20 into an integral. We begin by making use of the Euler-MacLaurin formula, finding outthat the corrections to the integral approximation are negligible ( supplementary material 5 ).Performing a change of variable to ν n = ν = ν T n , with ν T = 2 √ (cid:15) ξ T /c and ξ T = 2 πk B T / (cid:126) ,we find g ξ n > ( h ) = − c (cid:126) π h (cid:90) ∞ ν T (cid:101) R ( ν, x )(1 + hν ) e − hν dν (21)With (cid:101) R ( ν, x ) = R ( ν, x ) (cid:15) − / j − and j = (cid:16) d ln (cid:15) d ln ξ (cid:17) (this factor emerges from thechange of variable, as it is specified in supplementary material 5 ).The integral at the equation 21 is dominated by an exponential decay at large h , andby the algebraic decay of (cid:101) R ( ν, x ) when h tends to zero. Unfortunately, the crossover from10lgebraic to exponential decay is not amenable to analytical integration. Therefore, weintroduce an auxiliary exponential function , e − ν/ν ∞ , with the parameter ν ∞ chosen tomimic the scale at which the algebraic decay of (cid:101) R ( ν, x ) becomes significant. g ξ n > ( h ) = − c (cid:126) π h (cid:90) ∞ ν T (cid:101) Re ν/ν ∞ (cid:2) e − ν/ν ∞ (1 + hν ) e − hν (cid:3) dν (22)In this expression, the term in square brackets is an easily integrable weight function, whichconverges both at h → ∞ and at h →
0, as the original integrand. Accordingly, wecan exploit again the mean value theorem and perform a one point Gaussian quadrature( supplementary material 6 ), with f ( ν ) = (cid:101) R ( ν, x ) e ν/ν ∞ and w ( ν ) equal to the term insidethe brackets at 22. The outcome is the quadrature point ν ∗ = ν T + ν ∞ ξ , where ξ is anadimensional factor ξ = ( ν T h + 1)( ν ∞ h + 1) + 2 ν ∞ h ( ν ∞ h + 1) ( ν T h + 1) + ( ν ∞ h + 1) ν ∞ h (23)together with the approximate expression for the free energy as: g ξ n > ( h ) = − c (cid:126) ν ∞ π h (cid:101) R ∗ ξ (cid:101) F (24)where (cid:101) R ∗ ξ = (cid:101) R ( ν ∗ , x ) e ξ and (cid:101) F = ( ν T h + 1)( ν ∞ h + 1) + ν ∞ h ( ν ∞ h + 1) e − ν T h (25)Since (cid:101) R ∗ ξ depends on h weakly, the free energy is dominated by the function (cid:101) F . Inspectionof Eq.25 shows that g ξ n > depends on the two inverse length scales, ν ∞ and ν T . When h (cid:28) ν − ∞ , (cid:101) F becomes a constant, and g ξ n > ∝ h − . This is the Hamaker limit of non-retardedinteractions. For ν − ∞ (cid:28) h (cid:28) ν − T , (cid:101) F falls as 1 /h , so that g ξ n > ∝ h − , which correspondsto the Casimir regime of retarded interactions. Finally, when h (cid:29) ν − T , the evolution of theinteraction is dominated by the exponential e − ν T h . In this regime, g ξ n > vanishes altogether,and only the n = 0 term, Eq.19 of the van der Waals free energy survives. C. Analytical approximations for the four media correction
Consistent with the linearization approximation in Eq.16, the 4-media correction is givenas: ∆ g ( l, d ) = − k B T π ∞ (cid:88) n =0 (cid:48) (cid:90) ∞ ρ dρ R e − ρ l + ρ d ) (26)11here now R ( n, ρ ) = ∆ E ∆ E + ∆ M ∆ M and ∆ E,Mij as defined in Eq.9.In order to proceed, first notice from the definition of ρ i that one can write: ρ = ρ / −
12 ∆ (cid:15)c ξ n (27) ρ = ρ / + 12 ∆ (cid:15)c ξ n (28)where we have introduced ∆ (cid:15) = (cid:15) − (cid:15) , and ρ / = ρ +
12 ( (cid:15) + (cid:15) ) c ξ n . Then we assume small∆ (cid:15) and apply Taylor to get ρ ≈ ρ / − ξ n ∆ (cid:15)c ρ / (29) ρ ≈ ρ / + 14 ξ n ∆ (cid:15)c ρ / (30)Using these results, the exponential function in Eq.26 can be now factored into two simplerexponentials: e − ρ l + ρ d ) ≈ e − ρ / ( l + d ) e − ξ n ∆ (cid:15)c ρ / ( d − l ) (31)We call this the Similar Dielectric Function (SDF) approximation.Replacing this result in 26, we can now express ∆ g ( l, d ) as∆ g ( l, d ) = − k B T π (cid:80) ∞ n =0 (cid:48) (cid:82) ∞√ (cid:15) / ξnc dρ / ρ / R ( n, ρ / ) e − ρ / ( l + d ) e − ξ n ∆ (cid:15)c ρ / ( d − l ) (32)Where we have switched the integration variable to ρ / and define the mean dielectricresponse of the intervening media as: (cid:15) / = 12 ( (cid:15) + (cid:15) ) (33)With the purpose of following an analogous path to that previously described for the threemedia terms, we perform a second change of variable under the definition x = 2 ρ / ( l + d ),leading to ∆ g ( l, d ) = − k B T π ( l + d ) (cid:80) ∞ n =0 (cid:48) (cid:82) ∞ r n dx x R e ( n, x ) e − x (34)where R e ( n, x ) = R ( n, x ) e − ξ n ∆ (cid:15)c x ( d − l ) , and the lower integration limit is r n = 2( l + d ) √ (cid:15) / ξ n c .This result is now very similar to the three media potential of Eq.18, so we can findapproximate solutions along the same lines.12 . First Gaussian Quadrature Approximation First, we split the sum into the n = 0 term and the remaining contributions ( n > g ξ n =01234 ( l + d ) = − k B T π ( l + d ) (cid:18) (cid:15) − (cid:15) (cid:15) + (cid:15) (cid:19) (cid:18) (cid:15) − (cid:15) (cid:15) + (cid:15) (cid:19) (35)For the high frequency contributions, we proceed by performing a one-point Gaussianquadrature (again, supplementary material 4 ) under the definitions f ( x ) = R e ( n, x ) and w ( x ) = e − x , and get:∆ g ξ n > ( l, d ) = − k B T π ( l + d ) ∞ (cid:88) n =1 R e ( n, x )(1 + r n ) e − r n (36)where x = (2 + 2 r n + r n ) / (1 + r n ). To avoid any confusion, we highlight here that R e ( n, x )depends on ( d − l ), meaning that the dependence of g ξ n > on ( l, d ) is not as simple as ( l + d ).The equation 36 is the First Gaussian Quadrature Approximation under the Similar Dielec-tric Function Approximation (SDF-FGQA), and its validity will be proved quantitatively inthe results section.
2. Second Gaussian Quadrature Approximation
The same development of the previous SGQA is applied now to the equation 36. Byemploying the Euler-MacLaurin formula and operating the change of variable ν = ν T n =2 √ (cid:15) / ξ T n/c ( supplementary material 5 ) we reach∆ g ξ n > ( l, d ) = − c (cid:126) π ( l + d ) (cid:82) ∞ ν T (cid:101) R e ( ν, x )[1 + ( l + d ) ν ] e − ( l + d ) ν dν (37)With the definition (cid:101) R e ( ν, x ) = R e ( ν, x ) (cid:15) − / / j − / . Analogously, j / = (cid:16) d ln (cid:15) / d ln ξ n (cid:17) , andrecall here that ∆ (cid:15) = (cid:15) − (cid:15) , being 3 and 2 the subscripts of liquid water and ice. Thenwe introduce the auxiliary exponential, e − ν/ν ∞ , and proceed with the one point Gaussianquadrature over ν allowing f ( ν ) = (cid:101) R e ( ν, x ) e ν/ν ∞ and w ( ν ) = [1 + ( l + d ) ν ] e − ν/ν ∞ e − ν ( l + d ) ( supplementary material 6 ). This operation leads to∆ g ξ n > ( l + d ) = − c (cid:126) ν ∞ π ( l + d ) (cid:101) R e, ∗ ξ (cid:101) F ( l + d ) (38)13here (cid:101) R e, ∗ ξ = (cid:101) R e ( ν ∗ , x ) e ξ must be evaluated at the quadrature point ν ∗ = ν T + ν ∞ ξ , whilethe functions ξ and (cid:101) F , with h = l + d , are as Eq. 23 and Eq. 25, respectively.Since the leading order behavior of Eq.38 is given by (cid:101) F ( l + d ), the analogy of the correctionterm ∆ g ( l + d ) with the three media potential g ( h ) is made obvious. In practice, ∆ g ξ n > also depends on l − d , by virtue of the factor: R e ( ν ∗ , x ) = R ( ν ∗ , x ) e −
14 ∆ (cid:15)(cid:15) / ν ∗ x ( d − l ) (39)However, using the results for ν ∗ and x , one finds that R e ( ν ∗ , x ) provides only correctionsof order unity to the leading order dependence provided by the function (cid:101) F ( l + d ). D. Summary of results and outlook
In this section we have provided analytical expressions for the van der Waals free energyof two thick plates separated by two layers of variable thickness, l and d , g ( l, d ). Startingfrom the exact Lifshitz result (Eq. 7) with the appropriate dispersion relation for our system(Eq. 6), we perform the expansion of the logarithm and show that g ( l, d ) can be expressedin terms of the free energies for three media, g ( l ) and g ( d ), together with a correction∆ g ( l, d ).The terms g ( l ) and g ( d ) may be treated in a general form g RmL ( h ) (Eq. 18). Theirzero frequency contribution is taken independently (Eq. 19), while the rest is simplifiedthrough two Gaussian quadrature approximations and the Euler-MacLaurin formula (Eq.24).The four media term, g ( l, d ) can been worked out analogously, and yields, to leadingorder, exactly the same distance dependence than g RmL ( h ), with h = d + l .For qualitative purposes, this means that the complicated two variable dependence of g ( l, d ) can be described by a sum of one single variable functions, such that: g ( l, d ) ≈ C (cid:101) F ( l ) l + C (cid:101) F ( d ) d + C (cid:101) F ( l + d )( l + d ) (40)where the factors C ... are material parameters of the intervening media and (cid:101) F ( h ) is givenby Eq.25. A similar relation holds also under the approximation of purely pairwise additiveinteractions, with the function (cid:101) F ( h ) merely replaced by a constant factor .14 V. DESCRIPTION OF THE DIELECTRIC RESPONSE
So far we have described the different contributions to the surface van der Waals freeenergy, whose computation in the framework of the Lifshitz’s theory requires the knowledgeof the dielectrical properties of every substance implied, essentially through the ∆ ij functions(Eq. 9) that appear in the dispersion relations.In order to test our theory, we need to consider an explicit model for the dielectricproperties of the system. As a simple one single component test system, we consider theadsorption of water vapor on Silver Iodide just below water’s triple point. In this situation,we expect that a layer of ice of arbitrary thickness l can form, while, as the system approachesthe melting line, the ice surface can premelt and form a water layer of thickness d .Since the characterization of the dielectric properties is rather cumbersome, and the maingoal of this paper is to test the theory of the previous section, we provide here just a briefdescription. A complete bibliographic review of the extinction index of these substances,together with detailed explanation of the resulting parameterization has been presented asa part of a Master’s thesis , and will be published promptly in a forthcoming article. A. Dielectric response of AgI
The dielectric response of a substance evaluated at imaginary frequencies is a real andmonotonically decreasing function that drops at the frequencies at which that materialabsorbs. Every valid dielectric response must fulfill the limit (cid:15) ( iξ → ∞ ) = 1, meaningthat beyond the last absorption frequency only remains the dielectric response of the vac-uum. The dielectric function of AgI will be described following simple model called dampedoscillator , used when there are not many experimental optical properties available in thebibliography. This representation accounts for one absorption in the ultraviolet (UV) regionand another in the infrared (IR), and the parameterization requires only the knowledge ofthe static dielectric response, (cid:15) (0), the refractive index before the UV absorption, n UV , andthe absorption frequencies ω UV and ω IR . (cid:15) ( iξ ) = 1 + (cid:15) (0) − n UV ξ/ω IR ) + n UV −
11 + ( ξ/ω UV ) (41)The magnitudes corresponding to the characterization of the UV absorption, (cid:15) (0), n UV and ω IR , are reviewed from Bottger and Geddes . The remaining parameter, ω UV , is achieved15rom a calculation based on the evolution of the refractive index published by Cochrane and the Cauchy representation (see supplementary material 7 ). The complete parame-terization is displayed in table I. (cid:15) (0) n UV ω IR (eV) ω UV (eV)7.0 2.22 0.013 4.133 TABLE I: Parameterization of the dielectric response for AgI.
B. Dielectric responses of ice and liquid water
In this work we have employed a description of the dielectric functions of liquid waterand ice achieved from the numerical fit of experimental absorption spectra by means of theDrude model, also called Parsegian-Ninham model when it is employed in this framework .The complete bibliographic review of the extinction index of these substances, together withthe resulting parameterization, has been presented as a part of a Master’s thesis , and willbe the subject of a forthcoming article.For the case of water, we have selected extinction coefficients available in the literature,and choose those measured close to 0 degrees Celsius whenever possible . Relative tothe well known parameterization by Elbaum and Schick , our set of extinction coefficientsemploys measurements by Hayashi and Hiraoka which provide a complete description of thehigh energy band . The resulting parameterization is consistent with recent work whichuse the data of Ref. together with Infra Red absorption data at ambient temperature .For the dielectric response of ice, there appear to be far less recent measurements. Forthis reason, we have performed a parameterization largely based on the literature reviewby Warren . The resulting parameterization does not differ significantly from previouscalculations by Elbaum and Schick .As found by Fiedler et al. , the most significant feature in the novel parameterizationwith updated experimental data by Hayashi is that the dielectric constant of water remainshigher than that of ice at all relevant finite frequencies. As a result, the Hamaker functionfor the ice/water/air system is positive for all film thicknesses below the micron.16 . RESULTS AND DISCUSSIONA. Limiting cases First we will check the consistency of our exact solution of g ( l, d ) by computing theevolution of the function at the limits of l, d → ∞ and l, d →
0. These correspond to thecases displayed in Eq 12, 13, 14 and 15, and shown in the same order in Fig 2. The properfulfillment of these special conditions at infinite and zero thicknesses evinces the solidityof the exact Lifshitz result and its numerical solution in the present work. The continuityproperty of the surface potential g ( l, d ) is very convenient, since layers 2 and 3 need notadsorb preferentially onto the substrate 1. Whence, in the general case where one seeks anabsolute minimum of the surface potential, the case where either phase ’2’, phase ’3’ or bothare not favored thermodynamically is naturally built in. B. Numerical checks
Next we demonstrate the reliability of the First Gaussian Quadrature Approximation(FGQA) and the Second Gaussian Quadrature Approximation (SGQA) by comparing theoutcomes for the three media Hamaker function (Eq. 1). The exact result is given by Eq. 7with the dispersion relation expressed as in Eq. 8, while the FGQA calculation is performedthrough Eq. 20 and Eq. 19. The computation of the SGQA in Eq. 24 needs the knowledgeof the ν ∞ parameter. We achieve this by requiring the approximate expression, Eq.24 tomatch the exact free energy in the limit of vanishing film thickness. g SGQA ( h →
0) = g exact ( h →
0) (42)Essentially, this amounts to setting ν ∞ so as to match the exact Hamaker constant.Figure 3 presents the three media Hamaker function of the systems AgI/Ice/Water, Ice/Water/Air,AgI/Ice/Air and AgI/Water/Air, and illustrates the accuracy of the analytic approxima-tions that we have developed. The resulting Hamaker functions have been divided by k B T at this representation in order to make their values easier to handle, and offering also a ratioof their strength with respect to the thermal energy.With g AgI/Ice/W ater and g Ice/W ater/Air we have the first two terms required to describe com-pletely the free energy of the system (Eq. 16). The remaining contribution is given by17IG. 2: Limiting cases of g ( l, d ) at infinite and zero thicknesses. Top left exposes thecase contained in Eq. 12, top right represents Eq. 13, bottom left displays the limit in Eq.14, and bottom right the one in Eq. 15. Black lines correspond to exact calculations of thesurface potentials indicated on the y axis, while the symbols are the correspondingprediction obtained from g ( l, d ) as indicated in the cited equations. The labels 1, 2, 3, 4stand for AgI, ice, water and vapor, respectively. In practice for numerical purposes thelarge thickness limit is evaluated at 10 − m, while the vanishing thickness limit isevaluated at 10 − m. All surface energies are given in J/m .∆ g ( l, d ), whose values for the exact result and Second Gaussian Quadrature Approxima-tion under the Similar Dielectric Function Approximation (SDF-SGQA) are displayed in Fig4. Note that even if the dependence of ∆ g ( l, d ) on the thicknesses is more complicatedthan the simple ( l + d ) that appears in Eq 38 due to the ( d − l ) in Eq.39, as we canexpress this as ( l + d )( d − l ), it is still possible to employ the previous method to get ν ∞ for( l + d ) →
0. In the exact expression of ∆ g ( l, d ) in Eq. 26 we have to assume l = d → l + d ) → l = d .18IG. 3: Comparison of exact Hamaker functions with FGQA and from the SGQAapproximations used in this work. Here we show all relevant combinations of three mediaformed from AgI, ice, water and vapor: AgI/Ice/Water (top, left), Ice/Water/Air (top,right), AgI/Ice/Air (bottom, left) and AgI/Water/Air (bottom, right). Results are givenin units of k B T .The numerical solution of the exact formula (Eq. 7) with the dispersion relation of Eq. 11is presented in Fig. 5, left. The range is fixed from 0.3 nm to 30 nm, which is approximatelythe expected range of lengths of relevance of the Van der Waals forces, but the behavior isquite monotonous for wider thicknesses. This result is then taken as a reference to com-pare with energy curves emerging from the employment of the First Gaussian Quadratureunder the Similar Dielectric Function approximation (Fig. 5, right), whose results arisefrom solving equations 19, 20 and 36. The outcome never exceeds a relative error of 3%with respect to the exact result. Since we have proved how good this approximation works,we can now take advantage of its drastically lower computation time to perform a moreexhaustive calculation of g ( l, d ), using now a finer mesh that otherwise would require a19IG. 4: Correction ∆ g ( l, d ) obtained from the exact expression (left) and from theSDF-SGQA (right). It supplies a repulsive contribution when both thicknesses (’l’ for iceand ’d’ for liquid water) are low, preventing the total g ( l, d ) to present an absoluteminimum when both the liquid water layer and the ice layer disappear. Here thethicknesses are displayed in decimal logarithmic scale and surface energies are given in J/m .FIG. 5: Comparison of the full surface potential from numerical results with theoreticalcalculations. Exact g ( l, d ) (left), obtained from the Eq. 7 with the dispersion relation ofEq. 11. FGQA of g ( l, d ) (right) under the SDF approximation. In both cases, the axes’l’ and ’d’ are the thicknesses of the ice and the liquid water layers represented in decimallogarithmic scale, while surface free energies are given in J/m .huge amount of time to be completed.The FGQA approximation has also the advantage of being quite easier to compute interms of complexity of the algorithm. The computation of the SGQA is even faster, butrequires the calculation of the ν ∞ parameter.20 . Van der Waals forces of Water and ice adsorbed on AgI
1. Three media systems
In this section, we discuss the implications of every three media - Hamaker functiondisplayed in Fig. 3.Recall that from Eq. 1, the A Ham ( h ) has opposite sign to g ( h ). Thus the AgI/Ice/Waterand Ice/Water/Air functions (top in Fig. 3) provide a negative surface free energy, withabsolute minima at vanishing film thickness. This means that van der Waals forces donot favor the growth of ice at the AgI/water interface, and do not favor the growth of apremelting film of water at the ice/air interface either.On the other hand, both the AgI/Ice/Air and AgI/Water/Air systems (Fig. 3, bottom)present negative values of the Hamaker function, which implies positive and monotonousdecreasing surface free energies. Accordingly, van der Waals forces favor the growth of bothice and water thick films at the AgI/air interface.In practice, the ultimate behavior of growing films on the AgI surface is dictated by abalance of short range structural forces and long range van der Waals forces. Since theice/air interface is known to exhibit a significant amount of premelting , and van der Waalsforces appear to inhibit growth of a liquid film, we conclude that the existence of premeltingon ice is the result of short range structural forces, in agreement with recent findings fromcomputer simulation .
2. AgI/Ice/Water/Air
Once again, recall that according to Eq. 16, the four media Van der Waals free energy, g ( l, d ), is given as the sum of two three-media contributions and ∆ g ( l, d ). The resultsfor the four media surface free energy, displayed in the Fig. 5 (either left or right), illustratehow as ’l’ (the ice width) increases, the Van der Waals free energy in the Eq. 16 is completelygoverned by g ( d ), and analogously, if ’d’ (water layer thickness) becomes very large, thesystem behaves as if the air was not there, and the whole dispersive interaction comes nowby the hand of g ( l ). When ’l’ and ’d’ are negligible, the ∆ g ( l, d ) term in Eq. 16contributes with a positive energy, so that the total g ( l, d ) of the system does not havean absolute minimum at l = d = 0. 21ooking closely figure 5 at this scale, we appreciate how the lower set of minima appears forextremely low ’d’ ( d →
0, liquid water almost disappear) and for several - increasing valuesof ice width. Indeed, the Lifshitz theory extended to four media supports and confirms anintuitive result from the comparison between the Hamaker functions of three media systemsin the previous section: the Van der Waals interactions favor the growth of either ice orwater at the AgI/air interface only if the thickness of the other substance (water or ice)remains close to zero.This is confirmed by noticing that g ( l, d ) < γ AgI/air is smaller than the sum of γ AgI/ice , γ ice/water and γ water/air , implying that the adsorption of large ice and water films inbetween the AgI/air interface is unfavorable.
3. AgI/Water/Ice/Air
We have assumed all along that the nucleation occurs with the arrangement AgI/Ice/Water/Air.Nevertheless, this is not necessarily true, and it worth to check the behavior of the systemAgI/Water/Ice/Air. For this purpose we solve numerically the exact Eq. 7 with the disper-sion relation of Eq. 11, but this time with the following meaning of the indices: 1 = AgI, 2= Liquid Water, 3 = Ice, 4 = Air. The result is presented in Fig. 6,The resulting surface potential is a positive (and therefore repulsive) energy, whose zeros(minima) are placed at large values of both thicknesses, meaning that the system will tryto lower its energy by increasing the amount of both substances. Therefore, it appears thatdiverging layers of condensed water can grow on AgI if water is first adsorbed onto AgI andice grows in between water and air. On the contrary, growth of ice in between AgI andwater is not favored. The origin of this difference can be traced to the larger propensity ofice to grow in between water and air than that of water to form in between ice and air. Thiscan be checked by looking at Fig. 6 along the axis of large water thickness, where it is seenthat the free energy decreases by letting ice grow. This corresponds to the phenomenon ofsurface freezing, which here is seen to be favored by van der Waals forces, consistent withrecent calculations . 22IG. 6: Surface free energy for g ( l, d ) and the layer order AgI/Water/Ice/Air. The axes’l’ and ’d’ are now the thicknesses of the liquid water and the ice layers, in decimallogarithmic scale. Surface free energies in J/m .
4. Mechanism of ice nucleation
We will use this last part of the discussion to summarize some conclusions and consider-ations about the role of AgI as an ice nucleator.Our results show that van der Waals forces promote the condensation of both ice andwater on the AgI/air interface, a result which appears to be quite in agreement with theknown nucleation efficiency of AgI . In fact, we find that the Hamaker constant of theAgI/ice/air system is larger (in absolute value) than that of AgI/water/air. This means thatvan der Waals forces actually promote the freezing of water vapor over the condensationof liquid water onto AgI. This expectation from the Hamaker constants is confirmed byinspection of Fig.5, which shows that indeed, the free energies along the d → l → d → ∞ axis, or merely, by inspection of the Hamaker function of the AgI/ice/water system,Fig.3.From these observations, we can conclude that van der Waals forces actually favor adeposition mode of ice nucleation from the vapor phase. Of course, this is quite at oddswith experimental findings, which show that AgI can nucleate ice both from vapor or waterat a similar undercooling of about 4 ° C .The reason for this apparent discrepancy is that the ultimate behavior of the system isdictated by a balance of both structural and van der Waals forces. Since an undercooling ofat least 4 ° C is required for ice to grow from supercooled water vapor, there must be shortrange structural forces which oppose mildly to ice growth. Otherwise, in the absence of shortrange forces, van der Waals forces would favor nucleation of ice from the vapor without anyundercooling. This is very much consistent with computer simulations by Shevkunov, whichindicate that a one monolayer thick ice film can form at ice vapor saturation, but furthergrowth is mildly activated .On the contrary, van der Waals forces do not promote ice growth for AgI immersed inwater, but simulations consistently show that Ag + exposed surfaces readily nucleate ice atmild supercooling . This means that short range structural forces do favor ice growth atthe AgI/water surface, and a small activation is required because of the unfavorable van derWaals interactions.Interestingly, it has been suggested that AgI is actually most efficient in the contactmode, whereby ice nucleation is promoted for AgI particles in contact with condensed waterdroplets . An explanation for this effect can be provided by assuming that nucleationactually occurs at the three phase contact line formed between AgI, water and air , but thisrequires a favorable line tension. Our results lend support to such favorable phenomenon.Indeed, we find that van der Waals forces promote growth of a thin ice layer on AgI incontact with water vapor. Such growth is ultimately very slow, because of the low vaporpressure of ice. However, if such thin ice layer is formed at the contact line of an AgIparticle with the air/water interface, a large reservoir of undercooled bulk water would bemade available for the thin ice layer to continue growing at the water/air interface, where,as we have seen, van der Waals forces promote surface freezing.24 I. CONCLUSIONS
In this study we work out the exact Lifshitz theory of van der Waals forces for a substratein the neighborhood of three phase coexistence, where two adsorbed layers of variable thick-ness can form on the substrate. The theory goes well beyond a Hamaker theory of pairwiseadditive forces, and naturally builds in the role of polarization of the condensed phases. Bythis token, it provides the Hamaker regime for small layer thicknesses, but also incorporatesthe Casimir regime of retarded interactions.Accurate analytical approximations are provided which improve conventional treatmentsof non-retarded van der Waals forces and apply also to the regime of retarded interactions.This extends the validity of the calculations from film thicknesses barely decades of nanome-ters to arbitrary large values. Unlike the Gregory equation advocated by Israelachvili , thisis achieved without any ad-hoc parameter . By this token, we find one can now accuratelyevaluate the free energy in the Hamaker and Casimir regimes with the same data that isrequired to estimate Hamaker constants . The proper account of retardation does notonly provide better accuracy. It is a major qualitative improvement, as the crossover fromnon-retarded to retarded interactions is often accompanied by a sign reversal of the van derWaals forces.Our results are applied to the study of water vapor adsorption on AgI, where both layersof ice and water can form close to the triple point. Previous studies on the ice nucleationefficiency of AgI have provided insight into the short range interactions of either undercooledwater or vapor with the AgI surface . However, in the atmosphere both ice and watercompete simultaneously for the vapor phase and it is very difficult to assess the relativestability of thick water and ice films from simulation. Our results indicate that van derWaals forces stabilize the growth of ice films at the AgI/air interface, but on the contrary,inhibit the growth of ice in the immersion mode. Importantly, van der Waals forces alsopromote growth of thick ice films at the air/water interface. This explains the intriguingobservation of sub-surface nucleation in experiments and computer simulations , but alsohelps understand how the water/AgI/air contact line could promote ice nucleation . TheAgI/vapor surface provides the site for stabilized ice layers that serve as embryos for thegrowth of a stable ice film at the water/vapor interface, thus lending support to contactmode freezing of AgI particles found in experiments .25ur results provide a general framework to gauge the role of van der Waals forces inthe colloidal sciences, where multicomponent solutions display multiphase coexistence ubiq-uitously, and provide a solid background to assess how long range forces condition the icenucleation efficiency of atmospheric aerosols. ACKNOWLEDGMENTS
We thank F. Izquierdo-Ruiz and Pablo Llombart for helpful discussions and assistance.We gratefully acknowledge funds from the Spanish Agencia Estatal de Investigaci´on underGrant No. FIS2017-89361-C3-2-P.
AUTHORS CONTRIBUTIONS
Juan Luengo: Methodology, Formal analysis, Visualization, Software, Investigation,Writing-Original Draft; Luis G. MacDowell: Validation, Conceptualization, Methodology,Writing-Review and Editing. B. Derjaguin, Langmuir , 601 (1987). P. G. de Gennes, F. Brochard-Wyart, and D. Qu´er´e,
Capillarity and Wetting Phenomena (Springer, New York, 2004) pp. 1–292. J. N. Israelachvili,
Intermolecular and Surfaces Forces , 3rd ed. (Academic Press, London, 1991)pp. 1–674. M. Schick, in
Liquids at Interfaces , Les Houches Lecture Notes (Elsevier, Amsterdam, 1990) pp.1–89. I. E. Dzyaloshinskii, E. M. Lifshitz, and L. P. Pitaevskii, Soviet Physics Uspekhi , 153.175(1961). S. Takeya and R. Ohmura, J. Chem. Eng. Data , 1880 (2006). G. Aspenes, L. Dieker, Z. Aman, S. Høiland, A. Sum, C. Koh, and E. Sloan, J. Colloid. InterfaceSci. , 529 (2010). R. Nakane, E. Gima, R. Ohmura, I. Senaha, and K. Yasuda, J. Chem. Thermodyn. , 192(2019). M. Bostr¨om, R. W. Corkery, E. R. A. Lima, O. I. Malyi, S. Y. Buhmann, C. Persson, I. Brevik,D. F. Parsons, and F. Johannes, ACS Earth and Space Chem. , 1014 (2019). Y. Antonov and B. A. Wolf, Langmuir , 6508 (2014). S. Rafa¨ı, D. Bonn, and J. Meunier, Physica. A , 31 (2007). C. Hertlein, L. Helden, A. Gambassi, S. Dietrich, and C. Bechinger, Nature , 172 (2008). J. G. Dash, A. W. Rempel, and J. S. Wettlaufer, Rev. Mod. Phys. , 695 (2006). R. Rosenberg, Phys. Today , 50 (2005). B. Slater and A. Michaelides, Nat. Rev. Chem , 172 (2019). Y. Nagata, T. Hama, E. H. G. Backus, M. Mezger, D. Bonn, M. Bonn, and G. Sazaki, Acc.Chem. Res. , 1006 (2019). M. Elbaum and M. Schick, Phys. Rev. Lett. , 1713 (1991). J. Fiedler, M. Bostr¨om, C. Persson, I. Brevik, R. Corkery, S. Y. Buhmann, and D. F. Parsons, J.Phys. Chem. B , 3103 (2020), pMID: 32208624, https://doi.org/10.1021/acs.jpcb.0c00410. B. Vonnegut, Journal of applied physics , 593 (1947). C. Marcolli, B. Nagare, A. Welti, and U. Lohmann, Atmospheric Chemistry and Physics ,8915 (2016). G. Fraux and J. P. K. Doye, J. Chem. Phys. , 216101 (2014),https://doi.org/10.1063/1.4902382. S. A. Zielke, A. K. Bertram, and G. N. Patey, J. Phys. Chem. B , 9049 (2015). B. Glatz and S. Sarupria, J. Chem. Phys. , 211924 (2016),https://doi.org/10.1063/1.4966018. V. A. Parsegian,
Van der Waals Forces (Cambridge University Press, Cambridge, 2005) pp.1–311. B. W. Ninham and V. A. Parsegian, J. Chem. Phys. , 3398 (1970),https://doi.org/10.1063/1.1674507. R. Podgornik, P. L. Hansen, and V. A. Parsegian, J. Chem. Phys. , 1070 (2003),https://doi.org/10.1063/1.1578613. V. Esteso, S. Carretero-Palacios, L. G. MacDowell, J. Fiedler, D. F. Parsons, F. Spallek,H. M´ıguez, C. Persson, S. Y. Buhmann, I. Brevik, and M. Bostr¨om, Phys. Chem. Chem.Phys , 11362 (2020). N. G. van Kampen, B. R. A. Nijboer, and S. K., Phys. Lett. , 307 (1968). B. W. Ninham, V. A. Parsegian, and G. H. Weiss, J. Stat. Phys. , 323 (1970). H. B. Callen,
Thermodynamics and an Introduction to Thermostatistics (John Wiley & Sons,New York, 1985). L. G. MacDowell, J. Chem. Phys. , 081101 (2019). M. M¨uller, L. G. MacDowell, P. M¨uller-Buschbaum, O. Wunnike, and M. Stamm, J. Chem.Phys. , 9960 (2001). P. Loskill, H. H¨ahl, T. Faidt, S. Grandthyll, F. M¨uller, and K. Jacobs, Adv. Colloid InterfaceSci. , 107 (2012), interfaces, Wettability, Surface Forces and Applications: SpecialIssue in honour of the 65th Birthday of John Ralston. D. N. Simavilla, W. Huang, C. Housmans, M. Sferrazza, and S. Napolitano, ACS Cent. Sci ,755 (2018). J. Luengo and L. MacDowell,
Van der Waals Forces at Ice Surfaces with Atmospheric Interest ,Master’s thesis, Facultad de Ciencias (2020). D. B. Hough and L. R. White, Adv. Colloid Interface Sci. , 3 (1980). G. Bottger and A. Geddes, The Journal of Chemical Physics , 3000 (1967). G. Cochrane, Journal of Physics D: Applied Physics , 748 (1974). L. Bergstr¨om, Adv. Colloid Interface Sci. , 125 (1997). V. Parsegian and B. Ninham, Nature , 1197 (1969). H. R. Zelsmann, Journal of molecular structure , 95 (1995). D. J. Segelstein,
The complex refractive index of water , Ph.D. thesis, University of Missouri–Kansas City (1981). D. M. Wieliczka, S. Weng, and M. R. Querry, Applied optics , 1714 (1989). H. Hayashi and N. Hiraoka, The Journal of Physical Chemistry B , 5609 (2015). J. Wang and A. V. Nguyen, Adv. Colloid Interface Sci. , 54 (2017). S. G. Warren and R. E. Brandt, J. Geophys. Research , D14220 (2008). D. T. Limmer and D. Chandler, J. Chem. Phys. , 18C505 (2014). J. Benet, P. Llombart, E. Sanz, and L. G. MacDowell, Mol. Phys. , 2846 (2019). P. Llombart, E. G. Noya, D. N. Sibley, A. J. Archer, and L. G. MacDowell, Phys. Rev. Lett. , 065702 (2020). H. R. Pruppacher and J. D. Klett,
Microphysics of Clouds and Precipitation (Springer, Heidel-berg, 2010). S. V. Shevkunov, Colloid J. , 497–508 (2005). Y. S. Djikaev and E. Ruckenstein, J. Phys. Chem. A , 11677 (2008), pMID: 18925734,https://doi.org/10.1021/jp803155f. J. Gregory, J. Colloid. Interface Sci. , 138 (1981). H.-J. Butt and M. Kappl,
Surface and Interfacial Forces (Wiley-VCH, Weinheim, 2010). A. J. Durant and R. A. Shaw, Geo. Phys. Res. Lett. (2005), 10.1029/2005GL024175,https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2005GL024175. A. Haji-Akbari and P. G. Debenedetti, Proc. Natl. Acad. Sci. U.S.A. , 3316 (2017). uan Luengo and Luis G. MacDowell
1. Derivation of the dispersion relation of the exact Lifshitz formula
The dispersion relation of the system is an expression that all electromagnetic waves acrossthe system must hold. These waves have an electric, (cid:126)E ( t ), and magnetic, (cid:126)H ( t ), field functionof time, t , with the form (cid:126)E ( t ) = Re (cid:32)(cid:88) ω (cid:126)E ω e − iωt (cid:33) (43) (cid:126)H ( t ) = Re (cid:32)(cid:88) ω (cid:126)H ω e − iωt (cid:33) (44)Being i the imaginary unit, ω the frequency, and (cid:126)E ω , (cid:126)H ω the amplitudes of every field atthat frequency. Each must respect a certain wave equation ∇ (cid:126)E = (cid:15)µc ∂ (cid:126)E∂t (45) ∇ (cid:126)H = (cid:15)µc ∂ (cid:126)H∂t (46)With (cid:15) and µ the dielectric and magnetic permeability and c the velocity of light. We canturn Eq. 45 and 46 into differential equations of the components of the position solving thepartial derivative with time from Eq. 43 and 44, respectively. In our system, the vectors (cid:126)E and (cid:126)H have periodic x, y components with the general form f k ( z ) e i ( ux + vy ) . Then solve thederivatives in the resulting differential equation to get the general f (cid:48)(cid:48) ( z ) = ρ k f ( z ), with ρ k = u + v − (cid:15) k µ k ω c (47)Where we define later ρ = u + v . The general solution yields f k,α ( z ) = A k,α e ρ k z + B k,α e − ρ k z ,with α = x, y, z for every substance k at the system, here we name the layers generally asL/m/n/R. Finaly, we impose ∇ · (cid:126)E = 0 and ∇ · (cid:126)H = 0 to reach that A z = − iρ ( uA x + vA y ) (48) B z = iρ ( uB x + BA y ) (49)Which is true separately at every substance L, m, n, R. Next let us apply the boundaryconditions at every interface, setting the z axis perpendicular to the interfaces and the z = 0point at the interface L/m. Now the thickness of m will be ’l’ and that of the substance n30ill be ’d’. • z = 0, interface L/m. Here B Lx,y,z must be zero so that f L does not go to infinite when ztends to −∞ . Then we impose the conditions E Lx = E mx and E Ly = E my , sum the resultingequations and use Eq. 48 and 49 to get − ρ L A Lz + ρ m A mz − ρ m B mz = 0 (50)On the other hand, it must also be true that (cid:15) L E Lz = (cid:15) m E mz , that is (cid:15) L A Lz = (cid:15) m A mz + (cid:15) m B mz (51) • z = l , interface m/n. Here we apply basically the same procedure. Notice that now theexponentials do not vanish, so we reach from E mx = E nx and E my = E ny − ρ m A mz e ρ m l + ρ m B mz e − ρ m l + ρ n A nz e ρ n l − ρ n B nz e − ρ n l = 0 (52)And from the condition (cid:15) m E mz = (cid:15) n E nz (cid:15) m A mz e ρ m l + (cid:15) m B mz e − ρ m l = (cid:15) n A nz e ρ n l + (cid:15) n B nz e − ρ n l (53) • z = l + d , interface n/R. Here A Rx,y,z must be zero so that f R does not go to infinite as ztends to ∞ . From E nx = E Rx and E ny = E Ry now we have − ρ n A nz e ρ n ( l + d ) + ρ n B nz e − ρ n ( l + d ) − ρ R B Rz e − ρ R ( l + d ) = 0 (54)And finally, from (cid:15) n E nz = (cid:15) R E Rz we get (cid:15) n A nz e ρ n ( l + d ) + (cid:15) n B nz e − ρ n ( l + d ) = (cid:15) R B Rz e − ρ R ( l + d ) (55)So far we have used the boundaries for the electric fields. Solving the system of equationsformed by Eq. 50, 51, 52, 53, 54 and 55 for the six variables { A i } , { B i } , leads to thedispersion relation D M = 0, where D M is the determinant:31 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ρ L ρ m − ρ m (cid:15) L − (cid:15) m − (cid:15) m − ρ m e ρ m l ρ m e − ρ m l ρ n e ρ n l − ρ n e − ρ n l (cid:15) m e ρ m l (cid:15) m e − ρ m l − (cid:15) n e ρ n l − (cid:15) n e − ρ n l
00 0 0 − ρ n e ρ n ( l + d ) ρ n e − ρ n ( l + d ) − ρ R e − ρ R ( l + d ) (cid:15) n e ρ n ( l + d ) (cid:15) n e − ρ n ( l + d ) − (cid:15) R e − ρ R ( l + d ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = D M Solving for the determinant explicitly gives a sum of 32 terms wich is not particularlyinsightful. However, one notices that products of two matrix elements of the form d ij d kl have common factors with terms d il d kj , or alternatively, terms of the form ρ i (cid:15) j share commonfactors with ρ j (cid:15) i . Therefore, we organize the 32 terms as: D M =( ρ L (cid:15) m − ρ m (cid:15) L ) [( ρ R (cid:15) n − ρ n (cid:15) R )( ρ m (cid:15) n + ρ n (cid:15) m ) α − βγ − + ( ρ R (cid:15) n + ρ n (cid:15) R )( ρ n (cid:15) m − ρ m (cid:15) n ) α − β − γ ] δ − +( ρ m (cid:15) L + ρ L (cid:15) m ) [( ρ m (cid:15) n − ρ n (cid:15) m )( ρ n (cid:15) R − ρ R (cid:15) n ) αβγ − − ( ρ n (cid:15) m + ρ m (cid:15) n )( ρ R (cid:15) n + ρ n (cid:15) R ) αβ − γ ] δ − (56)where we have introduced the symbols α = e ρ m l , β = e ρ n l , γ = e ρ n ( l + d ) and δ = e ρ R ( l + d ) , forshort.Notice now that the determinat is a product of elements ( ρ i (cid:15) j + ρ j (cid:15) i ), which correspondto denominators of the function ∆ Mij in the main text (see also below, Eq.58); and elements( ρ i (cid:15) j − ρ j (cid:15) i ) which correspond to numerators of ∆ Mij . Since, from the dispersion relation, thisdeterminant must vanish, we can now multiply and divide by constants without changingthe result. Therefore, we divide by factors of the form ( ρ i (cid:15) j + ρ j (cid:15) i ), and further divide by αβ − γδ − to get: D M = 1 − ∆ MLm ∆ Mnm e − ρ m l − ∆ Mmn ∆ MRn e − ρ n d − ∆ MLm ∆ MRn e − ρ m l e − ρ n d (57)∆ Mij = ρ j (cid:15) i − ρ i (cid:15) j ρ j (cid:15) i + ρ i (cid:15) j (58)which is the sought result. Notice that here D M is symmetrical with respect to the inter-change of l and d .The result for a system of one single medium between two plates, with one plate coated bya layer of fixed thickness, d , can be obtained from this expression readily. In that situation32he physical requirement is that D M →
1, as l → ∞ at fixed d . From Eq.57, we find instead:lim l →∞ D M = 1 − ∆ Mmn ∆ MRn e − ρ n d (59)Therefore, a new dispersion relation, consistent with the mentioned physical requirementcan be obtained simply dividing D M by lim l →∞ D M . This yields: D (cid:48) M = 1 − ∆ MLm ∆ Mnm + ∆
MLm ∆ MRn e − ρ m d − ∆ Mmn ∆ MRn e − ρ n d e − ρ m l (60)which is the well known dispersion relation used for systems with one coated layer of fixedsize. Because of the choice of normalization condition, the dispersion relation is no longersymmetrical with respect to the interchange of l and d , since they now stand on a differentfooting.Notice that in practice, this normalization of the dispersion relation amounts to sub-stracting g mnR ( d ) to g LmnR ( l, d ) as discussed in Eq.14-15 of the main text for the particularchoice d → D E arises from a similar analysis of the boundaries for the magneticfields. The full dispersion relation emerges by setting D = D M D E , with D = 0.33 . Change of variable to ’x’ Our purpose here is to develop explicitly the route from the three media exact surfacefree energy to the equation over which the First Gaussian Quadrature Approximation isperformed. Then we start from the generalized form g LmR ( h ) = k B T π ∞ (cid:88) n =0 (cid:48) (cid:90) ∞ ρ dρ ln( D ELmR D MLmR ) (61)With D E,MLmR = 1 − ∆ E,MLm ∆ E,MRm e − ρ m h , being h here the thickness of the medium ’m’. Nextwe assume that the term after the 1 at D E,MLmR is small, and we expand the logarithm as ln (1 − x ) ≈ = − x . At the same time, let us change the variable to ρ m = ρ + (cid:15) m ξ n c . Weperform this transformation simply through ρ m dρ m = ρdρ , and the lower limit of the integralnow is ρ m ( ρ = 0) = (cid:113) (cid:15) m ξ n c g LmR ( h ) = − k B T π ∞ (cid:88) n =0 (cid:48) (cid:90) ∞ (cid:114) (cid:15)mξ nc ρ m dρ m (∆ MLm ∆ MRm + ∆
ELm ∆ ERm ) e − ρ m h (62)Being ∆ Mij defined as it is presented at the Eq. 58 and ∆
Eij ∆ Eij = ρ j − ρ i ρ j + ρ i (63)In the last step we perform a second change of variables to x = 2 ρ m h , so that x dx =4 h ρ m dρ m . Moreover, the lower limit is transformed again to 2 ρ (cid:113) (cid:15) m ξ n c , which we define as r n g LmR ( h ) = − k B T πh ∞ (cid:88) n =0 (cid:48) (cid:90) ∞ r n x dx (∆ MLm ∆ MRm + ∆
ELm ∆ ERm ) e − x (64)Both ∆ Mij and ∆
Eij may be expressed in terms of x i = (cid:112) x + ( (cid:15) i − (cid:15) m )(2 hξ n /c ) instead of ρ i , through the straightforward substitution x i = 2 ρ i h . The Eq. 64 corresponds to the stageimmediately before the FGQA at the main part of the work, and is the result to which wewanted to arrive here.
3. Gaussian Quadrature
The gaussian quadrature is an integration method where the integrand is separated betweena well shaped function, f ( x ), and a weight function, w ( x ). The generalized N points gaussianquadrature reads (cid:90) ba f ( x ) w ( x ) dx = N (cid:88) i =1 f ( x i ) m i (65)34ith f ( x ) evaluated at the nodes x i (also called quadrature points). The knowledge of thesets of { x i } and { m i } requires the capability of solving from the j = 0 up to the 2 N − I j = (cid:90) ba x j w ( x ) dx = N (cid:88) i =1 x ji m i (66)That leads to a system of 2N equations whose solution provides the quadrature points and all { m i } . In this work we have used the one point gaussian quadrature (N=1), so the integralsthat we have to solve are I = (cid:90) ba w ( x ) dx (67) I = (cid:90) ba x w ( x ) dx (68)Thus we have m = I , and x = I /I .
4. One point gaussian quadrature applied to the FGQA
We wish to specify here the quadrature performed over g ξ n > ( h ) = − k B T πh ∞ (cid:88) n =1 (cid:90) ∞ r n x dxR ( n, x ) e − x (69)Where we have generalized for both g LmR ( h ) and ∆ g ( l, d ). In the last case just considerthat h = ( l + d ), and employ R e ( n, x ) instead of R ( n, x ). We state then f ( x ) = R ( n, x ), and w ( x ) = xe − x . These are straightforward integrals of the form I = (cid:90) ∞ r n xe − x dx = e − r n (1 + r n ) (70) I = (cid:90) ∞ r n x e − x dx = e − r n (2 + 2 r n + r n ) (71)
5. Euler - MacLaurin formula
The Euler-MacLaurin formula allows the transformation of a summatory into an integralthrough an approximation, and reads b (cid:88) n = a f ( n ) = (cid:90) ba f ( n ) dn + 12 ( f ( a ) + f ( b )) + ∞ (cid:88) k =1 B k k ! ( f (2 k − ( b ) − f (2 k − ( a )) (72)Being { B i } the Bernoulli coefficients and f m the m-th derivative of f . From now on we willapply it up to the first corrective order, k = 1.35e start again from the generalization presented in the previous section. We have after theFGQA the expression g ξ n > ( h ) = − k B T πh ∞ (cid:88) n =1 R ( n, x )(1 + r n ) e − r n (73)And now we evaluate the corrective term of the Eq. 72, considering that R ( n, x ) is essen-tially constant with n compared to the exponential dependence. This is also true for the∆ g ( l, d ) particularization, since even if the exponential in R e ( n, x ) contains a n factor,the value of ∆ (cid:15) is extremelly close to zero for large n , and the value of ( d − l ) is quitesmall as well. These features make that exponential factor very close to one∆ g ξ n > = − k B T πh (cid:20) R (1 , x )(1 + r T ) e − r T + 112 R (1 , x ) r T e − r T (cid:21) (74)Where r T = r n =1 . This term is negligible for h (cid:28)
1, so we can state simply g ξ n > ( h ) = − k B T πh (cid:90) ∞ R ( n, x )(1 + r n ) e − r n dn (75)The following step is to change the variable to ν = (4 πk B T (cid:15) / m n ) / ( c (cid:126) ). Realize that (cid:15) m isalso a function of n, so we change the variable through dν = 4 πk B Tc (cid:126) (cid:15) / m dn (cid:20) d ln (cid:15) m d ln ξ n (cid:21) (76)The term inside the brackets is defined to be j m , and its value is approximately 1. Thenonce the variable is changed we achieve g ξ n > ( h ) = − c (cid:126) π h (cid:90) ∞ ν T (cid:101) R ( ν, x )(1 + hν ) e − hν dν (77)Where (cid:101) R ( ν, x ) = (cid:15) − / m j − m R ( ν, x ).
6. One point gaussian quadrature applied to the SGQA
Beginning at the Eq. 77, we multiply by e ν/ν ∞ e − ν/ν ∞ and perform the one point gaussianquadrature approximation with f ( ν ) = (cid:101) R ( ν, x ) e ν/ν ∞ , and w ( ν ) = e − ν/ν ∞ (1 + hν ) e − hν . Wesolve then I and I by parts I = (cid:90) ∞ ν T e − ν/ν ∞ (1 + hν ) e − hν dν (78) I = ν ∞ ( ν T h + 1)( ν ∞ h + 1) + ν ∞ h ( ν ∞ h + 1) e − ν T h − νTν ∞ (79)36 = (cid:90) ∞ ν T νe − ν/ν ∞ (1 + hν ) e − hν dν (80) I = ν ∞ ( ν T h + 1)( ν ∞ h + 1) ν T + (2 ν T h + 1)( ν ∞ h + 1) ν ∞ + 2 ν ∞ h ( ν ∞ h + 1) e − ν T h − νTν ∞ (81)
7. Parameterization of the damped oscillator model for AgI
The damped oscillator model employed for the AgI describes the dielectric function at imag-inary frequencies as (cid:15) ( iξ ) = 1 + (cid:15) (0) − n UV ξ/ω IR ) + n UV −
11 + ( ξ/ω UV ) (82)Notice here that it is constructed precisely to fulfill the properties associated to any validdielectric function: positive and decreasing function, at infinite frequencies it reaches theresponse of the vacuum ( (cid:15) ( iξ → ∞ ) = 1) and at zero frequencies it reaches the value of thestatic contribution ( (cid:15) ( iξ →