We demonstrate that combining an emerging approach to game theory with self-consistent mean field theory provides realistic treatments of diblock copolymer phase evolution. We especially examine order-order phase transformations upon quenched temperature change involving hexagonal cylinders, lamellae, and the gyroid. Our findings demonstrate that (i) the game theoretical dynamics produce realistic trajectories for the evolution of the local compositions, (ii) the predicted small-angle scattering follows experimentally observed trends, (iii) nucleation and growth is active when the system is quenched far from the critical point, and (iv) epitaxial growth is manifest. To our knowledge, the methodology presented provides the first merger of mean field game theory and statistical mechanics for soft matter systems, giving a new inroad to studying polymer dynamics.

Since the inception of the study of thermodynamics, nonequilibrium phenomena have both motivated and challenged much of the theoretical developments in this area dating all the way back to Boltzmann’s formulation of his H-theorem.1 Though much progress2 has been made on the theory over the past century and the advent of powerful computational hardware now permits the simulation of many complex systems, the subject of nonequilibrium statistical mechanics remains a subject under vigorous development.3 Elucidating the dynamical evolution of a system that has been perturbed (sometimes far) away from equilibrium often serves as the primary goal for those studying nonequilibrium phenomena. The concepts of entropy and free energy, in addition to enthalpy, facilitate the theoretical description of equilibrium systems, but quantification of these metrics in systems out of equilibrium is problematic. The latter class of systems also displays time asymmetry; one can observe a series of snapshots from the system and deduce their time-ordering.4 This characteristic must emerge from the complex dynamics of the system and special boundary conditions since, famously, the basic mechanical equations of motion do not encode such behavior. Moreover, many systems are not amenable to brute force computational interrogation because the relaxation times for dominant processes are very long and the number of objects needed to capture the relevant physics is large. Coarse-graining from particles to fields mitigates some of the latter difficulties, however. Thus, many formulations for addressing nonequilibrium phenomena produce estimates for the time evolution of the distribution of mass for systems prepared in an equilibrium state and subjected to a sudden change in environment that destabilizes that configuration. These formulations often require application of heuristics and/or recourse to approximate expansions of free energy expressions that are only well established at equilibrium. Along with others discussed briefly below, the time-dependent Ginzburg-Landau formulation (TDGL)5,6 stands prominent among those treatments. Though highly effective for many problems, that method and its existing counterparts typically utilize information from only one equilibrium configuration, the initial condition, and are prone to kinetic trapping in metastable states. Devising a treatment of nonequilibrium dynamics that self-consistently incorporates information about the final equilibrium state, as well as the initial condition, may provide a more robust computational approach to these complex systems and, perhaps, new insight into the statistical physics of this omnipresent class of phenomena without appealing to free energy functionals.

Thus motivated, we herein explore the application of a recent formulation of game theory to the study of nonequilibrium material dynamics. Mean Field Game Theory (MFGT)7–10 aims to address systems with large numbers of players and itself draws from concepts developed in statistical physics. As discussed in detail below, the components of this theory that make it appealing for application to nonequilibrium physical systems are its use of both past and future states to predict the path linking them and its construction from a principle akin to that of least action. To demonstrate the potential that this formulation of dynamical systems holds for application to the study of nonequilibrium material phenomena, we apply MFGT to order-order phase transitions in diblock copolymers. Nonequilibrium phase dynamics in heteropolymers are of interest because transient morphological defects may be trapped during materials manufacturing and potentially useful metastable phases may be achieved by exerting control over the phase separation process. Moreover, materials may be exposed to environmental changes during their application that effect modifications to their mesoscale organization. Block copolymers exhibit well-established equilibrium mesophases and provide model systems for exploring the dynamics of phase transitions. Much effort exploring order-order transitions, in particular, is reported in the literature. However, simulating nonequilibrium dynamics in concentrated high polymer systems remains a very challenging problem due to the slow dynamics associated with these materials. Modeling these systems with particle methods (e.g., molecular dynamics) for useful time scales remains beyond our computational grasp for realistic molecular weights, as the relaxation time varies as the length of the chains to the power 3.4.11 Nevertheless, these materials play major technological roles and are ubiquitous.

At equilibrium, semianalytical methods based on self-consistent mean field theory (SCFT)12 built from path integral formulations of coarse-grained models provide the principal computational approach to predicting the morphologies of these materials. Numerous approaches13–33 have been developed to address the dynamics of polymer ordering and fluctuations and remains a challenging area of study as interest in new mesophases and metastable phases grows.13 Exemplary of these approaches are methods drawn from SCFT, dynamical density functional theory, TDGL, complex Langevin formulations, and minimum energy pathway analysis. These methods capture much behavior that is suggested by experimental investigations, especially the role of nucleation and growth in the dynamics. Matsen14 applied a formulation of SCFT to study the transition pathway between the cylinder and gyroid phases. He argues that the free energy linking these two phases provides a pathway for the transition that supports nucleation and growth via epitaxy. Fraaije15 presented an early dynamic density functional theory based on a generalized Cahn-Hilliard dynamics. Fraaije noted transient defect states during the disorder-to-order phase transition. Qi and Wang19 used TDGL to study the transitions involving the hexagonal and body-centered-cubic phases, again finding evidence of epitaxy. Ganesan and Fredrickson20 applied the complex Langevin approach to study fluctuations about the SCFT solution and studied the impact of these fluctuations on the disorder-order phase transition. Ji et al.13 studied the kinetic pathways connecting lamellar and gyroid phases, using a combination of the “string method” and SCFT. Their results suggest a nucleation and growth mechanism and the formation of intermediate metastable states, including perforated lamellae, the O70 phase, and tetragonally perforated layers. Though successfully capturing much of the known phenomena in copolymer phase transitions, the complexity of these methods sometimes necessitates detailed model specialization for each transition studied and all rely, either explicitly or implicitly, on a free energy functional. We propose a simpler approach herein that should facilitate the rapid exploration of a wide range of transitions and formulate the models without recourse to a free energy functional.

