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.

## I. INTRODUCTION

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 progress^{2} 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 approaches^{13–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. Matsen^{14} 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. Fraaije^{15} 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 Wang^{19} used TDGL to study the transitions involving the hexagonal and body-centered-cubic phases, again finding evidence of epitaxy. Ganesan and Fredrickson^{20} 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 *O*^{70} 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.

## II. THEORETICAL APPROACH

We apply the Hamilton-Jacobi-Bellman (HJB)/Fokker-Planck (FP) formulation of MFGT as constructed by Huang, Malhamé, and Caines^{7} 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 *dR*_{i} = $vi$(*t*)*dt* + *σd*$w$. Here, *R*_{i} 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=\u222b0TM2vi2\u2212F[\rho ]d\tau +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/(2*M*)(*∂u*/*∂x*)^{2} + (*σ*^{2}/2)∇^{2}*u* = *F* with $v$(*t*) = −1/*M∂*_{x}*u*. One works backward from the target condition $CT$ in this formulation. The density *ρ*, in turn, obeys the Fokker-Planck (FP) equation^{34} 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 Swiecicki^{9} 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 $\Phi =exp\u2212u/(M\sigma 2)$ and Γ = *ρ*/Φ leads to the following equations:

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\u2192$ 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:

Read the initial and final density fields, $\rho (x\u2192,0)$ and $\rho (x\u2192,T)$, respectively.

Set the initial guess for $\rho (x\u2192,t\u2260T)=\rho (x\u2192,0)$.

Set the initial guess for $F\rho (x\u2192,t)$, described below.

Set the initial guess $\Phi (x\u2192,T)=exp\u2212F\rho (x\u2192,T)M\sigma 2$.

For *n* steps:

Solve Eq. (1) backwards to find $\Phi (x\u2192,t)$ for 0 ≤

*t*<*T*.Solve Eq. (2) forwards to find $\Gamma (x\u2192,t)$ for 0 <

*t*≤*T*with $\Gamma (x\u2192,0)=\rho (x\u2192,0)\Phi (x\u2192,0)$.Update the estimate of $\rho (x\u2192,t)\u2192(1\u2212\u03f5)\rho (x\u2192,t)+\u03f5\Phi (x\u2192,t)\Gamma (x\u2192,t)$ for 0 <

*t*<*T*.Normalize $\rho (x\u2192,t)$ for 0 <

*t*<*T*.Update the estimate of $F(x\u2192,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 $A\u2261\u2211t,x\u219212\Delta \rho \Delta t2\u2212F$ 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.

## III. RESULTS AND DISCUSSIONS

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$ = 64^{3} 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$ = 256^{3} simulation box. Both simulation sizes are used to eliminate any concerns over finite-size effects. Figure 1 illustrates the $Nv$ = 256^{3} 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$ = 64^{3} simulation boxes; we always use *T* = 50 for the $Nv$ = 256^{3} simulations. In order to establish that the dynamics are dominated by the distribution of mass $\rho (x\u2192,t)$ given by the boundary states at *t* = 0 and *t* = *T*, we apply a potential of the form $F\rho (x\u2192,t)=k(\rho (x\u2192,t)\u2212\rho (x\u2192,T))2\u2009+\u2009\rho (x\u2192,T)$ with values of *k* = 0, 0.2, 0.4, 0.6, 0.8, and 1 for the $Nv$ = 64^{3} simulations; *k* is always set to zero for the $Nv$ = 256^{3} systems. The quadratic term simply serves as a probe potential designed to establish that the linear term in $\rho (x\u2192,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.

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 Matsen^{14} 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 experiments^{46} 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.

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\u2261\chi N\u2211iNv\rho i,j1\u2212\rho 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 Ji^{13} 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 $\delta \rho \u2261\u2211iNv\u2211jT\rho i,j\u2212\rho i,T\u2212j\u20322$, where $\rho i,T\u2212j\u2032$ 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$ = 64^{3}, *T* = 100 and $Nv$ = 256^{3}, *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)=|\rho ^(q)|2$, where $\rho ^q$ is the Fourier transform of the density distribution *ρ* and the average is over all scattering vectors $q\u2192$ 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 2*q*^{*}, 3*q*^{*}, 4*q*^{*}, 5*q*^{*}, and 6*q*^{*}. 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 3*q*^{*}. 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 $dR\u2192i/dt=h(x\u2192,t)+g(x\u2192,t)\Gamma (x\u2192,t)$. The first term on the right-hand side, $h(x\u2192,t)$, is a velocity or drift term. The second term on the right, $g(x\u2192,t)\Gamma (x\u2192,t)$, mimics the thermal noise in the system. Typically, $g(x\u2192,t)$ is set to a constant that reflects the temperature and viscosity of the system, while $\Gamma (x\u2192,t)$ is assumed to average to zero and to be correlated locally in time and space ($<\Gamma (x\u2192,t)\Gamma (x0\u2192,t0)>\u2009\u221d\u2009\delta (t\u2212t0)\delta (x\u2212x0)$). As discussed in detail in the work of Risken,^{34} moving from discrete particle trajectories to probability distributions $W(x\u2192,t)$, and thus from Langevin to FP equations, involves identifying the Kramers-Moyal expansion coefficients with functions of $h(x\u2192,t)$ and $g(x\u2192,t)$. For a single particle, one has the FP equation $\u2202W(x\u2192,t)/\u2202t=\u2212\u2202/\u2202x\u2192\u2009D(1)(x\u2192,t)+\u22022/\u2202x\u21922\u2009D(2)W(x\u2192,t)$, where $D(1)=h(x\u2192,t)+\u2202g(x\u2192,t)/\u2202x\u2192\u2009g(x\u2192,t)$ and $D(2)=g2(x\u2192,t)$. We have taken the steps of identifying $W(x\u2192,t)$ with the density $\rho (x\u2192,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 Sekimoto^{23} provided the foundation for the former, elaborating on expressions analogous to *D*^{(2)} to account for chain dynamics. Fredrickson^{21} 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 Lubensky^{6}). 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.

## IV. CONCLUSIONS

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.

## V. SOFTWARE NOTES

In the preceding, all plots are constructed with the Matplotlib library.^{59} Simulation snapshots are prepared with ParaView.^{58} The smaller calculations with $Nv$ = 64^{3} grid elements are facilitated by the FFTW^{60} discrete Fourier transform library, while the larger $Nv$ = 256^{3} calculations rely on the NVIDIA® cuFFT library. The estimates of the small-angle scattering utilize the KernSmooth^{61} R library. All other software were written by us in Python 2 and C.

## ACKNOWLEDGMENTS

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.