We derive the effective Flory-Huggins parameter in polarizable polymeric systems, within a recently introduced polarizable field theory framework. The incorporation of bead polarizabilities in the model self-consistently embeds dielectric response, as well as van der Waals interactions. The latter generate a *χ* parameter (denoted $\chi \u0303$) between any two species with polarizability contrast. Using one-loop perturbation theory, we compute corrections to the structure factor $Sk$ and the dielectric function $\u03f5^(k)$ for a polarizable binary homopolymer blend in the one-phase region of the phase diagram. The electrostatic corrections to *S*(**k**) can be entirely accounted for by a renormalization of the excluded volume parameter *B* into three van der Waals-corrected parameters *B*_{AA}, *B*_{AB}, and *B*_{BB}, which then determine $\chi \u0303$. The one-loop theory not only enables the quantitative prediction of $\chi \u0303$ but also provides useful insight into the dependence of $\chi \u0303$ on the electrostatic environment (for example, its sensitivity to electrostatic screening). The unapproximated polarizable field theory is amenable to direct simulation via complex Langevin sampling, which we employ here to test the validity of the one-loop results. From simulations of *S*(**k**) and $\u03f5^(k)$ for a system of polarizable homopolymers, we find that the one-loop theory is best suited to high concentrations, where it performs very well. Finally, we measure $\chi \u0303N$ in simulations of a polarizable diblock copolymer melt and obtain excellent agreement with the one-loop theory. These constitute the first fully fluctuating simulations conducted within the polarizable field theory framework.

## I. INTRODUCTION

The electrostatic interactions between polymers and small molecules possessing polar, polarizable, and/or charged moieties are known to play a crucial role in determining structure and equilibrium properties. This interplay between electrostatics and structure in polymeric systems has become the subject of increased attention in recent years, due to the rapidly growing interest in charged polymers such as polyelectrolytes,^{1,2} ionomers,^{3,4} and polymerized ionic liquids,^{5,6} as well as salt-doped polymers.^{7–9} These systems are of interest for their potential applications in areas such as drug delivery and high-performance energy devices and as the basis for novel functional materials. Theoretical approaches to studying such materials must at least have the ability to incorporate ion-ion, ion-dipole, and dipole-dipole interactions into a molecular model, which for large, dense systems quickly becomes challenging to simulate by conventional methods.

For charge-neutral polymeric systems, phase behavior can be described using mean-field theories such as the random phase approximation (RPA),^{10,11} Flory-Huggins theory,^{12,13} and self-consistent field theory (SCFT).^{14,15} Typically, this involves the experimental determination of the Flory-Huggins interaction parameter *χ*, which characterizes the effective interaction between unlike chemical species. It is well known that the form of SCFT used in tandem with experimentally measured *χ* parameters is phenomenological^{16} since SCFT is mean-field yet the input “effective” *χ* parameter contains segment-scale correlation and fluctuation effects at the very least. The predictive power of phenomenological SCFT relies on the transferability of these measured *χ* parameters to other regions of the phase space, lest the experimentalist be required to measure *χ* everywhere. A great deal of effort has thus been made to understand the extent to which the effective *χ* parameter depends on effects such as segment- and *R*_{g}-scale fluctuations,^{17–21} solvent dilution,^{22,23} salt-doping,^{9,24–28} the incorporation of charged moieties in the chains,^{29,30} and other variables such as chain architecture, composition, and molecular weight.^{31–35}

Even in the simple context of a neat polymeric system with two monomeric species and no free or bound charges, the *χ* parameter contains both entropic and enthalpic contributions. The entropic contribution *χ*_{S} arises from differences in molecular structure and free volume between monomeric species, which can lead to segment packing effects that produce an excess entropy of mixing and must be absorbed into *χ*.^{32,36,37} The enthalpic contribution *χ*_{H}, on the other hand, is typically dominated by van der Waals (vdW) (or dispersion) forces between polarizable segments. For species *i* and *j* having segment polarizabilities *α*_{i} and *α*_{j}, these attractive interactions are characterized by energies that depend on the product of segment polarizabilities *α*_{i}*α*_{j} and lead to a *χ*_{H} which has a $(\alpha i\u2212\alpha j)2$ dependence.^{10}

Van der Waals interactions are rarely incorporated explicitly in either particle-based or field-based simulations of polymeric systems. Instead, their effects are commonly embedded implicitly via Lennard-Jones-type interaction potentials or by imposing a bare *χ* parameter, neither of which couple to the electrostatic environment. For problems involving the structure and thermodynamics of polymeric systems with charges or under applied electric fields, this may be undesirable since the very dipoles that produce a dielectric response to nearby ions or external electric fields are also the source of vdW attractions that contribute to *χ*_{H}. Indeed, the theory of vdW interactions, which accounts for their dependence on the dielectric and conducting properties of the medium through which they propagate, is well developed,^{38,39} yet few attempts seem to have been made to fold such theory into the existing frameworks for studying dense, inhomogeneous polymeric systems.

The neglect of explicit molecular polarizabilities in simulations is understandable. In the case of particle-based simulations, resolving the time scale associated with rapidly fluctuating induced dipoles, on top of already expensive-to-compute long-range Coulomb interactions, renders the problem highly resource-intensive. In the field-theoretic case, any fluctuation-induced electrostatic effect vanishes in the mean-field limit. This eliminates not just the vdW interactions but also charge correlation effects which produce coacervation and ion self-energies that serve to localize ions in high dielectric regions. This means that SCFT (the workhorse of the field-based simulation community) and other mean-field approaches (such as Poisson-Boltzmann theory) are by themselves often unsuitable for such problems.

Nonetheless, field-theoretic approaches are otherwise uniquely suited to the treatment of electrostatics since the particle-to-field transform renders the Coulomb interaction, which is long-ranged in the space of particle coordinates, short-ranged in the space of electrostatic field configurations.^{40} A variety of field-theoretic approaches have been proposed that either (a) supplement a mean-field framework with another method that is able to incorporate the effects of electrostatics,^{28–30,41–43} (b) relax the mean-field approximation more systematically, via Gaussian fluctuation approximations or loop expansions about the homogeneous state,^{24,44–50} or (c) relax the mean-field approximation entirely, which requires complex Langevin (CL) sampling of the fully fluctuating field theory.^{51–54}

Martin and co-workers^{55} recently introduced a polarizable field theory framework in which small molecules and polymers of arbitrary architecture can be constructed from elementary bead units that may be granted permanent or induced dipoles (and, optionally, monopole charges). In addition to the Coulomb interaction, the beads can also be subject to excluded volume interactions, Helfand compressibility constraints, and bare Flory-Huggins interactions. The field theory is rendered ultraviolet (UV) convergent by smearing the microscopic bead and charge densities into finite-ranged Gaussians.^{47} Martin *et al.* elucidated the structure of the field theory and demonstrated via a Gaussian fluctuation analysis that dielectric response, ion solvation energies, and vdW attractions emerge in a self-consistent manner as a result of the inclusion of bead polarizabilities. They showed that in a two-species system with polarizability contrast, the vdW attractions generate a *χ* parameter that has the expected $(\alpha i\u2212\alpha j)2$ dependence, although this calculation was restricted to the dilute limit. They also laid the groundwork for simulations of the unapproximated field theory via CL sampling, by presenting formulae for various observables and operators of interest, as well as “forces” that drive the sampling of field configurations in the complex Hamiltonian landscape.

In this paper, we explore the precise nature of the polarizability-generated vdW interactions and their contribution to the *χ* parameter, in the polymer context and within the polarizable field theory framework introduced by Martin and co-workers.^{55} Our analysis is restricted to non-polar but polarizable systems, in the absence of ions (ion effects and permanent dipoles will be considered in future work). Using one-loop perturbation theory, we have computed the structure factor and dielectric function for a binary homopolymer blend with polarizability contrast and identified corrections to the excluded volume parameter due to the vdW interactions, which define the vdW contribution to *χ* (here denoted $\chi \u0303$). Despite the limitations of perturbation theory, namely, its practical restriction here to the homogeneous state, the resulting expression for $\chi \u0303$ provides insight into the underlying charge physics and has a predictive capability that should prove useful when combined with direct simulations of the unapproximated field theory, allowing for the investigation of a wide range of inhomogeneous charge fluctuation phenomena. In that spirit and to test the predictions of the one-loop theory, we also present the first field-theoretic simulations (FTSs) of the unapproximated polarizable field theory, using complex Langevin sampling.^{56–59} We have conducted simulations of a single-species system of polarizable homopolymers and measured the structure factor and dielectric function for comparison with the one-loop theory. In addition, we measured the shift in the order-disorder transition (ODT) for a symmetric diblock copolymer melt with and without polarizability contrast, using an orientational persistence order parameter. This provides a direct measurement of $\chi \u0303N$, which we compare with the one-loop prediction.