We apply the Hamilton-Jacobi-Bellman (HJB)/Fokker-Planck (FP) formulation of MFGT as constructed by Huang, Malhamé, and Caines7 and by Lasry and Lions.8 More details on this branch of game theory may be found in Refs. 7–10. Imagine a set of players that are each obeying the stochastic differential equation dRi = vi(t)dt + σdw. Here, Ri is the position of player i, vi is the velocity choice at time t, and w is a Wiener process introducing noise into the system with standard deviation σ. The objective of the game is to determine the optimal set of velocities or policies vi such that a cost function C across the time frame T of the game is minimized and a terminal cost CT achieved. One assumes a kinetic energy penalty for players with mass M and a potential energy F that is a functional of the density of the players ρ. For a single class of players, from the Bellman principle of optimality, one has C=0TM2vi2F[ρ]dτ+CT. The brackets indicate averages over the noise. One aims to minimize the integral (the action, viz., principle of least action) across vi and achieve the terminal cost. From the Bellman principle, one may show that the stochastic HJB equation dictates the evolution of vi (indirectly) backward in time: ∂u/∂t − 1/(2M)(∂u/∂x)2 + (σ2/2)∇2u = F with v(t) = −1/M∂xu. One works backward from the target condition CT in this formulation. The density ρ, in turn, obeys the Fokker-Planck (FP) equation34 since the underlying positions are dictated by the stochastic differential equation above. This must be solved forward in time: ∂ρ/∂t + (vρ)/∂x − (σ2/2)∇2ρ = 0. Thus, finding solutions to this class of mean field games corresponds to finding the density distribution of players ρ by self-consistently solving both the FP and HJB equations for given starting and ending conditions, under the constraint of a time-dependent potential F. We propose that this formulation may prove effective in describing physical systems if one identifies the individual particles of a material with the players of a game, especially if only the probabilistic evolutionary dynamics are of interest.

The difficulties with the conflicting treatments of time have been partially mitigated by the recent work of Swiecicki9 and co-workers that showed a simple operation transformed both the HJB and FP equations into Schrödinger-like equations with imaginary time; this set of equations is very similar to that solved in the SCFT. Specifically, those authors showed that applying the transformations Φ=expu/(Mσ2) and Γ = ρ/Φ leads to the following equations:

Mσ2tΦ=Mσ422Φ+F[ρ]Φ,
(1)
Mσ2tΓ=Mσ422Γ+F[ρ]Γ.
(2)

Note that Eq. (1) must be solved backward in time, while Eq. (2) is solved forward in time. The two equations depend upon both space x and time t and are coupled through the density ρ. Because ρ changes with both space and time, the potential F is also space- and time-dependent. Noting the similarity of Eqs. (1) and (2) with those describing the propagators in the equilibrium polymer field theory, we solve the above equations with an Anderson update and Fourier method similar to that developed for SCFT.35 Specifically, we apply the following algorithm:

  1. Read the initial and final density fields, ρ(x,0) and ρ(x,T), respectively.

  2. Set the initial guess for ρ(x,tT)=ρ(x,0).

  3. Set the initial guess for Fρ(x,t), described below.

  4. Set the initial guess Φ(x,T)=expFρ(x,T)Mσ2.

For n steps:

  • Solve Eq. (1) backwards to find Φ(x,t) for 0 ≤ t < T.

  • Solve Eq. (2) forwards to find Γ(x,t) for 0 < tT with Γ(x,0)=ρ(x,0)Φ(x,0).

  • Update the estimate of ρ(x,t)(1ϵ)ρ(x,t)+ϵΦ(x,t)Γ(x,t) for 0 < t < T.

  • Normalize ρ(x,t) for 0 < t < T.

  • Update the estimate of F(x,t).

We apply this construction of MFGT to the simple case of a diblock copolymer. Though a straightforward extension of the scheme employed here can account for multiple chemical species, discussed briefly below, we only track the density of one block and enforce incompressibility by density normalization. For simplicity, we assume that M ≡ 1 and σ ≡ 1. We estimate the action as At,x12ΔρΔt2F and use it as a measure of convergence. Note that other metrics of numerical convergence could be used, such as total energy, and no special significance should be assigned to our use of A for this purpose. We appeal to SCFT,12 following the method outlined in Ref. 35, to set the equilibrium starting and ending states. However, one may also use density models reconstructed from experimental data, such as small-angle scattering measurements.

We especially focus on a diblock copolymer with fraction f = 0.385 for the tracked block because this division of the diblock yields a polymer that can access three different ordered phases by simply changing the value of χN, the product of the Flory-Huggins “χ-parameter” and the number of statistical segments in the polymer. The lamellae-to-hexagonal, the hexagonal-to-gyroid, the gyroid-to-lamellae, and their reverse transitions are considered. Guided by the mean field phase diagram for diblock copolymers,36 we use SCFT to calculate the lamellae phase (L) at χN = 40.0, the gyroid phase (G) at χN = 15.0, and the hexagonal phase (H) at χN = 13.0. For all phases, Nv = 643 grid elements separated by a spacing of 0.149 comprise the SCFT simulation lattice, with all lengths given in units of the radius of gyration of the polymer. Once converged, the SCFT results are replicated to form a Nv = 2563 simulation box. Both simulation sizes are used to eliminate any concerns over finite-size effects. Figure 1 illustrates the Nv = 2563 equilibrium boundary states for the dynamical calculations by presenting the color-coded density for all elements with ρ ≥ 0.5 (blue corresponds to ρ = 0.5 and red indicates ρ = 1.0). From the left, the simulation boxes represent the H, G, and L phases. In order to test the impact of the choice of T, we carry out a set of simulations with T = 30, 40, 50, 60, 70, 80, 90, and 100 in arbitrary units of time for the Nv = 643 simulation boxes; we always use T = 50 for the Nv = 2563 simulations. In order to establish that the dynamics are dominated by the distribution of mass ρ(x,t) given by the boundary states at t = 0 and t = T, we apply a potential of the form Fρ(x,t)=k(ρ(x,t)ρ(x,T))2+ρ(x,T) with values of k = 0, 0.2, 0.4, 0.6, 0.8, and 1 for the Nv = 643 simulations; k is always set to zero for the Nv = 2563 systems. The quadratic term simply serves as a probe potential designed to establish that the linear term in ρ(x,t) dominates the dynamics of the trajectory. That is, we use it to establish that the dynamics are dominated by the final equilibrium state and that the backward propagation of information from that state suffices to capture the phenomenology of these transitions.

FIG. 1.

