We present a new molecularly informed statistical field theory model of inhomogeneous polarizable soft matter. The model is based on fluid elements, referred to as beads, that can carry a net monopole of charge at their center of mass and a fixed or induced dipole through a Drude-type distributed charge approach. The beads are thus polarizable and naturally manifest attractive van der Waals interactions. Beyond electrostatic interactions, beads can be given soft repulsions to sustain fluid phases at arbitrary densities. Beads of different types can be mixed or linked into polymers with arbitrary chain models and sequences of charged and uncharged beads. By such an approach, it is possible to construct models suitable for describing a vast range of soft-matter systems including electrolyte and polyelectrolyte solutions, ionic liquids, polymerized ionic liquids, polymer blends, ionomers, and block copolymers, among others. These bead models can be constructed in virtually any ensemble and converted to complex-valued statistical field theories by Hubbard-Stratonovich transforms. One of the fields entering the resulting theories is a fluctuating electrostatic potential; other fields are necessary to decouple non-electrostatic interactions. We elucidate the structure of these field theories, their consistency with macroscopic electrostatic theory in the absence and presence of external electric fields, and the way in which they embed van der Waals interactions and non-uniform dielectric properties. Their suitability as a framework for computational studies of heterogeneous soft matter systems using field-theoretic simulation techniques is discussed.

## I. INTRODUCTION

There is considerable recent interest in heterogeneous soft materials that possess domains of varying dielectric properties and in which ions (either counter-ions or salt) are localized in particular regions, usually those with the highest polarizability. Specific classes of materials include polyelectrolyte block and graft copolymers,^{1–4} polymerized ionic liquids,^{5–9} ionomers,^{10–12} and polyelectrolyte complexes,^{13–16} among others. Such systems are of importance in energy devices such as battery electrolytes,^{17,18} separators,^{19} and fuel cell membranes,^{20} as well as media for separations,^{21–23} personal care products, adhesives, and drug and flavor/fragrance delivery systems.^{5,24,25} In designing materials for such applications, it is often critical to control the morphology of the domains possessing the ions for reasons of transport or mechanical properties. For example, in cases where macroscopic ionic conductivity is required, the ion-laden domains must percolate to create continuous ion channels.^{26}

The theoretical description of such systems is particularly challenging. Even “simple” polyelectrolyte solutions can demonstrate rich behavior such as counterion condensation^{27} and gelation arising from multivalent ion bridging^{28} or hydrophobic moieties. Polyelectrolyte block copolymers with charged and uncharged blocks in aqueous media are even more complex, producing mesophase structures with hydrated ion-containing domains and hydrophobic domains largely devoid of both water and ions. The local dielectric constant can vary as much as two orders of magnitude between such domains and if the hydrophobic “anchoring” domains are glassy or semi-crystalline, the morphology of the material is likely a non-equilibrium structure that reflects the processing history.^{1,29} Complexes formed by oppositely charged polyelectrolytes can be similarly complicated as they can “salt-out” as solid precipitates, or persist as homogeneous or structured dense fluid phases often referred to as complex coacervates.^{30}

Standard theoretical approaches follow either molecular simulations based on atomistic, or more frequently coarse-grained, particle models, field-theoretic methods such as self-consistent field theory and its extensions, or hybrids of particle and field descriptions.^{15,16,31–35} In most cases, the treatment of the electrostatic interactions is significantly oversimplified by assuming a uniform dielectric constant for the medium, which can be exceedingly inaccurate in a multi-domain system with large dielectric contrast. Moreover, the vast majority of studies neglect the “self-energies” of ions in the local dielectric environment. These reflect monopole-dipole interactions and are responsible for localizing ions in accessible domains with the largest dielectric constant. Investigators ignoring such effects are often confronted with the nonphysical behavior of ions leaking into low dielectric regions.

Other modeling challenges relate to external field effects. Several groups have reported on methods for extending self-consistent field theory (SCFT) to treat the self-assembly of polymer blends or block copolymers in static electric fields.^{36–39} In these approaches, linear constitutive laws^{33,36,40} are introduced relating the local dielectric constant to the local density of isotropically polarizable segments. While qualitatively successful, it is not obvious how to extend these theories beyond the mean-field approximation inherent in SCFT, nor in a consistent way to polymer chains possessing segments with anisotropic polarizability.^{41}

A recent series of papers by Zhen-Gang Wang and collaborators represent a significant advance towards addressing several of these challenges of modeling inhomogeneous charged soft materials.^{34,42,43} In the case of a simple electrolyte model with implicit solvent, Wang^{42} demonstrated a new field-theoretic approach that smears the microscopic density operator of the ions from a delta function to a Gaussian of finite width. This not only regularizes the field theory, rendering it free of ultraviolet divergences but also naturally introduces an ion self-energy of a form familiar from the Born solvation model^{44} and dependent on a locally averaged dielectric constant *ϵ*(**r**). It is important to note that the form of *ϵ*(**r**) must be externally supplied and is not self-consistently determined by the theory. With collaborators, Wang has used this powerful framework to address a range of interesting problems relating to ion distributions, thermodynamics, and self-assembly of inhomogeneous charged polymer systems.^{34,43}

A missing ingredient in Wang’s impressive formalism is the ability to self-consistently determine the dielectric function *ϵ*(**r**) within the confines of the theory, without the requirement of external specification or phenomenologically connecting it to the local species density. The reason that the dielectric properties are external to the theory is that the fundamental objects, namely ions and implicit solvent, were not considered to be *polarizable*, even though the ions were assigned monopoles. Indeed, it is the polarizability carried by these various species that contributes at a mesoscopic level to the dielectric contrast manifest in *ϵ*(**r**) as segregation into domains proceeds.

A first step towards such a more comprehensive theory was made by Nakamura, Shi, and Wang^{45} in considering the solvation of molecules with permanent and induced dipoles about a central fixed ion. However, the particular focus of this work did not reveal the full structure of a polarizable field theory and the analysis was restricted to a mean-field calculation about the fixed-ion broken symmetry state. There have been other attempts to include fixed and induced dipoles in molecular field theories, including “polarizable Poisson-Boltzmann” formalisms to treat dielectric inhomogeneity in polarizable solvent due to charges and charged surfaces,^{46,47} as well as electrolyte systems with polarizable counterions and a homogeneous dielectric background.^{48} These studies nevertheless similarly fall short of revealing an unapproximated field theory of explicitly enumerated molecular species possessing a polarizability and optionally a net free charge.

Here we show that Wang’s original field-theoretic formalism for inhomogeneous electrolytes^{42} can be straightforwardly extended to broad classes of heterogeneous dielectric or electrolyte fluids comprised of coarse-grained particles, or “beads,” that in addition to optionally carrying a charge (monopole), are intrinsically polarizable or carry a permanent dipole by means of a Drude distributed charge model. This approach, coupled with the aforementioned smearing strategy, yields a statistical field theory that not only self-consistently embeds spatially dependent dielectric constitutive behavior in inhomogeneous fluid or soft material structures, but also generates self-energy terms (one-body potentials manifesting monopole-dipole interactions) that serve to localize charges in high-dielectric domains. An elegant feature of the approach is that attractive van der Waals interactions among beads are automatically generated; in a two-component fluid this means that specifying distinct polarizabilities (*α _{A}*,

*α*) of the two bead types will uniquely determine an effective Flory

_{B}*χ*parameter proportional to (

*α*−

_{A}*α*)

_{B}^{2}addressing the tendency for phase separation, as well as fully determining the local dielectric constant

*ϵ*(

**r**) in relation to the local bead densities. Thus, unlike virtually all mesoscopic models of soft matter to date,

*χ*and

*ϵ*(

**r**) are not parameters to be independently or externally supplied.

We show that our molecularly informed field theory is fully consistent with continuum dielectric theory, including the effects of externally applied electric fields. Moreover, because the theory embeds exact fluctuation physics, it provides a framework for studying electric field effects on phase transitions and soft material morphology that transcend the mean-field approximation. It is further demonstrated that the approach can be extended to polymeric fluids of arbitrary architecture, including variable patterns of monopole charges and isotropic or anisotropic polarizability of chain segments. The implicit embedding of long-range van der Waals interactions has benefits beyond consistent treatment of *χ* and *ϵ*(**r**). For example, unlike field theories and coarse-grained particle models with purely local contact interactions, the present framework contains the mesoscopic physics necessary to study coalescence, break-up, wetting, and dewetting phenomena in drops and thin films of polymers and complex fluids.

The field theory models developed here are formally exact representations of underlying microscopic particle models. They can be approximately solved by analytical techniques such as loop expansions about homogeneous states, mean-field approximations such as SCFT (usually implemented numerically), or by approximation-free full numerical attack using the field-theoretic simulation (FTS) framework.^{49–51}

The structure of this paper is as follows. In Section II we discuss the general framework for building statistical field theory models of fluids composed of charged and polarizable beads. In Section III, we consider the special case of a homogeneous one-component fluid of uncharged, polarizable beads. Using a Gaussian approximation, we explore dielectric screening and constitutive behavior and elucidate the attractive van der Waals interactions among beads. In Section IV, we turn to a two-component dielectric bead model and use similar methods to identify and quantify the effective *χ* parameter and dielectric function implied by the difference in polarizabilities between the two bead species. Section V restores positive and negative monopole charges to the two-component model and, using a Gaussian approximation, recovers the limiting Debye-Hückel dilute behavior with an embedded effective dielectric constant. The self-energies of the ions are shown to have a form consistent with the theory of Wang for a homogeneous electrolyte.^{42} Section VI addresses the important extension to polymer chains built from polarizable and (optionally) charged beads. We illustrate the use of various chain models and cases of isotropic and anisotropic polarizability of polymer segments. Section VII discusses the application of static external electric fields and the connection to continuum dielectric theory. Finally, in Section VIII we discuss the prospects of this new polarizable field theory formalism for investigating a variety of inhomogeneous soft matter systems and the suitability of FTS techniques such as complex Langevin sampling for conducting efficient simulations without simplifying approximations.