This article is organized as follows: Sec. II outlines the microscopic ingredients for models of polarizable polymeric systems, as well as the field-theoretic representations of those models, which are used in the subsequent one-loop analysis and field-theoretic simulations. In Sec. III, the details of the one-loop calculation are presented, using a diagrammatic approach. The one-loop result clarifies the manner in which $\chi \u0303$ is sensitive to the electrostatic environment; implications for phase behavior are discussed. In Sec. IV, results of the field-theoretic simulations are presented and compared with the predictions of the one-loop theory, both for the polarizable homopolymer and the diblock copolymer melt. Section V concludes with a summary of the results and perspectives for future work. Extra details of the one-loop analysis can be found in Appendixes A–C.

## II. MICROSCOPIC MODEL AND FIELD THEORETIC REPRESENTATION

The general formalism, which allows for the elementary bead units to be granted monopole charges as well as permanent and induced dipoles, is described in detail in Ref. 55. Here we focus on the induced dipole case; thus in this article we are constructing models for non-polar but polarizable polymeric molecules. Each bead is granted polarizability via the attachment of a classical Drude oscillator to its center of mass. These Drude oscillators can be characterized by the magnitude of the partial charges *δq* and the spring constant *K* for the spring tethering them together; these parameters define the *bead polarizability* $\alpha \u030c=(\delta q)2/K$. The microscopic density for each bead, as well as the charge distribution for each partial charge, is smeared by a Gaussian $\Gamma (r)\u2009=\u2009(2\pi as2)\u22123/2\u2061exp\u2212r2/2as2$, where *a*_{s} characterizes the smearing range. The resulting smeared microscopic bead and charge densities are given by $\rho \xaf\u030c\u2009=\u2009\Gamma *\rho \u030c$ and $\rho \xaf\u030cc\u2009=\u2009\Gamma *\rho \u030cc$, respectively, where * indicates a spatial convolution. We note that the microscopic (as opposed to averaged field) variables are indicated by the check, and smeared quantities are indicated by the overbar. This smearing approach renders the eventual field theory UV convergent and has in recent years become a common regularization strategy.^{47,53,54,60,61}

Linear chains may be constructed by linking these beads together. Here we consider a system of volume $V$, containing *n* linear continuous Gaussian chains having equal degrees of polymerization *N* and statistical segment lengths *b* and composed of two species of beads that we shall denote by the labels *A* and *B*. For generality, at this point we do not specify the chain architecture or the number of polymeric species (we will consider specific architectures later). The bonded interaction energy per chain is given by

where *s* labels the position on the chain and $\beta =1kBT$. In addition to the bonded interaction, all pairs of beads interact via the familiar excluded-volume and Flory-Huggins interactions, which contribute interaction energies

where *u*_{0} and *χ*_{0} are the usual bare excluded volume and Flory-Huggins parameters, respectively, and *ρ*_{0} is the average bead density for the system. In addition, all partial charges in the system interact via the Coulomb interaction, contributing electrostatic energy

where *l*_{B} is a Bjerrum-like length and is given by $lB=\beta e24\pi \u03f50$, *e* and *ϵ*_{0} being the elementary charge and the vacuum permittivity, respectively.

### A. Auxiliary field representation

The microscopic model described above can be converted to a statistical field theory using standard techniques.^{40} The field-theoretic representation has the attractive feature that the chains are decoupled, as the non-bonded interactions (which couple the chains in the particle-based description) become mediated by a set of fluctuating fields. In the *auxiliary field representation*, the interaction matrix is diagonalized^{59} and the auxiliary fields are introduced via Hubbard-Stratonovich transformations, leading to the field-theoretic canonical partition function

where $w+$ and $w\u2212$ are auxiliary fields conjugate to the total density and the difference in the densities of the two species, respectively, and *ϕ* is a fluctuating electrostatic potential. The prefactor $Z0$ contains the normalization factors from the Hubbard-Stratonovich transforms, ideal gas contributions, and excluded-volume self-interaction corrections; it leads to a shift in the reference free energy that is inconsequential in this work and as such its form will be omitted here. The Hamiltonian *H*[$w+,w\u2212$, *ϕ*] has the form

where we have introduced a set of non-dimensionalized parameters for convenience. Here $C=nRg3/V$ is a dimensionless polymer chain number density, which corresponds to the average number of chains pervading the volume occupied by a given chain.^{40}^{,}$B=\beta u0N2Rg\u22123$ is a dimensionless excluded-volume parameter, describing the excluded volume interaction between statistical segments. $E=4\pi lBN2Rg\u22121$ is a dimensionless parameter setting the strength of electrostatic interactions in a vacuum, which is proportional to the dimensionless vacuum Bjerrum length *l*_{B}/*R*_{g}. Note that the actual strength of electrostatic interactions will depend on the explicit local dielectric environment as determined by the density, polarizability, etc., of material. The label *l* gives the polymer species (not to be confused with the bead species), and *ψ*_{l} is the volume fraction of polymer species *l*. Note that Eq. (6) can be straightforwardly generalized to the case where different polymeric species have different degrees of polymerization;^{59} however, in this work we are concerned only with the monodisperse case. We have also non-dimensionalized all length scales as well as the system volume $V$ using the unperturbed radius of gyration *R*_{g} as a reference length (defining the dimensionless $V=V/Rg3$) and absorbed a factor of *N* into the auxiliary fields. The single-chain partition function $Ql[w\xaf+,w\xaf\u2212,\varphi \xaf]$ accounts for the conformational statistics for a polymer chain of species *l*, subject to the smeared auxiliary fields $w\xaf+$, $w\xaf\u2212$, and $\varphi \xaf$. It is important to note that in the unapproximated theory, where each bead has embedded oscillating partial charges, $Ql[w\xaf+,w\xaf\u2212,\varphi \xaf]$ accounts for the degrees of freedom associated with the partial charges in addition to the usual chain conformation degrees of freedom. Here we use the dipole approximation to the oscillator charge distribution, which was also made in Ref. 55, to simplify *Q*_{l} to the form involving the familiar chain propagator *q*_{l}(**r**, *τ*),

where $ql(r,\tau ,[w\xaf+,w\xaf\u2212,\varphi \xaf])$ satisfies the modified diffusion equation for a linear polymer chain

with initial condition *q*_{l}(**r**, 0) = 1. Here we have introduced the scaled contour variable *τ* = *s*/*N*. A similar set of equations exists for the complementary propagator $ql\u2020(r,1\u2212\tau )$. Here Ω_{l}(**r**, *τ*) is the local chemical potential field experienced by a bead at contour position *τ* on a polymer of type *l*. For example, in the case of a binary *AB* homopolymer blend, Ω_{l}(**r**, *τ*) depends on the auxiliary fields according to

where *α*_{A} and *α*_{B} are scaled, dimensionless bead polarizabilities satisfying $\alpha A/B=\alpha \u030cA/B/(4\pi \u03f50lBNRg2)$. For an *AB* diblock copolymer with *A* block fraction *f*_{A}, Ω(**r**, *τ*) has the form

In the special case where *χ*_{0} = 0, the auxiliary field representation must be derived without the $w\u2212$ field, due to the breakdown of the Hubbard-Stratonovich transformation in the *χ*_{0} → 0 limit. In this case, the canonical partition function takes the form

where

The local chemical potential field for a binary homopolymer blend is given by

and for a diblock copolymer melt

In principle, thermodynamic observables are computed according to the appropriate ensemble average over field configurations. The ensemble average of an operator $O$ is given by

where we have returned to the description with non-zero *χ*_{0}. Due to the complex-valued nature of the Hamiltonian, the statistical weight that enters Eq. (15) has a phase that oscillates rapidly with the (real-valued) field configurations. Known as the sign problem, this renders traditional sampling methods, such as Monte Carlo, ineffective for sampling the unapproximated field theory. Here, we will employ the complex Langevin field-theoretic simulation (CL-FTS) approach, which solves the sign problem by promoting the fields to be complex variables. The CL equations of motion, for the models discussed above, have the form