Snapshots of the equilibrium boundary states as constructed from SCFT and subsequent replication to Nv = 2563. The density fields ρ are presented for all elements with ρ ≥ 0.5 and color coded with blue corresponding to ρ = 0.5 and red corresponding to ρ = 1.0. From the left, the hexagonal (H), the gyroid (G), and lamellae (L) are shown. All simulation snapshots are prepared with ParaView.58 

FIG. 1.

Snapshots of the equilibrium boundary states as constructed from SCFT and subsequent replication to Nv = 2563. The density fields ρ are presented for all elements with ρ ≥ 0.5 and color coded with blue corresponding to ρ = 0.5 and red corresponding to ρ = 1.0. From the left, the hexagonal (H), the gyroid (G), and lamellae (L) are shown. All simulation snapshots are prepared with ParaView.58 

Close modal

Figures 2–4(c) illustrate snapshots from intermediate time steps (t/T = 0.2, 0.4, 0.6, and 0.8) for calculations of the respective order-order phase transitions in both the forward and reverse directions. In particular, transitions between G and H phases are illustrated in Fig. 2. For the G-to-H transition illustrated in Fig. 2(a) (Multimedia view), one observes an initial reduction in ρ around the skeleton of the G phase, followed by an epitaxial formation of rough cylinders and an ultimate smoothing into the H phase. This epitaxial growth is consistent with the predictions of Matsen14 and has been observed experimentally.37–44 A nearly time-symmetric reversal of this phenomenon appears active in the H-to-G transition snapshots shown in Fig. 2(c) (Multimedia view). However, as discussed further below, this symmetry is not perfect. Similarly, Figs. 3(a) and 3(c) (Multimedia view) illustrate the transitions between the L and G phases. Figure 3(a) (Multimedia view) shows that during the L-to-G transition, ρ in the lamellae first thins as the lamellae begin to form undulations, which was also observed experimentally by Hajduk et al.45 Those undulations subsequently form protrusions that connect the lamellae in an epitaxial relation driven by the G phase. This leads to an enrichment of ρ in locations along the former lamellae planes that ultimately spreads out to establish the final G phase. Here again, at first glance, the reverse sequence appears active in the G-to-L transition illustrated in Fig. 3(c) (Multimedia view). The phenomena of epitaxial growth and morphological protrusion/pinching are most easily observed in the H and L transitions illustrated in Figs. 4(a) and 4(c) (Multimedia view). During the H-to-L transition illustrated by the simulation snapshots in Fig. 4(a) (Multimedia view), one observes an initial thinning of ρ and, curiously, a migration in cylindrical form of those features not in registry with the ultimate lamellae structure. Upon merging with cylinders closer to the final epitaxial location, undulating lamellae are again formed. The mass then redistributes within those sheets to form the final lamellar state. While this transition dynamic has not been observed experimentally, analogous dynamics in the reverse L-to-H transition were observed via scattering and microscopy experiments46 and are demonstrated in Fig. 4(c) (Multimedia view). The simulations presented here clearly show the mechanism proposed by Hajduk et al.,46 with budding and pinching off of cylinders from the lamellae that migrate to their ultimate location in a hexagonally packed structure.

FIG. 2.

Results for the order-order transitions involving the hexagonal (H) and gyroid (G) phases. (a) Snapshots of the G-to-H transition. (b) A snapshot of the barrier state in the G-to-H transition is shown on the left, while the barrier state for the H-to-G transition is shown on the right. (c) Snapshots of the H-to-G transition. All snapshots are thresholded and color coded as in Fig. 1. (d) The total pairwise interaction energy E normalized by the number of simulation grid elements Nv is presented as a function of time divided by the total time of the simulation T for the G-to-H transition; data for all combinations of Nv and T are shown. (e) The analogs of the data shown in (d) are presented for the H-to-G transition. (f) The differences between the densities of the trajectories δρ normalized by Nv as a function of the normalized time for the G-to-H and H-to-G transitions; data for Nv = 643 and T = 100 are shown in red, while the corresponding results for Nv = 2563 and T = 50 are shown in blue. (g) The estimated small-angle scattering profiles at different times are shown for the G-to-H transition. (h) The scattering profiles at selected times for the H-to-G transition. All plots are constructed with the Matplotlib library.59 Multimedia views: (a): https://doi.org/10.1063/1.5081829.1; (c): https://doi.org/10.1063/1.5081829.2

FIG. 2.

Results for the order-order transitions involving the hexagonal (H) and gyroid (G) phases. (a) Snapshots of the G-to-H transition. (b) A snapshot of the barrier state in the G-to-H transition is shown on the left, while the barrier state for the H-to-G transition is shown on the right. (c) Snapshots of the H-to-G transition. All snapshots are thresholded and color coded as in Fig. 1. (d) The total pairwise interaction energy E normalized by the number of simulation grid elements Nv is presented as a function of time divided by the total time of the simulation T for the G-to-H transition; data for all combinations of Nv and T are shown. (e) The analogs of the data shown in (d) are presented for the H-to-G transition. (f) The differences between the densities of the trajectories δρ normalized by Nv as a function of the normalized time for the G-to-H and H-to-G transitions; data for Nv = 643 and T = 100 are shown in red, while the corresponding results for Nv = 2563 and T = 50 are shown in blue. (g) The estimated small-angle scattering profiles at different times are shown for the G-to-H transition. (h) The scattering profiles at selected times for the H-to-G transition. All plots are constructed with the Matplotlib library.59 Multimedia views: (a): https://doi.org/10.1063/1.5081829.1; (c): https://doi.org/10.1063/1.5081829.2

Close modal
FIG. 3.

Results for the order-order transitions involving the lamellae (L) and gyroid (G) phases. (a) Snapshots of the L-to-G transition. (b) A snapshot of the barrier state in the L-to-G transition is shown on the left, while the barrier state for the G-to-L transition is shown on the right. (c) Snapshots of the G-to-L transition. All snapshots are thresholded and color coded as in Fig. 1. (d) The total pairwise interaction energy E normalized by the number of simulation grid elements Nv is presented as a function of time divided by the total time of the simulation T for the L-to-G transition; data for all combinations of Nv and T are shown. (e) The analogs of the data shown in (d) are presented for the G-to-L transition. (f) The differences between the densities of the trajectories δρ normalized by Nv as a function of the normalized time for the L-to-G and G-to-L transitions; data for Nv = 643 and T = 100 are shown in red, while the corresponding results for Nv = 2563 and T = 50 are shown in blue. (g) The estimated small-angle scattering profiles at different times are shown for the L-to-G transition. (h) The scattering profiles at selected times for the G-to-L transition. Multimedia views: (a): https://doi.org/10.1063/1.5081829.3; (c): https://doi.org/10.1063/1.5081829.4

