Field-Theoretic Simulations for Block Copolymer Melts Using the Partial Saddle-Point Approximation

Field-theoretic simulations (FTS) provide an efficient technique for investigating fluctuation effects in block copolymer melts with numerous advantages over traditional particle-based simulations. For systems involving two components (i.e., A and B), the field-based Hamiltonian, Hf[W−,W+], depends on a composition field, W−(r), that controls the segregation of the unlike components and a pressure field, W+(r), that enforces incompressibility. This review introduces researchers to a promising variant of FTS, in which W−(r) fluctuates while W+(r) tracks its mean-field value. The method is described in detail for melts of AB diblock copolymer, covering its theoretical foundation through to its numerical implementation. We then illustrate its application for neat AB diblock copolymer melts, as well as ternary blends of AB diblock copolymer with its A- and B-type parent homopolymers. The review concludes by discussing the future outlook. To help researchers adopt the method, open-source code is provided that can be run on either central processing units (CPUs) or graphics processing units (GPUs).


Introduction
Block copolymers refer to polymeric molecules composed of two chemically-distinct segments, generally denoted as A and B, that are grouped together into separate sections or rather blocks [1,2]. The simplest and most common is the AB diblock copolymer, a linear chain of N = N A + N B segments in which the first N A are of type A and the last N B are of type B. Here, we follow the common practice of defining A and B segments such that they occupy equal-sized volumes, ρ −1 0 . The segments will nevertheless have different statistical lengths, a A and a B , such that the natural end-to-end length of the entire molecule is R 0 = a A N 1/2 A + a B N 1/2 B . For convenience, we define an average segment length a such that R 0 = aN 1/2 . The segments will also differ in their interactions, usually resulting in an incompatibility characterized by a positive Flory-Huggins interaction parameter, χ. If the product χN is sufficiently large, then the A and B components will microphase segregate into a periodically ordered morphology with domain sizes of order R 0 .
The behavior of block copolymer melts is greatly simplified by the fact that it becomes universal in the limit of high molecular weight [3,4]. As a result, diblock copolymer melts are controlled by just four parameters: the segregation χN, the composition f ≡ N A /N, the ratio of segment lengths a A /a B , and the invariant polymerization indexN ≡ a 6 ρ 2 0 N. None of the other details of the system have any impact on the coarse-grained behavior. This is true for experimental systems as well as theoretical models so long as they include the essential physics. Consequently, they can all be mapped onto the standard Gaussian chain model (GCM) [5], which is a minimal model that treats block copolymer systems as incompressible melts of thin elastic threads that interact by pairwise contact forces. At present, simulations at a more realisticN ≈ 10 4 , the resulting phase diagram lacked both a direct gyroid-to-disorder transition and a stable Fddd region.
Near the same time in 2001, Reiter et al. [37] proposed a variant of FTS that avoids this issue by instead simulating H f [W − , w + ], where W + (r) is replaced by the saddle point w + (r). The advantage of this partial saddle-point approximation is that w + (r) is real valued, permitting the use of conventional simulation techniques. The justification for the approximation is that the fluctuations in W − (r) are far more important than those in W + (r), and indeed all evidence so far indicates that the approximation is accurate. Shortly after its introduction, Duchs et al. [38,39] used this variant of FTS to study ternary blends of diblock copolymer with its two parent homopolymers, but then the method received no further attention until 2013 when it was revived by Stasiak and Matsen [40]. The subsequent development, however, has resulted in FTS algorithms capable of problems beyond the reach of traditional particle-based simulations. As such, it is certain to become an important tool in the theoretical study of block copolymer systems. Here, we provide a tutorial on the method followed by demonstrations of its capabilities. To encourage its use, we also provide open-source code for neat diblock copolymer melts that can be readily modified to handle more complex systems.

Field-Theoretic Simulations
This section develops the FTS method starting from the underlying particle-based model, where the Hamiltonian, H p [{r α,i }], is a function of all the particle coordinates, {r α,i }. The first step involves a transformation to the field-based model, where the Hamiltonian H f [W − , W + ] is specified in terms of the fields, W − (r) and W + (r). The evaluation of H f [W − , W + ] requires the statistical mechanics of noninteracting polymers in the fields, which represents one of the two highly computational parts of the simulation. We explain how this is done numerically using fast Fourier transforms. Next, the pressure field, W + (r), is approximated by its saddle point, w + (r), which is the other highly computational part. We describe how to locate w + (r) iteratively using Anderson mixing [41]. The composition field, W − (r), is then evolved using conventional Langevin dynamics. The section concludes by discussing the Morse calibration, which converts the bare χ b parameter used in the simulations to an effective χ parameter corresponding to the standard GCM.

