Nonlinear susceptibilities of interacting polar molecules in the self-consistent field approximation P. M. Déjardin and F. Ladieu Citation: The Journal of Chemical Physics 140, 034506 (2014); doi: 10.1063/1.4855195 View online: http://dx.doi.org/10.1063/1.4855195 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/3?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Self-consistent field theory based molecular dynamics with linear system-size scaling J. Chem. Phys. 140, 134109 (2014); 10.1063/1.4869865 The multi-configuration self-consistent field method within a polarizable embedded framework J. Chem. Phys. 139, 044101 (2013); 10.1063/1.4811835 Methanol clusters (CH3OH) n , n = 3–6 in external electric fields: Density functional theory approach J. Chem. Phys. 135, 024307 (2011); 10.1063/1.3605630 An efficient, fragment-based electronic structure method for molecular systems: Self-consistent polarization with perturbative two-body exchange and dispersion J. Chem. Phys. 134, 094118 (2011); 10.1063/1.3560026 Finite-temperature infrared spectroscopy of polycyclic aromatic hydrocarbon molecules. II. Principal mode analysis and self-consistent phonons J. Chem. Phys. 133, 074303 (2010); 10.1063/1.3465554
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
THE JOURNAL OF CHEMICAL PHYSICS 140, 034506 (2014)
Nonlinear susceptibilities of interacting polar molecules in the self-consistent field approximation P. M. Déjardin1 and F. Ladieu2 1
Laboratoire de Mathématiques et Physique, Université de Perpignan Via Domitia, 52, Avenue Paul Alduy, F-66860 Perpignan Cedex, France 2 SPEC/SPHYNX(CNRS URA 2464), DSM/IRAMIS CEA Saclay, Bât. 772, F-91191 Gif-sur-Yvette, France
(Received 18 October 2013; accepted 10 December 2013; published online 16 January 2014) The nonlinear stationary response of an assembly of long range interacting electric dipoles is calculated via Berne’s nonlinear rotational diffusion equation [J. Chem. Phys. 62 1154 (1975)]. Analytical formulas are derived, showing that the behavior of ω and 3ω components of the nonlinear external field response spectra in such polar dielectrics strongly deviates from the Coffey-Paranjape formulas [Proc. R. Ir. Acad., Sect. A 78, 17 (1978)] as the long range dipole-dipole interactions increase, while the linear response remains qualitatively unaffected. By qualitatively comparing with recent experimental measurements on glycerol, it is further demonstrated that nonlinear dielectric response experiments provide a powerful tool for the characterization of short range intermolecular interactions in polar liquids. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4855195] I. INTRODUCTION
A system in thermal equilibrium at temperature T disturbed by an external stimulus evolves towards a new stationary state. Furthermore, if the energy stimulus is much lower than the thermal energy, deviations linear in the stimulus of the expectation value of the relevant dynamical variable are sufficient to calculate the generalized susceptibility using appropriate equilibrium correlation functions and linear response theory. The simplest illustration is the orientational electric polarization of noninteracting permanent dipoles in an alternative sinusoidal field E(t) treated by Debye.1 The calculation of the stationary nonlinear response is, however, much more involved because, for example, no connection between the transient and ac response exists. Nonlinear dielectric relaxation and the dynamic Kerr effect of permanent dipoles are naturally occurring examples. Now the first calculation concerning nonlinear dielectric relaxation of noninteracting dipoles was performed by Coffey and Paranjape.2 In effect, with ϑ the polar angle of the dipole moment vector μ, these authors show that in particular, amending cos ϑ(t) to terms of order E3 , the nonlinear response consists of a term oscillating at the fundamental superimposed on which is a term oscillating at the third harmonic. This prediction was experimentally verified 20 years later3 in a system for which the macromolecules have been considered noninteracting, putting in evidence the so-called nonlinear dynamical dielectric effect. All subsequent calculations existing in the literature4–6 systematically ignore the interactions between the dipoles, electric or magnetic, while the obtaining of a simple formula for the nonlinear response using the Debye-Fröhlich model, which includes intermolecular interactions via a double-well potential, is by no means easy.7, 8 Dipolar interactions in the orientational dynamics of an assembly of dipoles always lead to involved calculations as is apparent in dielectric relaxation of polar fluids.9, 10 For ex0021-9606/2014/140(3)/034506/13/$30.00
ample, the first Debye-Lorentz calculation of the static permittivity valid for polar gases at very low densities was improved upon by Onsager,10 and, as outlined by van Vleck,11 can be transposed to assemblies of Langevin paramagnets. In essence, Onsager, in order to improve Debye’s formula, accounted for the reaction of a given polar molecule on its surroundings by supposing that a typical dipole of an assembly resided at the center of a macroscopic spherical cavity in a dielectric medium, ultimately generalizing the Debye calculation of the static dielectric constant to higher densities. However, Onsager’s model was criticized by Kirkwood12 because it ignores the interactions between the dipole of the macroscopic cavity with its nearest neighbours. Consequently, Kirkwood, and later Fröhlich, have corrected Onsager’s formula by introducing the “Kirkwood correlation factor,” which accounts for the contribution of the neighbours of a molecule to the dielectric permittivity.9, 10 In spite of the largely successful static calculation, the determination of the dynamic permittivity of a dilute assembly of dipoles is far more involved than its static counterpart because, in general, the equation of motion of a typical dipole depends on the macroscopic properties that one is trying to calculate.9 In effect, this task leads to extremely difficult calculations because the time dependence of the field involved in the equation of motion of a typical dipole is unknown. Various attempts to generalize Onsager’s approach have been given leading to permittivity values (e.g., the Onsager-Cole formula) incompatible with any relaxation time distribution as the corresponding Cole-Cole plots may lie outside of the Debye semicircle.9 One approach to the problem is that of Nee and Zwanzig.13 By assuming a delayed response between the dipole moment and the Onsager reaction field, hence allowing for a special type of memory effect (dielectric friction), they were able to rederive the Fatuzzo-Mason equation9 for the dielectric permittivity ε(ω). A review of attempts to generalize the Kirkwood-Fröhlich theory to the dynamical case was given by Madden and Kivelson.14 Furthermore, a more
140, 034506-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-2
P. M. Déjardin and F. Ladieu
general extension to intermediate to large wavelengths was given by Bagchi and Chandra15 by adapting the equations of generalized molecular hydrodynamics to dielectric relaxation. Later, Zhou, and Bagchi16 have reconsidered Zwanzig’s model17 of dielectric relaxation of dipoles located at the sites of a cubic lattice. In particular, these authors pointed out the inadequacy of both the perturbation and the Nee-Zwanzig approaches to the collective dynamics, showing that these models for dielectric friction predict a too rapid decay of the dipole-dipole correlation function. Nevertheless, they still find Debye-like behavior for the permittivity for relatively dilute dipolar systems, in agreement with both the perturbation and continuum approaches. Yet another model that has been introduced for linear dielectric relaxation is Berne’s forced rotational diffusion model.18 Inspired by the literature on dielectric relaxation in electrolytes,19 Berne treated the rotational Brownian motion of interacting electric dipoles self-consistently so that the Fokker-Planck equation (FPE) describing the orientational dynamics of the polar molecules is nonlinear. The difficulty in handling this FPE is mainly a mathematical one, which Berne overcame by linearizing it about statistical equilibrium before taking its spatial Fourier transform. Then he took the zero wave-vector limit of the resulting equation in order to obtain the relevant correlation functions for the calculation of the dielectric constant of the dipolar assembly. This model in its nonlinear version has been also considered by Warchol and Vaughan.20 Berne’s model was furthermore recently extended to magnetic relaxation of single-domain ferromagnetic particles weakly coupled by magnetic dipole-dipole interactions.21 All the aforementioned works pertain to the static properties and to the linear response of permanent dipoles to alternating fields. Nonlinear response in interacting systems of electric dipoles is more involved, and calculations involving this dynamic response are hardly available.22 The purpose of this paper is therefore to derive formulas for the (lowfrequency) nonlinear response of electric dipoles to alternating electric fields using Berne’s model in its nonlinear version, amending terms up to O(E3 ).
J. Chem. Phys. 140, 034506 (2014)
W (ϑ, ϕ, t), viz.,
1 ∂ ∂V ∂W ∂ sin ϑ βW + 2τD W = ∂t sin ϑ ∂ϑ ∂ϑ ∂ϑ 1 ∂ ∂V ∂W + 2 βW + , ∂ϕ ∂ϕ ∂ϕ sin ϑ
where τ D = βζ /2, ζ is the rotational viscous coefficient, β = (kT)−1 , k is Boltzmann’s constant, T the absolute temperature, and V is a potential made of an effective single dipole free energy Us and a term Vint holding for the effect of interaction torques, viz., V (ϑ, ϕ, t) = Us (ϑ, ϕ, t) + Vint (ϑ, ϕ, t).
We consider in this paragraph independent groups of interacting molecules, each group being represented by a single dipole in a self-consistent potential. In each group there is no periodicity (neither in the positions nor in the orientations) and the positions of the molecules are drawn at random with an average distance a between the molecules. These groups are considered in equilibrium with a thermostat, so that the ensemble of independent groups constitute a canonical ensemble on which the methods of statistical mechanics may be applied. Thus, we may start our investigation with the Smoluchowski rotational diffusion equation for the distribution function of dipole orientations on the unit sphere
(2)
The term Us in Eq. (2) contains the orientational terms due to the externally applied fields, while Vint in Eq. (2) mimics the effects of long range Maxwell interaction fields. It has been shown21 that jk Vint (ϑ, ϕ, t) = UJ KLM (rj k , ϑj k , ϕj k ) J KLM j,k k=j
∗ × YJ K (ϑ, ϕ) YLM (ϑ, ϕ) (t),
(3)
where YJK is a spherical harmonic of order J and rank K, the star stands for the complex conjugate, the angular brackets denote an ensemble average over W, namely, 2π π YJ K (ϑ, ϕ)(t) =
YJ K (ϑ, ϕ)W (ϑ, ϕ, t) sin ϑdϑdϕ, 0 0
(4) (rjk , ϑ jk , ϕ jk ) is the set of spherical polar coordinates specifying the modulus and orientation of the vector rk − rj , rj being the location of dipole number j in the laboratory frame, the jk coefficients UJ KLM are given by21 jk
UJ KLM (rj k , ϑj k , ϕj k ) = (4π )4 i J −L+P IJ LP (rj k )KJ KLMP Q YP Q (ϑj k , ϕj k ), (5) P ,Q
where
KJ KLMP Q =
II. BASIC EQUATIONS AND EQUILIBRIUM SOLUTION OF THE SELF-CONSISTENT ROTATIONAL DIFFUSION MODEL
(1)
(2J + 1)(2P + 1) 4π (2L + 1)
× P 0J 0 |L0P QJ K |LM, P 0J 0 |L0 coefficients,
and
P QJ K |LM
are
(6)
Clebsch-Gordan
IJ LP (rj k ) =
[(L + J + P + 1)/2] +1 (J + 3/2) (L + 3/2) [1 − (L + J − P )/2] 8rjL+J k × dsZ(s)s J ds Z(s )s L π 3/2
×F4
J +L−P J +L+P +1 3 3 s 2 s 2 , ; J + , L+ ; 2 , 2 , 2 2 2 2 rj k rj k (7)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-3
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
and is the Euler gamma function,23 F4 is the two-variable hypergeometric series23 ∞ ∞ (γ )m+n (δ)m+n m n x y , F4 (γ , δ; ε, ζ ; x, y) = (ε) m (ζ )m m!n! m=0 n=0
(8)
Z(s) is the linear charge density, and the Pochhammer symbol is given by (a)n = (a + n)/ (a). We note in passing that the hypergeometric series in Eq. (7) converges only if s + s < rjk , which means that the theory has a meaning only if the mean interparticle distance is larger than the particle diameter, or, put differently, if the intermolecular distance is larger than the characteristic size of a typical molecule. Now, because of Eqs. (3) and (4), it is clearly seen that the Fokker-Planck Eq. (1) is nonlinear because V is a linear functional of W. Before starting to solve this equation, however, one must make some simplifying assumptions regarding the properties of the linear charge density Z(s). From now on, we assume that Z(s)ds = 0, sZ(s)ds = μ, (9) where μ is the dipole moment modulus, and s L Z(s)ds = 0, L > 1,
(10)
so that Appell’s hypergeometric series F4 in Eq. (7) may be replaced by unity. If an electric field is applied in some direction, this will act on the dipolar field in the same direction so that we may drop transverse components of the dipole field and consider only the longitudinal one. Hence, we assume that V = V (ϑ, t), so that the longitudinal and transverse components of the dipole are disentangled. Therefore, we may neglect the azimuthal dependence of the potential energy and also that of the distribution function. Thus, Eq. (1) becomes 1 ∂ ∂V ∂W ∂W = sin ϑ βW + . (11) 2τD ∂t sin ϑ ∂ϑ ∂ϑ ∂ϑ
tions are equivalent. In other words, it is thus implicitly assumed that the dielectric has a spherical shape. Now, on expanding the distribution function W as a series of Legendre polynomials, viz., W (ϑ, t) =
[ξ (t) − λf1 (t)] 2τD f˙n (t) + fn (t) = [fn−1 (t) − fn+1 (t)]. n(n + 1) (2n + 1) (17) This hierarchy is nonlinear in the expansion coefficients fn (t) and reflects the nonlinear integro-differential character of Eqs. (1) and (11). Because of this nonlinearity, one may solve Eq. (17) by numerical methods only. Fortunately, we only require the nonlinear response to a field such that ξ 0 1. Hence, we may use perturbation methods in order to obtain the various nonlinear responses. Now, we require the nonlinear response up to ξ03 , so that we set fn (t) = fn(0) + fn(1) (t) + fn(2) (t) + fn(3) (t) + · · · ,
P (t) = ρ0 μf1 (t). Next we combine Eqs. (17) and (18) to obtain the perturbation equations for n > 0 (f0(0) = 1), fn(0) =
λf1(0) (0) (0) fn+1 − fn−1 , 2n + 1
λf1(0) (1) (1) fn+1 (t) − fn−1 (t) , 2n + 1 (20)
2τD f˙(2) (t) + fn(2) (t) n(n + 1) n =
(13)
where ξ (t) = βμE0 cos ωt = ξ 0 cos ωt and λ = 4π βρ 0 μ2 /3. One may equivalently consider that the total field seen by a typical dipole is effectively F (t) = E(t) −
Hence by assuming that the field inside the system of dipoles is given by Eq. (14), one also assumes that all spatial direc-
(19)
ξ (t) − λf1(1) (t) (0) 2τD (0) f˙n(1) (t) + fn(1) (t) = fn−1 − fn+1 n(n + 1) 2n + 1
4πρ0 μ2 cos ϑcos ϑ(t), (12) 3 where ρ 0 is the number of molecules per unit volume. Next, we assume that the dipoles are subjected to an alternating field E(t) = E0 cos ωt, so that the potential energy reads
(15)
(18)
the superscript (i) meaning the order in the field strength, the polarization being given by Eq. (15) which we rewrite as follows:
+
4π P (t), (14) 3 where the polarization along the direction of the field P(t) is given by
(16)
and using this expansion in conjunction with Eqs. (11) and (13), we arrive at the infinite hierarchy of differential recurrence relations for the expansion coefficients fn (t), namely,
Vint (ϑ, t) =
P (t) = ρ0 μcos ϑ(t).
(n + 1/2)fn (t)Pn (cos ϑ)
m=0
We also have from Eqs. (3)–(10),18, 20, 21
βV (ϑ, t) = −ξ (t) cos ϑ + λcos ϑ(t) cos ϑ,
∞
ξ (t) − λf1(1) (t) (1) (1) fn−1 (t) − fn+1 (t) (2n + 1) +
λf1(2) (t) (0) (0) f − fn−1 (2n + 1) n+1
+
λf1(0) (2) (2) f (t) − fn−1 (t) , (2n + 1) n+1
(21)
and 2τD f˙(3) (t) + fn(3) (t) n(n + 1) n =
ξ (t) − λf1(1) (t) (2) (2) fn−1 (t) − fn+1 (t) (2n + 1)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-4
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
+
λf1(3) (t) (0) (0) fn+1 − fn−1 (2n + 1)
+
λf1(2) (t) (1) (1) fn+1 (t) − fn−1 (t) (2n + 1)
which steady-state solution yields the linear dynamical polarization P1 (t) as
P1 (t) = ρ0 μf1(1) (t) = ρ0 E0 α1(1) (ω) cos ωt + α1(1) (ω) sin ωt ,
+
λf1(0) (2n + 1)
(3) (3) fn+1 (t) − fn−1 (t) .
(30) (22)
Since in Eq. (19), the functions fn(0) are not driven by the field, these functions coincide with the moments of the equilibrium distribution function W0 (ϑ), namely, fn(0) = Pn (cos ϑ)0 ,
π Pn (cos ϑ)W0 (ϑ) sin ϑdϑ,
(24)
0
and the equilibrium distribution function W0 (ϑ) is given by W0 (ϑ) =
1 e−λcos ϑ0 cos ϑ , Z(cos ϑ0 )
α1(1) (ω) =
(31) and where τ1 =
Z(cos ϑ0 ) =
(25)
(32)
2τD f˙(2) (t) + fn(2) (t) n(n + 1) n =
e−λcos ϑ0 cos ϑ sin ϑdϑ.
ξ (t) − λf1(1) (t) (1) (1) fn−1 (t) − fn+1 (t) δn,2 , (2n + 1)
so that
0
From Eq. (25), it appears that the equilibrium moments (24) are completely determined by the specification of f1(0) = cos ϑ0 . Instead of evaluating Eq. (24), we alternatively use Eq. (19) with n = 1. This leads to the solution f1(0)
τD . 1 + λ/3
Thus, a Debye spectrum with shorter relaxation time is obtained for the complex linear polarizability, as may be expected from such a model. Now we may proceed by calculating f1(3) (t). On writing down Eq. (21) and accounting for the aforementioned simplifications, we have
with π
βμ2 βμ2 ωτ1 , α1(1) (ω) = 2 2 (3 + λ) 1 + ω τ1 (3 + λ) 1 + ω2 τ12
(23)
where the angular brackets denote an equilibrium ensemble average given by Pn (cos ϑ)0 =
where the in-phase and out of phase “polarizabilities” α1(1) (ω) and α1(1) (ω) are given by
= 0,
which means that dipole-dipole interactions cannot spontaneously orient the dipolar assembly due to their frustrating character. Hence, from Eq. (19), for n > 0, fn(0) = 0.
(26)
We now may proceed to the evaluation of the linear and nonlinear polarizabilities.
fn(2) (t) = 0, n = 2,
(33)
and 3 (2) 3ξ (t) (1) 3λ (1) 2 f˙2(2) (t) + f2 (t) = f1 (t) − f (t) . (34) τD 5τD 5τD 1 This equation may be integrated by quadratures with f2(2) (−∞) = 0 to yield f2(2) (t)
3ξ0 = 5τD
t
e−3(t−t )/τD f1(1) (t ) cos ωt dt −
−∞
t ×
2 e−3(t−t )/τD f1(1) t dt ,
3λ 5τD
(35)
−∞
III. CALCULATION OF THE LINEAR AND NONLINEAR RESPONSES TO THE APPLIED FIELD
Here we give the solution of Eqs. (20)–(22) which will allow us to calculate the linear and nonlinear polarizabilities. First, Eq. (20) reads 2τD f˙(1) (t) + fn(1) (t) = 0, n > 1, n(n + 1) n
(27)
meaning that for n > 1, we have fn(1) (t) = 0. For n = 1, we have the differential equation: λ ξ (t) (1) ˙ τD f1 (t) + 1 + f1(1) (t) = , 3 3
(28)
(29)
which pertains, e.g., to the dynamic Kerr effect. As we do not focus our attention on the dynamic birefringence, we do not explicitly give the result of the integrals here. For λ = 0, the result of the calculation agrees with that of Benoit24 and that of Coffey and Paranjape.2 Owing to Eqs. (26) and (28) and (22) and (33) reads 2τD f˙(3) (t) + fn(3) (t) n(n + 1) n ξ (t) − λf1(1) (t) (2) (2) fn−1 (t) − fn+1 (t) = (2n + 1) λf1(3) (t) (0) (0) f + − fn−1 δn,1 , (2n + 1) n+1
(36)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-5
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
so that for n = 1 we obtain the equation λ λf (1) (t) − ξ (t) (2) (3) ˙ f1(3) (t) = 1 f2 (t), τD f1 (t) + 1 + 3 3 (37) which again may be integrated by quadratures with f1(3) (−∞) = 0. The result of the calculation yields the
nonlinear polarization as
P3 (t) = ρ0 E03 α3(1) (ω) cos ωt + α3(1) (ω) sin ωt + α3(3) (ω) cos 3ωt + α3(3) (ω) sin 3ωt ,
where the in-phase and out-of-phase components of the nonlinear polarizabilities α3(s) (ω) = α3(s) (ω) − iα3(s) (ω), s = 1, 3 are provided by
3β 3 μ4 (3 + λ)2 (4λ − 39) ω4 τ14 + 378 + 522λ + 51λ2 ω2 τ12 + 729 =− ,
3 20 (3 + λ)4 81 + 4 (3 + λ)2 ω2 τ 2 1 + ω2 τ 2
(39)
β 3 μ4 ωτ1 81 (2λ − 21) + 3 (3 + λ) (λ (λ − 45) − 198) − (3 + λ)2 ω2 τ12 ω2 τ12 = , 3 10 (3 + λ)4 81 + 4 (3 + λ)2 ω2 τ12 1 + ω2 τ12
(40)
α3(1) (ω)
1
α3(1) (ω)
(38)
1
9β 3 μ4 (3 + λ)2 (51 + 20λ) ω6 τ16 + 3 (279 − 4λ (λ (6 + λ) − 21)) ω4 τ14 − 3 (λ (78 + λ) − 99) ω2 τ12 − 81 = , (41)
3 20 (3 + λ)4 81 + 4 (3 + λ)2 ω2 τ12 1 + ω2 τ12 1 + 9ω2 τ12
9β 3 μ4 ωτ1 9 (2λ − 21) + λ3 − 243λ − 297 ω2 τ12 − 3 (3 + λ) (3 + λ (19 + 4λ)) − (3 + λ)2 ω2 τ12 ω4 τ14 (3) α3 (ω) = . (42)
3 10 (3 + λ)4 81 + 4 (3 + λ)2 ω2 τ12 1 + ω2 τ12 1 + 9ω2 τ12 α3(3) (ω)
For λ = 0, these formulas reduce to the Coffey-Paranjape formulas2 for the nonlinear response of non-interacting dipoles to alternating electric fields, viz., α3(1) (ω) =
13ω2 τD2 − 27 β 3 μ4 , 180 1 + ω2 τ 2 2 9 + 4ω2 τ 2 D
α3(1) (ω) = −
(43)
D
β 3 μ4 ωτD 21 + ω2 τD2 2 , 9 + 4ω2 τD2 90 1 + ω2 τD2
(44)
17ω2 τD2 − 3 β 3 μ4 , = 60 1 + ω2 τD2 9 + 4ω2 τD2 1 + 9ω2 τD2 (45) β 3 μ4 ωτD 3ω2 τD2 − 7 , α3(3) (ω) = 30 1 + ω2 τD2 9 + 4ω2 τD2 1 + 9ω2 τD2 (46) where we have used τ 1 = τ D for λ = 0. Thus, it appears that in order to reproduce the Coffey-Paranjape results, the formal limit λ → 0 is taken in the expressions for the complex polarizabilities, and this limit refers to the very diluted limit, where intermolecular interactions are negligible. Our equations (39)–(42) also formally agree with the results of Felderhof and Jones,22 who considered nonlinear response of magnetic dipoles, provided their interaction parameter C is replaced by −λ/3. We emphasize that the nonlinear susceptibilities provided by Eqs. (39)–(42) are (mathematically) valid for an arbitrary positive value of λ, and constitute the central result of our work. By inspection of Eqs. (39)–(42), we remark that the time scales of the nonlinear response are α3(3) (ω)
τ1 =
τD , 1 + λ/3
τ2 =
2 (3 + λ) τ1 , 9
τ3 = 3τ1 .
If τ 1 represents the fundamental time scale, then the time scales of the nonlinear response occurring in Eqs. (39)–(42) behave in the same fashion as in the Coffey-Paranjape formulas. However, the frequency behavior of Eqs. (39)–(42) is more complicated. In order to complete our calculation, we need now to introduce Maxwell field corrections to the nonlinear response.
(47)
IV. MAXWELL FIELD CORRECTIONS
As alluded to above, we need to account for Maxwell field corrections to the nonlinear susceptibilities in order to obtain quantities which can be compared with experiment. Here, we have in mind the experimental work recently published by Brun et al.25–29 on glycerol, where in particular, the authors seek to obtain universal susceptibility curves by representing the results as a function of ωτ , where τ is a certain relaxation time. In order to accomplish this, we may proceed in two ways. In the following, the Roman superscripts I and II will refer to the two ways. We first start with the linear response.
A. Determination of the linear response
First way: here we account for static effects of the Maxwell field only. Then we remark that the in-phase and out-of-phase molecular responses to E0 are given by Eq. (31), which means that the complex linear response per molecule is given by α1(1) (ω) = α1(1) (ω) − iα1(1) (ω). The linear polarization frequency response P˜1(1) (ω) is therefore given by P˜1(1) (ω) = ρ0 α1(1) (ω)E0 . Since we have assumed that the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-6
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
molecules have only a permanent dipole moment, the static polarizability of an isolated molecule is therefore given by βμ2 /3. We therefore may write the linear polarization response as E0 1 λ , (48) P˜1(1) (ω) = 4π 1 + iωτ1 1 + λ/3 which means that the Maxwell field amplitude El is given by El = (1 + λ/3)−1 E0 . If we write P˜1(1) (ω) = χ1(1) (ω)El then, the linear dielectric constant ε1(1) (0) is given by ε1(1) (0) − 1 = 4π χ1(1) (0) = λ,
(49)
which means that the Maxwell field amplitude is given by El = 3[2 + ε1(1) (0)]−1 E0 , e.g., by the Lorentz static formula. The susceptibility spectrum, given by χ1(1),I (ω) =
1 λ 4π 1 + iωτ1
(50)
is then Lorentzian with the relaxation time τ 1 . Second way: Here, we account for the dynamical effects of the Maxwell field on the linear response. In order to accomplish this, we follow the spirit of the work of Madden and Kivelson.14 As our calculations are valid for a dielectric sphere only, the linear polarization frequency response is given in terms of the linear complex permittivity by 3 ε1(1),I I (ω) − 1 P˜1(1) (ω) = E0 . 4π ε1(1),I I (ω) + 2
(51)
This expression should be equal to the linear polarization spectrum expressed in terms of molecular quantities, P˜1(1) (ω) = ρ0 α1(1) (ω)E0 . By accomplishing this we obtain the following expression for the dynamical susceptibility, namely, 4π χ1(1),I I (ω) = ε1(1),I I (ω) − 1 =
ε(1) (0) − 1 λ = 1 , 1 + iωτD 1 + iωτD (52)
namely, again a Lorentzian spectrum, with now a relaxation time τ D . Here, we note in passing that Jadzyn et al.30 have shown that the relaxation time of the linear dielectric response is the same in two nematogenic compounds, where the molecular moment differs by a factor as large as 2: even if our calculations do not apply to liquid crystals, this experimental result qualitatively supports our previous finding, namely, the fact that the linear response does not “feel” qualitatively the existence of long range dipolar interactions. Now, since the qualitative frequency behavior of Eqs. (50) and (52) is the same, so that we have no clear indication regarding the way to proceed in order to account for the Maxwell field in the determination of the nonlinear response. This is why we now turn to the calculation of the nonlinear susceptibilities using the two ways in section B of this paragraph. B. Determination of the nonlinear response
First way: Here we account for the static Maxwell field effects only. Thus we first write the nonlinear steady polariza-
tion response as follows: P3 (t) =
3El3 (1),I χ3 (ω) cos ωt + χ3(1),I (ω) sin ωt 4 E3 + l χ3(3),I (ω) cos 3ωt + χ3(3),I (ω) sin 3ωt , (53) 4
where again El = (1 + λ/3)−1 E0 . Next, we identify Eq. (53) with Eq. (38). We have λ 3 ρ0 (s) (s),I α3 (ω), s = 1, 3, (54) χ3 (ω) = 1 + 3 n(i) 3 (3) where n(1) 3 = 3/4 and n3 = 1/4. Second way: Here we account on the dynamical Maxwell field effects on the nonlinear response. Because the derivations are more involved than for the static field corrections, we only give the results and leave the details of calculation to Appendix A. The expressions for the nonlinear susceptibilities are provided by the following equations, viz. 2 ε1(1),I I (ω)+2 4ρ0 ε1(1)∗,I I (ω)+2 (1),I I (ω) = α3(1) (ω), χ3 3 3 3 (55)
χ3(3),I I (ω)
= 4ρ0
ε1(1),I I (ω) + 2 3
3 α3(3) (ω),
(56)
(s) where α3(s) (ω) = α(s) 3 (ω) − iα3 (ω) and the in-phase and out-of-phase parts of this expression are given by Eqs. (39)– (42). One notes that if ω is set equal to zero in the complex permittivity factors, one indeed recovers the static results Eq. (54). Hence, accounting for dynamical effects of the Maxwell field, the nonlinear susceptibilities χ3(1) (ω) and χ3(3) (ω) are given by Eqs. (55) and (56), respectively, where the molecular nonlinear responses to E0 are corrected by the dynamical Lorentz cavity field factor. The quantities of experimental interest are, for the two ways: (s),j ρ0 (s),j X3 (ω) = (1) 2 χ3 (ω), s = 1, 3 ; j = I, I I. β χ1 (0) (57) The quantities (57) are dimensionless and render a static value of 0.2 in the interaction-free limit, so facilitating the comparison with the work of Brun et al.25–29 Notice also that for both ways, the nonlinear susceptibilities agree with the CoffeyParanjape ones if the Maxwell field factors are replaced by unity and the nonlinear polarizabilities by Eqs. (43)–(46).
V. RESULTS AND DISCUSSION
The expressions we have obtained for the nonlinear response, Eqs. (39)–(42) are easily computed and therefore numerical analysis of them is straightforward. Furthermore, as alluded to in the text, they can be used for any value of the interaction parameter λ, which can in practice be considered as a fitting parameter with reasonable values (as demonstrated by Eq. (49) λ is of the order of the linear static dielectric constant ε1(1) (0)). We first focus onto the frequency dependence of the (s) (s) nonlinear response by studying A(s) 3 (ω) = |α3 (ω)/α3 (0)|,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-7
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
FIG. 3. Dimensionless modulus of the ω nonlinear susceptibility component (static Maxwell field corrections) as a function of ωτ 1 for various values of λ. The value λ = 0 is the Coffey-Paranjape result. FIG. 1. Normalized modulus of the ω nonlinear polarizability component as a function of ωτ 1 for various values of λ. The value λ = 0 is the CoffeyParanjape result.
(s),j
then we discuss the behavior of X3 (ω) which are the quantities to be compared with experimental data. Figures 1 and 2 illustrate the behavior of the ω and 3ω normalized nonlinear responses to the externally applied field E0 (t), without Maxwell field corrections, where the quantitiesA(s) 3 (ω) are plotted against ωτ 1 . Note that the (possibly strong) temperature dependence of the rotational friction coefficient is totally absorbed by the chosen ωτ 1 horizontal scale. Figures 1 and 2 allow to see how A(s) 3 (ω) deviate from the Coffey-Paranjape formulas. In effect, the nonlinear response to the applied ac field exhibits a humped shape as λ increases, and therefore our results clearly demonstrate that the nonlinear response is affected by intermolecular interactions. This contrasts with the linear response which in our model keeps a Debye shape for all the values of λ. Now, Figures 3–6 show the frequency behavior of the modulus of the dimensionless nonlinear susceptibilities (s),j X3 (ω), where Maxwell field corrections are accounted for. Therefore, these quantities can be directly compared with experimental data. Furthermore, these quantities are of di-
FIG. 2. Normalized modulus of the 3ω nonlinear polarizability component as a function of ωτ 1 for various values of λ. The value λ = 0 is the CoffeyParanjape result.
rect interest with respect to experiments such as those of Refs. 25–29 in which one of us is directly involved. Figures 3 and 4 show the behavior of the nonlinear responses X3(s),I (ω), e.g., when only static effects due to the Maxwell field are accounted for. One may see that the static values of the nonlinear responses decrease, and a humped shape develops while λ increases, the peak value being lower than the value for noninteracting dipoles, while the high-frequency behavior is the same for all λ values. Figures 5 and 6 show the behavior of the nonlinear responses X3(s),I I (ω), when dynamical effects due to the Maxwell field are accounted for. Indeed the behavior of the static value is the same as those for X3(s),I (ω), while the humped shape has now disappeared completely whatever the λ values. At all frequencies, the value of the nonlinear susceptibilities are all less than the those provided by the nonlinear susceptibilities of the ideal gas and in contrary with X3(s),I (ω) the high frequency values depend on λ. In both cases I and (s),j II, the values of the nonlinear susceptibilities X3 (ω) are less than the ones returned by the Coffey-Paranjape formulas. Thus, we may infer that if the experimental values are in excess with respect to those returned by the Coffey-Paranjape formulas, then this excess is due to short range intermolecular interactions. Comparison of our results with previous ones can only be achieved in the context of linear response. First, we
FIG. 4. Dimensionless modulus of the 3ω nonlinear susceptibility component (static Maxwell field corrections) as a function of ωτ 1 for various values of λ. The value λ = 0 is the Coffey-Paranjape result.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-8
P. M. Déjardin and F. Ladieu
FIG. 5. Dimensionless modulus of the ω nonlinear susceptibility component (dynamic Maxwell field corrections) as a function of ωτ 1 for various values of λ. The value λ = 0 is the Coffey-Paranjape result.
should compare our results with those obtained by Berne.18 As demonstrated in Appendix B, Berne’s results are valid for weak densities, and for a sample which shape is similar with an infinite slab. In effect, Berne (and also Warchol and Vaughan for the nonlinear version20 ) uses the field F(t) = E(t) − 4π PZ (t)eZ , which is the macroscopic field in a uniformly polarized infinite slab and which differs from that in Eq. (14). The latter is isotropic and for a single microscopic relaxation time, we have a single macroscopic relaxation time. In contrast, Berne uses a priori an anisotropic field, which is compatible with the Fatuzzo-Mason theory for purely polar media, while we use the Kubo formula,9 which is compatible with Eq. (14). Nevertheless, as we demonstrate in Appendix B, the linear complex permittivity (e.g., the macroscopic response) for both models is the same. Now, strictly speaking, Berne’s calculation is valid only for weak dipolar concentrations, as we demonstrate in Appendix B. The calculation presented there clearly demonstrates that (a) Onsager’s equation for polar liquids is recovered at weak densities and (b) the Fatuzzo-Mason theory for purely polar liquids renders the Debye formula for the complex permittivity, just as in our work. In effect, the range of validity of Berne’s calculation may fairly well be extended, and that both Berne’s and Onsager’s formulations are valid at weak densities only, left apart short range interactions. In contrast, the (extended) Berne’s theory may now have a range of validity from an experimental point
FIG. 6. Dimensionless modulus of the 3ω nonlinear susceptibility component (dynamic Maxwell field corrections) as a function of ωτ 1 for various values of λ. The value λ = 0 is the Coffey-Paranjape result.
J. Chem. Phys. 140, 034506 (2014)
of view, which is probably the same as ours (see the next paragraph) in its present stage, simply because when the concentration increases (or the temperature decreases), short range intermolecular interactions have to come into play for a realistic description of the statics and the dynamics of the system. The comparison of our results with those of Chandra and Bagchi31 can also be made. For linear response, we obviously obtain the same answer in the zero wave vector limit and very weak densities, as it must be. However, we have to notice that the calculation of the nonlinear response is not possible from the Chandra-Bagchi Fokker-Planck equation (FPE). Here, we have to stress that Berne’s FPE used by us, and that of Chandra and Bagchi are similar, but not identical, because they are not derived from the same starting point and the drift coefficients are different in these equations. Therefore, even if there are some similarities, they definitely produce different quantitative results. In particular, and this is reckoned by Chandra and Bagchi themselves,31 the direct correlation function in their formulation is time-independent, which means that the pair correlation function and the pair distribution function are also time-independent, therefore not allowing one to study interacting systems far from equilibrium. Furthermore, the identification of Berne’s drift coefficients with theirs is only formal. Also, strictly speaking, the Bagchi-Chandra nonlinear FPE is valid for diluted polar systems and this is also reckoned by them.31 Now, it is very difficult to work out the nonlinear response from their equation and compare this eventually feasible calculation with our present results because in order to achieve such a task, the knowledge of the time dependence of the direct correlation function is required. In effect, this task goes beyond the Ornstein-Zernike formalism which is established for systems in equilibrium only. We remark in passing that it is in principle possible to calculate the wave vectordependent complex nonlinear polarizabilities from our model. This calculation is, however, far more involved than what we present in our work and is also beyond our scope. (s) (ω) Now, in Refs. 25–29 the measured values of X3,gly in supercooled glycerol close to its glass temperature transition temperature Tg were contrasted with the “trivial” values (s) (ω) obtained by Coffey and Paranjape2 for noninteractX3,triv (s) (ω) has a humped shape in freing dipoles. The fact that X3,gly quency with an amplitude which depends significantly on the (s) (ω) does temperature T was contrasted with the fact that X3,triv not depend on temperature and is monotonous in frequency. The conclusion of Refs. 25–29 was that the temperature de(2) (ω) gave some information pendence of the hump of X3,gly about the glassy correlations—as predicted in32 —i.e., about the temperature dependence of the number Ncorr of dynamically correlated molecules in supercooled glycerol. In this context, the model developed in the present paper gives a re(s) (ω) since it goes beyond the ideal fined calculation of X3,triv case of infinitely diluted dipoles and takes into account the finite concentration effects of the dipoles and the resulting finite long range dipole-dipole interactions. This is why we shall now contrast the calculation displayed in Figures 3–6 with the experimental data of Refs. 25–28. In case II (Figs. 5 and 6), X3(s),I I (ω) has the same qualitative frequency dependence as the Coffey-Paranjape formulae, and X3(s),I I (ω) is always decreasing when λ
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-9
P. M. Déjardin and F. Ladieu
increases. As a result, in case II, the conclusions of Refs. 24–28 are qualitatively reinforced by the present refined (s) (ω). calculation of X3,triv In the case I (Figs. 3 and 4), a more detailed discussion must be done: (i) The values computed in the present paper have a humped shape in frequency: this hump is absent at λ = 0, where our Eqs. (39)–(42) reduce to the CoffeyParanjape result. This is qualitatively important since this seems to indicate that a model without glassy correlations, such as the present one, can yield a humped frequency behavior (as recently found by Diezemann33 ). However, we emphasize that, even if our calculation is technically valid for any value of λ, our model cannot capture the “dipolar glassy phase” that one expects when λ 1. This comes from the mean field approximation used in our derivation, where we assume that the environment is strictly the same for all the dipoles. We thus cannot hope to describe accurately the low temperature phase where the dipole-dipole interaction supersedes the thermal energy. Little is known about this dipolar glassy phase. To the best of our knowledge, the physical situation which is the closest to our calculation is that of concentrated ferrofluids where single domain anisotropic magnetic nanoparticles interact through dipolar couplings. In that case, a transition, at low temperature, towards a superspin glass phase has been evidenced and the cubic susceptibility is singular,34 revealing the amorphous long range order which sets in at this transition. This is why, coming back to Figures 3 and 4 and focusing on the range 0
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
THE JOURNAL OF CHEMICAL PHYSICS 140, 034506 (2014)
Nonlinear susceptibilities of interacting polar molecules in the self-consistent field approximation P. M. Déjardin1 and F. Ladieu2 1
Laboratoire de Mathématiques et Physique, Université de Perpignan Via Domitia, 52, Avenue Paul Alduy, F-66860 Perpignan Cedex, France 2 SPEC/SPHYNX(CNRS URA 2464), DSM/IRAMIS CEA Saclay, Bât. 772, F-91191 Gif-sur-Yvette, France
(Received 18 October 2013; accepted 10 December 2013; published online 16 January 2014) The nonlinear stationary response of an assembly of long range interacting electric dipoles is calculated via Berne’s nonlinear rotational diffusion equation [J. Chem. Phys. 62 1154 (1975)]. Analytical formulas are derived, showing that the behavior of ω and 3ω components of the nonlinear external field response spectra in such polar dielectrics strongly deviates from the Coffey-Paranjape formulas [Proc. R. Ir. Acad., Sect. A 78, 17 (1978)] as the long range dipole-dipole interactions increase, while the linear response remains qualitatively unaffected. By qualitatively comparing with recent experimental measurements on glycerol, it is further demonstrated that nonlinear dielectric response experiments provide a powerful tool for the characterization of short range intermolecular interactions in polar liquids. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4855195] I. INTRODUCTION
A system in thermal equilibrium at temperature T disturbed by an external stimulus evolves towards a new stationary state. Furthermore, if the energy stimulus is much lower than the thermal energy, deviations linear in the stimulus of the expectation value of the relevant dynamical variable are sufficient to calculate the generalized susceptibility using appropriate equilibrium correlation functions and linear response theory. The simplest illustration is the orientational electric polarization of noninteracting permanent dipoles in an alternative sinusoidal field E(t) treated by Debye.1 The calculation of the stationary nonlinear response is, however, much more involved because, for example, no connection between the transient and ac response exists. Nonlinear dielectric relaxation and the dynamic Kerr effect of permanent dipoles are naturally occurring examples. Now the first calculation concerning nonlinear dielectric relaxation of noninteracting dipoles was performed by Coffey and Paranjape.2 In effect, with ϑ the polar angle of the dipole moment vector μ, these authors show that in particular, amending cos ϑ(t) to terms of order E3 , the nonlinear response consists of a term oscillating at the fundamental superimposed on which is a term oscillating at the third harmonic. This prediction was experimentally verified 20 years later3 in a system for which the macromolecules have been considered noninteracting, putting in evidence the so-called nonlinear dynamical dielectric effect. All subsequent calculations existing in the literature4–6 systematically ignore the interactions between the dipoles, electric or magnetic, while the obtaining of a simple formula for the nonlinear response using the Debye-Fröhlich model, which includes intermolecular interactions via a double-well potential, is by no means easy.7, 8 Dipolar interactions in the orientational dynamics of an assembly of dipoles always lead to involved calculations as is apparent in dielectric relaxation of polar fluids.9, 10 For ex0021-9606/2014/140(3)/034506/13/$30.00
ample, the first Debye-Lorentz calculation of the static permittivity valid for polar gases at very low densities was improved upon by Onsager,10 and, as outlined by van Vleck,11 can be transposed to assemblies of Langevin paramagnets. In essence, Onsager, in order to improve Debye’s formula, accounted for the reaction of a given polar molecule on its surroundings by supposing that a typical dipole of an assembly resided at the center of a macroscopic spherical cavity in a dielectric medium, ultimately generalizing the Debye calculation of the static dielectric constant to higher densities. However, Onsager’s model was criticized by Kirkwood12 because it ignores the interactions between the dipole of the macroscopic cavity with its nearest neighbours. Consequently, Kirkwood, and later Fröhlich, have corrected Onsager’s formula by introducing the “Kirkwood correlation factor,” which accounts for the contribution of the neighbours of a molecule to the dielectric permittivity.9, 10 In spite of the largely successful static calculation, the determination of the dynamic permittivity of a dilute assembly of dipoles is far more involved than its static counterpart because, in general, the equation of motion of a typical dipole depends on the macroscopic properties that one is trying to calculate.9 In effect, this task leads to extremely difficult calculations because the time dependence of the field involved in the equation of motion of a typical dipole is unknown. Various attempts to generalize Onsager’s approach have been given leading to permittivity values (e.g., the Onsager-Cole formula) incompatible with any relaxation time distribution as the corresponding Cole-Cole plots may lie outside of the Debye semicircle.9 One approach to the problem is that of Nee and Zwanzig.13 By assuming a delayed response between the dipole moment and the Onsager reaction field, hence allowing for a special type of memory effect (dielectric friction), they were able to rederive the Fatuzzo-Mason equation9 for the dielectric permittivity ε(ω). A review of attempts to generalize the Kirkwood-Fröhlich theory to the dynamical case was given by Madden and Kivelson.14 Furthermore, a more
140, 034506-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-2
P. M. Déjardin and F. Ladieu
general extension to intermediate to large wavelengths was given by Bagchi and Chandra15 by adapting the equations of generalized molecular hydrodynamics to dielectric relaxation. Later, Zhou, and Bagchi16 have reconsidered Zwanzig’s model17 of dielectric relaxation of dipoles located at the sites of a cubic lattice. In particular, these authors pointed out the inadequacy of both the perturbation and the Nee-Zwanzig approaches to the collective dynamics, showing that these models for dielectric friction predict a too rapid decay of the dipole-dipole correlation function. Nevertheless, they still find Debye-like behavior for the permittivity for relatively dilute dipolar systems, in agreement with both the perturbation and continuum approaches. Yet another model that has been introduced for linear dielectric relaxation is Berne’s forced rotational diffusion model.18 Inspired by the literature on dielectric relaxation in electrolytes,19 Berne treated the rotational Brownian motion of interacting electric dipoles self-consistently so that the Fokker-Planck equation (FPE) describing the orientational dynamics of the polar molecules is nonlinear. The difficulty in handling this FPE is mainly a mathematical one, which Berne overcame by linearizing it about statistical equilibrium before taking its spatial Fourier transform. Then he took the zero wave-vector limit of the resulting equation in order to obtain the relevant correlation functions for the calculation of the dielectric constant of the dipolar assembly. This model in its nonlinear version has been also considered by Warchol and Vaughan.20 Berne’s model was furthermore recently extended to magnetic relaxation of single-domain ferromagnetic particles weakly coupled by magnetic dipole-dipole interactions.21 All the aforementioned works pertain to the static properties and to the linear response of permanent dipoles to alternating fields. Nonlinear response in interacting systems of electric dipoles is more involved, and calculations involving this dynamic response are hardly available.22 The purpose of this paper is therefore to derive formulas for the (lowfrequency) nonlinear response of electric dipoles to alternating electric fields using Berne’s model in its nonlinear version, amending terms up to O(E3 ).
J. Chem. Phys. 140, 034506 (2014)
W (ϑ, ϕ, t), viz.,
1 ∂ ∂V ∂W ∂ sin ϑ βW + 2τD W = ∂t sin ϑ ∂ϑ ∂ϑ ∂ϑ 1 ∂ ∂V ∂W + 2 βW + , ∂ϕ ∂ϕ ∂ϕ sin ϑ
where τ D = βζ /2, ζ is the rotational viscous coefficient, β = (kT)−1 , k is Boltzmann’s constant, T the absolute temperature, and V is a potential made of an effective single dipole free energy Us and a term Vint holding for the effect of interaction torques, viz., V (ϑ, ϕ, t) = Us (ϑ, ϕ, t) + Vint (ϑ, ϕ, t).
We consider in this paragraph independent groups of interacting molecules, each group being represented by a single dipole in a self-consistent potential. In each group there is no periodicity (neither in the positions nor in the orientations) and the positions of the molecules are drawn at random with an average distance a between the molecules. These groups are considered in equilibrium with a thermostat, so that the ensemble of independent groups constitute a canonical ensemble on which the methods of statistical mechanics may be applied. Thus, we may start our investigation with the Smoluchowski rotational diffusion equation for the distribution function of dipole orientations on the unit sphere
(2)
The term Us in Eq. (2) contains the orientational terms due to the externally applied fields, while Vint in Eq. (2) mimics the effects of long range Maxwell interaction fields. It has been shown21 that jk Vint (ϑ, ϕ, t) = UJ KLM (rj k , ϑj k , ϕj k ) J KLM j,k k=j
∗ × YJ K (ϑ, ϕ) YLM (ϑ, ϕ) (t),
(3)
where YJK is a spherical harmonic of order J and rank K, the star stands for the complex conjugate, the angular brackets denote an ensemble average over W, namely, 2π π YJ K (ϑ, ϕ)(t) =
YJ K (ϑ, ϕ)W (ϑ, ϕ, t) sin ϑdϑdϕ, 0 0
(4) (rjk , ϑ jk , ϕ jk ) is the set of spherical polar coordinates specifying the modulus and orientation of the vector rk − rj , rj being the location of dipole number j in the laboratory frame, the jk coefficients UJ KLM are given by21 jk
UJ KLM (rj k , ϑj k , ϕj k ) = (4π )4 i J −L+P IJ LP (rj k )KJ KLMP Q YP Q (ϑj k , ϕj k ), (5) P ,Q
where
KJ KLMP Q =
II. BASIC EQUATIONS AND EQUILIBRIUM SOLUTION OF THE SELF-CONSISTENT ROTATIONAL DIFFUSION MODEL
(1)
(2J + 1)(2P + 1) 4π (2L + 1)
× P 0J 0 |L0P QJ K |LM, P 0J 0 |L0 coefficients,
and
P QJ K |LM
are
(6)
Clebsch-Gordan
IJ LP (rj k ) =
[(L + J + P + 1)/2] +1 (J + 3/2) (L + 3/2) [1 − (L + J − P )/2] 8rjL+J k × dsZ(s)s J ds Z(s )s L π 3/2
×F4
J +L−P J +L+P +1 3 3 s 2 s 2 , ; J + , L+ ; 2 , 2 , 2 2 2 2 rj k rj k (7)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-3
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
and is the Euler gamma function,23 F4 is the two-variable hypergeometric series23 ∞ ∞ (γ )m+n (δ)m+n m n x y , F4 (γ , δ; ε, ζ ; x, y) = (ε) m (ζ )m m!n! m=0 n=0
(8)
Z(s) is the linear charge density, and the Pochhammer symbol is given by (a)n = (a + n)/ (a). We note in passing that the hypergeometric series in Eq. (7) converges only if s + s < rjk , which means that the theory has a meaning only if the mean interparticle distance is larger than the particle diameter, or, put differently, if the intermolecular distance is larger than the characteristic size of a typical molecule. Now, because of Eqs. (3) and (4), it is clearly seen that the Fokker-Planck Eq. (1) is nonlinear because V is a linear functional of W. Before starting to solve this equation, however, one must make some simplifying assumptions regarding the properties of the linear charge density Z(s). From now on, we assume that Z(s)ds = 0, sZ(s)ds = μ, (9) where μ is the dipole moment modulus, and s L Z(s)ds = 0, L > 1,
(10)
so that Appell’s hypergeometric series F4 in Eq. (7) may be replaced by unity. If an electric field is applied in some direction, this will act on the dipolar field in the same direction so that we may drop transverse components of the dipole field and consider only the longitudinal one. Hence, we assume that V = V (ϑ, t), so that the longitudinal and transverse components of the dipole are disentangled. Therefore, we may neglect the azimuthal dependence of the potential energy and also that of the distribution function. Thus, Eq. (1) becomes 1 ∂ ∂V ∂W ∂W = sin ϑ βW + . (11) 2τD ∂t sin ϑ ∂ϑ ∂ϑ ∂ϑ
tions are equivalent. In other words, it is thus implicitly assumed that the dielectric has a spherical shape. Now, on expanding the distribution function W as a series of Legendre polynomials, viz., W (ϑ, t) =
[ξ (t) − λf1 (t)] 2τD f˙n (t) + fn (t) = [fn−1 (t) − fn+1 (t)]. n(n + 1) (2n + 1) (17) This hierarchy is nonlinear in the expansion coefficients fn (t) and reflects the nonlinear integro-differential character of Eqs. (1) and (11). Because of this nonlinearity, one may solve Eq. (17) by numerical methods only. Fortunately, we only require the nonlinear response to a field such that ξ 0 1. Hence, we may use perturbation methods in order to obtain the various nonlinear responses. Now, we require the nonlinear response up to ξ03 , so that we set fn (t) = fn(0) + fn(1) (t) + fn(2) (t) + fn(3) (t) + · · · ,
P (t) = ρ0 μf1 (t). Next we combine Eqs. (17) and (18) to obtain the perturbation equations for n > 0 (f0(0) = 1), fn(0) =
λf1(0) (0) (0) fn+1 − fn−1 , 2n + 1
λf1(0) (1) (1) fn+1 (t) − fn−1 (t) , 2n + 1 (20)
2τD f˙(2) (t) + fn(2) (t) n(n + 1) n =
(13)
where ξ (t) = βμE0 cos ωt = ξ 0 cos ωt and λ = 4π βρ 0 μ2 /3. One may equivalently consider that the total field seen by a typical dipole is effectively F (t) = E(t) −
Hence by assuming that the field inside the system of dipoles is given by Eq. (14), one also assumes that all spatial direc-
(19)
ξ (t) − λf1(1) (t) (0) 2τD (0) f˙n(1) (t) + fn(1) (t) = fn−1 − fn+1 n(n + 1) 2n + 1
4πρ0 μ2 cos ϑcos ϑ(t), (12) 3 where ρ 0 is the number of molecules per unit volume. Next, we assume that the dipoles are subjected to an alternating field E(t) = E0 cos ωt, so that the potential energy reads
(15)
(18)
the superscript (i) meaning the order in the field strength, the polarization being given by Eq. (15) which we rewrite as follows:
+
4π P (t), (14) 3 where the polarization along the direction of the field P(t) is given by
(16)
and using this expansion in conjunction with Eqs. (11) and (13), we arrive at the infinite hierarchy of differential recurrence relations for the expansion coefficients fn (t), namely,
Vint (ϑ, t) =
P (t) = ρ0 μcos ϑ(t).
(n + 1/2)fn (t)Pn (cos ϑ)
m=0
We also have from Eqs. (3)–(10),18, 20, 21
βV (ϑ, t) = −ξ (t) cos ϑ + λcos ϑ(t) cos ϑ,
∞
ξ (t) − λf1(1) (t) (1) (1) fn−1 (t) − fn+1 (t) (2n + 1) +
λf1(2) (t) (0) (0) f − fn−1 (2n + 1) n+1
+
λf1(0) (2) (2) f (t) − fn−1 (t) , (2n + 1) n+1
(21)
and 2τD f˙(3) (t) + fn(3) (t) n(n + 1) n =
ξ (t) − λf1(1) (t) (2) (2) fn−1 (t) − fn+1 (t) (2n + 1)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-4
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
+
λf1(3) (t) (0) (0) fn+1 − fn−1 (2n + 1)
+
λf1(2) (t) (1) (1) fn+1 (t) − fn−1 (t) (2n + 1)
which steady-state solution yields the linear dynamical polarization P1 (t) as
P1 (t) = ρ0 μf1(1) (t) = ρ0 E0 α1(1) (ω) cos ωt + α1(1) (ω) sin ωt ,
+
λf1(0) (2n + 1)
(3) (3) fn+1 (t) − fn−1 (t) .
(30) (22)
Since in Eq. (19), the functions fn(0) are not driven by the field, these functions coincide with the moments of the equilibrium distribution function W0 (ϑ), namely, fn(0) = Pn (cos ϑ)0 ,
π Pn (cos ϑ)W0 (ϑ) sin ϑdϑ,
(24)
0
and the equilibrium distribution function W0 (ϑ) is given by W0 (ϑ) =
1 e−λcos ϑ0 cos ϑ , Z(cos ϑ0 )
α1(1) (ω) =
(31) and where τ1 =
Z(cos ϑ0 ) =
(25)
(32)
2τD f˙(2) (t) + fn(2) (t) n(n + 1) n =
e−λcos ϑ0 cos ϑ sin ϑdϑ.
ξ (t) − λf1(1) (t) (1) (1) fn−1 (t) − fn+1 (t) δn,2 , (2n + 1)
so that
0
From Eq. (25), it appears that the equilibrium moments (24) are completely determined by the specification of f1(0) = cos ϑ0 . Instead of evaluating Eq. (24), we alternatively use Eq. (19) with n = 1. This leads to the solution f1(0)
τD . 1 + λ/3
Thus, a Debye spectrum with shorter relaxation time is obtained for the complex linear polarizability, as may be expected from such a model. Now we may proceed by calculating f1(3) (t). On writing down Eq. (21) and accounting for the aforementioned simplifications, we have
with π
βμ2 βμ2 ωτ1 , α1(1) (ω) = 2 2 (3 + λ) 1 + ω τ1 (3 + λ) 1 + ω2 τ12
(23)
where the angular brackets denote an equilibrium ensemble average given by Pn (cos ϑ)0 =
where the in-phase and out of phase “polarizabilities” α1(1) (ω) and α1(1) (ω) are given by
= 0,
which means that dipole-dipole interactions cannot spontaneously orient the dipolar assembly due to their frustrating character. Hence, from Eq. (19), for n > 0, fn(0) = 0.
(26)
We now may proceed to the evaluation of the linear and nonlinear polarizabilities.
fn(2) (t) = 0, n = 2,
(33)
and 3 (2) 3ξ (t) (1) 3λ (1) 2 f˙2(2) (t) + f2 (t) = f1 (t) − f (t) . (34) τD 5τD 5τD 1 This equation may be integrated by quadratures with f2(2) (−∞) = 0 to yield f2(2) (t)
3ξ0 = 5τD
t
e−3(t−t )/τD f1(1) (t ) cos ωt dt −
−∞
t ×
2 e−3(t−t )/τD f1(1) t dt ,
3λ 5τD
(35)
−∞
III. CALCULATION OF THE LINEAR AND NONLINEAR RESPONSES TO THE APPLIED FIELD
Here we give the solution of Eqs. (20)–(22) which will allow us to calculate the linear and nonlinear polarizabilities. First, Eq. (20) reads 2τD f˙(1) (t) + fn(1) (t) = 0, n > 1, n(n + 1) n
(27)
meaning that for n > 1, we have fn(1) (t) = 0. For n = 1, we have the differential equation: λ ξ (t) (1) ˙ τD f1 (t) + 1 + f1(1) (t) = , 3 3
(28)
(29)
which pertains, e.g., to the dynamic Kerr effect. As we do not focus our attention on the dynamic birefringence, we do not explicitly give the result of the integrals here. For λ = 0, the result of the calculation agrees with that of Benoit24 and that of Coffey and Paranjape.2 Owing to Eqs. (26) and (28) and (22) and (33) reads 2τD f˙(3) (t) + fn(3) (t) n(n + 1) n ξ (t) − λf1(1) (t) (2) (2) fn−1 (t) − fn+1 (t) = (2n + 1) λf1(3) (t) (0) (0) f + − fn−1 δn,1 , (2n + 1) n+1
(36)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-5
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
so that for n = 1 we obtain the equation λ λf (1) (t) − ξ (t) (2) (3) ˙ f1(3) (t) = 1 f2 (t), τD f1 (t) + 1 + 3 3 (37) which again may be integrated by quadratures with f1(3) (−∞) = 0. The result of the calculation yields the
nonlinear polarization as
P3 (t) = ρ0 E03 α3(1) (ω) cos ωt + α3(1) (ω) sin ωt + α3(3) (ω) cos 3ωt + α3(3) (ω) sin 3ωt ,
where the in-phase and out-of-phase components of the nonlinear polarizabilities α3(s) (ω) = α3(s) (ω) − iα3(s) (ω), s = 1, 3 are provided by
3β 3 μ4 (3 + λ)2 (4λ − 39) ω4 τ14 + 378 + 522λ + 51λ2 ω2 τ12 + 729 =− ,
3 20 (3 + λ)4 81 + 4 (3 + λ)2 ω2 τ 2 1 + ω2 τ 2
(39)
β 3 μ4 ωτ1 81 (2λ − 21) + 3 (3 + λ) (λ (λ − 45) − 198) − (3 + λ)2 ω2 τ12 ω2 τ12 = , 3 10 (3 + λ)4 81 + 4 (3 + λ)2 ω2 τ12 1 + ω2 τ12
(40)
α3(1) (ω)
1
α3(1) (ω)
(38)
1
9β 3 μ4 (3 + λ)2 (51 + 20λ) ω6 τ16 + 3 (279 − 4λ (λ (6 + λ) − 21)) ω4 τ14 − 3 (λ (78 + λ) − 99) ω2 τ12 − 81 = , (41)
3 20 (3 + λ)4 81 + 4 (3 + λ)2 ω2 τ12 1 + ω2 τ12 1 + 9ω2 τ12
9β 3 μ4 ωτ1 9 (2λ − 21) + λ3 − 243λ − 297 ω2 τ12 − 3 (3 + λ) (3 + λ (19 + 4λ)) − (3 + λ)2 ω2 τ12 ω4 τ14 (3) α3 (ω) = . (42)
3 10 (3 + λ)4 81 + 4 (3 + λ)2 ω2 τ12 1 + ω2 τ12 1 + 9ω2 τ12 α3(3) (ω)
For λ = 0, these formulas reduce to the Coffey-Paranjape formulas2 for the nonlinear response of non-interacting dipoles to alternating electric fields, viz., α3(1) (ω) =
13ω2 τD2 − 27 β 3 μ4 , 180 1 + ω2 τ 2 2 9 + 4ω2 τ 2 D
α3(1) (ω) = −
(43)
D
β 3 μ4 ωτD 21 + ω2 τD2 2 , 9 + 4ω2 τD2 90 1 + ω2 τD2
(44)
17ω2 τD2 − 3 β 3 μ4 , = 60 1 + ω2 τD2 9 + 4ω2 τD2 1 + 9ω2 τD2 (45) β 3 μ4 ωτD 3ω2 τD2 − 7 , α3(3) (ω) = 30 1 + ω2 τD2 9 + 4ω2 τD2 1 + 9ω2 τD2 (46) where we have used τ 1 = τ D for λ = 0. Thus, it appears that in order to reproduce the Coffey-Paranjape results, the formal limit λ → 0 is taken in the expressions for the complex polarizabilities, and this limit refers to the very diluted limit, where intermolecular interactions are negligible. Our equations (39)–(42) also formally agree with the results of Felderhof and Jones,22 who considered nonlinear response of magnetic dipoles, provided their interaction parameter C is replaced by −λ/3. We emphasize that the nonlinear susceptibilities provided by Eqs. (39)–(42) are (mathematically) valid for an arbitrary positive value of λ, and constitute the central result of our work. By inspection of Eqs. (39)–(42), we remark that the time scales of the nonlinear response are α3(3) (ω)
τ1 =
τD , 1 + λ/3
τ2 =
2 (3 + λ) τ1 , 9
τ3 = 3τ1 .
If τ 1 represents the fundamental time scale, then the time scales of the nonlinear response occurring in Eqs. (39)–(42) behave in the same fashion as in the Coffey-Paranjape formulas. However, the frequency behavior of Eqs. (39)–(42) is more complicated. In order to complete our calculation, we need now to introduce Maxwell field corrections to the nonlinear response.
(47)
IV. MAXWELL FIELD CORRECTIONS
As alluded to above, we need to account for Maxwell field corrections to the nonlinear susceptibilities in order to obtain quantities which can be compared with experiment. Here, we have in mind the experimental work recently published by Brun et al.25–29 on glycerol, where in particular, the authors seek to obtain universal susceptibility curves by representing the results as a function of ωτ , where τ is a certain relaxation time. In order to accomplish this, we may proceed in two ways. In the following, the Roman superscripts I and II will refer to the two ways. We first start with the linear response.
A. Determination of the linear response
First way: here we account for static effects of the Maxwell field only. Then we remark that the in-phase and out-of-phase molecular responses to E0 are given by Eq. (31), which means that the complex linear response per molecule is given by α1(1) (ω) = α1(1) (ω) − iα1(1) (ω). The linear polarization frequency response P˜1(1) (ω) is therefore given by P˜1(1) (ω) = ρ0 α1(1) (ω)E0 . Since we have assumed that the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-6
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
molecules have only a permanent dipole moment, the static polarizability of an isolated molecule is therefore given by βμ2 /3. We therefore may write the linear polarization response as E0 1 λ , (48) P˜1(1) (ω) = 4π 1 + iωτ1 1 + λ/3 which means that the Maxwell field amplitude El is given by El = (1 + λ/3)−1 E0 . If we write P˜1(1) (ω) = χ1(1) (ω)El then, the linear dielectric constant ε1(1) (0) is given by ε1(1) (0) − 1 = 4π χ1(1) (0) = λ,
(49)
which means that the Maxwell field amplitude is given by El = 3[2 + ε1(1) (0)]−1 E0 , e.g., by the Lorentz static formula. The susceptibility spectrum, given by χ1(1),I (ω) =
1 λ 4π 1 + iωτ1
(50)
is then Lorentzian with the relaxation time τ 1 . Second way: Here, we account for the dynamical effects of the Maxwell field on the linear response. In order to accomplish this, we follow the spirit of the work of Madden and Kivelson.14 As our calculations are valid for a dielectric sphere only, the linear polarization frequency response is given in terms of the linear complex permittivity by 3 ε1(1),I I (ω) − 1 P˜1(1) (ω) = E0 . 4π ε1(1),I I (ω) + 2
(51)
This expression should be equal to the linear polarization spectrum expressed in terms of molecular quantities, P˜1(1) (ω) = ρ0 α1(1) (ω)E0 . By accomplishing this we obtain the following expression for the dynamical susceptibility, namely, 4π χ1(1),I I (ω) = ε1(1),I I (ω) − 1 =
ε(1) (0) − 1 λ = 1 , 1 + iωτD 1 + iωτD (52)
namely, again a Lorentzian spectrum, with now a relaxation time τ D . Here, we note in passing that Jadzyn et al.30 have shown that the relaxation time of the linear dielectric response is the same in two nematogenic compounds, where the molecular moment differs by a factor as large as 2: even if our calculations do not apply to liquid crystals, this experimental result qualitatively supports our previous finding, namely, the fact that the linear response does not “feel” qualitatively the existence of long range dipolar interactions. Now, since the qualitative frequency behavior of Eqs. (50) and (52) is the same, so that we have no clear indication regarding the way to proceed in order to account for the Maxwell field in the determination of the nonlinear response. This is why we now turn to the calculation of the nonlinear susceptibilities using the two ways in section B of this paragraph. B. Determination of the nonlinear response
First way: Here we account for the static Maxwell field effects only. Thus we first write the nonlinear steady polariza-
tion response as follows: P3 (t) =
3El3 (1),I χ3 (ω) cos ωt + χ3(1),I (ω) sin ωt 4 E3 + l χ3(3),I (ω) cos 3ωt + χ3(3),I (ω) sin 3ωt , (53) 4
where again El = (1 + λ/3)−1 E0 . Next, we identify Eq. (53) with Eq. (38). We have λ 3 ρ0 (s) (s),I α3 (ω), s = 1, 3, (54) χ3 (ω) = 1 + 3 n(i) 3 (3) where n(1) 3 = 3/4 and n3 = 1/4. Second way: Here we account on the dynamical Maxwell field effects on the nonlinear response. Because the derivations are more involved than for the static field corrections, we only give the results and leave the details of calculation to Appendix A. The expressions for the nonlinear susceptibilities are provided by the following equations, viz. 2 ε1(1),I I (ω)+2 4ρ0 ε1(1)∗,I I (ω)+2 (1),I I (ω) = α3(1) (ω), χ3 3 3 3 (55)
χ3(3),I I (ω)
= 4ρ0
ε1(1),I I (ω) + 2 3
3 α3(3) (ω),
(56)
(s) where α3(s) (ω) = α(s) 3 (ω) − iα3 (ω) and the in-phase and out-of-phase parts of this expression are given by Eqs. (39)– (42). One notes that if ω is set equal to zero in the complex permittivity factors, one indeed recovers the static results Eq. (54). Hence, accounting for dynamical effects of the Maxwell field, the nonlinear susceptibilities χ3(1) (ω) and χ3(3) (ω) are given by Eqs. (55) and (56), respectively, where the molecular nonlinear responses to E0 are corrected by the dynamical Lorentz cavity field factor. The quantities of experimental interest are, for the two ways: (s),j ρ0 (s),j X3 (ω) = (1) 2 χ3 (ω), s = 1, 3 ; j = I, I I. β χ1 (0) (57) The quantities (57) are dimensionless and render a static value of 0.2 in the interaction-free limit, so facilitating the comparison with the work of Brun et al.25–29 Notice also that for both ways, the nonlinear susceptibilities agree with the CoffeyParanjape ones if the Maxwell field factors are replaced by unity and the nonlinear polarizabilities by Eqs. (43)–(46).
V. RESULTS AND DISCUSSION
The expressions we have obtained for the nonlinear response, Eqs. (39)–(42) are easily computed and therefore numerical analysis of them is straightforward. Furthermore, as alluded to in the text, they can be used for any value of the interaction parameter λ, which can in practice be considered as a fitting parameter with reasonable values (as demonstrated by Eq. (49) λ is of the order of the linear static dielectric constant ε1(1) (0)). We first focus onto the frequency dependence of the (s) (s) nonlinear response by studying A(s) 3 (ω) = |α3 (ω)/α3 (0)|,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-7
P. M. Déjardin and F. Ladieu
J. Chem. Phys. 140, 034506 (2014)
FIG. 3. Dimensionless modulus of the ω nonlinear susceptibility component (static Maxwell field corrections) as a function of ωτ 1 for various values of λ. The value λ = 0 is the Coffey-Paranjape result. FIG. 1. Normalized modulus of the ω nonlinear polarizability component as a function of ωτ 1 for various values of λ. The value λ = 0 is the CoffeyParanjape result.
(s),j
then we discuss the behavior of X3 (ω) which are the quantities to be compared with experimental data. Figures 1 and 2 illustrate the behavior of the ω and 3ω normalized nonlinear responses to the externally applied field E0 (t), without Maxwell field corrections, where the quantitiesA(s) 3 (ω) are plotted against ωτ 1 . Note that the (possibly strong) temperature dependence of the rotational friction coefficient is totally absorbed by the chosen ωτ 1 horizontal scale. Figures 1 and 2 allow to see how A(s) 3 (ω) deviate from the Coffey-Paranjape formulas. In effect, the nonlinear response to the applied ac field exhibits a humped shape as λ increases, and therefore our results clearly demonstrate that the nonlinear response is affected by intermolecular interactions. This contrasts with the linear response which in our model keeps a Debye shape for all the values of λ. Now, Figures 3–6 show the frequency behavior of the modulus of the dimensionless nonlinear susceptibilities (s),j X3 (ω), where Maxwell field corrections are accounted for. Therefore, these quantities can be directly compared with experimental data. Furthermore, these quantities are of di-
FIG. 2. Normalized modulus of the 3ω nonlinear polarizability component as a function of ωτ 1 for various values of λ. The value λ = 0 is the CoffeyParanjape result.
rect interest with respect to experiments such as those of Refs. 25–29 in which one of us is directly involved. Figures 3 and 4 show the behavior of the nonlinear responses X3(s),I (ω), e.g., when only static effects due to the Maxwell field are accounted for. One may see that the static values of the nonlinear responses decrease, and a humped shape develops while λ increases, the peak value being lower than the value for noninteracting dipoles, while the high-frequency behavior is the same for all λ values. Figures 5 and 6 show the behavior of the nonlinear responses X3(s),I I (ω), when dynamical effects due to the Maxwell field are accounted for. Indeed the behavior of the static value is the same as those for X3(s),I (ω), while the humped shape has now disappeared completely whatever the λ values. At all frequencies, the value of the nonlinear susceptibilities are all less than the those provided by the nonlinear susceptibilities of the ideal gas and in contrary with X3(s),I (ω) the high frequency values depend on λ. In both cases I and (s),j II, the values of the nonlinear susceptibilities X3 (ω) are less than the ones returned by the Coffey-Paranjape formulas. Thus, we may infer that if the experimental values are in excess with respect to those returned by the Coffey-Paranjape formulas, then this excess is due to short range intermolecular interactions. Comparison of our results with previous ones can only be achieved in the context of linear response. First, we
FIG. 4. Dimensionless modulus of the 3ω nonlinear susceptibility component (static Maxwell field corrections) as a function of ωτ 1 for various values of λ. The value λ = 0 is the Coffey-Paranjape result.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-8
P. M. Déjardin and F. Ladieu
FIG. 5. Dimensionless modulus of the ω nonlinear susceptibility component (dynamic Maxwell field corrections) as a function of ωτ 1 for various values of λ. The value λ = 0 is the Coffey-Paranjape result.
should compare our results with those obtained by Berne.18 As demonstrated in Appendix B, Berne’s results are valid for weak densities, and for a sample which shape is similar with an infinite slab. In effect, Berne (and also Warchol and Vaughan for the nonlinear version20 ) uses the field F(t) = E(t) − 4π PZ (t)eZ , which is the macroscopic field in a uniformly polarized infinite slab and which differs from that in Eq. (14). The latter is isotropic and for a single microscopic relaxation time, we have a single macroscopic relaxation time. In contrast, Berne uses a priori an anisotropic field, which is compatible with the Fatuzzo-Mason theory for purely polar media, while we use the Kubo formula,9 which is compatible with Eq. (14). Nevertheless, as we demonstrate in Appendix B, the linear complex permittivity (e.g., the macroscopic response) for both models is the same. Now, strictly speaking, Berne’s calculation is valid only for weak dipolar concentrations, as we demonstrate in Appendix B. The calculation presented there clearly demonstrates that (a) Onsager’s equation for polar liquids is recovered at weak densities and (b) the Fatuzzo-Mason theory for purely polar liquids renders the Debye formula for the complex permittivity, just as in our work. In effect, the range of validity of Berne’s calculation may fairly well be extended, and that both Berne’s and Onsager’s formulations are valid at weak densities only, left apart short range interactions. In contrast, the (extended) Berne’s theory may now have a range of validity from an experimental point
FIG. 6. Dimensionless modulus of the 3ω nonlinear susceptibility component (dynamic Maxwell field corrections) as a function of ωτ 1 for various values of λ. The value λ = 0 is the Coffey-Paranjape result.
J. Chem. Phys. 140, 034506 (2014)
of view, which is probably the same as ours (see the next paragraph) in its present stage, simply because when the concentration increases (or the temperature decreases), short range intermolecular interactions have to come into play for a realistic description of the statics and the dynamics of the system. The comparison of our results with those of Chandra and Bagchi31 can also be made. For linear response, we obviously obtain the same answer in the zero wave vector limit and very weak densities, as it must be. However, we have to notice that the calculation of the nonlinear response is not possible from the Chandra-Bagchi Fokker-Planck equation (FPE). Here, we have to stress that Berne’s FPE used by us, and that of Chandra and Bagchi are similar, but not identical, because they are not derived from the same starting point and the drift coefficients are different in these equations. Therefore, even if there are some similarities, they definitely produce different quantitative results. In particular, and this is reckoned by Chandra and Bagchi themselves,31 the direct correlation function in their formulation is time-independent, which means that the pair correlation function and the pair distribution function are also time-independent, therefore not allowing one to study interacting systems far from equilibrium. Furthermore, the identification of Berne’s drift coefficients with theirs is only formal. Also, strictly speaking, the Bagchi-Chandra nonlinear FPE is valid for diluted polar systems and this is also reckoned by them.31 Now, it is very difficult to work out the nonlinear response from their equation and compare this eventually feasible calculation with our present results because in order to achieve such a task, the knowledge of the time dependence of the direct correlation function is required. In effect, this task goes beyond the Ornstein-Zernike formalism which is established for systems in equilibrium only. We remark in passing that it is in principle possible to calculate the wave vectordependent complex nonlinear polarizabilities from our model. This calculation is, however, far more involved than what we present in our work and is also beyond our scope. (s) (ω) Now, in Refs. 25–29 the measured values of X3,gly in supercooled glycerol close to its glass temperature transition temperature Tg were contrasted with the “trivial” values (s) (ω) obtained by Coffey and Paranjape2 for noninteractX3,triv (s) (ω) has a humped shape in freing dipoles. The fact that X3,gly quency with an amplitude which depends significantly on the (s) (ω) does temperature T was contrasted with the fact that X3,triv not depend on temperature and is monotonous in frequency. The conclusion of Refs. 25–29 was that the temperature de(2) (ω) gave some information pendence of the hump of X3,gly about the glassy correlations—as predicted in32 —i.e., about the temperature dependence of the number Ncorr of dynamically correlated molecules in supercooled glycerol. In this context, the model developed in the present paper gives a re(s) (ω) since it goes beyond the ideal fined calculation of X3,triv case of infinitely diluted dipoles and takes into account the finite concentration effects of the dipoles and the resulting finite long range dipole-dipole interactions. This is why we shall now contrast the calculation displayed in Figures 3–6 with the experimental data of Refs. 25–28. In case II (Figs. 5 and 6), X3(s),I I (ω) has the same qualitative frequency dependence as the Coffey-Paranjape formulae, and X3(s),I I (ω) is always decreasing when λ
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.248.155.225 On: Tue, 25 Nov 2014 06:16:38
034506-9
P. M. Déjardin and F. Ladieu
increases. As a result, in case II, the conclusions of Refs. 24–28 are qualitatively reinforced by the present refined (s) (ω). calculation of X3,triv In the case I (Figs. 3 and 4), a more detailed discussion must be done: (i) The values computed in the present paper have a humped shape in frequency: this hump is absent at λ = 0, where our Eqs. (39)–(42) reduce to the CoffeyParanjape result. This is qualitatively important since this seems to indicate that a model without glassy correlations, such as the present one, can yield a humped frequency behavior (as recently found by Diezemann33 ). However, we emphasize that, even if our calculation is technically valid for any value of λ, our model cannot capture the “dipolar glassy phase” that one expects when λ 1. This comes from the mean field approximation used in our derivation, where we assume that the environment is strictly the same for all the dipoles. We thus cannot hope to describe accurately the low temperature phase where the dipole-dipole interaction supersedes the thermal energy. Little is known about this dipolar glassy phase. To the best of our knowledge, the physical situation which is the closest to our calculation is that of concentrated ferrofluids where single domain anisotropic magnetic nanoparticles interact through dipolar couplings. In that case, a transition, at low temperature, towards a superspin glass phase has been evidenced and the cubic susceptibility is singular,34 revealing the amorphous long range order which sets in at this transition. This is why, coming back to Figures 3 and 4 and focusing on the range 0