FIG. 3.

Results for the order-order transitions involving the lamellae (L) and gyroid (G) phases. (a) Snapshots of the L-to-G transition. (b) A snapshot of the barrier state in the L-to-G transition is shown on the left, while the barrier state for the G-to-L transition is shown on the right. (c) Snapshots of the G-to-L transition. All snapshots are thresholded and color coded as in Fig. 1. (d) The total pairwise interaction energy E normalized by the number of simulation grid elements Nv is presented as a function of time divided by the total time of the simulation T for the L-to-G transition; data for all combinations of Nv and T are shown. (e) The analogs of the data shown in (d) are presented for the G-to-L transition. (f) The differences between the densities of the trajectories δρ normalized by Nv as a function of the normalized time for the L-to-G and G-to-L transitions; data for Nv = 643 and T = 100 are shown in red, while the corresponding results for Nv = 2563 and T = 50 are shown in blue. (g) The estimated small-angle scattering profiles at different times are shown for the L-to-G transition. (h) The scattering profiles at selected times for the G-to-L transition. Multimedia views: (a): https://doi.org/10.1063/1.5081829.3; (c): https://doi.org/10.1063/1.5081829.4

Close modal
FIG. 4.

Results for the order-order transitions involving the hexagonal (H) and lamellae (L) phases. (a) Snapshots of the H-to-L transition. (b) A snapshot of the barrier state in the H-to-L transition is shown on the left, while the barrier state for the L-to-H transition is shown on the right. (c) Snapshots of the L-to-H transition. All snapshots are thresholded and color coded as in Fig. 1. (d) The total pairwise interaction energy E normalized by the number of simulation grid elements Nv is presented as a function of time divided by the total time of the simulation T for the H-to-L transition; data for all combinations of Nv and T are shown. (e) The analogs of the data shown in (d) are presented for the L-to-H transition. (f) The differences between the densities of the trajectories δρ normalized by Nv as a function of the normalized time for the H-to-L and L-to-H transitions; data for Nv = 643 and T = 100 are shown in red, while the corresponding results for Nv = 2563 and T = 50 are shown in blue. (g) The estimated small-angle scattering profiles at different times are shown for the H-to-L transition. (h) The scattering profiles at selected times for the L-to-H transition. Multimedia views: (a): https://doi.org/10.1063/1.5081829.5; (c): https://doi.org/10.1063/1.5081829.6

FIG. 4.

Results for the order-order transitions involving the hexagonal (H) and lamellae (L) phases. (a) Snapshots of the H-to-L transition. (b) A snapshot of the barrier state in the H-to-L transition is shown on the left, while the barrier state for the L-to-H transition is shown on the right. (c) Snapshots of the L-to-H transition. All snapshots are thresholded and color coded as in Fig. 1. (d) The total pairwise interaction energy E normalized by the number of simulation grid elements Nv is presented as a function of time divided by the total time of the simulation T for the H-to-L transition; data for all combinations of Nv and T are shown. (e) The analogs of the data shown in (d) are presented for the L-to-H transition. (f) The differences between the densities of the trajectories δρ normalized by Nv as a function of the normalized time for the H-to-L and L-to-H transitions; data for Nv = 643 and T = 100 are shown in red, while the corresponding results for Nv = 2563 and T = 50 are shown in blue. (g) The estimated small-angle scattering profiles at different times are shown for the H-to-L transition. (h) The scattering profiles at selected times for the L-to-H transition. Multimedia views: (a): https://doi.org/10.1063/1.5081829.5; (c): https://doi.org/10.1063/1.5081829.6

Close modal

A specific intermediate state for each transition merits closer scrutiny. Several experimental and theoretical studies of order-order phase transitions in copolymers have concluded that a nucleation and growth mechanism is active in the dynamics.13,14,33,40–43,45,47,48 This implies the existence of an energy barrier along the transition pathway. To evaluate whether this is true in our proposed model, we calculate postsimulation the field theoretical pairwise energy for the system as EjχNiNvρi,j1ρi,j, where i corresponds to different spatial locations within the simulation box and j corresponds to time along the trajectory. From this calculation, one may readily identify a maximum value for E at a given time during the transition for all transitions studied here, a state henceforth referred to as the “barrier state.” Figures 2–4(b) contain images from the corresponding barrier states, while Figs. 2–4(e) contain plots of E normalized by Nv against time normalized by T. Note that in those plots, data from all combinations of Nv, T, and k are shown; once normalized by the system size and the time frame for the transition, no discernible dependence on any of these parameters is detected. The left-hand side of Fig. 2(b) presents the snapshot of the barrier state in the G-to-H transition. Compared to the corresponding image of the barrier state for the H-to-G transition shown on the right-hand side of Fig. 2(b), little difference can be detected between the two transitions. However, examining the plots of E/Nv vs t/T contained in Figs. 2(d) and 2(e) reveals asymmetry in both the time of the transition and the barrier height. The barrier state arrives later and the barrier energy is lower for the G-to-H transition compared to the H-to-G transition. Analogous data are presented in Figs. 3(b), 3(d), and 3(e) for the L-to-G and G-to-L transitions. Here, however, notable differences are visible between the image of the L-to-G barrier state, shown left, and the G-to-L barrier state, shown right. While the former shows clear signatures of the L state with a superposed G morphology, the latter shows a less fully developed analog and thinner mass draped on the morphological scaffold. As above, one also notes asymmetry in both the barrier height and times; the barrier is much lower and arrives later in the L-to-G transition than the G-to-L. The barrier arrives later in the L-to-G transition compared to that calculated by the string method and reported by Ji13 and co-workers. The string method approach of that study, however, requires insertion of a preformed nucleus, making a one-to-one trajectory comparison impossible. Figures 4(b), 4(d), and 4(e) present similar results for the H-to-L and L-to-H transitions. Such asymmetry between forward and reverse H-to-G, L-to-G, and H-to-L transitions has been observed experimentally, especially in the works of Hajduk et al.45 and Floudas et al.41 Comparison of the barrier state for the H-to-L process, shown left, and the L-to-H, shown right, reveals significant differences in the concentrations of those cylinders that fall close to registry with the lamellae sheets. We also note that the profiles for E/Nv vs t/T are of the same form as that we observe for the other order-order transitions. One may have expected that passage through the G state that connects the L and H morphologies on the phase diagram should be captured by the dynamics. However, the simulation method in its current form does not explicitly consider temperature. Moreover, we are unable to find any experimental support that such an intermediate exists.