## II. GENERAL FORMALISM: POLARIZABLE BEAD FLUID MODELS

### A. Model and field theory representation

As described previously, we construct fluid models from coarse-grained soft particles that will be referred to as “beads.” As depicted in Figure 1, the *i*th bead has a characteristic diameter *a* and three sites at positions **r**_{ip} indexed by *p* = 1, 2, and 3. The bead center of mass is located at position **r**_{i1} ≡ **r**_{i} and has associated charge *q*_{i1} ≡ *q _{i}*, while sites 2 and 3 correspond, respectively, to positive and negatively charged sites with partial charges

*q*

_{i2}=

*δq*and

_{i}*q*

_{i3}= −

*δq*that are connected by a harmonic spring with spring constant

_{i}*K*to form a Drude oscillator. As

**r**

_{i1}= (

**r**

_{i2}+

**r**

_{i3})/2 is the center of mass of the oscillator, only two of the three site positions are independent; for convenience we choose the center of mass

**r**

_{i}and the spring displacement,

**s**

_{i}≡

**r**

_{i2}−

**r**

_{i3}, as the two independent coordinates to specify the configuration of the bead and partial charges. The net charge on the

*i*th bead

*q*is the monopole charge placed at

_{i}**r**

_{i}. Evidently the polarization of a bead results from the response of the bead’s Drude oscillator to an electric field applied externally or created internally by the charge distribution of other beads in the fluid.

For a collection of *n* such beads, we follow Wang’s approach of smeared charges and define a smeared microscopic charge density operator by the expression

where Γ(**r**) = (2*πa*^{2})^{−3/2}exp(−*r*^{2}/2*a*^{2}) is a Gaussian smearing function of characteristic range *a*, which we choose arbitrarily to be the size of a bead. The function is normalized so that Γ(**r**) approaches a three-dimensional Dirac delta function for *a* → 0^{+}. Thus, Eq. (1) reflects a microscopic density that is smeared or convolved with the Gaussian Γ(**r**).

The Coulomb electrostatic energy associated with the 3*n* charges of the *n* beads in a particular configuration can be written as

where we employ CGS electrostatic units and energy units with *k _{B}T* = 1. It is important to note that Eq. (2) includes site self-interaction terms and intra-bead site-site interactions in addition to inter-bead electrostatic interactions. As will be demonstrated below, and by analogy with Wang’s formalism, the former two types of interactions contribute to inconsequential shifts in reference chemical potential in a homogeneous system and therefore have no thermodynamic significance.

In contrast, in an *inhomogeneous* system, these single-molecule contributions produce physically important self-energy terms with the form of Born solvation energies, which (as discussed by Wang) are important to localizing charges in high dielectric regions and positioning them non-uniformly close to external surfaces. Beyond rendering these self-energy terms finite, the smearing of the charge density has the added benefit of making the resulting field theory ultraviolet convergent, which eliminates complications in both analytical calculations and FTS simulations.^{52,53}

We shall see that the electrostatically interacting Drude oscillators in the ensemble of beads will naturally create *attractive* van der Waals interactions among beads. The only additional interaction thus required to construct a suitable model of a bead liquid is a repulsive interaction among bead centers to prevent collapse and stabilize a uniform fluid. There are various possible choices, including a general pair interaction,

where $ \rho \u0304 ( r ) = \u2211 i = 1 n \Gamma ( r \u2212 r i ) $ is the smeared density operator of the bead centers of mass. The simplest choice for the pair potential *u*(*r*) is a repulsive contact potential, *u*(*r*) = *uδ*(**r**) with *u* > 0 an excluded volume parameter. It is easily shown that Eq. (3) with the contact potential for *u*(*r*) is equivalent to a soft Gaussian repulsion of strength *u* and range $ 2 a$ between bead centers.^{52} Here we have chosen to smear the density operator with the same Γ function used in smearing the charge density operator; other choices are evidently possible. Finite self-interaction terms of no thermodynamic consequence are included in Eq. (3) and the smearing of both microscopic charge and bead densities is required to render the corresponding field theories ultraviolet convergent.

Having specified the form of the potential energy *U* = *U _{e}* +

*U*describing electrostatic and non-electrostatic interactions, one can construct partition functions, expressed as integrals over bead and Drude oscillator coordinates, in any desired ensemble. These particle-based partition functions can be exactly transformed to statistical field theories by following the well-worn path of Hubbard-Stratonovich (HS) transforms.

_{r}^{49,54}For example, the canonical (

*nVT*) partition function for the bead model described above, assuming identical beads, can be expressed in field-theoretic form as

where *H*[*ϕ*, *w*] is an effective Hamiltonian (or action) functional given by

In these expressions, *ϕ*(**r**) and *w*(**r**) are real-valued fields that were introduced to decouple electrostatic and non-electrostatic interactions, respectively, and $\u222bD\varphi \u222bDw$ represents functional integrals over the two fields. Smeared versions of these two fields, reflecting convolutions with the function Γ(**r**), are denoted by $ \varphi \u0304 ( r ) =\u222bd r \u2032 \Gamma ( r \u2212 r \u2032 ) \varphi ( r \u2032 ) $ and $ w \u0304 ( r ) =\u222bd r \u2032 \Gamma ( r \u2212 r \u2032 ) w ( r \u2032 ) $, respectively. The object *Z*_{0} = *V ^{n}*/

*n*! is the ideal gas translational partition function of the bead fluid, while

*Z*

_{ϕ}and

*Z*are normalizing denominators associated with the Hubbard-Stratonovich transforms defined by the following functional integrals:

_{w}*n*-dependence, but do depend on volume

*V*, and so are important in computing pressure and stress operators. Most significantly, the HS transforms have decoupled the beads so that they interact only with the two fields. The remaining intra-bead degrees of freedom are contained in the

*single-bead partition function*$Q [ i \varphi \u0304 , i w \u0304 ] $, whose form will be discussed below. Crucially,

*Q*depends on purely imaginary ($i\u2261 \u2212 1 $) smeared fields. The imaginary character is an artifact of the HS transforms and renders the resulting Hamiltonian

*H*complex-valued.

In the grand-canonical (*zVT*) ensemble, the relevant partition function and Hamiltonian are similar

where *z* ∝ exp(*βμ _{b}*) is the activity of beads and

*μ*the bead chemical potential.

_{b}We next turn to the form of $Q [ i \varphi \u0304 , i w \u0304 ] $, distinguishing cases where the distributed charge on a bead reflects induced or permanent electric dipoles.

### B. Single bead partition function: Induced dipole case

The induced dipole case corresponds to the Drude oscillator model described previously. Denoting the displacement of the two partial charges by **s** = **r**_{+} − **r**_{−} and the bead center of mass by **r** = (**r**_{+} + **r**_{−})/2, the single-bead partition function can be written as

where *q* is the monopole charge and ±*δq* are the partial distributed charges connected by a spring with spring constant *K*. We have normalized *Q* such that *Q*[0, 0] = 1. While the above expression is exact, a more useful approximate form is obtained by assuming that $ \varphi \u0304 ( r ) $ is slowly varying across the bead and carrying out a multipole expansion of the last two factors in the integrand. The resulting expression is accurate to second order in gradients of $ \varphi \u0304 $,

where

is the *polarizability* of a bead. We shall employ the simpler approximate form of Eq. (10) in all subsequent analyses that relate to beads with induced, rather than fixed, dipoles.

### C. Single bead partition function: Permanent dipole case

To treat the case of a bead with a permanent dipole moment, we replace the spring connecting the two partial charges in Figure 1 with a rigid link of length *d*. The partial charge separation vector can then be expressed as **s** = *d***n**, where **n** is a unit vector. Invoking the same multipole expansion to second order, we have

where

is the *dipole moment* of the bead. The integral over **n** in Eq. (12) is taken over the surface of the unit sphere and can be carried out explicitly leading to the simpler expression

where *j*_{0}(*x*) = (sin*x*)/*x* is the familiar spherical Bessel function.

In comparing the two cases of induced and permanent dipoles, we see from Eqs. (10) and (14) that *Q* is determined in both cases by a simple volume integral over integrands depending on both $ \varphi \u0304 $ and $ w \u0304 $, as well as the gradient of $ \varphi \u0304 $. For the special case of no net charge (monopole) on the bead (*q* = 0), *Q* depends only on $ w \u0304 $ and $\u2207 \varphi \u0304 $. In the induced case, the polarizability *α* is a free parameter of the model to regulate the strength of bead polarization, while the dipole moment *μ* plays a similar role in the permanent dipole case.

### D. Observables and operators

The field theories summarized by Eqs. (4), (5), (7), and (8) are equivalent representations of the corresponding interacting bead models. Beyond the usual thermodynamic connection formulas between the partition functions *Z* and Ξ and macroscopic thermodynamic properties, it is useful to define *field operators* that can inform about structure and other observable characteristics of the fluid. For a generic operator $ O \u0303 [ \varphi , w ] $, which is a functional of the two fields, the average value can be expressed (in the canonical ensemble) as

Because of the complex nature of the field theory, $ O \u0303 $ may be a complex number for instantaneous field values, but if it reflects a physical observable, its imaginary part must vanish on average. As a specific example, we consider an operator for the density of bead centers at some location **r** in the fluid, $ \rho \u0303 ( r ; [ \varphi , w ] ) $. By inserting a source term *U _{h}* = − ∫

*d*

**r**

*h*(

**r**)

*ρ*(

_{m}**r**) in the original particle model, with $ \rho m ( r ) = \u2211 i = 1 n \delta ( r \u2212 r i ) $ the (unsmeared) microscopic bead-center density, the average bead-center density can be generated from

Tracking the *h*(**r**) source through the HS transformation into the field theory and taking the functional derivative indicated in Eq. (16) leads to the following expression for the field operator corresponding to the local bead density:

where *ρ*_{0} = *n*/*V* is the average bead number density, and we have used the induced-dipole expression of Eq. (10) for *Q*. Similarly, an operator for the excess chemical potential (relative to the ideal gas) of the fluid in the canonical ensemble can be written as