Particle-Based Model
Most FTS have been based on continuous Gaussian chains. However, this requires numerical integration along the chain contours, which then leads to numerical inaccuracies. To avoid this, we will instead model the n polymers of the system by discrete bead-spring chains each with N monomers [42]. The position of monomer i of molecule α will be denoted by the vector r α,i . In this particle-based representation, the Hamiltonian is divided into bonded and nonbonded interactions each expressed in terms of the particle coordinates. The bonded interactions take the simple form which treats them as harmonic springs. Here, the spring constant is expressed in terms of the natural bond length, a. For simplicity, we assume conformational symmetry, but it is straightforward to assign different statistical lengths, a A and a B , to the A-A and B-B bonds, respectively. The nonbonded interactions are restricted to contact forces between the A and B monomers. As such, their energy can be expressed as where χ b is the bare interaction parameter and are dimensionless concentrations of the A and B monomers, respectively. It will become convenient to reexpress the A and B concentrations in terms of the composition,φ − (r), and total concentration,φ + (r), defined bŷ In terms of these new concentrations, the energy of the nonbonded interactions becomes As a result of incompressibility (i.e.,φ + (r) = 1), the first term reduces to a constant, ρ 0 Vχ b /4, and therefore is typically dropped. However, we retain it because of its dependence on χ b , which will have implications regarding an ultraviolet divergence discussed later [42].
To determine the equilibrium behavior of this model, we need to perform statistical mechanics. This requires calculation of the partition function Note that the simple particle-based Hamiltonian in Equation (1) lacks any energy penalty for deviations from melt density, and therefore a Dirac delta functional has to be inserted to enforce incompressibility. Although the form of the partition function remains relatively simple, the integrations over the 3nN particle coordinates are nevertheless impossible to perform, which is the motivation for switching from particle coordinates to fields.

Field-Based Model
The transformation [14][15][16]18] to a field-based representation employs a Hubbard-Stratonovich identity for the Boltzmann weight involving a functional integral over a composition field, W − (r), and a Fourier representation of the Dirac delta functional involving a functional integral over a pressure field, W + (r). Once the identities are inserted into Equation (8), the 3nN integrations over the particle coordinates can be performed analytically leaving just the two integrations over W ± (r). The partition function can then be recast as where the field-based Hamiltonian takes the form involving a relatively simple single-chain partition function for an individual molecule in an equivalent system of noninteracting polymers acted upon by only the fields. The Hamiltonian for Q[W − , W + ] is given by where the bond energy and concentrations are for a single molecule (i.e., n = 1), which allows us to drop the α index on the particle coordinates, r i .

System of Noninteracting Polymers
The single-chain partition function can be conveniently reexpressed as where the Boltzmann weight in the integrand of Equation (13) is written as a product of separate factors for the N − 1 bonds and the N monomers. Here, is the Boltzmann weight for an individual bond and is the Boltzmann weight for the field acting on monomer i. To distinguish the two components of the copolymer, we define γ i = 1 for A monomers (i.e., i ≤ N A ) and γ i = −1 for B monomers (i.e., i > N A ).
To evaluate Equation (15), we define a partial partition function, q i (r), for the first i monomers of the chain with its i'th monomer constrained to position r. It is obtained recursively using starting from q 1 (r) = h 1 (r) [16,43]. Similarly, we define an analogous partial partition function, q † i (r), for the last N + 1 − i monomers, which is obtained by iterating starting from q † N (r) = h N (r). Once both partial partition functions have been calculated, the single-chain partition function is given by The integral for Q can be evaluated using any value of 1 ≤ i ≤ N, but one typically sets i = N. We will also require ensemble averages of the composition and total concentration in the system of noninteracting polymers, which are given by respectively.

Numerical Method
The recursion relation in Equation (18) is calculated using the convolution theorem. This is done by taking a Fourier transform of the propagator (23) and multiplying by The inverse Fourier transform of the product provides the integral which is then multiplied by h i+1 (r) to give q i+1 (r). The second recursion relation in Equation (19) is calculated in an analogous manner. The above equations are generally solved in an orthorhombic L x × L y × L z simulation box of volume V = L x L y L z with periodic boundary conditions. The box is overlaid with an m x × m y × m z grid of uniform spacing (i.e., ∆ α = L α /m α for α = x, y, and z). As such, each grid point corresponds to a volume of V cell = ∆ x ∆ y ∆ z . Spatially-dependent quantities are then converted into arrays over the M = m x m y m z vertices of the grid. Likewise, we define an analogous grid in Fourier space that extends from −π/∆ α to π/∆ α with a spacing of 2π/L α in each direction α. The Fourier-space grid also contains M vertices, but approximately half of them are redundant for real-valued functions, f (r), due to the fact that f (−k) equals the complex conjugate of f (k).