In order to further probe the time asymmetry suggested by the barrier states, we calculate the differences between the dynamical states for the pairs of transitions involving the same ordered states. To quantify this, we calculate δρiNvjTρi,jρi,Tj2, where ρi,Tj indicates the time-reversed states for the corresponding forward order-order transition. For example, the density for the L-to-H transition in element i = 2 at time step j = T − 1 is subtracted from the density for the H-to-L transition in element i = 2 at time step j = 1. Thus, we simply compare the densities for the two transitions as if we reversed the clock for one of them. The results of the calculations are presented in Figs. 2–4(f) for all transitions studied herein. Note that we calculate δρ for both the Nv = 643, T = 100 and Nv = 2563, T = 50 simulations. Once δρ is normalized by Nv and T, as shown in the plots, no differences can be detected between the simulations with differing values of Nv and T, suggesting that the variations observed cannot be attributed to simple numerical errors. Nevertheless, all values obtained for δρ/Nv are small. Figure 2(f) displays the results for the trajectories for the transitions involving the H and G morphologies. As suggested by the images in Fig. 2(b), very little variation is observed between the H-to-G and G-to-H transitions, though some variation is noted near the time boundaries, and the variation is itself nearly symmetric along the transition pathway. In contrast, the transitions involving the L phase display more than an order-of-magnitude higher variation, and the variation along the pathway is asymmetric. These results strongly suggest that nonequilibrium processes involving the L phase display time asymmetry; one should be able to view a series of snapshots from these transitions and discern the time ordering of the observations.

To facilitate further comparison with experimental findings, we estimate the small-angle radiation scattering signatures for the transitions by I(q)=|ρ^(q)|2, where ρ^q is the Fourier transform of the density distribution ρ and the average is over all scattering vectors q with the same magnitude. For clarity, these signatures are only presented for t/T = 0.02 and 0.98, as well as the barrier state for each transition, in Figs. 2–4(h). Each scattering curve possesses multiple peaks that reflect the ordered morphology. For example, in the L-to-G transition shown in Fig. 3(g), the t/T = 0.02 curve is consistent with the scattering signature for the L phase; the primary peak q* = 0.25 gives higher order reflections at 2q*, 3q*, 4q*, 5q*, and 6q*. In the same plot, the t/T = 0.98 curve reflects the G phase; the primary peak is located at q* = 0.28, and higher order reflections are observed at 3q*,2q*,7q*, and 13q*. The scattering curve for the barrier state at t/T = 0.65 has features that reflect the mixed lamellae-gyroid morphology, with higher order peaks at values of 3q*,2q*, and 3q*. Analogous trends are observed for the reverse transition shown in Fig. 3(h), as well as in the G-H [Figs. 2(g) and 2(h)] and H-L [Figs. 4(g) and 4(h)] transitions. These mixed-phase morphology results are consistent with scattering signatures observed experimentally during the order-order transitions studied here,41,45,46,49,50 as well as in other order-order transitions.51–53 

Though these results are promising, several extensions may be made to the formulation presented herein to improve the fidelity of the dynamics to actual physical systems, as well as accommodate external processes coupled to the phase evolution. As noted above and reflected in plots (d)–(f) of Figs. 2–4, the time scales for the dynamics prove to vary linearly with the artificially imposed total time of the phase transition T in the current form of the model. This property derives from our simple application of the FP and HJB equations [as transformed into the form of Eqs. (1) and (2)] with a constant diffusion coefficient. Real polymeric materials display richer dynamics due to the local bonding of the individual chain segments, chain entanglement, and many-body interactions. These constraints result in correlated motions along the chains and hydrodynamic coupling between spatially distant portions of the material. This physics is lost in our formulation, and we must interpret the dynamics presented herein as a zeroth order approximation dominated by the longest relaxation processes active in the material, specifically long-time diffusive motion. This is a reasonable treatment, if the focus is on the coarse-grained phase evolution.12 However, as noted by Fredrickson,21 momentum relaxation (associated with hydrodynamics) also evolves slowly. This latter phenomenon may be expected to become important in cases involving applied fluid forces, such as in rheological models, but is less of a concern for the simple quiescent phase separation studied here. Moreover, we have herein only studied the evolution of one order parameter, the density of one monomer type in the copolymer. We justify this by noting the incompressible nature of the material during the process examined. However, many interesting polymeric processes require tracking more than a single order parameter. The route to expanding the model presented here to facilitate the study of these more complicated processes should prove straightforward, even if somewhat technically challenging to implement.

Here, we briefly outline a path for elaborating on the presented formulation. The underlying Langevin equation of motion for the players given above may be written in the general form dRi/dt=h(x,t)+g(x,t)Γ(x,t). The first term on the right-hand side, h(x,t), is a velocity or drift term. The second term on the right, g(x,t)Γ(x,t), mimics the thermal noise in the system. Typically, g(x,t) is set to a constant that reflects the temperature and viscosity of the system, while Γ(x,t) is assumed to average to zero and to be correlated locally in time and space (<Γ(x,t)Γ(x0,t0)>δ(tt0)δ(xx0)). As discussed in detail in the work of Risken,34 moving from discrete particle trajectories to probability distributions W(x,t), and thus from Langevin to FP equations, involves identifying the Kramers-Moyal expansion coefficients with functions of h(x,t) and g(x,t). For a single particle, one has the FP equation W(x,t)/t=/xD(1)(x,t)+2/x2D(2)W(x,t), where D(1)=h(x,t)+g(x,t)/xg(x,t) and D(2)=g2(x,t). We have taken the steps of identifying W(x,t) with the density ρ(x,t), D(1) with the velocity to be learned by the MFGT algorithm, and D(2) with a simple diffusion constant. Similarly, these choices inform the coefficients and role of the velocity in the HJB equation. Modifying these terms may account for polymer processes that relax more quickly than the coarse density evolution studied here and provide a link to dynamically imposed fluid flow. Kawasaki and Sekimoto23 provided the foundation for the former, elaborating on expressions analogous to D(2) to account for chain dynamics. Fredrickson21 developed a theory based on the complex Langevin model to account for momentum transfer during hydrodynamic flow processes. We propose that replacing the continuity equation of motion for the density in the formulation with an appropriate FP equation coupled to a control expression may provide a path to introducing this physics into the model proposed here. Several approaches to studying system dynamics dependent upon multiple order parameters may extend MFGT to address more complicated problems; the most straightforward is simply appealing to a heuristic coupling through the potential F as is constructed in the Cahn-Hilliard model (“model B” in the nomenclature of Chaikin and Lubensky6). In summary, the key idea for applying MFGT lies in coupling a control expression to an equation of motion in which the position choices are to be learned by back-propagating the terminal state while obeying the optimality principle; depending upon the complexity of the physics of interest, this may entail elaborating both the FP and HJB expressions by informing them with the appropriate functions for D(1) and D(2), or coupling sets of many forward motion and control equations.