A variety of other field operators, including operators for density-density correlation functions and stress operators, can be derived following similar approaches, as outlined in Ref. 49.

A field operator unique to the present theory relates to the local *polarization* within the fluid. We can define an unsmeared microscopic polarization (dipole) density by the expression

Applying a source term to the electrostatic energy of the form −∫*d***rP**_{m} ⋅ **E**_{a}, with an externally applied electric field **E**_{a}(**r**), the average polarization density can be written as

Propagating the source term into the field theory representation of the model and assuming the induced dipole form of *Q* from Eq. (10), we find the following expression for the *polarization density field operator*:

where

is a field operator for the *local internal electric field* created by the charges in the fluid. Again we emphasize that operators such as $ \rho \u0303 $, $ P \u0303 $, and $ E \u0303 $ are complex-valued in specific field configurations, but their expectation values over all field configurations are real.

Equation (21) has an appealing physical interpretation—the local polarization in the fluid is locally proportional to the internal electric field, with an isotropic dielectric susceptibility operator given by $ \chi \u0303 ( r ) =\alpha \rho \u0303 ( r ) $, i.e., the segment polarizability times the density of segments. The isotropy of the local response evidently emerges from the isotropy of the polarizable beads. Nonetheless, the appearance of such operator products can be deceiving; the individual elements are complex and correlated objects prior to averaging. Only in a mean-field approximation, where one breaks averages to write

does this physical interpretation based in continuum macroscopic dielectric theory rigorously hold. Indeed, beyond mean-field theory the dielectric susceptibility is not strictly linear in either *α* or the fluid density.

Another field operator specific to the present theory relates to the density of charge in the fluid. We can define an unsmeared microscopic charge density by the expression

If we introduce an externally applied potential field *ϕ _{a}* that interacts with the fluid according to

*U*= − ∫

_{a}*d*

**r**

*ρ*

_{e,m}

*ϕ*, the thermally averaged charge density is given by

_{a}with

a *charge density field operator* that corresponds to the total charge density arising from the three charges on each bead unit. As indicated by the second line of Eq. (26), this operator can be further decomposed into a *free charge density operator*, $ \rho \u0303 f $, which reflects the contribution from the monopole charge carried by each bead, and a *bound charge density operator*, $\u2212\u2207\u22c5 P \u0303 $, which gives the contribution from the polarization associated with the paired partial charges (dipole) of each bead. The former is given by

which has an unsurprising dependence on the product of the local bead density and the net free charge carried by the beads.

Other field operators of importance are the *forces* required for finding mean-field solutions or conducting field-theoretic simulations by the complex Langevin technique. These are the functional derivatives of the Hamiltonian with respect to the two fields. For the canonical model Hamiltonian of Eq. (5), the *w*-force is easily seen to be

For this force to be zero, we see that the smeared bead density must balance the molecular field *w* created by the bead-bead repulsions. Similarly, the *ϕ*-force is given by

This expression becomes more illuminating if we smear both sides of the equation and multiply by 4*πi*,

Here $ D \u0303 $ is a *dielectric displacement* (or induction) field operator defined by

We note that for vanishing *ϕ*-force, Eq. (30) recovers a smeared analog of the Maxwell equation ∇ ⋅ **D** = 4*πρ _{f}*. Similarly, Eq. (31) is readily recognized as a fluctuating and smeared analog of the familiar relation

**D**=

**E**+ 4

*π*

**P**relating displacement, field, and polarization from continuum dielectric theory.

A further expression follows in the case of induced dipoles by combining Eqs. (21) and (31) to write

where $ \u03f5 \u0303 $ is a *nonlocal dielectric permittivity operator* given by

We thus see nonlocal dielectric constitutive behavior involving the local bead density and bead polarizability arising from the smearing of the charge density in the model. Again we caution the reader that while such expressions are physically suggestive in the context of the constitutive behavior, **D** = *ϵ***E**, of continuum dielectric theory, the operators remain fluctuating and complex-valued objects, so no rigorous interpretation should be ascribed.

We shall return to consider in more detail the connection of these expressions to macroscopic dielectric theory in Section VII.

## III. ONE-COMPONENT DIELECTRIC BEAD FLUID

In illustrating the construction of a polarizable field theory in Sec. II, we were a bit cavalier about setting the monopole charge on all the *n* beads to the same value *q _{i}* =

*q*. Unless

*q*= 0, which corresponds to the case of a one-component fluid of dielectric beads, the fluid carries a net charge of

*nq*, which violates charge neutrality. Thus, to build a sensible electrolyte model, one needs a minimum of two types of beads carrying opposite charges or a neutralizing background. We shall defer that situation to Section V and consider in more detail here the

*q*= 0 case of a purely dielectric one-component bead fluid possessing induced dipoles.

The model of interest in this section corresponds to a homogeneous fluid in the canonical ensemble with the Hamiltonian given in Eq. (5) and single bead partition function given by

Within the *mean-field approximation*, which is obtained by finding the saddle point (*ϕ*^{∗}, *w*^{∗}) configurations of the fields from

the theory admits a homogeneous, purely imaginary solution in both fields. The free energy has no contributions from a spatially uniform *ϕ*^{∗}, so we choose *ϕ*^{∗} = 0 for convenience. An expression for *w*^{∗} follows from Eq. (28) and the property that smearing a constant with a Γ convolution recovers the same constant, resulting in

Within the mean-field approximation, the Helmholtz free energy is readily evaluated as *A* = − ln*Z*, and the pressure follows from Π = − ∂*A*/∂*V*|_{T,n}, leading to

The two terms in this mean-field equation of state are, respectively, an ideal gas term and a mean-field pairwise interaction term.

Beyond the mean-field approximation, one can conduct systematic loop expansions^{55} to assess the role of field fluctuations about the homogeneous saddle point. Here, we pursue the simple Gaussian approximation, which is obtained by expanding *H* to quadratic order in *ϕ* and *ω* = *w* − *w*^{∗}. As detailed elsewhere,^{52} the functional integral defined by Eq. (4) can be reexpressed as a Gaussian integral over each Fourier mode of the field fluctuations

where $ \varphi \u02c6 ( k ) $ and $ \omega \u02c6 ( k ) $ are the discrete Fourier modes of the field fluctuations, and

is a Hamiltonian associated with *Gaussian-order* field fluctuations. The objects $ G \u02c6 \varphi \varphi $ and $ G \u02c6 w w $ give the Fourier representation of pair correlations (Green’s functions) associated with the two fields; their purely local wavevector dependence implies *independent contributions from each fluctuation mode*. For the various models treated in the present paper using the Gaussian approximation, the partition function is given in each case by Eqs. (38) and (39), and differs only in the definition of the average bead density, *ρ*_{0}, and in the form of the Green’s functions, $ G \u02c6 \varphi \varphi $ and $ G \u02c6 w w $. For the immediate case of a single-component fluid of dielectric bead molecules, the latter are given by

and

with $ \u03f5 \u02c6 ( k ) $ a *wavevector-dependent dielectric constant*

and $ \gamma \u02c6 ( k ) $ a wavevector-dependent interaction kernel

Here, $ \Gamma \u02c6 ( k ) =exp [ \u2212 ( k a ) 2 / 2 ] $ is the three-dimensional Fourier transform of the smearing function Γ(**r**) with the property $ \Gamma \u02c6 ( 0 ) =1$. Carets throughout denote Fourier-transformed functions.

For *k* → 0, or long wavelengths, we see that $ G \u02c6 \varphi \varphi \u223c1/ ( \u03f5 k 2 ) $, which is the Green’s function of a dielectric medium with dielectric constant *ϵ* given by