where the set of coefficients ${\lambda j}$ are real-valued mobility coefficients (*j* runs over $w+,w\u2212$, and *ϕ*) and the ${\eta j}$ are real, Gaussian-distributed white noise variables that satisfy ⟨*η*_{j}(**r**, *t*)⟩ = 0 and the fluctuation-dissipation theorem ⟨*η*_{j}(**r**, *t*)*η*_{k}(**r**′, *t*′)⟩ = 2*λδ*(**r** − **r**′)*δ*(*t* − *t*′)*δ*_{jk}. Comprehensive discussions of the technical details of the CL approach, including the CL forces for models that incorporate electrostatics and polarizability, can be found elsewhere.^{40,54,55,59,61} We will thus omit these details for brevity.

## III. ONE-LOOP PERTURBATION THEORY

The exact field-theoretic partition functions presented in Sec. II contain contributions to the Hamiltonian from all orders of the fluctuating fields, due to the complicated functional dependence of the single-chain partition functions *Q*_{l}[Ω_{l}] on Ω_{l}. In order to render the calculation of thermodynamic observables analytically tractable, an approximation must be made. Here we expand *Q*_{l} in powers of the field fluctuations about the homogeneous saddle-point, truncating at fourth order. The standard techniques of diagrammatic perturbation theory^{62,63} provide a prescription for how to treat the third- and fourth-order terms in the expansion. To see how the incorporation of bead polarizabilities affects structure and dielectric response, we will compute the structure factor and dielectric function in the one-loop approximation.

For simplicity, we will focus initially on the binary homopolymer blend, for the case *χ*_{0}*N* = 0. In this case, the auxiliary field representation of the field theory has the Hamiltonian

with species fields $\Omega A(r)=iw\xaf(r)+\alpha A2|\u2207\varphi \xaf|2$ and $\Omega B(r)=iw\xaf(r)+\alpha B2|\u2207\varphi \xaf|2$. The field $w$(**r**) may be written as $w$(**r**) = *ω*(**r**) + $w*$, where $w*$ = −*iBC* is the homogeneous saddle-point for the $w$ field. The insensitivity of the Hamiltonian to a spatially constant *ϕ* means that we can take *ϕ*^{*} = 0. The reference free energy shift due to $w*$ may be absorbed into $Z0$ and the resulting Hamiltonian written as a sum over Fourier modes

where we are using the Fourier conventions $f^k=\u222bdr\u2009f(r)e\u2212ik\u22c5r$ and $f(r)=1V\u2211kf^keik\u22c5r$. The Fourier-transformed species fields are given by $\Omega ^k=i\omega \xaf^k\u2212\alpha 2V\u2211k\u2032k\u2032\u22c5(k\u2212k\u2032)\varphi \xaf^k\u2032\varphi \xaf^k\u2212k\u2032$. Note that although the auxiliary fields $\omega ^k$ and $\varphi ^k$ vanish at **k** = 0, the species fields do *not* vanish at **k** = 0, but rather satisfy $\Omega ^0=\alpha 2V\u2211k\u2032k\u20322\varphi \xaf^k\u2032\varphi \xaf^\u2212k\u2032$. The expansion of the homopolymer single-chain partition function *Q*[Ω] to fourth-order in the field Ω is given by^{64}

where $\Lambda ^(m)(k1,\u2026,km)$ is the *m*-point segment density correlation function of an ideal (non-interacting) homopolymer,^{64}

where *i* and *j* sum from 1 to *m* and *χ*(*τ*_{i}, *τ*_{j}) is a function of *τ*_{i}, *τ*_{j} defined in Ref. 64, which is not to be confused with the Flory-Huggins interaction parameter.

A product of $\Omega ^$ fields, such as those in Eq. (21), will introduce mode-coupling terms in powers of the $\omega ^$ and $\varphi ^$ fields. Truncating these terms at fourth-order, the Hamiltonian is obtained by combining (20) and (21),

where *H*_{0} is the Gaussian Hamiltonian and *U* contains the higher-order terms. The field vertex functions $\Lambda ^(\omega 3)$, $\Lambda ^(\omega \varphi 2)$, $\Lambda ^(\omega 4)$, $\Lambda ^(\omega 2\varphi 2)$, and $\Lambda ^(\varphi 4)$, which depend on the homopolymer vertex functions, are defined in Appendix A.

Using the Hamiltonian of Eq. (23), we seek to generate corrections to the bare (Gaussian) propagators $\u27e8\varphi ^k\varphi ^\u2212k\u27e90$ and $\u27e8\omega ^k\omega ^\u2212k\u27e90$, which appear in the Gaussian Hamiltonian *H*_{0}[*ω*, *ϕ*],