We herein present a new method for studying the nonequilibrium dynamics operative in copolymer order-order phase transitions. This method is motivated by computational expediency and supported by an underpinning optimality principle similar to the principle of least action. As demonstrated above, this simple approach surprisingly reproduces much of the known phenomenology of order-order phase transitions in block copolymer systems. These include mechanistic details such as epitaxial growth of the morphological features, the presence of energy barriers that give rise to the nucleation and growth dynamics, and the asymmetry of the transition barriers separating the ordered states. However, the physical notions driving this approach require comment. The procedure self-consistently calculates the trajectories of mass distributions using information about both the initial and final states of the system; knowledge of the future state of the material informs the path followed by the material throughout the timeframe of the transition. In this sense, the method draws from a weak form of retrocausality, a concept that has been explored to both clarify the physical meaning of calculational procedures in other domains of physics and address paradoxes that arise at the foundation of quantum mechanics.54 For example, as noted by Eddington,55 the Born rule for calculating probabilities from the underlying probability amplitudes can be interpreted as combining information from waves propagating both forward and backward in time. However, as noted by Wharton,56 the idea that one may use knowledge of the future to calculate the present and past seems to fall into the category of metaphysics, even if it may be used to resolve apparent inconsistencies in otherwise sound theories; if we require knowledge of the future to apply the concept, we cannot use it to make meaningful predictions. However, there are systems for which we do know, at least probabilistically, the future state. As noted by Prigogine,57 equilibrium states act as an attractor for nonequilibrium dynamics. Therefore, equilibrium statistical mechanics permits the prediction of the future with high precision in many systems. Our scheme aims to predict the nonequilibrium dynamics for the approach to the equilibrium state. Thus, the use of two bounding equilibrium conditions in the approach outlined herein utilizes information from both the past and the future to provide a computational procedure for predicting the intervening dynamics by anchoring the calculation with two well-understood physical states. The self-consistent use of two equilibrium points distinguishes the method from other treatments of nonequilibrium dynamics that are often informed by a single equilibrium point and rely upon a free energy functional. Though one may use symmetry heuristics, TDGL, for example, is strictly formulated around a Landau expansion of the equilibrium theory’s free energy functional.5 The method proposed herein does not appeal to any such free energy functional. Moreover, this approach need not be bounded by equilibrium states. Rather, the dynamics of the formation of experimentally observed metastable configurations may also be probed.

The principal advantages offered by this theoretical approach are that (i) information from two equilibrium states may be used to inform the calculation of the nonequilibrium dynamics; (ii) implementation of a computational solution is relatively simple and here we use a straightforward extension of a numerical scheme used to solve SCFT, the same theory that provides the equilibrium states bounding the dynamical trajectory; (iii) the method does not require theoretical input for the bounding states, but may also take advantage of information gleaned from experimental measurements; and (iv) transitions between metastable states may also be studied by explicitly specifying them rather than hunting for them in the dynamics.

In addition to establishing the method, our primary finding is that an elaborate theoretical framework is not needed to predict many of the qualitative features of order-order phase transitions. Nevertheless, along with the extensions discussed above, the highly successful machinery of polymer field theory may be incorporated into the MFGT approach to improve the method’s fidelity to experimental observations. Greater detail may be built into the model, for example, by appealing to the full conventional field theoretical description of diblock copolymers and studying the evolution of the chemical potential and pressure fields, similar to the external potential dynamics model of Fraaije and co-workers.15–18 Finally, the simplicity of the method presented here suggests that it could find application to other classes of nonequilibrium phenomena manifest by polymers and soft matter, in general.

In the preceding, all plots are constructed with the Matplotlib library.59 Simulation snapshots are prepared with ParaView.58 The smaller calculations with Nv = 643 grid elements are facilitated by the FFTW60 discrete Fourier transform library, while the larger Nv = 2563 calculations rely on the NVIDIA® cuFFT library. The estimates of the small-angle scattering utilize the KernSmooth61 R library. All other software were written by us in Python 2 and C.

This work was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at the Los Alamos National Laboratory, operated by Los Alamos National Security, LLC, under Contract No. DE-AC52-06NA25396. This research was supported by the Laboratory Directed Research and Development program under Project No. 20170001DR.