Thus, in the long-wavelength limit we see an enhancement of the vacuum permittivity, 1, by a term proportional to the product of the fluid bead density and the bead polarizability. Correspondingly, at length scales less than the bead diameter *a*, i.e., *ka* ≫ 1, $ \u03f5 \u02c6 ( k ) \u22481$, which recovers the pure vacuum dielectric behavior. Our model therefore exhibits *dielectric screening*; in real-space at scales less than *a* we observe the vacuum Coulomb response of the potential away from a point charge at the origin, i.e., *G*_{ϕϕ}(*r*) ∼ 1/*r*. At scales much larger than *a*, the decay of *G*_{ϕϕ}(*r*) ∼ 1/(*ϵr*) remains algebraic, but is reduced in amplitude by a factor of 1/*ϵ* due to screening by the dielectric medium. Finally, we note that the expression for $ \u03f5 \u02c6 ( k ) $ given in Eq. (42) is easily recovered for the special case of a homogeneous fluid by averaging and Fourier transforming the more general expression Eq. (33).

The integrals in Eq. (38) can be evaluated exactly, and the thermodynamic connection formula to the pressure Π applied, leading to the following Gaussian-order equation of state:

The two integrals appearing in Eq. (45) represent excess contributions (beyond the ideal gas and mean-field terms) to the pressure associated with a Gaussian treatment of the *w* and *ϕ* fluctuations. The integrals are straightforward to evaluate numerically for specified values of *uρ*_{0} and *αρ*_{0}, but here we restrict our analysis to an expansion in powers of *ρ*_{0}, retaining only the leading $O ( \rho 0 2 ) $ terms. At this order, the integrals can be performed analytically; the first integral providing a fluctuation correction to the mean-field interaction term, while the second integral results in new physics not present in the mean-field approximation—the emergence of attractive *fluctuating dipole-induced dipole interactions*. We obtain

The final term on the first line is a first correction to the mean-field excluded volume interaction associated with *w* field fluctuations. It is small in relative importance for *u*/*a*^{3} ≲ 1. Of more interest is the first term on the second line, which reflects attractive dipole-dipole interactions created by *ϕ* fluctuations at the Gaussian level. This has the form of a second virial contribution with virial coefficient *B*_{2} scaling as *B*_{2} ∼ − *α*^{2}/*a*^{3}. The scaling can be understood by regularizing the usual formula for *B*_{2} with a cutoff at *r* = *a*,

and inserting the London formula^{56} *w _{DD}*(

*r*) ∼ −

*α*

^{2}/

*r*

^{6}for the attractive VDW interactions between a pair of orientationally averaged fluctuating induced dipoles. The latter formula is derived in the Appendix using the Drude model adopted here.

It is important to note that the VDW interactions emerging from our theory arise from thermal fluctuations of the classical Drude dipoles and thus must vanish as *T* → 0. In contrast, the observed London interaction for nonpolar liquids and gases is a *quantum effect* that is non-vanishing for *T* → 0.^{56,57} Nevertheless, as demonstrated by the derivation in the Appendix, our simple Drude model correctly infers the scaling of the London interaction and thus can be imposed phenomenologically, provided its temperature dependence is suppressed and *α* is specified to enforce consistency with the temperature-independent *C*_{6} coefficient of the desired chemistry.^{58} Such an approach is akin to the procedure employed in a recent set of papers,^{59,60} wherein a “quantum Drude oscillator” model of polarizable molecules is permitted to undergo fluctuations within a prescribed inhomogeneous electric potential (informed by an instantaneous molecular dynamics configuration), recovering London interactions due to dipole and higher-order correlations. We note that the latter can be recovered in the present theory by reverting to the unapproximated form of the single-molecule induced-dipole partition function, Eq. (9), which retains all multipole moments. For the permanent dipole case, we arrive at an analogous electrostatic contribution to the pressure, with *α* replaced by *μ*^{2}/3. Evidently, this is the Keesom interaction, which arises from thermal fluctuations in molecular species possessing a permanent dipole.^{56,61}

In summary, as previously advertised, our field theory naturally embeds attractive dipole-dipole interactions among beads created by fluctuations in the electrostatic potential field *ϕ*.

## IV. TWO-COMPONENT DIELECTRIC BEAD FLUID

Here we generalize the theory of Sec. III to the case of a binary mixture of type A and B polarizable bead molecules with respective polarizabilities *α _{A}* and

*α*and no monopoles (

_{B}*q*=

_{A}*q*= 0). For simplicity, we work in the canonical ensemble and assume that all molecules have the same soft mutual excluded volume interaction and smearing function Γ(

_{B}**r**) as in Sec. III. The Hamiltonian is given by

where *n _{A}* and

*n*are the numbers of molecules of each component and the single-molecule partition functions

_{B}*Q*are given by Eq. (34) with

_{j}*α*replaced by

*α*.

_{j}We now repeat the Gaussian fluctuation analysis of Sec. III, assuming a homogeneous fluid. At the mean-field level, again only the *w* field has a thermodynamically relevant saddle point value, which is given by *iw*^{∗} = *u*(*ρ _{A}* +

*ρ*), where

_{B}*ρ*=

_{j}*n*/

_{j}*V*denotes the average number density of the

*j*th species. The corresponding mean-field pressure of the mixture is given by Eq. (37) with

*ρ*

_{0}interpreted here as the total bead density

*ρ*+

_{A}*ρ*. At the Gaussian level of fluctuations, the

_{B}*ϕ*and

*w*modes are again independently distributed and characterized by two Green’s functions $ G \u02c6 w w ( k ) $ and $ G \u02c6 \varphi \varphi ( k ) $. The

*w*-propagator again has the form of Eq. (41) where $ \gamma \u02c6 ( k ) $ is given by Eq. (43), but with

*ρ*

_{0}=

*ρ*+

_{A}*ρ*. Likewise, $ G \u02c6 \varphi \varphi ( k ) $ for the binary mixture can be expressed as Eq. (40), but with a modified expression for the

_{B}*k*-dependent dielectric constant

Thus, as in the one-component case, the electrostatic correlations are described by a screened Coulomb interaction, but here the dielectric profile is determined by a density-weighted contribution from the polarizability of each component. This is analogous to the linear constitutive model used by Refs. 36–39 to model non-uniform dielectric profiles in microphase separated block copolymers. We shall explore this analogy in more detail in Section VII B.

It follows that to second order in the species densities, the pressure is generalized from Eq. (46) to the two-component expression

where the attractive term on the second line is the contribution from VDW interactions. This formula can be evidently generalized to a *K*-component dielectric fluid by simply extending the sums from 2 to *K* components. In the two-component case, it is instructive to re-express Eq. (50) in terms of the total density *ρ*_{0} = *ρ _{A}* +

*ρ*and the density difference

_{B}*ρ*

_{−}=

*ρ*−

_{A}*ρ*,

_{B} The final term of the first line of this expression is a VDW attraction proportional to $ ( \alpha A + \alpha B ) 2 \rho 0 2 $ that does not distinguish species and competes against the repulsive excluded volume interaction. The remaining attractive VDW interaction terms appear on the second line; the first of these, proportional to $ ( \alpha A 2 \u2212 \alpha B 2 ) \rho 0 \rho \u2212 $, would not cause a demixing transition of A and B since it is first order in *ρ*_{−}. However, the final VDW term in the second line proportional to $ ( \alpha A \u2212 \alpha B ) 2 \rho \u2212 2 $ can produce phase separation of A and B for sufficiently large |*α _{A}* −

*α*|. Equating this pressure contribution to the corresponding interaction term $\u2212 ( 1 / 4 ) v\chi \rho \u2212 2 $ in a Flory-Huggins or regular solution model

_{B}^{62,63}leads to the following expression for an

*effective Flory interaction parameter*

*χ*between A and B beads:

Here *v* is a reference lattice site volume that one might choose to be *O*(1/*ρ*_{0}).

In summary, we find that Flory-Huggins type interaction terms are naturally generated from *ϕ* field fluctuations in a multi-component polarizable fluid. For a two-component, purely dielectric bead fluid, a demixing phase transition can occur when *χ* ≳ 1 or |*α _{A}* −

*α*| ≳ (

_{B}*a*

^{3}/

*ρ*

_{0})

^{1/2}. In physical terms, phase separation creates favorable VDW interactions among the more polarizable bead molecules that more than compensates the loss of translational mixing entropy.

## V. TWO-COMPONENT ELECTROLYTE BEAD FLUID

Here we further generalize the model of Sec. IV by permitting the A and B beads to carry net charges of *q _{A}* and

*q*, respectively. We again work in the canonical ensemble and impose the constraint of electroneutrality,

_{B}*q*+

_{A}n_{A}*q*= 0. The Hamiltonian of Eq. (48) is unchanged, but the monopole charges

_{B}n_{B}*q*are restored in the single-bead partition functions

_{j}Repeating the Gaussian analysis of Secs. III and IV with the assumption of a homogeneous fluid, we find that the *ϕ* and *w* fluctuations remain uncorrelated at Gaussian order and that $ G \u02c6 w w ( k ) $ is unchanged from the charge-free case of Eq. (41). However, the monopole terms included in the present electrolyte model fundamentally modify the electrostatic Green’s function to a Debye-Hückel form^{64}

where $ \kappa \u02c6 ( k ) $ is a wavevector-dependent inverse Debye screening length defined by

The wavevector-dependent dielectric permittivity $ \u03f5 \u02c6 ( k ) $ appearing in these expressions is unchanged from Eq. (49).

For large wave vectors, $ \kappa \u02c6 ( k ) $ vanishes exponentially and $ \u03f5 \u02c6 ( k ) $ approaches 1, so $ G \u02c6 \varphi \varphi ( k ) $ reduces to the bare Coulomb Green’s function ∼1/*k*^{2}. In real space, this corresponds to *G*_{ϕϕ}(**r**) ∼ 1/*r*, indicating no dielectric interference or Debye-Hückel screening of electrostatic interactions on small length-scales below *a*. In the opposite limit of small wave vectors, $ \kappa \u02c6 ( k ) $ tends to the bulk inverse Debye screening length of the fluid: $\kappa \u2261 \kappa \u02c6 ( 0 ) = [ 4 \pi ( q A 2 \rho A + q B 2 \rho B ) / \u03f5 ] 1 / 2 $, with $\u03f5\u2261 \u03f5 \u02c6 ( 0 ) =1+4\pi ( \alpha A \rho A + \alpha B \rho B ) $ the macroscopic dielectric constant. In this regime, $ G \u02c6 \varphi \varphi ( k ) \u223c1/ [ \u03f5 ( k 2 + \kappa 2 ) ] $, which in real-space corresponds to long wavelength electrostatic screening both by the polarizable medium and correlated “Debye clouds” of oppositely charged ions, *G*_{ϕϕ}(**r**) ∼ exp(−*κr*)/(*ϵr*).

The equation of state of the polarizable binary electrolyte model has the same ideal gas and VDW terms as appearing in Eq. (50) for the dielectric fluid mixture. However, the modified form of $ G \u02c6 \varphi \varphi ( k ) $ produces additional electrostatic contributions to the pressure in the electrolyte case. These can be written as

This expression can be re-expressed as a sum of two terms, Π_{q} = Π_{DH} + Π_{ID}, with

and

For dilute solutions, the first integral, Π_{DH}, is dominated by small wavevectors *k* ∼ *κ* much less than 1/*a*. In this case, $ \kappa \u02c6 ( k ) $ can be approximated by $\kappa = \kappa \u02c6 ( 0 ) $ in the integrand and the integral performed exactly

Equation (59) is the classic *Debye-Hückel* result for a dilute electrolyte,^{64} representing the contribution to the pressure from screened monopole-monopole interactions. The expression is $O ( \rho i 3 / 2 ) $ and is thus the dominant correction to the ideal gas equation of state in the limit of very dilute fluids. The VDW and excluded volume terms in Eq. (50) are smaller by $O ( \rho i 1 / 2 ) $ in that regime. It is also important to note that *κ* depends on the macroscopic dielectric constant *ϵ*, so the polarizable nature of the charged molecules is implicitly included in this Debye-Hückel formula. Finally, we remark that a more realistic model of an electrolyte solution would include a polarizable (but uncharged) third “solvent” component that would further contribute to *ϵ*.

The second pressure term arising from monopole charges, Π_{ID}, is defined by an integral whose dominant contribution is from large wavevectors with *k* ∼ 1/*a*. As a consequence, the denominator of the integrand can be approximated by $ \u03f5 \u02c6 ( k ) k 2 $, and the integral performed analytically

where we only report the leading order $ \rho i 2 $ terms. The higher order contributions arise from the density dependence of $ \u03f5 \u02c6 ( k ) $.

This *ion-dipole* pressure contribution evidently arises from the net attraction between the monopoles on bead molecules and the fluctuating induced dipoles of surrounding molecules in the fluid. Like the VDW term in Eq. (50), Π_{ID} depends on the smearing scale *a*, is quadratic in species density, and contains a factor of $ \u2211 j \alpha j \rho j $. However, Π_{ID} is distinguished by its linear dependence on the ionic strength $ \u2211 j q j 2 \rho j $. We show in the Appendix that these scalings are consistent with the ion-dipole contributions to the second virial coefficients for a binary gas of polarizable and charged hard spheres with diameter *a*.

### A. Electrostatic self-energies

The analysis of this section has been restricted to *homogeneous* fluids of charged and polarizable molecules. In such a system, the electrostatic self-energies of the monopoles and dipoles contribute reference terms to the chemical potentials of each species, but do not influence thermodynamic properties such as the pressure. In contrast, for *inhomogeneous* fluids, the self-energy terms are critical in establishing the equilibrium spatial distribution of ions as well as dictating properties such as interfacial tension. While the self-energy terms of the present homogeneous fluid are not of thermodynamic significance, their form is helpful for understanding the structure of the present theory. Here we show that the two-component electrolyte model has electrostatic self-energy contributions similar to those deduced by Wang,^{42} but generalized to the case of polarizable molecules.

We can evaluate the electrostatic self-energy of a bead molecule of type *j* by considering its chemical potential in the Gaussian approximation

Here, the first two terms give the ideal gas and mean-field excluded volume contributions, respectively, to the bead chemical potential. The last term of the expression, *v _{j}*, gives the fluctuation contribution of the excluded volume interaction; its form is omitted for brevity. The remaining term,

*u*, is then the fluctuation contribution of electrostatics to the chemical potential, which we define to be the

_{j}*electrostatic self-energy*. For a homogeneous system, it is given by the integral

This integral can be immediately decomposed into a sum of self-energy contributions *u*_{j,M} + *u*_{j,D} arising, respectively, from the monopole and induced dipole residing on the molecule. The monopole contribution can be conveniently expressed as

As with the analysis of the pressure, the first integral in this equation is dominated by small wavevectors *k* ∼ *κ* such that $ \u03f5 \u02c6 $ and $ \kappa \u02c6 $ can be replaced by their macroscopic values. In contrast, the second integral is dominated by large wavevectors, *k* ∼ 1/*a*. Performing the resulting integrals leads to an expression for the monopole contribution to the self-energy

where *ϵ _{M}* is an average dielectric constant defined by

The parameter *ϵ _{M}* differs from the bulk dielectric constant

*ϵ*because of the loss of dielectric shielding below the

*a*scale; to first order in

*ρ*it can be approximated by

_{i} the second term differing from the corresponding term in *ϵ* by a factor of $ 2 $. Overall, Eq. (64) has the same structure as revealed by Wang:^{42} the self-energy of a monopole in an electrolyte solution can be expressed as the sum of a *Born solvation energy*,^{44} proportional to $ q j 2 / ( \u03f5 M a ) $ and a term proportional to $ q j 2 \kappa /\u03f5$ describing the average electrostatic interaction of the ion with the Debye cloud surrounding it. In Wang’s work, *ϵ _{M}* was replaced by

*ϵ*because the macroscopic dielectric permittivity was imposed at all scales, rather than self-consistently computed.

We now turn to the self-energy contribution from the induced dipole of a type-*j* molecule,

This integral has a complicated dependence on the concentration of molecules through the *ρ _{i}* dependence of both $ \u03f5 \u02c6 $ and $ \kappa \u02c6 $, reflecting the average interaction of the dipole with the monopoles and induced dipoles surrounding it. At leading order for

*ρ*→ 0, the dominant contribution comes from large wavevectors

_{i}*k*∼ 1/

*a*, in which case we can neglect $ \kappa \u02c6 2 $ and approximate $ \u03f5 \u02c6 \u22481$ in the denominator of the integrand. This leads to the expression

In summary, in the dilute regime, the dominant contributions to the electrostatic self-energy will be the Born solvation contribution to *u*_{j,M} and the leading order dipole term (Eq. (68)), as they are both $O ( \rho i 0 ) $. Consequently, for an *inhomogeneous* phase-separated system, we can expect charged beads to accumulate in regions of locally higher dielectric constant, so as to minimize the sum of the two terms.

## VI. EXTENSION TO POLYMERIC FLUIDS

The models presented in Secs. II–V are readily extended to fluids containing polarizable polymer molecules of a broad range of architectures. The approach is simply to connect beads with induced or permanent dipoles (and optional monopoles) into chains using rigid or spring linkages. Here we illustrate the approach with a one-component fluid of *n* induced-dipole bead chains, each chain containing *N* + 1 beads with bead polarizability *α*, a monopole of valence *q*, and a soft excluded volume interaction of strength *u*. In the canonical ensemble, the Hamiltonian is given by Eq. (5), where $Q [ i \varphi \u0304 , i w \u0304 ] $ is now a *single-chain partition function* given by

with *N* the label given to the last [(*N* + 1) th] bead in the chain, the first bead being labeled 0. The object $ q p ( r , j ; [ i \varphi \u0304 , i w \u0304 ] ) $ is a “chain propagator,” which represents the spatially varying statistical weight that the end of a chain containing *j* + 1 beads is at position **r**. It is defined recursively by

with “initial condition”

Equation (70) reflects the Markov property of chain molecules in an external field; the statistical weight of a chain with *j* + 1 beads can be built up from the corresponding weight for a chain with one fewer bead. Φ(Δ**r**; **r**′) is a bond transition probability that defines the type of bead linker used in the model.^{49} It specifies the probability density that a bead at position **r**′ is linked by a bond vector Δ**r** = **r** − **r**′ to a subsequent bead at position **r**. The model specified by Eq. (70) is quite general, although it evidently describes chains with *isotropic polarizability* in that the dielectric properties of each bead are independent of the adjacent bond orientations.

One such model is the “freely-jointed chain,” where the bond length between two adjacent beads is fixed to a value *b*; for such a model, the bond transition probability function is given by

Another popular model for bead connectivity assumes a harmonic potential between adjacent beads, referred to as the “discrete Gaussian chain.” The bond transition probability in this case is

Crucially, Φ(Δ**r**; **r**′) for both models is a function only of Δ**r** and is independent of **r**′. The integral appearing in Eq. (70) is thus of convolution form and can be efficiently computed using Fourier methods.

Finally, a “continuous Gaussian chain” model can be obtained by defining a chain of *N _{s}* + 1 beads, replacing

*b*

^{2}in the Gaussian linker expression above by

*b*

^{2}(

*N*/

*N*), and taking the limit of

_{s}*N*→ ∞.

_{s}^{49}The overall result is to replace Eq. (70) with a diffusion-like partial differential equation

where we have suppressed functional dependence on the fields *ϕ* and *w* in the notation for *q _{p}*. Here

*s*∈ [0,

*N*] is a continuous segment index that substitutes for the discrete index

*j*of the previous model. The initial condition of Eq. (71) is further simplified to $ q p ( r , 0 ; [ i \varphi \u0304 , i w \u0304 ] ) =1$. Eq. (74) has the virtue of being linear and is amenable to asymptotic analysis or numerical solution using a variety of techniques.

^{49}

### A. Binary dielectric polymer blend

To illustrate the types and behaviors of polarizable polymer models that can be constructed using the above formalism, we consider a binary blend of *n _{A}* chains of type

*A*and

*n*chains of type

_{B}*B*. We adopt the continuous Gaussian chain model for both species and assign them isotropic segmental polarizabilities

*α*,

_{A}*α*, no monopoles (

_{B}*q*=

_{A}*q*= 0), and polymerization indices

_{B}*N*,

_{A}*N*, respectively. Although we omit monopoles in the present analysis, it is straightforward to extend the model to include segment monopoles; indeed, such systems have previously been considered (using both Gaussian fluctuation approximate and fully fluctuating FTS techniques) in the absence of segment polarizability.

_{B}^{15,16,65}If we further assume that segments of all types have a mutual soft excluded volume repulsion of strength

*u*, the canonical ensemble Hamiltonian of the present blend is identical in form to that of the two-component small molecule system, Eq. (48). Here the single-chain partition function for species