Here we have introduced the average polarizability *α*_{avg} = *ψ*_{A}*α*_{A} + *ψ*_{B}*α*_{B}, the Debye function $\u011dD(k)=2k\u22124(e\u2212k2\u22121+k2)$, and the Fourier-transformed smearing function which is given by $\Gamma ^k=e\u2212k2a2/2$, where *a* = *a*_{s}/*R*_{g}. The Gaussian propagators were also given by Martin and co-workers.^{55} Equation (24) indicates screening of the bare Coulomb interaction, due to the dielectric response of the beads, and defines the dielectric function $\u03f5^0(k)$ in the Gaussian approximation. We note that over short length scales (**k** → *∞*) the vacuum Coulomb interaction is recovered as $limk\u2192\u221e\u03f5^0(k)=1$, and over long length scales (**k** → 0), $\u03f5^0(k)$ approaches the bulk dielectric constant $\u03f5^0(0)=1+EC\alpha avg$. Equation (25) is the familiar screened excluded volume potential^{65,66} for the regularized Edwards solution model.

### A. Diagrammatic summations

Corrections to the bare propagators due to the higher-order terms in *U*[*ω*, *ϕ*] may be calculated using the diagrammatic technique, via the partial summation of connected diagrams. Figures 1 and 2 present the irreducible one-loop diagrams that contribute to $\u27e8\varphi ^k\varphi ^\u2212k\u27e9$ and $\u27e8\omega ^k\omega ^\u2212k\u27e9$, respectively. We initially compute the sum of the connected diagrams that produce corrections to $\u27e8\varphi ^k\varphi ^\u2212k\u27e90$ and $\u27e8\omega ^k\omega ^\u2212k\u27e90$, which may be constructed out of the irreducible diagrams given in Figs. 1 and 2, respectively. In the thermodynamic limit, the result in both cases is a geometric series with the sum given by the Dyson equation,^{67,68}

where Σ_{ϕ} and Σ_{ω} are the self-energies, given in Appendix B. The contributions to Σ_{ϕ} are from Figs. 1(a) and 1(c), and the contributions to Σ_{ω} are from Figs. 2(b)–2(d). Figures 1(b) and 2(a), which are derived from the *ω*^{2}*ϕ*^{2} vertex, evaluate identically to zero (this is demonstrated in Appendix B).

The dielectric function $\u03f5^(k)$ in the one-loop approximation may be computed from Eq. (26), inserting the expression for Σ_{ϕ}, and defining $\u27e8\varphi ^k\varphi ^\u2212k\u27e9\u22121\u2261k2\u03f5^(k)EV$,

where we have introduced the average squared polarizability $(\alpha 2)avg=\psi A\alpha A2+\psi B\alpha B2$. The second and third terms in Eq. (28) correspond to corrections to the dielectric function due to inter- and intra-chain correlations, respectively. We note that the inter-chain correlation term introduces a weak dependence of $\u03f5^(k)$ on the system compressibility, which we explore further in the supplementary material. The incompressible limit of Eq. (28) may be obtained by taking *B* → *∞*; the result is

where we have used the identity $(\alpha 2)avg=\alpha avg2+\psi A\psi B\alpha A\u2212\alpha B2$. Equation (29) indicates that in the absence of polarizability contrast, the one-loop correction to the dielectric function vanishes in the incompressible limit.

The one-loop corrections to the dielectric function also have an *N*-dependence that vanishes in the large-*N* limit at constant segment density $\rho 0=(nA+nB)N/V$. This is consistent with our expectation that the dielectric response primarily depends upon monomer-scale physics. As *N* becomes large, the loop integrals in Eqs. (28) and (29) become dominated by wave-vectors that are large compared with $Rg\u22121$, permitting the approximation $\u011dD(k)\u22482k2$ in their integrands. Noting that *C* ∼ *N*^{1/2}, *E* ∼ *N*^{3/2}, *B* ∼ *N*^{1/2}, *α* ∼ *N*^{−2}, and *a* ∼ *N*^{−1/2}, it is straightforward to see that in this limit, Eqs. (28) and (29) become *N*-independent. We point the reader to the supplementary material for further discussion of the *N*-dependence of $\u03f5^(k)$.

The particle-center structure factor *S*(**k**) (corresponding to the structure factor of unsmeared densities) can be related to the propagator $\u27e8\omega ^k\omega ^\u2212k\u27e9$ by the following relationship:

which, invoking Eq. (27) and the expression for Σ_{ω} given in Appendix B, takes the form

*Ĵ*(**k**) is a structure function containing the one-loop excluded volume and electrostatic fluctuation corrections and can be written as *Ĵ*(**k**) = *ĝ*_{D}(**k**) + *Ĵ*_{1}(**k**) + *Ĵ*_{2}(**k**) + *Ĵ*_{3}(**k**), where

The terms *Ĵ*_{1}(**k**) and *Ĵ*_{2}(**k**) are regularized versions of well-known excluded volume fluctuation corrections in the Edwards solution model, due to intra-chain and inter-chain correlations, respectively.^{64,65} The term *Ĵ*_{3}(**k**) is an inter-chain correlation term due to polarizability effects and is of particular interest here. The diagrammatic representation of the corrections to $\u27e8\omega ^k\omega ^\u2212k\u27e9$ which generate *Ĵ*_{3}(**k**) [those due to (c) from Fig. 2] provides some useful insight into the underlying physics.

Figure 3 shows two diagrammatic representations of one of the electrostatic corrections to $\u27e8\omega ^k\omega ^\u2212k\u27e9$ appearing in the geometric series given by Eq. (27). (a) uses the diagrammatic representation introduced earlier, whereas (b) shows an alternate representation, in which the homopolymer segment density correlation functions are represented explicitly. Here the thick solid vertical lines represent individual polymers, and it should be understood that the contour positions for insertions of the $\omega ^$ and $\varphi ^$ fields are to be integrated over, as in Eq. (22). Momentum flows are also indicated in Fig. 3(b). Similar diagrammatic representations can be found in the literature.^{67,69,70} The electrostatic interaction between polymers in Fig. 3, which involves a loop formed by two $\u27e8\varphi ^\varphi ^\u27e90$ propagators, is evidently the vdW interaction between polarizable beads on different polymers, as indicated by Fig. 3(b). This is in contrast to the screened excluded volume interaction between beads on different polymers, which involves a single $\u27e8\omega ^\omega ^\u27e90$ propagator and does not form a loop.

The series of corrections of the type of Fig. 3(a) is included in Eq. (27) and produces the term *Ĵ*_{3}(**k**). These diagrams involve a “chain” of polymers interacting via *alternating* vdW and screened excluded volume interactions. Any diagrams involving a chain of polymers interacting via *consecutive bare excluded volume interactions* are also included since they emerge at the Gaussian level and give rise to the screened excluded volume interaction. However, any diagrams involving a chain of polymers interacting via *consecutive vdW interactions* have not been accounted for. Such a diagram involves the *ϕ*^{4} vertex, and an example is shown in Fig. 4, in both diagrammatic representations. The alternate representation shows that this term contains two contributions, corresponding to Figs. 4(b) and 4(c). The former and its higher order analogs are precisely the diagrams that should be accounted for, involving consecutive vdW interactions. We will ignore diagrams of the form of Fig. 4(c) since they involve many-body induced dipole interactions rather than the simple two-body interactions.

The required series of diagrams may be included via a renormalized *ωϕ*^{2} vertex. Figure 5 shows the diagrammatic representation of this vertex renormalization, the calculation of which is described in Appendix C. The result may be incorporated into $\u27e8\omega ^k\omega ^\u2212k\u27e9$ by replacing *one* of the bare *ωϕ*^{2} vertices in the self-energy part of Fig. 2(c) with the renormalized vertex given by Fig. 5. The result modifies the correction *Ĵ*_{3}(**k**),

where for convenience we have introduced $I(k)$, which is an electrostatic fluctuation loop integral,

In order to understand the nature of the electrostatic fluctuation corrections to *S*(**k**) with *Ĵ*_{3}(**k**) given by Eq. (35), let us momentarily consider *S*(**k**) in the absence of the one-loop excluded volume corrections. In this case, we may write *Ĵ*(**k**) = *ĝ*_{D}(**k**) + *Ĵ*_{3}(**k**) in Eq. (31). After rearrangement and again making use of the identity $(\alpha 2)avg\u2009=\u2009\alpha avg2\u2009+\u2009\psi A\psi B\alpha A\u2009\u2212\u2009\alpha B2$, we obtain for *S*(**k**)

The electrostatic corrections in Eq. (37) can be understood entirely as a set of vdW-corrected excluded volume parameters *B*_{AA}(**k**), *B*_{AB}(**k**), and *B*_{BB}(**k**), given by

To see this, one can compare Eq. (37) to *S*(**k**) within the RPA for an unsmeared binary blend of homopolymers with excluded volume parameters *B*_{AA}, *B*_{AB}, and *B*_{BB} and no polarizability effects. The RPA *S*(**k**) is

### B. The effective $\chi $ parameter

As indicated by Eqs. (38)–(40), the van der Waals attraction between two polarizable beads reduces the effective excluded volume parameter between those two beads. In the two-species case, where *α*_{A} ≠ *α*_{B}, the set of *B*_{AA}(**k**), *B*_{AB}(**k**), and *B*_{BB}(**k**) define a wave-vector-dependent effective Flory-Huggins parameter $\chi \u0303(k)$, according to

where $v=\rho 0\u22121$ is the volume of a statistical segment. Equation (42) may be written in terms of the dimensionless parameters as

where $C=(nA+nB)Rg3V\u22121$ is the total dimensionless chain concentration. Equations (38)–(40) and (43) reveal the sensitivity of the vdW attractions to dielectric screening, contained in the dependence of the electrostatic loop integral $I(k)$ on $\u03f5^0(k)$. This effect was not captured by the dilute-limit estimate of $\chi \u0303$ given in Ref. 55. Indeed, for small *C* and at **k** = 0 and in the continuum limit ($1V\u2211k\u2032\u2192\u222bdk\u2032(2\pi )3$), Eq. (43) reduces to the result of Ref. 55,

which can also be obtained simply by writing $\u03f5^0=1$ in $I(0)$.

At constant segment density $\rho 0=(nA+nB)N/V$, the vdW-generated $\chi \u0303$ is independent of the degree of polymerization *N*, which can be seen most easily in Eq. (44) by accounting for the *N*-dependence of the dimensionless parameters *C* ∼ *N*^{1/2} and *E* ∼ *N*^{3/2}, the dimensionless polarizabilities *α* ∼ *N*^{−2}, and the dimensionless smearing range *a* ∼ *N*^{−1/2}. This is preserved in the full expression for $\chi \u0303$ given by Eq. (43), where we note that the *a*^{−3} dependence is contained in the loop integral $I(k)$. This molecular weight independence of $\chi \u0303$ is consistent with our expectations since contributions to *χ* should involve local structure and correlations beneath the *R*_{g} scale.

Figures 6(a)–6(c) show the dependence of the continuum limit of $(\chi \u0303N)(k)$ on the parameters *k*, *C*, and *E*, respectively, according to Eq. (43). Typical values^{71} for the statistical segment length (*b* ≈ 0.5 nm), monomeric volume ($v$ ≈ 0.1 nm^{3}), temperature (*T* ≈ 300 K), and degree of polymerization (*N* ≈ 100–1000) can be used to estimate reasonable values of *C* ≈ 1–5 and *E* ≈ 10^{6}–10^{8}, corresponding roughly to the regime of experimental interest for a polymer melt. Note that large values of *N* result in large values of *E* as a result of the *N*^{3/2} dependence of *E*, as well as the dielectric background being that of vacuum with corresponding Bjerrum length *l*_{B} ≈ 56 nm at room temperature. In all three figures, values of the smearing range *a* = *a*_{s}/*R*_{g} = 0.1, 0.15, and 0.2 are considered. The magnitude of $(\chi \u0303N)(k)$ depends strongly on *a*, which is unsurprising since vdW attractions should be strongly sensitive to the short-wavelength fluctuations in the *ϕ* field that manifest from fluctuating induced dipoles. The decrease in $(\chi \u0303N)(k)$ as *a* is increased naturally reflects the weakening of vdW attractions when these short-wavelength fluctuations are suppressed. The sources of the wave-vector-dependence of $(\chi \u0303N)(k)$ are twofold: from the fluctuation spectrum of the *ϕ* field itself (i.e., the propagator $\u27e8\varphi ^k\varphi ^\u2212k\u27e9$) and the smearing function $\Gamma ^k$. Figure 6(a) shows that the cumulative effect of these wave-vector dependences is a decrease in $(\chi \u0303N)(k)$ as *k* increases and that $(\chi \u0303N)(k)\u21920$ as *k* → *∞*. We are primarily interested in $(\chi \u0303N)(k)$ to the extent that it controls phase behaviour; that is, we are concerned with its value at a critical wave-vector *k*^{*} for phase separation in a given system. For most polymeric systems of interest, *k*^{*} will not be significantly larger than unity (for example, *k*^{*} = 0 for a binary homopolymer blend and *k*^{*} ≈ 2 for a symmetric diblock copolymer). Given that the *k*-dependence of $(\chi \u0303N)(k)$ at small *k* is relatively weak, it is sufficient for our purposes to focus on $(\chi \u0303N)(k=0)$. For brevity, we will henceforth refer to $(\chi \u0303N)(0)$ simply as $\chi \u0303N$.

Figure 6(b) shows the *C* dependence of $\chi \u0303N$, with all other parameters fixed. Note that $\chi \u0303N\u21920$ both for *C* → 0 (the dilute limit) and *C* → *∞* (the mean-field limit). In the dilute limit, chains are isolated and thus interact only with themselves, precluding the possibility of a polarizability contribution to *χN* (although the intra-chain vdW attractions could still drive, for instance, a coil to globule transition^{72}). As *C* → *∞*, on the other hand, fluctuations are completely suppressed and the vdW attractions themselves vanish. Figure 6(c) indicates that $\chi \u0303N$ increases rapidly for small values of *E*, but the growth slows at larger *E* due to the increasing effect of dielectric screening on the vdW attractions.

The system composition also has an important effect on $\chi \u0303N$. Figure 7(a) shows that $\chi \u0303N$ increases as the volume fraction of the more polarizable component (in this case, species *A*) decreases. The source of this effect is the dielectric function in the electrostatic fluctuation loop integral: $\u03f5^0(k)$ *decreases* as the system becomes less rich in the more polarizable component, thus reducing the screening effect on the van der Waals attractions and producing a larger $\chi \u0303N$. This effect will produce an asymmetry in the phase diagram which can becomedramatic in cases where the polarizability effects produce a large contribution to the total *χ* parameter. To see this, consider a binary homopolymer blend, with *N*_{A} = *N*_{B} and a non-zero bare *χ*_{0} parameter, in which one species is more polarizable than the other. The phase behavior with polarizability effects can be understood approximately from the Flory-Huggins spinodal for the incompressible binary blend,

where the total *χ* parameter has been enhanced by polarizability effects due to Eq. (43). Figure 7(b) shows the spinodal, given by Eq. (45), in terms of *χ*_{0}*N* for the case where species *A* is polarizable (*α*_{A} = 2 × 10^{−7}) and species *B* is not (*α*_{B} = 0). For *a* = 0.2, the polarizability effects are sufficiently weak that the spinodal does not take on a significant asymmetry. However, for both *a* = 0.15 and *a* = 0.1, the spinodal is altered dramatically, exhibiting windows of *ψ*_{A} for which a bare *χ*_{0} parameter is not required to cross into the unstable region. The binodal will exhibit a similar behavior.

A series of recent papers examining the effect of dielectric mismatch on the phase behavior of charge-containing systems have predicted asymmetries in the spinodal that qualitatively resemble the spinodals of Fig. 7. Nakamura^{49} studied spinodal decomposition in a mixture of a polymer and an ionic liquid, via a Gaussian fluctuation analysis of the free energy in a field-theoretic model. Rumyantsev and Kramarenko^{50} also computed spinodals, in solutions of polyelectrolytes and of ionomers, using RPA. In both cases, the dielectric contrast led to an asymmetry in the spinodal which indicates enhanced phase separation for systems rich in the component with lower dielectric constant. Sing and co-workers^{29,30} have also predicted similar trends in spinodals for blends of charged and neutral polymers (with counterions), using a hybrid self-consistent field/liquid state theory (SCFT-LS) approach. Although in their case the components had the same dielectric constant, the presence of Debye-Hückel screening for the polyelectrolyte and counterions (absent for the neutral polymer) evidently produces a similar effect as dielectric mismatch. In all of these studies, charge correlations between ions are the driving force for phase separation, whereas in our case the vdW attractions are the driving force. Nevertheless, the common feature of spinodal asymmetry is produced by a contrast in the electrostatic screening properties of the components of the system.

### C. Self-consistent one-loop theory

An obvious extension of the one-loop theory presented above consists of computing *Ĵ*(**k**) and $\u03f5^(k)$ self-consistently. Such a procedure is similar in spirit to the so-called “Hartree” approximation which has been applied in various forms to canonical field theory models.^{17,19–21,67,73} In our case, this procedure consists of a renormalization of the self-energies Σ_{ω} and Σ_{ϕ} by replacing the bare propagators $\u27e8\omega ^k\omega ^\u2212k\u27e90$ and $\u27e8\varphi ^k\varphi ^\u2212k\u27e90$ inside the loop integrals with $\u27e8\omega ^k\omega ^\u2212k\u27e9$ and $\u27e8\varphi ^k\varphi ^\u2212k\u27e9$. A set of coupled equations for *Ĵ*(**k**) and $\u03f5^(k)$ is generated, given by

and

where $\u27e8\varphi ^k\varphi ^\u2212k\u27e9=EVk2\u03f5^(k)$, $\u27e8\omega ^k\omega ^\u2212k\u27e9=BV1+BC\u0134(k)\Gamma ^k2$, and $I(k)$ contains the corrected $\u03f5^(k)$ in place of $\u03f5^0(k)$. Here we solve Eqs. (46) and (47) self-consistently in the continuum limit using an iterative numerical procedure. An initial guess is made for *Ĵ*(**k**) and $\u03f5^(k)$: we choose the obvious initial guesses *Ĵ*(**k**) = *ĝ*_{D}(**k**) and $\u03f5^(k)=\u03f5^0(k)$. Equations (46) and (47) are then evaluated, and their results are compared to the guesses. The difference between the guesses and the results are then used to inform the guesses for the next iteration. This process is then continued until a fixed point is found, to within a small tolerance ($guess\u2212result<10\u22126$). In this self-consistent variant of the one-loop theory, the expression for $\chi \u0303N=(\chi \u0303N)(0)$ is obtained by replacing the Gaussian dielectric function $\u03f5^0$ in $I(0)$ with the self-consistently determined $\u03f5^$. In the continuum limit, $\chi \u0303N$ is given by

As discussed in Sec. III B, the bare one-loop expression for $\chi \u0303$ is independent of *N*. We also note that there is no explicit dependence of $\chi \u0303$ on the compressibility of the system [which would require *B* appearing in Eq. (43) or Eq. (48)]. Thus, in the self-consistent one-loop theory, $\chi \u0303$ will only depend on *B* or *N* if the one-loop corrections to the dielectric function $\u03f5^(k)$ depend on *B* or *N*. As discussed in Sec. III A and the supplementary material, the one-loop $\u03f5^(k)$ increases with *B* but becomes *B*-independent for large *B* (incompressible limit) and decreases with *N* but becomes *N*-independent for large *N*. In both cases, the dependences are relatively weak, producing a change in $\u03f5^(k)$ of less than 10% over the entire range of *B* and *N* for the parameters that we have tested.

In Sec. IV, results of the bare and self-consistent one-loop theories are presented, alongside field-theoretic simulations of the exact field theory model using CL sampling. In the one-loop theory, the integrals over **k**^{′} are evaluated numerically over the discretized (*k*′, *u*) plane, where *k*′ is the magnitude of **k**′ and *u* = cos *θ*, with *θ* being the angle between **k** and **k**′. The three- and four-point homopolymer vertex functions are pre-calculated in the discretized (*k*, *k*′, *u*) space, according to Eq. (22). In the calculation of the vertex functions, the chain contour is discretized into steps with Δ*τ* = 0.05, where *τ* ∈ [0, 1]. For consistency, in the CL-FTS results in Sec. IV A, we use the same chain contour step size Δ*τ* = 0.05 for the modified diffusion equation solver.

## IV. COMPARISON OF THE ONE-LOOP THEORY WITH CL-FTS

In this section, we present a comparison of the one-loop predictions for $\u03f5^(k)$, *S*(**k**), and $\chi \u0303N$ from Sec. III to the results from CL-FTS. Initially, we will present measurements of the dielectric function and structure factor for a system containing a single species of polarizable homopolymers (chosen for simplicity), alongside the one-loop predictions for the same model. To obtain a quantitative measure of $\chi \u0303N$, we have also characterized the order-disorder transition (ODT) in field-theoretic simulations of a symmetric polarizable diblock copolymer melt. Such a system is particularly convenient for the quantitative measurement of the ODT. Although the one-loop theory in Sec. III has been developed in the context of a binary system of homopolymers, the one-loop expression for $\chi \u0303N$ indicates that it is dominated by sub-*R*_{g} scale physics that is not explicitly sensitive to chain architecture. It therefore seems that any chain architecture dependence of $\chi \u0303N$ would enter through the one-loop corrections to the dielectric function $\u03f5^(k)$. We expect such a dependence to be a secondary effect at the worst, and therefore Eq. (48) should be applicable for our purposes to the polarizable diblock copolymer.

The properties of interest can be measured without approximation using CL-FTS; the only sources of error in these measurements are the sampling error associated with computing operator averages over a finite simulation interval and errors due to time and space discretization. Here, the structure factor and dielectric function are obtained directly from their relationships to $\u27e8\u0175k\u0175\u2212k\u27e9$ and $\u27e8\varphi ^k\varphi ^\u2212k\u27e9$, which have already been described in Sec. III. Note that in this section we consider the structure factor for the smeared density correlations rather than the particle-center structure factor which was considered in Sec. III A. The relevant CL operators are simply given by $\u0175k\u0175\u2212k$ and $\varphi ^k\varphi ^\u2212k$. To characterize the diblock ODT, we use an orientational persistence order parameter *λ* (see Appendix D).^{74}

The CL-FTS simulations presented here were performed in a cubic cell with periodic boundary conditions. Throughout, a spatial collocation mesh resolution of Δ*x* = *a* = 0.2 was used (all lengths are in units of *R*_{g}). The resolution Δ*x* = *a* has been shown to be sufficient to resolve the fluctuation spectrum in smeared models.^{61,75} The modified diffusion equations for the chain propagators were solved using a pseudospectral operator-splitting approach,^{76,77} and the CL equations of motion were propagated using a time step Δ*t* = 0.05 and the exponential time difference predictor-corrector (ETDPEC) algorithm.^{78} A large mobility coefficient *λ*_{ϕ} = 100 is used to accelerate thermalization of the *ϕ* field due to the large values of *E* used in this work; mobility coefficients of unity are used for the $w$ fields. In the polarizable homopolymer results of Sec. IV A, a cell with sides of length *L* = 12.8 was used, and the chain propagators were solved with chain contour resolution Δ*τ* = 0.05, for consistency with the chain contour resolution of homopolymer vertex functions in the one-loop theory. In the polarizable diblock copolymer results of Sec. IV B, *L* = 6.4 and a chain contour resolution Δ*τ* = 0.01 were used. All simulations were performed on NVIDIA Tesla M2075, C2070, or K80 graphics processing units (GPUs).^{79} Results for each state point correspond to a simulation on a single GPU, with typical runtimes of 2-4 days for roughly 2 × 10^{6}–5 × 10^{6} time steps.

### A. Polarizable homopolymer: Structure factor and dielectric function

The one-loop expressions for the dielectric function and structure factor of the polarizable homopolymer solution model can be obtained from Eqs. (28) and (31), by setting either *ψ*_{A} = 1 and *ψ*_{B} = 0 or *α*_{A} = *α*_{B}. The self-consistent one-loop equations are straightforwardly related to these equations according to the description in Sec. III C. Figure 8(a) shows the structure factor *S*(**k**) from the one-loop (1L) theory, self-consistent one-loop (SC-1L) theory, and CL-FTS, for *C* = 1, *B* = 1, *α* = 1 × 10^{−4}, and for various values of *E* ranging from *E* = 5 × 10^{3} to 3 × 10^{4}. The Gaussian fluctuation structure factor is plotted as a dotted line in Fig. 8(a); we remind the reader that since the *w* and *ϕ* fields are uncoupled at the level of Gaussian fluctuations, this structure factor is unresponsive to polarizability effects (and is therefore independent of *E*). The Gaussian fluctuation structure factor provides a useful reference, showing the degree to which the density fluctuations in the other *S*(**k**) curves are enhanced by van der Waals attractions. This effect becomes more pronounced as *E* is increased, as expected.

The one-loop and self-consistent one-loop predictions for *S*(**k**) are in excellent agreement with CL-FTS for the smaller values of *E* in Fig. 8(a). However, a discrepancy emerges as *E* increases, wherein both variations of the one-loop theory overestimate the enhancement of density fluctuations seen in the simulations. Perhaps surprisingly, the bare one-loop theory outperforms its self-consistent counterpart in the comparison with CL. The expression for the one-loop polarizability correction to the excluded volume parameter [which can be generalized to a single-species system from Eqs. (38)–(40)] reveals that for these larger values of *E*, *B*_{eff}(0) is approaching zero. This indicates that the system is likely not far from its stability limit, beyond which spontaneous phase separation, into high- and low-density phases, would occur. Indeed, as *E* is increased further the self-consistent one-loop theory eventually fails to converge (for the parameters in Fig. 8, around *E* ≈ 5 × 10^{4}), indicating that the fluctuation theory is no longer permitting a homogeneous solution. It is clear from Fig. 8(a) that both one-loop theories are likely to reach their stability limits before the fully fluctuating system.

Evidently, fluctuations in the vicinity of the spinodal are engaging higher-order terms in the loop expansion that are not included in the self-consistent or bare one-loop theories. Such a term, which is notably absent in the one-loop theory, is an explicit correction to the *single-chain structure* due to the polarizability effects. The absence of such a correction can be attributed to the fact that in the diagrams (see Fig. 2) for one-loop corrections to $\u27e8\omega ^k\omega ^\u2212k\u27e9$ from Sec. III, Fig. 2(a) evaluates to zero. This particular diagram would seem to be the electrostatic analog of the one-loop excluded volume correction to the single-chain structure, which is given by Fig. 2(b). In fact, the lowest-order loop diagram that produces an explicit, non-zero polarizability correction to the single-chain structure is of two-loop order. It is possible that such a correction could be required for a more robust theory of density fluctuations near the spinodal for this system and also for low concentrations where corrections to the single-chain structure cannot be ignored. However, our interest in this work is neither in characterizing the spinodal nor studying the dilute or semi-dilute regimes: we are interested in the dielectric properties of the dense homogeneous phase and in the $\chi \u0303$ parameter that is generated by polarizability contrast.

Figure 8(b) shows the dielectric functions from the Gaussian fluctuation theory, one-loop theory, and CL-FTS, for the same systems as in Fig. 8(a). The Gaussian approximation $\u03f5^0(k)$ overestimates the dielectric function significantly for the larger values of *E* in Fig. 8. The bare one-loop theory in this case is virtually indistinguishable from CL-FTS, except at small *k* and there presumably only because the CL data become noisy. For larger values of *E*, the self-consistent one-loop theory begins to slightly underestimate $\u03f5^(k)$. Since the van der Waals attractions are subject to dielectric screening, the fact that the self-consistent one-loop theory underestimates the dielectric function is consistent with its simultaneous overestimation of the density fluctuations shown in Fig. 8(a).

If we instead consider a system at higher concentration and far from the spinodal, the accuracy of the bare and self-consistent one-loop theories becomes more in line with expectations. Figure 9(a) shows the structure factor for*C* = 5, *B* = 1, *α* = 10^{−4}, and *a* = 0.2, with *E* ranging from 5 × 10^{4} to 2 × 10^{5}. Here, the self-consistent one-loop theory clearly does a better job than the bare one-loop theory of capturing the density fluctuations in the fully fluctuating system. As *E* increases, the self-consistent one-loop theory overestimates the density fluctuations; this trend should increase as *E* is increased further. In Fig. 9(b), the dielectric functions are shown for the same parameters as in Fig. 9(a). Both the bare and self-consistent one-loop theories have excellent agreement with CL-FTS; it is difficult to tell which is superior in this case.

The results of this section indicate that the one-loop theory is best suited to high concentrations, far from the spinodal. This is particularly the case when comparing the structure factors from the bare and self-consistent one-loop theories to fully fluctuating simulations and can be seen in the contrast between Figs. 8(a) and 9(a). On the other hand, the dielectric function in the one-loop theory fares better across the range of concentrations that have been tested and provides a significantly better estimate of $\u03f5^(k)$ in all cases than does the Gaussian approximation. This is despite the fact that it is coupled to the density fluctuations which, particularly in the case of the *C* = 1 data, are becoming less well described by the one-loop theory as *E* is increased. Evidently, the enhancement of density fluctuations due to vdW interactions is strongly sensitive to changes in the dielectric function, whereas the dielectric function is only weakly sensitive to the density fluctuations. In Sec. IV B, we will move to a highly concentrated (*C* ≥ 5) system of diblock copolymers with polarizability contrast and compare the measured $\chi \u0303N$ from CL-FTS with the predictions of the self-consistent one-loop theory.

### B. Symmetric polarizable diblock copolymer melt: Determination of $\chi \u0303N$

Here we consider a dense, weakly compressible system containing symmetric *AB* diblock copolymers with polarizabilities *α*_{A} and *α*_{B}, supplemented with a bare *χ*_{0} parameter. We have measured the contribution to *χN* from polarizability contrast for this system, for various values of *C*, by computing in CL-FTS the orientational persistence order parameter *λ* (defined in Appendix D) as a function of *χ*_{0}*N* in the vicinity of the ODT. As the species polarizability contrast $\alpha A\u2212\alpha B$ is increased from zero, the ODT shifts to smaller values of *χ*_{0}*N* by an amount which we define to be $\chi \u0303N$. This manifests as a shift in the *λ* versus *χ*_{0}*N* curve, from which $\chi \u0303N$ can be measured. We fix *E* = 4 × 10^{6} and vary *B* with *C* such that the compressibility is fixed according to *BC* = 100. Maintaining constant *BC* in this manner ensures that the system behaves as a weakly compressible melt with the same compressibility for all concentrations.

In Fig. 10, the *λ* versus *χ*_{0}*N* data are plotted for simulations with *C* = 5 (black), *C* = 10 (blue), *C* = 15 (red), *C* = 20 (green), and *C* = 25 (brown). For each value of *C*, two sets of simulations have been conducted over a range of *χ*_{0}*N*: a set with *α*_{A} = *α*_{B} = 5 × 10^{−8} (zero polarizability contrast, dotted lines) and a set with *α*_{A} = 1 × 10^{−6} and *α*_{B} = 5 × 10^{−8} (solid lines). All simulations shown in Fig. 10 were initialized with a disordered configuration. We note that in our relatively small cells, no hysteresis is observed when comparing with simulations initialized with a lamellar configuration (3 lamellar periods are formed along the box diagonal when ordered). The relatively small system size lowers the free energy barrier for the system to hop between disordered and ordered states when near the ODT, enabling the relevant configuration space to be explored within reasonable simulation time. However, the precision with which the *λ* versus *χ*_{0}*N* curve is able to delineate the ODT is lessened as the system size is decreased (and also as the fluctuation strength is increased, by lowering *C*, an effect which can be seen in Fig. 10). Although such effects are normally a concern for locating the ODT, in our case we are interested in the *relative shift* for two *λ* versus *χ*_{0}*N* curves at the same value of *C*. In that case, the precise delineation of $(\chi 0N)ODT$ is not required as long as the *λ* versus *χ*_{0}*N* curves, for a given *C* but different polarizability contrast, are not altered qualitatively in any way other than by a horizontal shift. Figure 10 indicates that this is indeed approximately the case.

In order to compute the shift in *χ*_{0}*N* between two *λ* versus *χ*_{0}*N* curves, we choose one data set to be the reference and apply a uniform horizontal shift Δ*χ*_{0}*N* to the non-reference data set. We numerically locate the value of Δ*χ*_{0}*N* that minimizes the (appropriately normalized) area bounded vertically by the shifted curve and the reference curve and bounded horizontally (as well as normalized) by the range of *χ*_{0}*N* which the two curves share. We linearly interpolate between data points for the determination of the bounded area, as is illustrated in Fig. 11. Note that the area of the shaded region in Fig. 11 is the unnormalized bounded area. In order to determine the uncertainty in Δ*χ*_{0}*N*, we allow the values of *λ* for the data points in both sets to fluctuate, each according to a Gaussian distribution with a variance given by the standard error *δλ* for that data point. For each realization of the {*λ*} for the two data sets, the shift Δ*χ*_{0}*N* is computed and the accumulated statistics are then used to compute the uncertainty *δ*Δ*χ*_{0}*N*.

This procedure, carried out for the data sets in Fig. 10, yields the shift in the ODT $\Delta (\chi 0N)ODT$, as a function of concentration *C*, between symmetric diblock copolymer melts with and without polarizability contrast. The ODT shift gives $\chi \u0303N$ via $\chi \u0303N=\u2212\Delta (\chi 0N)ODT$. Figure 12 plots $\chi \u0303N$ versus concentration *C* from CL-FTS, along with the bare and self-consistent one-loop predictions. The SC-1L prediction for $\chi \u0303N$ is in excellent agreement with CL-FTS and is a significant improvement over the bare 1L theory, except at *C* = 5 where both one-loop theories are statistically indistinguishable in their comparison with CL-FTS. This is consistent with the observation in Sec. IV A that the SC-1L theory is best suited to high concentrations. The bare one-loop prediction (dashed line) systematically underestimates $\chi \u0303N$ due to the fact that the Gaussian approximation $\u03f5^0(k)$, which in Sec. IV A was seen to overestimate the dielectric function, is used in the electrostatic loop integral. Nevertheless, the bare one-loop theory should provide a useful estimate of $\chi \u0303N$ in cases where computing the self-consistent one-loop correction to $\u03f5^0(k)$ is deemed more effort than it is worth.

Figure 12 indicates that the one-loop prediction for $\chi \u0303N$ does an extremely good job of capturing the enhancement of *χN* due to the van der Waals attractions in a system with species polarizability contrast. Evidently, neither $\chi \u0303N$ nor $\u03f5^(k)$ are especially sensitive to chain architecture, given that the one-loop equations were formulated for a binary homopolymer blend in Sec. III yet were successfully applied here to a diblock copolymer melt. It was suggested in Sec. IV A that the self-consistent one-loop theory should be best suited to systems at high concentrations, perhaps in part due to the absence of electrostatic corrections to single-chain structure at the one-loop level; here, we see that SC-1L indeed performs well at melt concentrations. We emphasize also that we have performed simulations at values of *C* low enough to enter the regime of experimentally relevant molecular weights (*C* = 5 corresponds to an invariant degree of polymerization of $N\xaf=5.4\xd7103$).

## V. CONCLUSIONS

We have addressed the effect of the incorporation of bead polarizabilities on the density fluctuations, dielectric response, and phase behavior of neutral, non-polar polymeric systems, using the polarizable field theory framework introduced by Martin and co-workers.^{55} In this theory, polymers are constructed out of beads that have been made polarizable by the attachment of classical Drude oscillators to their centers of mass. The resulting molecules exhibit dielectric response and van der Waals attractions in a manner that is self-consistently determined by their polarizabilities. The vdW interactions generate a $\chi \u0303$ parameter between species with polarizability contrast, making the approach well suited to address phase separation in charged systems or in systems with applied electric fields since the vdW interactions are generated consistently within the electrostatic environment rather than being imposed phenomenologically.

The framework is amenable to standard analytical approximations (such as loop expansions) as well as direct simulation via complex Langevin sampling, both of which we have employed in this paper. Using a one-loop approximation, we computed corrections to the structure factor and dielectric function for a binary homopolymer blend with polarizability contrast. The electrostatic one-loop corrections to the structure factor can be entirely accounted for by a renormalization of the excluded volume parameter *B*, into three effective parameters *B*_{AA}(**k**), *B*_{AB}(**k**), and *B*_{BB}(**k**), which in turn define $\chi \u0303N$. The diagrammatic approach from Sec. III provides some insight into the nature of the vdW interaction within perturbation theory. In contrast to the excluded volume interaction, which is propagated between two beads by a single $\u27e8\omega ^\omega ^\u27e9$ line, the vdW interaction is propagated between two polarizable beads by two $\u27e8\varphi ^\varphi ^\u27e9$ lines that form a closed loop. The electrostatic one-loop correction to *S*(**k**) thus sums only those diagrams where there is at most a single vdW interaction between any pair of chains and *no* intra-chain vdW interactions (which are of at least two-loop order). The result of the correction is more closely analogous to the random phase approximation^{69,70,80} [as indicated by the comparison of Eqs. (37) and (41)] than it is, for instance, to the excluded volume one-loop corrections, which include both intra-chain interactions and polymer pairs interacting at two points simultaneously.^{65}

In order to investigate the validity of the one-loop theory, we also performed simulations of the fully fluctuating theory using CL-FTS. These constitute the first direct simulations conducted within the polarizable field theory framework. We measured the structure factor and dielectric function for a system containing a single species of polarizable homopolymers and compared them with self-consistent and bare versions of the one-loop theory. We found that the one-loop theory performs best at high concentrations, where the agreement with CL-FTS is excellent. We then measured the shift in the value of *χ*_{0}*N* at the ODT due to polarizability contrast in a compositionally symmetric polarizable diblock copolymer melt, which allowed us to extract the vdW-generated $\chi \u0303N$ from simulations. The measured $\chi \u0303N$ is in excellent agreement with the prediction of the self-consistent one-loop theory, particularly at high concentrations.

We suggest several straightforward extensions of the work presented here. Given that the polarizability-generated $\chi \u0303$ is coupled to the electrostatic environment, it will be interesting to investigate the effects on the phase behavior of polarizable polymeric systems due to the incorporation of charged components or application of an electric field. The addition of salts to block copolymers in which one block preferentially solvates the ions, for instance, is known to produce an additional driving force for phase separation^{7} that is often interpreted as an increase in *χ*_{eff}.^{8,9,27,28} This effect should be captured by the polarizable field theory; however, as we have demonstrated in this paper, the vdW-generated $\chi \u0303$ is also sensitive to the electrostatic screening properties of the medium, which are altered upon salt-doping. Another possible direction for future work is in applying the loop expansion technique used here to account for corrections to the single-chain structure due to vdW interactions, which would require a two-loop expansion. The vdW interactions should be able to drive a coil-to-globule transition, which might then be induced or suppressed by electrostatic stimulus.

A few comments are in order regarding the mapping of polymer chemistries onto the polarizable model. Since the dielectric constants for the components, and the $\chi \u0303$ parameter between them, are self-consistently determined by the polarizabilities rather than being independent inputs to the theory, the electrostatic parameters should be chosen judiciously. As discussed in Sec. III A and the supplementary material, neither $\u03f5^(k)$ nor $\chi \u0303$ depend strongly on the compressibility for a concentrated system. Thus, the mapping of the excluded volume parameter *B*, which here serves mainly to set the compressibility, is not as crucial as long as the system compressibility is appropriately weak for a concentrated solution or melt. For a given degree of polymerization *N*, reference volume *v*, and temperature *T*, the parameters *C* and *E* may be estimated using tabulated values of statistical segment lengths.^{71} The bulk dielectric constants *ϵ*_{A} and *ϵ*_{B} can then be tuned roughly to the desired values by adjusting *α*_{A} and *α*_{B}, guided by the Gaussian approximation $\u03f5^A/B(0)\u22481+CE\alpha A/B$. As discussed in Sec. III B, $\chi \u0303$ is strongly sensitive to the smearing range *a*, whereas the bulk dielectric constant $\u03f5^(0)$ depends on *a* only through the one-loop corrections, which are small (although not negligible) in comparison with $\u03f5^0(0)$. With *α*_{A} and *α*_{B} fixed, *a* can therefore be tuned as a model parameter to produce the desired $\chi \u0303$ to match experimental data.

As an example of the mapping heuristic, consider a symmetric poly(ethylene oxide)-*b*-polystyrene (PEO-*b*-PS) diblock copolymer. Measurements at *T* = 210–230 °C have indicated that $\chi PEO/PS=\u22121.73\xd710\u22122+23.7T$, for a reference volume $v$ = 0.1 nm^{3}.^{71,81} Assuming a degree of polymerization of roughly *N* ≈ 200 and measured statistical segment lengths *b*_{PEO} = 0.72 nm and *b*_{PS} = 0.50 nm,^{71} we can estimate *C* ≈ 1–2 and *E* ≈ 10^{6}–10^{7} at *T* ≈ 200 °C. Setting, for example, *C* = 1.5 and *E* = 10^{7}, dielectric constants close to *ϵ*_{PEO} = 10 and *ϵ*_{PS} = 2.5 can be achieved for *α*_{PEO} = 7.5 × 10^{−7} and *α*_{PS} = 1.0 × 10^{−7}. We must now adjust the smearing range *a* in order to produce a $\chi \u0303N$ that is consistent with the enthalpic contribution to *χ*_{PEO/PS} ($\chi HN=23.7NT\u224810.0$ for *N* = 200 and at *T* = 200 °C), assuming that *χ*_{H} is entirely due to vdW interactions. The one-loop theory indicates that for these parameters, $\chi \u0303N\u224810.0$ can be achieved if *a* ≈ 0.115. The remaining (entropic) contribution to *χ* can be accounted for with a bare *χ*_{0} parameter. Including the one-loop corrections for these parameters with *a* = 0.115 gives dielectric constants of *ϵ*_{PEO} ≈ 10.4 and *ϵ*_{PS} ≈ 2.4.

## SUPPLEMENTARY MATERIAL

See supplementary material for a discussion of the dependence of the one-loop dielectric function on the excluded volume parameter *B* and the degree of polymerization *N*.

## ACKNOWLEDGMENTS

This research was partially supported by the NSF DMR-CMMT Program under Award No. DMR-1506008, by the NSF MRSEC Program under Award No. DMR 1720256, and by the Institute for Collaborative Biotechnologies through Grant No. W911NF-09-001 from the U.S. Army Research Office. The content of the information does not necessarily reflect the position or the policy of the Government, and no official endorsement should be inferred. Use was made of the computational facilities of the Center for Scientific Computing at the CNSI and MRL: an NSF MRSEC (No. DMR-1720256) and No. NSF CNS-1725797.

### APPENDIX A: FIELD VERTEX FUNCTIONS

Here we present the formulae for the field vertex functions referred to in Eq. (23). Combining Eqs. (20) and (21) and the expressions for the Fourier-transformed species fields $\Omega ^A$ and $\Omega ^B$ leads to the following vertex functions:

and

These field vertex functions are related to the *m*-point homopolymer segment density correlation functions $\Lambda ^(m)(k1,\u2026,km)$, which are defined in Eq. (22). The vertex functions $\Lambda ^(\omega 3)$ and $\Lambda ^(\omega 4)$, given by Eqs. (A1) and (A3), are simply the smeared analogs of the vertex functions that are well known from the Edwards solution model.^{16,64}

### APPENDIX B: ONE-LOOP SELF-ENERGIES AND CANCELLATION OF *ω*^{2}$\varphi 2$ ONE-LOOP DIAGRAMS

The self-energies Σ_{ϕ}(**k**) and Σ_{ω}(**k**) that appear in Eqs. (26) and (27), respectively, are given by

and

The contribution to Σ_{ϕ} from Fig. 1(b) evaluates to zero, as does the contribution to Σ_{ω} from Fig. 2(a). Let us denote these contributions as $I1b(k)$ and $I2a(k)$, respectively. Both of these diagrams are generated by the *ω*^{2}*ϕ*^{2} vertex. We note the following property of the three-point homopolymer segment density correlation function $\Lambda ^(3)$, which is a consequence of Eq. (22):

Equation (B3) also holds for any permutation of the wave-vector arguments. The contribution of Fig. 2(a) to Σ_{ω} is given by

### APPENDIX C: $\omega \varphi 2$ VERTEX RENORMALIZATION

Figure 5 illustrates diagrammatically how to sum the desired series of corrections to $\Lambda ^(\omega \varphi 2)$ due to $\Lambda ^(\varphi 4)$. Denoting the renormalized *ωϕ*^{2} vertex as $\Lambda ^r(\omega \varphi 2)$, we may write

where we have chosen to absorb the renormalization into $\Lambda ^r(2)$. The second equation in Fig. 5 can be written, after simplification, in terms of $\Lambda ^r(2)$ as

which is easily solved for $\Lambda ^r(2)$. Note that $I(k)$ is defined by Eq. (36), and we have taken only the corrections of the form of Fig. 4(b) [ignoring Fig. 4(c) and its higher-order analogs]. The contribution to the self-energy Σ_{ω} from Fig. 2(c), which we denote $I2c(k)$, must be modified to include the correction by replacing one of the bare *ωϕ*^{2} vertices in Fig. 2(c) with the renormalized $\Lambda ^r(\omega \varphi 2)$. The result is

which should replace the third term in Eq. (B2).

### APPENDIX D: ORIENTATIONAL PERSISTENCE ORDER PARAMETER *λ*

We use an orientational persistence order parameter *λ* to characterize the presence of long-range order in the diblock copolymer system of Sec. IV B. A local director field operator *ñ*_{i}(**r**; [{$w$}]) is defined, which describes the orientation of domains,

where it is assumed to be understood that the functional dependence of *ñ*_{i} on {$w$} is a dependence on *all* of the auxiliary fields, including *ϕ*. Here the density operator $\rho \u0303A(r)$ for a diblock copolymer is given by

The local director defines an operator for the nematic order parameter $Q\u0303ij(r;[{w}])$ according to

where *D* = 3 is the dimensionality of the system. The global order parameter is based on an orientational persistence correlation function,

The order parameter *λ* is non-zero in the presence of long-range order and is thus well suited to characterize the ODT in a symmetric diblock copolymer.^{61}