1.
For an interesting historical perspective, see:
I.
Prigogine
,
The End of Certainty
(
The Free Press
,
New York
,
1996
).
2.
I.
Prigogine
,
Non-Equilibrium Statistical Mechanics
(
Interscience Publishers
,
New York
,
1962
).
3.
P.
Attard
,
Non-Equilibrium Thermodynamics and Statistical Mechanics: Foundations and Applications
(
Oxford University Press
,
Oxford
,
2012
).
4.
P. C. W.
Davies
,
The Physics of Time Asymmetry
(
University of California Press
,
Berkeley
,
1977
).
5.
L. D.
Landau
and
E. M.
Lifshitz
,
Statistical Physics
(
Addison-Wesley Publishing, Inc.
,
Reading
,
1958
).
6.
P. M.
Chaikin
and
T. C.
Lubensky
,
Principles of Condensed Matter Physics
(
Cambridge University Press
,
Cambridge
,
2000
).
7.
M.
Huang
,
R. P.
Malhamé
, and
P. E.
Caines
, “
Large population stochastic dynamic games: Closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle
,”
Commun. Inf. Syst.
6
,
221
252
(
2006
).
8.
J.-M.
Lasry
and
P.-L.
Lions
, “
Mean field games
,”
Jpn. J. Math.
2
,
229
260
(
2007
).
9.
I.
Swiecicki
,
T.
Gobron
, and
D.
Ullmo
, “
Schrödinger approach to mean field games
,”
Phys. Rev. Lett.
116
,
128701
(
2016
).
10.
A.
Bensoussan
,
J.
Frehse
, and
P.
Yam
,
Mean Field Games and Mean Field Type Control Theory
(
Springer
,
New York
,
2013
).
11.
P. G.
de Gennes
,
Scaling Concepts in Polymer Physics
(
Cornell University Press
,
New York
,
1979
).
12.
G. H.
Fredrickson
,
The Equilibrium Theory of Inhomogeneous Polymers
(
Oxford University Press, Inc.
,
New York
,
2006
).
13.
N.
Ji
,
P.
Tang
,
F.
Qiu
, and
A.-C.
Shi
, “
Kinetic pathways of lamellae to gyroid transition in weakly segregated diblock copolymers
,”
Macromolecules
48
,
8681
8693
(
2015
).
14.
M. W.
Matsen
, “
Cylinder-gyroid epitaxial transitions in complex polymeric liquids
,”
Phys. Rev. Lett.
80
,
4470
4473
(
1998
).
15.
J. G. E. M.
Fraaije
, “
Dynamic density functional theory for microphase separation kinetics of block copolymer melts
,”
J. Chem. Phys.
99
,
9202
9212
(
1993
).
16.
N. M.
Maurits
and
J. G. E. M.
Fraaije
, “
Mesoscopic dynamics of copolymer melts: From density dynamics to external potential dynamics using nonlocal kinetic coupling
,”
J. Chem. Phys.
107
,
5879
5889
(
1997
).
17.
J. G. E. M.
Fraaije
,
B. A. C.
van Vlimmeren
,
N. M.
Maurits
,
M.
Postma
,
O. A.
Evers
,
C.
Hoffmann
, and
G.
Goldbeck-Wood
, “
The dynamic mean-field density functional method and its application to the mesoscopic dynamics of quenched block copolymer melts
,”
J. Chem. Phys.
106
,
4260
4269
(
1997
).
18.
J. G. E. M.
Fraaije
and
J. A.
Sevink
, “
Model for pattern formation in polymer surfactant nanodroplets
,”
Macromolecules
36
,
7891
7893
(
2003
).
19.
S.
Qi
and
Z.-G.
Wang
, “
Kinetic pathways of order-disorder and order-order transitions in weakly segregated microstructured systems
,”
Phys. Rev. Lett.
76
,
1679
1682
(
1996
).
20.
V.
Ganesan
and
G. H.
Fredrickson
, “
Field-theoretic polymer simulations
,”
Europhys. Lett.
55
,
814
820
(
2001
).
21.
G. F.
Fredrickson
, “
Dynamics and rheology of inhomogeneous polymeric fluids: A complex Langevin approach
,”
J. Chem. Phys.
117
,
6810
6820
(
2002
).
22.
G. H.
Fredrickson
,
V.
Ganesan
, and
F.
Drolet
, “
Field-theoretic computer simulation methods for polymers and complex fluids
,”
Macromolecules
35
,
16
39
(
2002
).
23.
K.
Kawasaki
and
K.
Sekimoto
, “
Dynamical theory of polymer melt morphology
,”
Physica A
143
,
349
413
(
1987
).
24.
G. H.
Fredrickson
and
K.
Binder
, “
Kinetics of metastable states in block copolymer melts
,”
J. Chem. Phys.
91
,
7265
7275
(
1989
).
25.
R.
Hasegawa
and
M.
Doi
, “
Adsorption dynamics. Extension of self-consistent field theory to dynamical problems
,”
Macromolecules
30
,
3086
3089
(
1997
).
26.
G. J. A.
Sevink
,
A. V.
Zvelindovsky
,
B. A. C.
van Vlimmeren
,
N. M.
Maurits
, and
J. G. E. M.
Fraaije
, “
Dynamics of surface directed mesophase formation in block copolymer melts
,”
J. Chem. Phys.
110
,
2250
2256
(
1999
).
27.
C.
Yueng
and
A. C.
Shi
, “
Formation of interfaces in incompatible polymer blends: A dynamical mean field study
,”
Macromolecules
32
,
3637
3642
(
1999
).
28.
E.
Reister
,
M.
Müller
, and
K.
Binder
, “
Spinodal decomposition in a binary polymer mixture: Dynamic self-consistent-field theory and Monte Carlo simulations
,”
Phys. Rev. E
64
,
041804
(
2001
).
29.
T.
Araki
and
H.
Tanaka
, “
Three-dimensional numerical simulations of viscoelastic phase separation: Morphological characteristics
,”
Macromolecules
34
,
1953
1963
(
2001
).
30.
M.
Nonomura
and
T.
Ohta
, “
Kinetics of morphological transitions between mesophases
,”
J. Phys.: Condens. Matter
13
,
9089
9112
(
2001
).
31.
K.
Yamada
,
M.
Nonomura
, and
T.
Ohta
, “
Kinetics of morphological transitions in microphase-separated diblock copolymers
,”
Macromolecules
37
,
5762
5777
(
2004
).
32.
D. Q.
Ly
,
T.
Honda
,
T.
Kawakatsu
, and
A. V.
Zvelindovsky
, “
Hexagonally perforated lamella-to-cylinder transition in a diblock copolymer thin film under and electric field
,”
Macromolecules
41
,
4501
4505
(
2008
).
33.
X.
Cheng
,
L.
Lin
,
W.
E
,
P.
Zhang
, and
A.-C.
Shi
, “
Nucleation of ordered phases in block copolymers
,”
Phys. Rev. Lett.
104
,
148301
(
2010
).
34.
H.
Risken
,
The Fokker-Planck Equation
(
Springer
,
New York
,
1996
).
35.
K. Ø.
Rasmussen
and
G.
Kalosakas
, “
Improved numerical algorithm for exploring block copolymer mesophases
,”
J. Polym. Sci., Part B: Polym. Phys.
40
,
1777
1783
(
2002
).
36.
M. W.
Matsen
and
F. S.
Bates
, “
Unifying weak- and strong-segregation block copolymer theories
,”
Macromolecules
29
,
1091
1098
(
1996
).
37.
M. F.
Schulz
,
F. S.
Bates
,
K.
Almdal
, and
K.
Mortensen
, “
Epitaxial relationship for hexagonal-to-cubic phase transition in block copolymer mixture
,”
Phys. Rev. Lett.
73
,
86
89
(
1994
).
38.
S.
Förster
,
A. K.
Khandpur
,
J.
Zhao
,
F. S.
Bates
,
I. W.
Hamley
,
A. J.
Ryan
, and
W.
Bras
, “
Complex phase behavior of polyisoprene-polystyrene diblock copolymers near the order-disorder transition
,”
Macromolecules
27
,
6922
6935
(
1994
).
39.
M. E.
Vigild
,
K.
Almdal
,
K.
Mortensen
,
I. W.
Hamley
,
J. P. A.
Fairclough
, and
A. J.
Ryan
, “
Transformations to and from the gyroid phase in a diblock copolymer
,”
Macromolecules
31
,
5702
5716
(
1998
).
40.
G.
Floudas
,
R.
Ulrich
, and
U.
Wiesner
, “
Microphase separation in poly(isoprene-b-ethylene oxide) diblock copolymer melts. I. Phase state and kinetics of the order-to-order transitions
,”
J. Chem. Phys.
110
,
652
663
(
1999
).
41.
G.
Floudas
,
R.
Ulrich
,
U.
Wiesner
, and
B.
Chu
, “
Nucleation and growth in order-to-order transitions of a block copolymer
,”
Europhys. Lett.
50
,
182
188
(
2000
).
42.
C.-Y.
Wang
and
T. P.
Lodge
, “
Unexpected intermediate state for the cylinder-to-gyroid transition in a block copolymer solution
,”
Macromol. Rapid Commun.
23
,
49
54
(
2002
).
43.
C.-Y.
Wang
and
T. P.
Lodge
, “
Kinetics and mechanisms for the cylinder-to-gyroid transition in a block copolymer solution
,”
Macromolecules
35
,
6997
7006
(
2002
).
44.
H.-W.
Park
,
J.
Jung
,
T.
Chang
,
K.
Matsunaga
, and
H.
Jinnai
, “
New epitaxial phase transition between DG and HEX in PS-b-PI
,”
J. Am. Chem. Soc.
131
,
46
47
(
2009
).
45.
D. A.
Hajduk
,
R.-M.
Ho
,
M. A.
Hillmyer
,
F. S.
Bates
, and
K.
Almdal
, “
Transition mechanisms for complex ordered phases in block copolymer melts
,”
J. Phys. Chem. B
102
,
1356
1363
(
1998
).
46.
D. A.
Hajduk
,
S. M.
Gruner
,
P.
Rangarajan
,
R. A.
Register
,
L. J.
Fetters
,
C.
Honeker
,
R. J.
Albalak
, and
E. L.
Thomas
, “
Observation of a reversible thermotropic order-order transition in a diblock copolymer
,”
Macromolecules
27
,
490
501
(
1994
).
47.
V.
Abetz
and
P. F. W.
Simon
, “
Phase behaviour and morphologies of block copolymers
,”
Adv. Polym. Sci.
189
,
125
212
(
2005
).
48.
C.
Wang
,
K.
Jiang
,
P.
Zhang
, and
A.-C.
Shi
, “
Origin of epitaxies between ordered phases of block copolymers
,”
Soft Matter
7
,
10552
10555
(
2011
).
49.
I. W.
Hamley
,
K. A.
Koppi
,
J. H.
Rosedale
,
F. S.
Bates
,
K.
Almdal
, and
K.
Mortensen
, “
Hexagonal mesophases between lamellae and cylinders in a diblock copolymer melt
,”
Macromolecules
26
,
5959
5970
(
1993
).
50.
I. W.
Hamley
,
V.
Castelletto
,
G.
Floudas
, and
F.
Schipper
, “
Templated crystallization from oriented gyroid and hexagonal melt phases in a diblock copolymer
,”
Macromolecules
35
,
8839
8845
(
2002
).
51.
C. Y.
Ryu
and
T. P.
Lodge
, “
Thermodynamic stability and anisotropic fluctuations in the cylinder-to-sphere transition of a block copolymer
,”
Macromolecules
32
,
7190
7201
(
1999
).
52.
K.
Kimishima
,
T.
Koga
, and
T.
Hashimoto
, “
Order-order phase transition between spherical and cylindrical microdomain structures of block copolymer. I. Mechanism of the transition
,”
Macromolecules
33
,
968
977
(
2000
).
53.
R.
Krishnamoorti
,
A. S.
Silva
,
M. A.
Modi
, and
B.
Hammouda
, “
Small-angle neutron scattering study of a cylinder-to-sphere order-order transition in block copolymers
,”
Macromolecules
33
,
3803
3809
(
2000
).
54.
See for example:
M. S.
Leifer
and
M. F.
Pusey
, “
Is a time symmetric interpretation of quantum theory possible without retrocausality?
,”
Proc. R. Soc. A
473
,
20160607
(
2017
).
55.
A. S.
Eddington
,
The Nature of the Physical World
(
Cambridge at the University Press
,
Cambridge
,
1929
).
56.
W. R.
Wharton
, “
Backward causation and the EPR paradox
,” e-print arXiv:quant-ph/9810060v1 (
1998
).
57.
I.
Prigogine
, “Time, Structure and Fluctuations,” Nobel Lecture,
1977
.
58.
J.
Ahrens
,
B.
Geveci
, and
C.
Law
, “
ParaView: An end-user tool for large data visualization
,” in
The Visualization Handbook
, edited by
C. D.
Hansen
and
C. R.
Johnson
(
Butterworth-Heinemann
,
Burlington
,
2005
), pp.
717
731
.
59.
J. D.
Hunter
, “
Matplotlib: A 2D graphics environment
,”
Comput. Sci. Eng.
9
,
90
95
(
2007
).
60.
M.
Frigo
and
S. G.
Johnson
, “
The design and implementation of FFTW3
,”
Proc. IEEE
93
,
216
231
(
2005
).
61.
M. P.
Wand
and
M. C.
Jones
,
Kernel Smoothing
(
Chapman and Hall
,
London
,
1995
).