*j*is given by

where *q*_{p,j} satisfies Eq. (74) with *α* replaced by *α _{j}* and

*q*= 0.

To explore the properties of the dielectric blend, we again perform a Gaussian fluctuation analysis about the homogeneous state, expanding *Q _{j}* to quadratic order in

*w*and

*ϕ*, such that fluctuations are approximated as Gaussian. The Green’s function associated with the excluded volume interaction is given by

where $ \gamma \u02c6 c ( k ) $ is a wavevector-dependent interaction kernel analogous to the one given by Eq. (43), with additional contributions arising from the conformational entropy of the chains. It is defined as

where *ρ _{j}* ≡

*n*/

_{j}N_{j}*V*is the

*monomer density*of the

*j*th component,

*R*

_{g,j}≡ (

*N*/6)

_{j}^{1/2}

*b*is the bare radius of gyration of chains of type

*j*, and $ g \u02c6 D ( x ) \u2261 ( 2 / x 2 ) ( e \u2212 x + x \u2212 1 ) $ is the Debye function.

The contribution of the excluded volume interaction to the pressure of the blend is thus given by the expected mean-field term and an integral determining the Gaussian fluctuation contribution

In examining the asymptotics of the integrand of the fluctuation correction term, it is clear that the wavevector-dependence of $ \gamma \u02c6 c $ indicates three competing length scales: *a*, *R*_{g,A}, and *R*_{g,B}. Fluctuation modes with wavevectors that are very large (*k* ≳ 1/*a*) are exponentially suppressed in the integrand. Assuming that *R*_{g,A} ∼ *R*_{g,B} ≫ *a*, we find that the integral is dominated by modes with wavevectors in the intermediate range 1/*R*_{g,j} ≪ *k* ≪ 1/*a*, where the approximation $ g \u02c6 D ( x ) \u22482/x$ applies. With this replacement, the integral can be evaluated analytically

This is the two-component version of a well-known result due to Edwards and Olvera de la Cruz.^{66,67}

The corresponding electrostatic Green’s function for the blend is given by Eq. (40), with a wavevector-dependent dielectric constant given by Eq. (49), where *ρ _{A}* and

*ρ*are now monomer densities. It follows that to Gaussian order, the electrostatic contribution to the pressure of a homogeneous binary mixture of polarizable chains is

_{B}*insensitive to chain connectivity*and is given by the VDW term appearing in Eq. (50) for a mixture of polarizable small molecules. Moreover, the driving force for phase separation is the same in both cases, with an effective Flory parameter given by Eq. (52). An important distinction, however, is that the ideal gas terms resisting phase separation are smaller in the polymer blend system, Π

_{ig}=

*ρ*/

_{A}*N*+

_{A}*ρ*/

_{B}*N*, reflecting a reduced translational entropy due to bead/monomer connectivity.

_{B}### B. Anisotropically polarizable polymer chains

While the above models of isotropically polarizable chains have the benefit of simplicity, real polymers can respond anisotropically to an applied electric field. The most straightforward way to account for such a possibility is to consider the connecting *bonds* to be polarizable, in addition to the beads. For a bond oriented along Δ**r**, one can define a polarizability tensor ** α**(Δ

**r**), which describes the linear response of the segment polarization to the local electric field. We assume a form for bonds possessing cylindrical symmetry

where **I** is the unit tensor, *α*_{∥} describes the polarizability of the segment along its bond vector, and *α*_{⊥} is the polarizability perpendicular to the bond vector. Equivalently, it is convenient to introduce two linear combinations *α* = *α*_{⊥} and *ν* = *α*_{∥} − *α*_{⊥} that specify the isotropic and anisotropic contributions, respectively, to the polarizability tensor.

In assigning polarizability to the bonds of a polymer, we assume that the electrostatic potential energy of a discrete chain, subject to the effective electrostatic potential field $i \varphi \u0304 $, is given by

If we group the isotropic bond contributions with those of the beads, the resulting chain statistics are identical to Eq. (70), with the exception that Φ(Δ**r**; **r**′) is replaced by an anisotropic bond transition probability, given by

where Φ is the bond transition probability of the chain model in the absence of electrostatics.

An unfortunate feature of Eq. (82) is its dependence on the local electric field, resulting in the integral in Eq. (70) no longer having the form of a convolution. The usual computational approach of attacking the integral pseudo-spectrally and exploiting fast Fourier transforms would thus fail, as the integral is non-local in both real and Fourier spaces. Instead, we assume weak anisotropy of ** α** and expand Φ

_{e}to first order in

*ν*,

If Φ possesses translational invariance, i.e., is a function of Δ**r** but is independent of **r**′, as in the freely jointed and discrete Gaussian chain models, this expression renders Eq. (70) a convolution form. The result, valid to first order in *ν*, is

Here $ q \u02c6 p ( k , j ) $ is the Fourier transform of *q _{p}*(

**r**,

*j*), defined by ∫

*d*

**r**exp(−

*i*

**k**⋅

**r**)

*q*(

_{p}**r**,

*j*), and $ F \u2212 1 $ denotes the inverse Fourier transform. Provided a closed form exists for $ \Phi \u02c6 ( k ) $ and its dyadic $ \u2207 k \u2207 k \Phi \u02c6 $, both terms in Eq. (84) can be efficiently evaluated by a Fourier transform of

*q*, a local multiplication in reciprocal space, an inverse Fourier transform, and a local multiplication in real space.

_{p}The introduction of molecular species with anisotropic polarizability implies an extension of the framework introduced in Section II to fluids for which the polarization operator is not necessarily parallel to the local electric field. We can define an unsmeared, tensor-valued *microscopic polarizability density* for a specified configuration of an anisotropically polarizable chain

This object gives the contribution of a representative single chain to the local polarizability of space near **r**, given chain configuration **r**^{N+1}. Equation (81) for the electrostatic potential energy of the chain can be re-expressed in terms of this microscopic density

Evidently, the microscopic polarizability density is thermodynamically conjugate to the dyadic product $ ( 1 / 2 ) \u2207 \varphi \u0304 \u2207 \varphi \u0304 $. It follows that we can define a field-theoretic *polarizability density operator* for a fluid of *n* anisotropically polarizable chains by direct functional differentiation of the single chain partition function

To efficiently evaluate this object, the weak anisotropy approximation (Eq. (83)) produces an expression that contains convolution integrals suitable for pseudo-spectral evaluation

where

is the bead density operator for a collection of *n* discrete chains, derived elsewhere.^{49}

For polymer systems possessing anisotropic polarizability, a generalization of the linear dependence between the polarization operator $ P \u0303 ( r ; [ \varphi , w ] ) $ and the local internal electric field $ E \u0303 ( r ; [ \varphi ] ) \u2261\u2212i\u2207 \varphi \u0304 ( r ) $ (Eq. (21)) emerges

In the limit of vanishing *ν*, $ \alpha \u0303 $ is purely isotropic and Eq. (21) is immediately recovered.

The expressions following Eq. (22) are similarly altered for an anisotropic polymer fluid, such that all instances of $\alpha \rho \u0303 $ are replaced by the tensor $ \alpha \u0303 $. The *ϕ*-force for a fluid of anisotropic chains is again given by Eqs. (30) and (31), where $ P \u0303 $ is now the anisotropic expression given in Eq. (90). Evidently, such a fluid will possess a nonlocal, *anisotropic dielectric permittivity* operator defined by

## VII. EXTERNAL ELECTRIC FIELDS

The formalism developed in this paper is also a suitable framework for studying the effect of externally applied electric fields on self-assembly behavior, thermodynamics, and structure of soft polarizable materials. A variety of field-theoretic methods have been used in the past to study such phenomena,^{36–39,68,69} but they generally lack internal consistency because the dielectric constitutive behavior is externally supplied, rather than a consequence of the intrinsic polarizability of the molecules and particular configuration of the system. Moreover, many approaches, such as those used to describe electric-field effects on block copolymer self-assembly, invoke mean-field treatments of electrostatic and/or compositional fields. As a result, the structure of a more comprehensive framework allowing for full treatment of thermal fluctuations is masked and not self-evident.

For the sake of conciseness, we illustrate the approach to external fields for the one-component, induced-dipole bead fluid model of Section III, but the method can be easily generalized to any of the systems discussed above. An applied electric field **E**_{a}(**r**) is introduced as in the derivation of Eq. (20), namely by including a source term in the potential energy of the form

where **P**_{m} is the microscopic polarization density defined in Eq. (19). Tracing this term through the Hubbard-Stratonovich transform leads to a field theory of the same form as Eqs. (4) and (5), except that the single bead partition function *Q* now has explicit dependence on the applied field

where |⋯| denotes the magnitude of the vector, not the complex modulus. Recalling the definition, Eq. (22), of the local internal electric field operator, $ E \u0303 $, this expression can be rewritten as

It is important to understand the distinction between the two electric field objects appearing in this expression. **E**_{a}(**r**) is an applied field that would normally be generated by some charges external to the system, such as the charge maintained on the conductive plates of a simple capacitor. In a planar capacitor geometry, **E**_{a} is the field that would be present if there were no dielectric fluid between the plates; this is the **D** field of continuum dielectric theory. As **E**_{a} is determined externally, it does not depend on the fluid degrees of freedom manifest in the configurations of the *ϕ* and *w* fields. In contrast, $ E \u0303 ( r ; [ \varphi ] ) $ is a fluctuating internal electric field that reflects the thermal motion of the molecules/beads in the fluid and their fluctuating induced dipoles, as well as the electrostatic correlations present in the system. Complicating this interpretation, however, is the fact that $ E \u0303 $ is a complex-valued object for *ϕ*(**r**) evaluated on the real path of functional integration. Nonetheless, its average value is real, both in the mean-field approximation and the exact theory.

In the presence of an applied field, we can generalize Eq. (21) for the polarization density field operator, $ P \u0303 $, by means of the following relationships:

Explicit evaluation of the functional derivative in the middle expression using the field-theoretic representation of the model leads to

where $ \rho \u0303 $ is a generalization of the bead density operator defined in Eq. (17) to include an applied field

Equation (96) has an appealing, although non-rigorous, physical interpretation. The local polarization density in a bead fluid subject to an applied field is proportional to the *total* local electric field, i.e., the sum of the internal and external fields, $ E \u0303 + E a $, multiplied by a local isotropic dielectric susceptibility given by $\alpha \rho \u0303 ( r ) $. Nonetheless, the apparent linear dependence on particle density *ρ*_{0} and applied field **E**_{a}, as well as the locality of Eq. (96), can be misleading because the fluctuation spectrum of the *ϕ* and *w* fields is dependent on *ρ*_{0} and **E**_{a} and the fields are spatially correlated.

### A. Continuum dielectric theory and the Clausius-Mossotti (CM) formula

The operator and force expressions provided in the present section are *exact* relations that can be generalized to provide a broad foundation for studying external electric field effects in soft matter systems ranging from block copolymers and blends to polymerized ionic liquids. Here we show that a mean-field approximation to the molecularly founded field theory can recover a classic expression relating molecular polarizability to the continuum dielectric permittivity—the so-called *Clausius-Mossotti formula*.^{70}

Our starting point is Eq. (96), which is averaged to obtain

Following the spirit of Kirkwood’s molecular derivation of the Clausius-Mossotti formula,^{71} we make a mean-field approximation by breaking the average in the first term on the right hand side of Eq. (98). Further noting that in a homogeneous fluid $\u3008 \rho \u0303 ( r ) \u3009= \rho 0 $, we obtain

Next, an expression is required for $\u3008 E \u0303 ( r ) \u3009$. Thermodynamic perturbation theory can be used to derive the following linear response formula, rigorous to first order in the amplitude of the applied field **E**_{a}:

where $ P \u0303 0 $ is the polarization density operator in the absence of an external field, i.e., Eq. (21), and 〈⋯〉_{0} denotes an average given by Eq. (15) with **E**_{a} = 0. Specializing to the case of a uniform applied field, this expression can be rewritten as

where the tensor **T** is defined by

and we have used the identity ∇Γ(**r** − **r**′) = − ∇′Γ(**r** − **r**′), which follows from the spherical symmetry of the charge smearing function. The formula for **T** can be approximately evaluated by replacing the field-free Green’s function by the long-wavelength expression

where *ϵ* is the macroscopic dielectric constant. In the limit of vanishing smearing range, *a* → 0, **T** simplifies considerably to

an integral that has previously been considered by Kirkwood^{71} describing the average interaction between a pair of point dipoles separated by vector **r**. Here the domain of integration is critical; we consider a plane condenser configuration with conductor surfaces oriented normal to $ z \u02c6 $ confining the dielectric. In Kirkwood’s molecular theory, he argues that a sphere of infinitesimal radius can be excluded about the origin from the integration domain, since hard-core interactions between molecules would give such a separation on statistical weight. A version of Green’s theorem can then be used to replace the volume integral in Eq. (104) by a sum of surface integrals over the two faces of the dielectric slab and the infinitesimal sphere around the origin. This leads to

It is important to note that neither the radius of the excluding sphere nor the thickness of the dielectric slab appear in this expression. Returning to Eq. (102) at finite *a*, and once the integrals over **r**_{1} and **r**_{2} have been performed, the smearing renders the integrand finite at **r**′ = 0. As a consequence, we can again exclude an infinitesimal sphere about the origin from the domain of integration without changing the value of the integral. While we have no closed form solution in this case, the unsmeared result in Eq. (105) should be a good approximation for small values of *a*. Combining Eqs. (101) and (105) thus leads to

The first term in this expression can be physically interpreted as the field produced by the polarization charge created on the two planar surfaces of the dielectric, while the second is the familiar Lorentz field^{72} that a polarized dielectric exerts on a cavity. As a next step, we identify **E**_{a} with the continuum dielectric displacement **D** and $\u3008 P \u0303 \u3009$ with the continuum polarization density **P** and use the relation **D** = *ϵ***E** to combine Eqs. (99) and (106) to obtain

Finally, the continuum relations *D _{z}* =

*ϵE*and

_{z}*P*= (

_{z}*ϵ*− 1)

*E*/(4

_{z}*π*) are used to eliminate an overall factor of

*E*, recovering the familiar Clausius-Mossotti (CM) formula relating the macroscopic permittivity to the molecular polarizability and density of the fluid

_{z}^{71}

with *ζ* ≡ 4*παρ*_{0}/3.

To first order in *ζ* ∝ *αρ*_{0}, the CM formula agrees with the expression *ϵ* = 1 + 4*παρ*_{0} derived in Section III based on the Gaussian approximation. At the same order, it is also consistent with the average of Eq. (33) in a homogeneous fluid. Taylor expansion of the CM formula, however, predicts universal corrections to this result at *O*(*ζ*^{2}) and higher. These corrections are inexact due to the mean-field approximation employed; a full analysis would render the *O*(*ζ*^{2}) coefficients non-universal—dependent on both the smearing parameter *a* and the excluded volume parameter *u*.

### B. Polymer self-consistent field theory in an electric field

As a final exercise, we show the consistency of the present formalism in the mean-field approximation with the standard self-consistent field theory (SCFT) framework that is used to study the behavior of polymer blends and block copolymers in the presence of an external electric field.

For simplicity, we consider a binary blend of homopolymers with *isotropically polarizable segments* in the canonical ensemble, identical to the system considered in Section VI A. From the discussion of Sec. VII A, it is apparent that in the presence of an external electric field **E**_{a} the effective local field experienced by beads of chains of species *j* is given by

where we again interpret $\u2212i\u2207 \varphi \u0304 $ as the local internal electric field $ E \u0303 $. As a consequence, in the presence of an external field, $i w \u0304 + ( \alpha j / 2 ) |\u2207 \varphi \u0304 | 2 $ should be replaced by *ψ _{j}* in expressions such as Eqs. (70) and (71) determining the chain propagators and single-chain partition functions

*Q*.

_{j}In the mean-field or SCFT approximation, the stationary configurations of *ϕ* and *w* are obtained by setting the *w* and *ϕ* forces, analogous to Eqs. (28) and (29), to zero. In the case of the *w* force, one obtains

where $ \rho \u0303 j $ is the local density operator of monomer species *j* given by a suitably modified version of Eq. (89) for a discrete chain,

or an obvious generalization for a continuous chain.^{49}

For stationary *ϕ* force, we have correspondingly

where $ \u03f5 \u0303 ( r , r \u2032 ; [ \varphi , w , E a ] ) $ is a generalization of the non-local dielectric kernel given in Eq. (33), with parametric dependence on the applied electric field **E**_{a} carried through to all operators

Equations (110)-(113) represent the coupled SCFT equations that must be solved simultaneously for *w* and *ϕ* (or alternatively $ E \u0303 $) at a specified external field **E**_{a}. However, they simplify considerably in the unsmeared (*a* → 0), incompressible ($ \u2211 j \rho \u0303 j = \rho 0 $) polymer melt case, as is typical in the literature.^{36–39} Specifically, a component contribution to the dielectric permittivity, *ϵ _{j}*, can be identified as proportional to the component polarizability,

*ϵ*≡ (1/

_{j}*ρ*

_{0}) + 4

*πα*, such that the constitutive model of previous studies of dielectric behavior with SCFT is recovered. This simplification amounts to

_{j}where $\u03f5 ( r ; [ \varphi , w , E a ] ) =\u222bd r \u2032 \u03f5 \u0303 ( r , r \u2032 ; [ \varphi , w , E a ] ) $ is a nonuniform dielectric profile associated with the local density profile in the melt, given by

and the mean-field *iϕ* is identified with the real-valued, macroscopic electrostatic potential field. Equation (115) is the linear dielectric constitutive law normally employed in studies of external electric field effects on the microphase separation of block copolymers.^{36–39}

Finally, we note that it is straightforward to generalize such expressions to incorporate anisotropic polarizability of polymer chain segments as discussed in Section VI B.

## VIII. DISCUSSION AND OUTLOOK

The formalism described in Secs. II–VII is evidently a versatile framework for building statistical field theory models that can serve as the basis for analyzing the structure and thermodynamics of a broad range of polymer and soft matter systems, both at equilibrium and under stationary conditions in the presence of an external electric field. Our exposition was primarily aimed at demonstrating how microscopic models of polarizable simple and complex fluids can be constructed by embedding a Drude oscillator and, optionally, a monopole within an elementary bead unit. Charged and polarizable polymers, including those with permanent or induced dipoles and anisotropic segment polarizabilities, are then built up by connecting beads with rigid or spring linkers according to standard chain models. Beyond electrostatic interactions, excluded volume and other types of bead-bead attractions or repulsions are easily incorporated. Systems as varied as polymer alloys, polyelectrolytes, polyelectrolyte block and graft copolymers, polymerized ionic liquids, and ionomers are all within the scope of the framework.

Once a complete microscopic model has been posed in bead coordinates, its transformation to a field-theory proceeds by standard Hubbard-Stratonovich (HS) transforms. The auxiliary field *ϕ* introduced to decouple the electrostatic interactions can be viewed as a fluctuating electrostatic potential, and this is the field whose structure and correlations impart polarization effects, VDW interactions, and macroscopic dielectric and electrostatic behaviors on the medium. Other HS auxiliary fields, such as the *w* field introduced here, are responsible for dictating the compressibility of the fluid and the density patterns that emerge from confinement or phase separation.

It is worth commenting that in both the binary dielectric fluid model of Section IV and the dielectric polymer blend model of Section VI A, we did not explicitly include a net repulsive contact interaction between dissimilar bead species, i.e., a “Flory *χ* parameter.” Instead, we argued that the VDW interactions generated by the electrostatic correlations of the *ϕ* field automatically provide such a term of the form *χ* ∼ (*α _{A}* −

*α*)

_{B}^{2}. This is indeed the molecular origin of

*χ*(apart from packing entropy) in small molecule fluid and polymer mixtures, recognizing, however, that the true quantum nature of the

*α*discussed in Section III requires a phenomenological treatment. A disadvantage of such an approach is that no phase separation can take place within a mean-field level of treatment, since the VDW interactions only emerge when

_{j}*ϕ*fluctuations are included. In contrast, an artificially imposed

*χ*can produce phase separation in the mean-field approximation, e.g., within SCFT. We view this not as an argument against the formalism, but as providing more flexibility in modeling the approach. If, for example, one intends to simulate the full theory without approximation, the present approach of not specifying

*χ*, but only the segmental polarizabilities, has the benefit of elegance and parameter efficiency. The

*α*automatically determine an effective

_{j}*χ*, as well as effective permittivities of phase separated domains. Alternatively, if one intends an SCFT study, an external

*χ*parameter can be used to induce phase separation, but the dielectric properties of the domains will then be uniquely determined by the specified

*α*and the equilibrium compositions of the domains.

_{j}The field theories described here are exact representations of the underlying microscopic bead models and preserve information about molecular architecture and interactions. They are in general characterized by a complex-valued Hamiltonian functional, so the statistical weight, exp(−*H*[*ϕ*, *w*]), is a complex number and not a true Boltzmann weight. Such theories are said to possess a “sign problem.” For approximate analytical calculations, such as the Gaussian fluctuation analysis (about a homogeneous state) used to elucidate the behavior of various models throughout this manuscript, the sign problem is of no consequence. Nor is it a difficulty for analytical or numerical treatment of the field theories in the mean-field or SCFT approximation. However, in direct numerical attack of the exact field theories by computer simulation, one does encounter the sign problem. While it is nearly insurmountable in other contexts, such as fermionic quantum field theories, we expect that the sign problem for the present theories can be effectively overcome by the same complex Langevin (CL) simulation technique^{49–51} that has enabled a broad range of studies of strongly fluctuating polymer field theories, including models of polyelectrolyte complexation.^{15,16} Improvements in model formulation and the introduction of advanced numerical algorithms for stably and accurately time-stepping the CL equations^{52,53,73,74} suggest that the present polarizable field theory models will succumb to similar treatment.

We envision several branches of soft matter research that should benefit from the present model construction. One relates to the ordering or phase separation behavior of heterogeneous dielectric soft matter systems in external electric fields. Prior work on systems such as binary fluids, blends, and block copolymers has largely applied the mean-field approximation to molecularly based models^{36–39,68} or used approximate analytical techniques to analyze correlations between structure and electrostatics in phenomenological continuum models.^{69} Heretofore, the framework for building molecularly-faithful field-theoretic models beyond the mean-field approximation and enabling their direct, unapproximated numerical simulation did not exist. We anticipate reporting on progress in this line of research in the near future.

A second promising direction involves application of the formalism to ion-containing heterogeneous systems such as polyelectrolyte block copolymers, ionomers, and polymerized ionic liquids. In such systems, the equilibrium distribution of ions is intricately linked to the local dielectric environment; ions will tend to migrate to domains with the highest permittivity caused by local enrichment by highly polarizable segments or solvent molecules. The present framework is ideal for studying such phenomena, since the molecular entities that enter the theory, e.g., segments, solvents, and ions, are specified (in part) by prescribing their microscopic polarizabilities and charges. The collective influence of these parameters is felt in terms of the VDW interactions that are generated and serve to drive phase separation, as well as the dielectric contrast that emerges naturally as the molecules partition in domains. Within the present formalism, we have the convenience of a continuum description and access to parameters that specify molecular details, yet with exact embedded physics of coupled electrostatic, density, and compositional fluctuations in a heterogeneous soft material. Computationally demanding techniques for dealing with long-range Coulomb interactions in particle simulations of charged systems (e.g., Ewald sums) are further avoided, because only the inverse of the Coulomb operator appears in the present field theories.

A third area for future research relates to studies of mesoscopic phenomena such as coalescence, breakup, wetting, and dewetting of films, drops, and other objects in structured soft matter systems. Such phenomena are often dominated by long-range VDW interactions, yet the models and methods used for their study are limited. Atomistic or coarse-grained particle models that lack polarizability, yet include Lennard-Jones type interactions with a ∼1/*r*^{6} attractive tail, can in principle capture the phenomena, but simulations based on such models are often challenged by large systems, especially if they are mesostructured and contain polymers. SCFT models, in turn, lack VDW interactions completely and are restricted to the mean-field approximation. Continuum models can easily incorporate VDW interactions, but lack the ability to connect with molecular details and predict/embed mesostructure. Evidently, the present framework has none of these deficiencies, other than the sign problem that can in principle be effectively addressed with the complex Langevin simulation technique.

## Acknowledgments

J.M.M. was supported by the Institute for Collaborative Biotechnologies through Grant No. W911NF-09-D-0001-0042 from the U.S. Army Research Office. K.T.D. and G.H.F. were partially supported by the NSF DMR-CMMT Program under Award No. DMR-1506008 and by the Institute for Collaborative Biotechnologies through Grant No. W911NF-09-0001 from the U.S. Army Research Office. The content of the information presented here does not necessarily reflect the position or the policy of the U.S. Government, and no official endorsement shall be inferred. Extensive use was made of the computational facilities of the Center for Scientific Computing at the CNSI and MRL: an NSF MRSEC (Grant No. DMR-1121053) and NSF No. CNS-0960316.

### APPENDIX: SECOND VIRIAL COEFFICIENTS OF A POLARIZABLE IONIC HARD SPHERE GAS

In this appendix we investigate the contributions from ion-dipole and dipole-dipole interactions to the second virial coefficients of a dilute gas of charged and polarizable atoms with hard sphere repulsions. We demonstrate consistency in the scaling behavior of the virial coefficients with the $O ( \rho i 2 ) $ coefficients in the expansion of the equation of state for the field-theoretic model considered in Section V.

We consider a pair of polarizable and charged spherically symmetric molecules with hard sphere repulsions at diameter *a* in vacuum. Molecules of species *i* are assigned a polarizability *α _{i}* and an unbalanced (monopole) charge of

*q*. The ion-ion electrostatic interaction between the monopole charges on the two molecules is neglected, as this decays as 1/

_{i}*r*and makes a divergent contribution to the virial coefficients (but is accounted for in the Debye-Hückel contribution to the pressure). Bare dipole-dipole (van der Waals) and ion-dipole interactions, however, are sufficiently short-ranged for finite virial contributions.

The polarizability of each molecule is implemented with a Drude oscillator model, assigning two paired opposite partial charge values ±*δq _{i}* and their connecting spring constant $\delta q i 2 / \alpha i $, such that the internal energy of a molecule of type

*i*is given by

where **d** is the separation vector between the two partial charges.

The ion-dipole interaction between molecule 1, of type *i*, and molecule 2, of type *j*, can be determined by summing the Coulomb interactions of each molecule’s monopole with the opposing molecule’s paired partial charges

where **r**_{1} and **r**_{2} are the centers of mass of the two molecules and **d**_{1} and **d**_{2} are the separation vectors between the partial charges of each molecule, representing two dipoles.

Next, we invoke a multipole expansion truncated after *O*(**d**_{i}), assuming |**d**_{i}| ≪ |**r**_{1} − **r**_{2}|, to derive the instantaneous ion-dipole interaction

Averaging with a Boltzmann weight based on the sum of the Drude model internal energies given in Eq. (A1) over all possible polarizations (vectors **d**_{1}, **d**_{2}) for the two molecules, and shifting molecule 1 to the origin leads to the following orientationally averaged ion-dipole interaction:

with *r* = |**r**_{1} − **r**_{2}|. Unsurprisingly, the average ion-dipole interaction is attractive when the dipoles are allowed to fluctuate, as configurations with the dipole moment of each molecule aligned along the electric field generated by the monopole of the other molecule are highly preferred.

In a similar fashion, we can derive the orientationally averaged dipole-dipole interaction by starting with the Coulomb interactions of all of the paired partial charges

Again we employ a multipole expansion, recovering the instantaneous dipole-dipole interaction

Averaging over the polarization degrees of freedom and retaining only the most slowly decaying term recovers the familiar London dispersion interaction^{56} scaling as 1/*r*^{6},

Next, we consider the contribution of each orientationally averaged interaction to the second virial coefficient associated with the *ρ _{i}ρ_{j}* term of the pressure. For this purpose, we integrate (−1/2) times the Mayer-

*f*function, which we approximate as (

*e*

^{−wij}− 1 ≈ −

*w*), over all vector separations of the two molecules, with the cutoff |

^{ij}**r**

_{1}−

**r**

_{2}| ≥

*a*imposed by the assumption of hard sphere repulsions. This leads to the following ion-dipole and dipole-dipole contributions to the virial coefficients:

The contribution of each interaction to the pressure is quadratic in component densities, and given by $ \Pi e = \u2211 i \u2211 j B e i j \rho i \rho j $. For the case of a two-component fluid of polarizable ionic molecules, performing this sum explicitly recovers terms that differ from their analogs in the Gaussian-approximated field theory of Sections IV and V by simple prefactors; the latter arising from the difference in regularization procedures,

## REFERENCES

_{2}/light gas separations