Partial Saddle-Point Approximation
Although the two forms of the partition function in Equations (8) and (11) are mathematically equivalent, the Boltzmann weight in the field-based representation cannot be interpreted in terms of probability thus precluding the use of conventional simulation techniques. This is because the integration of the pressure field is performed along the imaginary axis from W + (r) = −i∞ to +i∞, which makes the Boltzmann weight complex valued. However, this problem can be overcome by evaluating the integral over W + (r) using the saddle-point approximation. Given that the Boltzmann weight is an analytic function of W + (r), the integration path can be deformed in the complex plane so as to pass through the saddle point of the Hamiltonian along a trajectory of constant phase, which concentrates the integration to the saddle-point region. Ignoring irrelevant prefactors, the partition function reduces to where w + (r) denotes the value of W + (r) at the saddle point. As such, w + (r) is determined by solving DH f /DW + = 0, which equates to the mean-field incompressibility condition φ + (r) = 1 (27) Not only does the partial saddle-point approximation reduce Z to a single functional integration, w + (r) is real valued and thus the statistical mechanics can then be performed using standard techniques. The approximation does, however, require w + (r) to be reevaluated every time W − (r) changes.

Anderson Mixing
The saddle point, w + (r), is updated iteratively using the Anderson-mixing scheme [41,[44][45][46]. To facilitate this, we label the saddle-point values by w (k) +,j , where j = 1, 2, . . . , M is the grid index and k = 1, 2, 3, . . . is the iteration index. The iterations begin anew from k = 1 after each change in the composition field, W − (r), generally starting with w (1) +,j equal to the saddle-point solution corresponding to the preceding composition field.
The first step of the k'th iteration is to evaluate the deviation from incompressibility at each grid point, j. The overall error is then quantified by If the error exceeds some given tolerance, then an improved estimate, w (k+1) +,j , is obtained from the preceding n r iterations. This is done by evaluating the symmetric matrix and vector where m, n = 1, 2, . . . , n r . From these, we calculate coefficients used to combine the fields and deviations from past histories as respectively. The next iteration of the field is then obtained by mixing the fields and deviations together as follows where λ > 0 is referred to as the mixing parameter. In the absence of any histories (i.e., n r = 0), Anderson mixing reduces to what is generally referred to as simple mixing.
The convergence tends to be faster for larger λ, but it becomes unstable if λ is too large. Increasing the number of histories improves the convergence as well as the stability, but only up to a point. Therefore, based on previous experience [45], we increase the number of histories up to a maximum of n h = 20 while ramping up the mixing parameter to a maximum of N. This is done by setting n r = min(k − 1, n h ) and λ = N(1.0 − 0.9 k ). We typically continue the iterations until ε ≤ 10 −4 .

Langevin Dynamics
A Markov sequence of configurations for the partition function in Equation (26) is generated by evolving the composition field, W − (r; τ), using Langevin dynamics, where thermal noise is provided by random numbers, N (0, σ), generated from a normal distribution of zero mean and σ 2 = 2δτ/V cell ρ 0 variance. To improve accuracy, we apply the predictor-corrector algorithm [47,48]. A predicted field at τ + δτ is first obtained using Equation (36) with and then a corrected field is obtained with evaluated using the predicted field and corresponding composition. Note that the predictor and corrector steps use the same set of random numbers. As usual, the Langevin dynamics are applied for a sufficient amount of time, τ eq , to equilibrate the system before observables are sampled. There is generally significant correlation between successive time steps and so it is common to only sample periodically, such as once every 10δτ. Naturally, one is interested in the physical quantities of the particle-based model. It is therefore necessary to relate these quantities to those of the field-based model. The average composition, for example, is given by [37] φ − (r) = − 2 and the two-point correlation function is given by [37] φ − (r)φ − (r ) = 4 It can be advantageous to evaluate the average composition using φ − (r) = φ − (r) , because the composition in the system of noninteracting polymers fluctuates less than the composition field. However, be aware that φ Another quantity of interest is the structure function, which is proportional to the Fourier transform of the two-point correlation function in Equation (40).

Effective Interaction Parameter
Because we switched from continuous to discrete chains, our results need to be mapped back onto the standard GCM. Furthermore, there is an ultraviolet (UV) divergence that occurs as ∆ α → 0, which generally requires renormalization of χ and a [29,49]. In this case, however, the divergence in the segment length is avoided by our use of the partial saddle-point approximation. As it turns out, the change in the model and the UV divergence can be dealt with simultaneously by employing the Morse calibration [4,6].
For large N, the interactions are generally weak and therefore the dependence of the effective χ on the bare χ b can be approximated by the linear relationship [6] where the proportionality factor, z ∞ , is given by the relative number of intermolecular contacts in the limit of χ b → 0 and N → ∞. This factor is conveniently expressed as where V cell ρ 0 is the total number of contacts experienced by a given monomer and P i is the probability that two monomers of an infinite chain, separated by i bonds along the contour, occupy the same cell of the grid. For orthorhombic simulation boxes [42], Note that P 0 = 1. In the case of cubic grids with a resolution less than the bond length (i.e., ∆ a), Equation (43) reduces to which is the proportionality factor first used by Stasiak and Matsen [40]. Figure 1 illustrates the quality of the linear calibration by comparing the structure function, S(k), at different grid resolutions forN = 10 4 . Plot (a) demonstrates what happens if the UV divergence is ignored. When χ b N is held constant, the segregation decreases as ∆ = L/m → 0, which is evident by a reduction in the peak height of S(k). Plot (b) corrects for this by comparing S(q) at a fixed z ∞ χ b N, which results in good agreement over the full range of resolutions. Furthermore, the resulting peak height is similar to that of the ROL prediction, which implies that the linear approximation, χ ≈ z ∞ χ b , is reasonably accurate forN = 10 4 . However, the peak is slightly higher than the ROL prediction, which means that the true value of χ is somewhat larger.
For smallerN, it becomes necessary to go beyond the linear approximation of χ. In the Morse calibration, the nonlinear correction is calculated by fitting the peak of the structure function, S(k * ), from FTS to ROL predictions for symmetric diblocks. Beardsley and Matsen [13,48] demonstrated the calibration for ∆ = a and ρ 0 = 8/a 3 , assuming the empirical relationship [4] between the bare χ b of the FTS and the effective χ of the ROL theory. The resulting calibration, z ∞ = 0.7084, c 1 = 1.246, and c 2 = 1.367, is plotted in the inset of Figure 2.
The main plot of Figure 2 illustrates that the quality of the fit is generally good. There is, however, some significant deviation at weak segregations, where the FTS results coincide with the mean-field or rather random-phase approximation (RPA) [22] as opposed to the more accurate ROL theory. This is a consequence of the partial saddle-point approximation, and so it could potentially be remedied by performing FTS with fluctuations in both fields. Nevertheless, the inaccuracy remains acceptable forN 10 3 .    (42) and (46) are plotted in the inset with dashed and solid lines, respectively. Reproduced from Ref. [48].

Applications
The initial applications of FTS have largely focused on two test cases. The first one is neat diblock copolymer melts, which is a natural choice given its relative simplicity and the fact that it is the most thoroughly understood [13,33,34,40,42,48,50]. As mentioned, the greatest impact of fluctuations is on the disordered phase, and thus particular attention is paid to the order-disorder transition (ODT). The second one is ternary blends of diblock copolymer with its two parent homopolymers [38,39,[51][52][53]. Naturally, this system has a far larger parameter space, and so studies have concentrated on symmetric blends where the diblock has a composition of f = 0.5, the two homopolymers have equal polymerizations N h , and the homopolymers have either equal concentrations or equal chemical potentials. For experimental convenience [54][55][56][57][58][59][60][61][62][63][64], the diblock polymerization is generally set to N c = N h /α with α = 0.2, such that the ODT occurs at χN h ≈ 2 over the full range of diblock volume fractions,φ c . In this system, particular attention is paid to a bicontinuous microemulsion (BµE), where the two homopolymers segregate into interweaving microdomains separated by a monolayer of diblock copolymer. The BµE is regarded as a pivotal test of the FTS, as it is a fluctuation-induced phase that is entirely absent from the SCFT phase diagram [38,39,52,65].

Diblock Copolymer Melts
The classical example of a fluctuation effect is the shift in the ODT of symmetric diblock copolymers from the SCFT prediction, (χN) ODT = 10.495. The earliest treatment by Fredrickson and Helfand [24] predicted a shift of 41.0N −1/3 , whereN ≡ a 6 ρ 2 0 N is the invariant polymerization index. Particle-based simulations of different models [8][9][10][11] have since provided the more accurate prediction The ODT is generally found by matching the free energies of the ordered and disordered phases. Although FTS do not provide direct access to the free energy, ensemble averaging a derivative of the Hamiltonian provides the corresponding derivative of the free energy. These can then be integrated to obtain changes in free energy, so long as a phase transition does not occur along the integration path. Lennon et al. [33] were the first to apply thermodynamic integration to FTS. They integrated a composite Hamiltonian, H = λH f + (1 − λ)H ec , from λ = 0 to 1 in order to bridge between the known free energy of an Einstein crystal and the free energy of the polymeric system for a specified value of χ. One integration was performed to evaluate the free energy of the disordered phase below the transition and another to evaluate the free energy of the ordered phase above the transition. The ODT was then located by integrating with respect to the χ parameter.
Beardsley and Matsen [48] located the ODT of symmetric diblocks using a slight variation of the technique where the ensemble average and integration are performed simultaneously [66]. Figure 3 shows the resulting free energy difference between the lamellar and disordered phases, F L − F dis , for diblocks of N = 28 and ρ 0 a 3 = 8 in a cubic simulation box of size L = 24a. The free energy comparison is performed for lamellar periods of D = 8a and D = 24a/ √ 8 ≈ 8.5a, which are the periods that spontaneously form when the disordered phase is quenched into the ordered region. In this case, the thermodynamic integration implies that the shorter period is favored, given that it results in a lower ODT of (χN) ODT = 16.25.   (42) and (46), respectively. The linear χ results in considerable inaccuracy forN 10 4 , but this is largely corrected for by the nonlinear χ. Still, there is a slight deviation for the shorter polymers, which could be due to inaccuracy in the partial saddle-point approximation. It could also be due to a breakdown in universality, which will ultimately occur if N becomes too small.   (46), interaction parameters, respectively. The curve compares the universal prediction from Equation (47) [8]. Adapted from Ref. [48].
A more recent study [42] mapped the ODTs over a range of compositions, f , for diblock copolymers of a fixed N = 90. In that study, the FTS were performed at a sufficiently largeN = 10 4 to justify the use of the linear χ. However, at this relatively largeN, the thermodynamic integration became problematic for the nonlamellar phases. The relative ease of forming defects created noisy free energy curves and the resulting lack of metastability prevented a clear crossing of the curves. Nevertheless, because of the lack of metastability, the ODT could instead be located by monitoring the disappearance of Bragg peaks in the structure function of the ordered phase, S(k), as χN was slowly decreased. Figure 5 illustrates the procedure for diblocks of N A = 36 and N B = 54. A single unit cell of the gyroid phase [67][68][69] was first obtained in a series of different-sized simulation boxes by quenching the disordered phase. The resulting morphology was then periodically repeated creating eight unit cells, and S(k) was evaluated revealing Bragg peaks consistent with the Ia3d symmetry of the gyroid phase. The box corresponding to the strongest peaks was assumed to be commensurate with the preferred period, and the segregation of that box was then reduced in small steps of 0.1 until the peaks disappeared at χN = 13.7. The simulations at χN = 13.8 were then repeated for a prolonged interval to ensure that the gyroid phase was indeed stable at the higher segregation. The fact that it was implies an ODT of (χN) ODT = 13.75 ± 0.05. Notably, these are the first results to show that even weak fluctuations atN = 10 4 are sufficient to produce a direct gyroid-disorder transition, which is consistent with experiments on PE-PEP diblock copolymers [23]. Although the experiments actually observed a transition to perforated lamellae [70], this is understood to be an intermediate state that eventually converts to gyroid [71][72][73].    [74] becomes a potential candidate. Indeed, while the gyroid phase disordered at χN = 13.7, the Fddd phase remained stable down to χN = 13.4 as shown in Figure 6. Furthermore, we have subsequently confirmed that when the segregation of the disordered melt at χN = 13.4 is increased back to χN = 13.5, the first three peaks in S(k) reappear. Although the defects are too long lived for the higher order peaks to return, it is clear that the Fddd phase is stable at χN = 13.5, which then implies an ODT of (χN) ODT = 13.45 ± 0.05. Previous calculations [26] and simulations [34] have suggested that fluctuations destroy the Fddd phase, contrary to its observation in PS-PI diblock copolymers melts [27,28]. These new FTS results are the first to corroborate the experiments. Figure 7 compares the above ODTs and others to the SCFT phase diagram. Consistent with experiments [23], the ODT is shifted upwards relative to the SCFT prediction producing direct transitions between the disordered phase and the ordered lamellar, Fddd, gyroid, cylindrical and spherical phases. The preferred symmetries of the ordered phases were obvious for most of the compositions. The only exception was at N A = 34 and N B = 56 (denoted by a cross in Figure 7), where the cylindrical and gyroid phases both remained stable down to χN = 14.2. Nevertheless, this is consistent with the fact that this composition coincides closely with the cylindrical-gyroid boundary in the SCFT phase diagram.   Figure 7. Order-disorder transitions for diblock copolymer melts atN = 10 4 . The different symmetries of the ordered phases are denoted by circles (lamellar), stars (Fddd), diamonds (gyroid), squares (cylindrical), and hexagons (spherical). The cross marks an ODT for which the relative stability of the cylindrical and gyroid phases was indistinguishable. For comparison purposes, the SCFT phase diagram [21] is overlaid with dashed curves. Adapted from Ref. [42].
Interestingly, the fluctuation-induced shift of the ODT in Figure 7 is reasonably uniform across the range of compositions. This, in fact, agrees with particle-based simulations, which observe shifts of 2.61 and 2.65 at f = 0.25 and 0.5, respectively, [8,12]. However, the shift predicted by FTS is somewhat smaller by about 0.4. Based on the results in Figure 4, this is readily attributed to the inaccuracy of the linear χ = z ∞ χ b .

Ternary Diblock-Homopolymer Blends
The addition of extra molecular species, A-and B-type homopolymers in this case, leads to the possibility of macrophase separation. To deal with this, it is convenient to work in the grand-canonical ensemble, where the concentrations of the different species are controlled by chemical potentials. This generally poses a problem for particle-based simulations because of the difficulty of inserting large macromolecules into a dense melt [75,76]. However, there is no such problem in FTS, and in fact only minor modifications are required to switch between different systems and ensembles [51]. The first term in the Hamiltonian of Equation (12) just needs to be substituted by the free energy of noninteracting molecules corresponding to the system of interest in the ensemble of interest. Similarly, the expressions for φ − (r) and φ + (r) in Equations (21) and (22) are replaced by the composition and total concentration for the same noninteracting system. Figure 8 shows results from a grand-canonical simulation performed in a cubic simulation box of size L = 72R h [53]. In particular, it shows the average copolymer concentration, φ c , as a function of the copolymer chemical potential, µ c , for a blend withN h = 10 4 , N c = 5 × 10 4 , and χN h = 2.31. Here, the linear χ = z ∞ χ b is used and lengths are expressed in terms of the unperturbed end-to-end length of a homopolymer molecule, R h = aN 1/2 h . Interestingly, the plot exhibits a sudden jump in copolymer concentration near µ c ≈ 0.1k B T. Inspection of the configurations (see insets) reveals that this is due to a transition from a homopolymer rich phase at low µ c to a bicontinuous microemulsion (BµE) at high µ c . As a result of the symmetry, the A-homopolymer and B-homopolymer rich phases are identical in free energy. Thus, the low-µ c region corresponds to two-phase coexistence (A+B) and the transition corresponds to three-phase coexistence (A+B+BµE). Based on the jump in φ c , the binodals of the three-phase region are approximatelyφ c ≈ 0.072 and 0.082. However, there is some inaccuracy in these values due to the uncertainty regarding the chemical potential of coexistence. With a sufficiently large simulation box, it is possible to simulate the three-phase coexistence in the canonical ensemble. Spencer and Matsen [53] did so by splicing together configurations of the three phases from the grand-canonical simulation (e.g., the insets of Figure 8), creating an orthorhombic simulation box of size 216R h × 72R h × 72R h . The image at the top of Figure 9 shows an equilibrated configuration with well-defined interfaces separating the three coexisting phases. Note that the simulation box used reflecting boundaries on all sides, and thus there is an A/B interface at x ≈ 0 as well as two homopolymer/microemulsion interfaces at x ≈ 72R h and 144R h . Figure 9a,b plot the composition and copolymer concentration averaged over the y-z plane as a function of the long dimension x. The dashed curves show fits using the standard hyperbolic profile for the interfaces [77]. The interfacial widths can be extracted directly from the fits and the interfacial tensions can be obtained from ensemble averages of appropriate derivatives of the Hamiltonian [51,53,78]. The fits also provide accurate binodals,φ c ≈ 0.0734 and 0.0823, due to the fact that the system is able to equilibrate the concentrations by adjusting the volumes of the coexisting phases.
The boundaries of the A+B+BµE region at χN h = 2.31 are included in the FTS phase diagram shown in Figure 10 [52,53]. As the diagram illustrates, the A+B+BµE coexistence switches to A+B+L coexistence for χN h 2.38 and it narrows in the opposite direction resulting in a tricritical point at χN h ≈ 2.22. Beyond that is a critical line separating A+B coexistence from disordered melts, referred to as the Scott line [79]. The Scott line was calculated using a finite-size scaling analysis [51,80,81], and the L/dis boundary was calculated in much the same way as it was for neat diblock copolymer melts. Note that the lamellar and disordered phases must, in principle, be separated by L+dis coexistence; however, the width of this coexistence was too narrow to resolve.
The SCFT phase diagram for this system [38], corresponding toN h → ∞, is shown in Figure 10 with dashed curves. In this limit, the L/dis boundary and Scott line meet at a Lifshitz point, above which there is A+B+L coexistence. The effect of fluctuations is quite simple. As expected, fluctuations shift both the L/dis boundary and the Scott line to higher χN h . However, the shift of the L/dis boundary is larger, and consequently it intersects the three-phase coexisting region splitting it into A+B+BµE below the point of intersection and A+B+L above. Experiments [54][55][56][57][58][59][60][61][62][63][64] have proposed a somewhat different phase diagram, where the A+B and L regions are separated by a channel of BµE rather than three-phase coexistence. It is difficult to fathom how this proposed phase diagram could possibly converge to the SCFT diagram, which it must, in principle, do asN h → ∞. It is entirely possible that the experiments overlooked the three-phase coexistence due to insufficient annealing. Indeed, a recent experiment by Xie et al. [82] has, in fact, observed macrophase separation. Thus, we believe that the topology of the FTS phase diagram in Figure 10 represents the true equilibrium behavior.

Discussion and Future Outlook
Even with relatively few applications, FTS have demonstrated remarkable potential. For instance, the simulations of three-phase coexistence in Figure 9 involve approximately 10 8 molecules, which is orders of magnitude beyond what is possible in particle-based simulations [75]. This capability has been facilitated by a number of significant advances. For instance, the transition from CPUs to GPUs has allowed the system size to be scaled up considerably [83,84]. A similar improvement was also achieved by switching from Monte Carlo dynamics used in the early studies [13,[38][39][40]50,51,80] to the Langevin dynamics in Equation (36) [48].
There will undoubtedly be further opportunities to improve the FTS method described here. One possibility is replacing the predictor-corrector algorithm with a better scheme [47]. There is also scope for tuning the Anderson-mixing algorithm, such as how the mixing parameter, λ, is adjusted. Alternatively, there may be a better numerical method for iterating w + (r) [85]. Naturally, the speed of the simulations would be enhanced by reducing the polymerization, N, although one must be careful not to go too far or otherwise universality will be lost. Likewise, it might be possible to increase the grid spacing, ∆ α , but it should nevertheless remain fine enough to resolve the relevant coarse-grained details such as the width of the internal A/B interfaces.
One of the main strengths of the field-theoretic approach is its incredible versatility. In particular, it is capable of handling complicated polymeric architectures with a minimal increase in computational effort relative to that of the simple diblock [86,87], which is most certainly not the case for traditional particle-based simulations. Furthermore, it can be adapted to a variety of ensembles [51,88,89], which is very useful when dealing with blends. For the case of AB-type systems, the only change involves the statistical mechanics for the noninteracting polymers. This just requires simple modifications to the first term in the field-theoretic Hamiltonian, Equation (12), and the expressions for φ ± (r), Equations (21) and (22). The extension to three or more chemically-distinct components (e.g., ABC-type systems) is also possible, although the modifications are less trivial [47].
It is important to remember that the version of FTS described here uses a partial saddlepoint approximation, which equates to a mean-field enforcement of incompressibility. Past studies [38,90] have shown that the approximation is accurate, but it becomes less so asN is reduced. The clearest evidence for this is the inability to capture the departure of S(k * ) from the RPA prediction at small χN, as seen in Figure 2. Although complex-Langevin simulations could, in principle, capture this effect, it is uncertain whether, in practice, they can; the compressibility required to access the smaller values ofN in Figure 2 may actually destroy the universality. In any case, the effect remains relatively minor forN 10 3 , and thus the partial saddle-point approximation appears justified. For systems with more than two components, there will be additional fields that may also be imaginary, depending on the relative size of the different interaction parameters [47]. It is likely that they can also be treated accurately with the partial saddle-point approximation, although this remains to be seen.
Naturally, there will be some systems for which it is important to treat incompressibility rigorously, and thus the partial saddle-point approximation would be inappropriate or at least significantly inaccurate. An example is bottlebrush architectures involving backbones with densely-grafted side-chains. In these systems, steric interactions cause the side-chains to form a brush exerting a tension on the backbone, which reduces its flexibility. This steric stiffening of the backbone is completely neglected by mean-field treatments of incompressibility, but the effect does appear to be captured by complex-Langevin FTS [91]. Another example is block copolymer nanocomposites. Without a proper treatment, the nanoparticles can overlap and polymers can penetrate their interiors [92,93]. Although complex-Langevin FTS have been applied to these nanocomposites [35,88,94], they may not adequately account for the more significant steric interactions of bulky nanoparticles. Such steric interactions should remain in the limit of infiniteN, but this cannot possibly be the case in complex-Langevin FTS as they reduce to SCFT in that limit. On the other hand, complex-Langevin FTS do appear capable of handling the steric interactions in polymer solutions involving small solvent molecules [95].
Despite considerable progress, there remain some significant issues to resolve. An important one is determining the exact nature of the UV divergence, or in other words exactly how various quantities depend on the grid resolution. The calibration of the interaction parameter, χ, ensures that melts behave equivalently on different sized grids. Nevertheless, the free energy, for example, still retains a dependence on the grid resolution. This can be calculated analytically for homopolymer melts. In this special case, the fieldbased Hamiltonian reduces to where the plus and minus signs correspond to N A = N and N A = 0, respectively. In both cases, the Hamiltonian is just a system of M independent harmonic oscillators, and therefore its free energy is 1 2 Mk B T ln(χ b M/ρ 0 V) to within an irrelevant constant. The issue is that the free energy diverges as M increases or equivalently as ∆ α → 0. Interestingly, we find that the UV divergence for diblock copolymer melts appears to be independent of composition (i.e., N A ) to within the numerical inaccuracy of our thermodynamic integration. This suggests that the UV divergence remains identical to that of simple homopolymer melts, but this still needs to be confirmed.
The ODT boundaries in Figures 4, 7, and 10 were all determined by comparing ordered and disordered phases in identical simulation boxes, in which case the exact form of the UV divergence becomes irrelevant. This also remains true for comparisons between different sized simulation boxes provided their grid resolutions, ∆ α , are the same. However, this still represents a serious constraint, in large part because the efficiency of the Fourier transforms relies on m α factorizing into small prime numbers [46]. We note that Delaney and Fredrickson [34] have avoided the UV divergence altogether by smearing the molecular interactions, which could likewise be implemented in FTS with the partial saddle-point approximation. In order to obtain universal results, the range of the interactions should be small relative to the width of the internal A/B interfaces. In addition though, the grid resolution needs to smaller than the range of the interactions, probably by a factor of at least four in order to adequately resolve the interaction profile. This implies that the 64 × 64 × 64 grid used to simulate the gyroid phase in Figure 5 would have to be switched to a 256 × 256 × 256 grid, which is currently unfeasible even on GPUs.
Another related issue is how to adjust the size of the simulation box to be commensurate with the equilibrium period of an ordered phase. For the lamellar and cylinder phases, the period can determined by minimizing the free energy with respect to the aspect ratio of the simulation box subject to a constant volume constraint [50,66]. Unfortunately, this strategy cannot be extended to triply-periodic phases. For the gyroid and Fddd phases in Figures 5 and 6, respectively, the box size was adjusted to maximize the Bragg peaks, but this procedure is computationally expensive. The difficulty with using free energy is that changing the volume of the box affects the UV divergence. However, solving the previous problem of how to compare phases in boxes of different grid resolution would provide the means of solving this problem as well.
Now that field-theoretic simulations can handle large numbers of polymers with complicated architectures and multiple species [51,[86][87][88][89], we expect them to gain widespread use. Hopefully, the open-source code provided in the Supplementary Materials will help facilitate this. The increased activity will undoubtedly spawn further developments that will expand its range of applications. As such, FTS is certain to become an invaluable complement to SCFT, which itself is one of the most successful theories in soft condensed matter physics [18].

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/polym13152437/s1, Open-source code for field-theoretic simulations of diblock copolymer melts (read SM.pdf for a description of the files). Funding: This work was supported by NSERC of Canada, and computer resources were provided by Compute Canada.