Evolution in off-critical diblock copolymer melts

We study the evolution of diblock copolymer melts in which one component has small volume fraction. In this case one observes phase morphologies which consist of small spheres of the minority component embedded in the other component. Based on the Ohta-Kawasaki free energy one can set up an evolution equation which has the interpretation of a gradient flow. We restrict this gradient flow to morphologies in which the minority phase consists of spheres and derive monopole approximations for different parameter regimes. We use these approximations for simulations of large particle systems.


Diblock copolymers.
A diblock copolymer molecule is a linear sub-chain of N A A-monomers grafted covalently to another sub-chain of N B B-monomers.Because of the repulsion between the unlike monomers, the different type subchains tend to segregate, but as they are chemically bonded in chain molecules, segregation of sub-chains cannot lead to a macroscopic phase separation.Only a local micro-phase separation occurs: micro-domains rich in A and B emerge.These micro-domains form morphology patterns on a larger scale.In this paper we are mainly interested in the case that the fraction of, say, A-monomers in a chain is small.Then we observe patterns which consist of many small separated spheres rich in A-monomers.1.2.The Ohta-Kawasaki free energy.The Ohta-Kawasaki [11] free energy of an incompressible diblock copolymer melt is a functional of the A-monomer density field.Let u(x) be the relative A-monomer number density at a point x in the sample.When there is high A-monomer concentration at x, u(x) is close to 1; when there is high concentration of B-monomers at x, u(x) is close to 0. A value of u(x) between 0 and 1 means that a mixture of A-and B-monomers occupies x.The free energy of the system is written as defined in X a = {u ∈ W 1,2 (D) : 1 |D| D u = ρ}.When the free energy I η is minimized, the three terms of the integrand (the energy density of the material) have different preferences.The first term likes large blocks of monomers, thereby reducing the total size of the interface between the two monomers.The function W in the second term is a double well potential with two global minima at 0 and 1, reflecting its preference for segregated monomers over mixtures.The third nonlocal term is minimized if u ≡ ρ.However, this configuration of u makes the second term large.The second best choice for the third term is to have u oscillate rapidly around ρ.This way, because of the compact nature of (−∆) −1/2 (when D is bounded), (−∆) −1/2 (u − ρ) becomes close to 0 and hence the third term becomes small.As these preferences compete and compromise, small blocks rich in A-monomers and B-monomers appear.This phenomenon is known as micro-phase separation.For a further analysis of the scaling of the energy we refer to [4].
1.3.The sharp interface limit.If η is small the interfacial regions become smaller and one can replace I η by its sharp interface limit.In [10] the Ohta-Kawasaki theory is formulated on a bounded domain as a singularly perturbed variational problem with a nonlocal term and the limiting sharp interface problem is identified.The latter is rigorously derived in [12] as a Γ-limit of the singularly perturbed variational problem.
In this paper we consider an idealized situation, in which D = R 3 , so we formulate the energy directly for this case.The corresponding energy is defined for all Ω ∈ R 3 with |Ω| = ρ and is given by where χ is the characteristic function of Ω and

1.4.
A time dependent model.We are interested in how an initial configuration of a diblock copolymer melt evolves towards a state of minimal energy.An appropriate evolutionary model set in whole space which reduces E and keeps the volume fraction of both phases conserved is the following Mullins-Sekerka type free boundary problem.In this model the normal velocity v of the interface where [∇u with no-flux boundary conditions for u at infinity.Aspects of a local well-posedness theory for the Mullins-Sekerka-model (that is the case γ = 0) can be found for example in [2], [3] and [6].
1.5.The gradient flow structure of the evolutionary model.The evolution defined by ( 2)-( 4) has an interpretation as a gradient flow on a Riemannian manifold, more precisely it is the gradient flow of the energy (1) with respect to the H −1 norm in the bulk.To reveal this structure consider the manifold M := {Ω ⊂ R 3 | |Ω| = ρ} of all subsets of R 3 with fixed volume.The tangent space T Ω M at an element Ω ∈ M is described by all kinematically admissible normal velocities of ∂Ω, that is, The metric tensor on the tangent space is given by Note that u α is well-defined up to an additive constant (see [5] for a rigorous set-up).
The gradient flow of the energy ( 1) is now the dynamical system where at each time the velocity is the element of the tangent space which points in the direction of steepest descent of the energy.In other words, v is such that for all ṽ ∈ T Ω(t) M. Choosing ṽ = v we immediately obtain the energy estimate associated with each gradient flow, which is for all T > 0.
It is well known that the first variation of the surface energy, the first summand of E, is the mean curvature.The second nonlocal part is given by R 3 |∇µ| 2 dx, where µ solves −∆µ = χ.Hence R 3 |∇µ| 2 dx = R 3 µχ dx = Ω µ dx and the first variation in direction ṽ ∈ T Ω M is given by ∂Ω µṽ dS.Thus (5) reads after an integration by parts in the metric tensor ∂Ω uṽ dS = ∂Ω (κ + γµ)ṽ dS for all ṽ ∈ T Ω(t) M. This is up to an irrelevant constant just (4).
The advantage of the gradient flow perspective is that it can be restricted in a natural way to a lower dimensional submanifold.We will make use of this idea in the next section.

2.1.
Off-critical mixtures and restriction to spherical particles.In the following we are interested in the regime where the fraction of A-monomers is much smaller than the one of B-monomers.In this case the A-phase consists of a set of many small disconnected approximately spherical particles.In view of this fact it seems natural to restrict the evolution (2)-( 4) to spherical particles by restricting the gradient flow to such morphologies.For the case γ = 0 it has been shown rigorously in [1] that such an approximation is valid.We also refer to [13], in which stable states of the energy are established which consist of a collection of spheres.
We define the submanifold N ⊂ M consisting of all sets Ω which are the union of disjoint balls Ω = ∪ i B R i (X i ), where the centers {X i } i and the radii {R i } i are variable.Hence N can be identified with an open subspace of the hypersurface {Y , where N is the number and i = 1, • • • , N an enumeration of the particles.The tangent space can be identified with the hyperplane such that V i describes the rate of change of the radius of particle i and ξ i the rate of change of its center.The restriction of the metric tensor can be expressed as where the function w α solves For the following it will be convenient to split the metric tensor into the radial and shift part respectively.We write w = u + φ where u and φ are harmonic in-and outside the particles and where [∇u It turns out to be notationally convenient to consider the normalized energy E = E surf + γE nl , where We obtain the differentials of the energies in the direction of a tangent vector Z = { Ṽi , ξi } i as We will see in Section 2.2 that both, g Y (Z, Z) and DE nl , Z , can be expressed explicitly in terms of {R i , X i } i .For the moment we just state once more the gradient flow equation.For any t ≥ 0 we chose Z(t) such that for all Z ∈ T Y N .We can expect that if initial data {R i (0), X i (0)} are given such that particles do not overlap, a smooth solution to (7) exists at least for short times.If a particle disappears the evolution is not smooth; we can, however, extend the solution continuously by just starting again with the new configuration.Our evolution cannot be further extended when two or more particles collide.We cannot completely exclude such an event a priori.However, as we will see below, one part of the nonlocal energy is given by i j =i , which again shows that this part of the energy prefers uniformly distributed particles.If the energy is initially bounded, it remains bounded for all times as long as a solution exists.Hence, centers X i and X j cannot come arbitrarily close unless the particles become very small at the same time.
In the next Section 2.2 we solve the gradient flow equation.The main idea is that we can solve the equations for the potentials by superposition of suitable monopoles.We will see that the leading order terms in the equations differ, depending on whether the domain covered by the particles is larger or smaller than the well-known screening length that describes the effective interaction range of particles.In Section 2.3 we state the approximate equations for the gradient flow, first in the simplest case where the screening length is much larger than the system size, and second in the case in which the system size is of the order of the screening length.Those equations are the starting point for our numerical simulations.In Section 2.4 we derive the corresponding mean-field models which arise if one passes to a description with densities and in Section 3 we identify stationary states.Finally, in Section 4 we present results of numerical simulations of the monopole approximation.

Solution by monopoles.
In this section we explicitly give the solution of the gradient flow equation based on the fact that we can solve the equations for the respective potentials in the metric tensor and for the nonlocal part of the energy explicitly by the superposition of monopoles.
The relevant length scales and parameters in our system are the following.We assume that initially we have a uniform distribution of particles contained in the box (0, L) 3 .This property will not be conserved by the evolution but within the typical time scale particles remain in a box of order of size L.
We denote by 1 d 3 the number density of particles, such that d is the typical distance of one particle to its nearest neighbor, and we call R the typical size of the radii of the particles.We are interested in the case that the A-phase has small volume fraction and hence we always assume that It is well-known and it will also become apparent in the computations below that the crucial intrinsic length scale in diffusional interactions is the so-called screening length which describes the effective range of particle interactions.The following heuristic analysis shows that the leading order dynamics significantly differ, depending on whether L ≪ L sc or whether L ∼ L sc .

The energy.
In a first step we express the energy in terms of {R i , X i } i , which is trivial for the surface energy but also simple for the nonlocal energy.Indeed, the solution of −∆µ = χ is given by µ = i µ i , where Thus we find, using the mean-value theorem for Hence, the differential of the energy We can now compare the order of size of the different terms.First notice that in order that surface and nonlocal energy balance, the parameter γ has to be of size of order R −3 which we assume from now on.Furthermore, in what follows, we approximate sums by the corresponding integral, i. e.
, and we find that j =i Equivalently, we can neglect the last two terms on the right hand side of (11) compared to the first two if L ≪ L sc .If, however, L ∼ L sc , all terms are of the same order of size.

The metric tensor.
In this section we compute the metric tensor for given Z = {V i , ξ i } i solving the equations for the potential w = u + φ respectively by superposition of single particle solutions.Recall that g , where u and φ are given by (6).
One easily verifies that in case of a single particle B R i (X i ) the solutions of problem (6) are given by Hence, for a system with N particles the potentials are given by We first observe that Now Hence, the mean value theorem for harmonic functions implies and we arrive at We compare the order of the first and second term in (13) respectively, as we did for the energy in Section 2.2.1.Keeping in mind that the number of particles is of order d 3 , we find that the first term is of size of order , where V R denotes the order of size of the velocities of the radii.Again approximating sums by integrals, the order of size of the second term is Hence, the first term dominates if and only if L 2 ≪ d 3 R , that is, again, if the box size is much smaller than the screening length.
Similarly as in (12) we compute Here and in the following we use the notation for mean values.It is easily calculated, for instance by using polar coordinates, that and thus The second term is computed with the help of Gauss' theorem and again the mean value theorem, It turns out to be of higher order than the diagonal term.Indeed, the latter is of order , where V X is the order of size of the midpoint velocities, and the order of the ratio of the former and this term can be estimated as before by Finally, for the mixed term of the metric tensor we find In the last equalities we used that • n dS = 0 and, as in the derivation of (13), that φ j is harmonic in B R i (x i ) for j = i.
To summarize we obtain and the diagonal terms dominate if L ≪ L sc .
2.3.Leading order approximations.We now set up an approximation of the gradient flow, keeping only the leading order terms in the energy and the metric tensor respectively.This will be the starting point for numerical simulations whose results we present in Section 4 below.As we have seen, there are two regimes of interest.The first one is the dilute case, where the system size is much smaller than the screening length.
2.3.1.Case I: L ≪ L sc .In case L ≪ L sc we obtain, in view of (10) and the subsequent discussion, that whereas the metric tensor can be approximated due to (15) by The gradient flow equation , where λ is a Lagrange multiplier to ensure i R 2 i V i = 0. Consequently, we obtain for the direction of steepest descent that where 2.3.2.Case II: L ∼ L sc .With the same procedure as above we infer from ( 11) and (15) that V i solves the linear system of equations where λ is such that i R 2 i V i = 0, and we get the nontrivial evolution equation for the centers.Notice that (20) implies that V X is of order L/d 3 .Using this in (19), we see that the last term on the right hand side is of order R 4 L 2 /d 6 = ε(L/L sc ) 2  and can be neglected compared to the other terms, which are of order 1.We finally obtain that the evolution of the radii is to leading order governed by 2.4.Mean-field models.In this section we derive mean-field models which arise if one passes from the discrete setting of finitely many particles to a description with densities.Our derivation is purely formal and starts from ( 17) and ( 20)-( 21) respectively.A rigorous derivation from the full model in the spirit of [9] seems feasible under some assumptions on the distribution of particles.This is, however, not within the scope of this paper.
For the derivation of mean-field models we need to go over to suitably rescaled variables.First we recall the relevant parameters and length scales in our system.We denote by N the total number of particles, by 1 d 3 the number density, and by R the typical particle radius.Our system size, that is the area covered with particles, is then L 3 ∼ N d 3 .The crucial intrinsic length scale is the screening length L sc ∼ (d 3 /R) 1/2 .Note that these quantities are in general time dependent.As introduced here they refer to the initial configuration and we derive the meanfield models for finite times in which these length scales remain of the same order.
Recall also that in order to balance surface and nonlocal energy our parameter γ must be of order R −3 .Hence, we set from now on γ := γR 3 .

2.4.1.
Case I: L ≪ L sc .We start with the simple case L ≪ L sc .In view of (17) we measure the time in which radii change in units of R 3 , that is, we introduce the new time scale τ = t/R 3 .Since ξ i = 0 in the dilute regime, it suffices to consider the distribution of particle radii.Let and define where ζ ∈ C ∞ 0 ((0, ∞)).Then, using (17), we find where λ = λR is such that r 3 ν t (dr) ≡ 3ρ 4πN =: ρ.Equation ( 23) just means that ν t satisfies in the sense of distributions.
For γ = 0 we recover the LSW-theory for coarsening of particles.For γ > 0 the additional term γ i R 5 i 5 in the energy prevents coarsening and hence the long-time behavior of (24) is different from the one for γ = 0. Indeed, while the energy (16) together with the volume constraint has no stationary point for γ = 0, we will see in Section 3 below, that (16) has a unique global minimizer for γ > 0 which is also the only stationary point with respect to an appropriate topology.

2.4.2.
Case II: L ∼ L sc .We assume now that L ∼ L sc = d 3 R 1/2 .In addition to r i and v i as in (22) we also introduce Now we define the joint distribution of particle centers and radii via We are going to derive that ν t satisfies the following system of equations in a distributional sense: where ε is the volume fraction (recall (8)), (notice that the integral is finite since ν τ has finite support in the x-variable) and ū = ū(x, τ ) satisfies Indeed, we first recall (25), to find Equation ( 20) gives where ū(x) (neglecting in the notation the dependence on t) is given by Then (21) reads where we used (27) in the last step as well as the fact RN = L sc , which is due to L = L sc .Inserting (32) into (31), we obtain If we take the Laplacian in the last equation, we obtain (28).Finally, (26) follows from ( 29), ( 30) and (32).
3. Stationary points of the energy.In this section we investigate the stationary states of the energies.Since we are interested in configurations with a large number of particles, we consider the continuous case.As we will argue later, in our setting stationary points only exist in the dilute case, so we first concentrate on this.For notational convenience we rename ρ and γ by ρ and γ again.Thus, let ν be a measure on (0, ∞) such that is fixed.The energy of ν is Proof.Let µ be a measure so that µ(dr) = r 3 ν(dr). Then Since the integrand is convex, by Jensen's inequality we have where r = 4π 3ρ r µ(dr) is the average of r under µ.Note that equality in (33) holds only when µ = 3ρ 4π δ r.In this case ν = 3ρ 4πr 3 δ r.So it suffices to minimize E among such measures and, since .
Next we investigate stationary points of E. The notion of stationary points depends on the variations and hence the topology one uses in the space of measures.The coarsest topology is the one induced by the Wasserstein distance, which metrizes the weak topology of measures.In one dimension the Wasserstein distance is easy to compute.Let be the radius distributions of the measures ν 1 and ν 2 and denote by r 1 (z), r 2 (z) their right-continuous inverses.The L p -Wasserstein distance between ν 1 and ν 2 is given by In the next lemma we show that the only stationary point of the energy with respect to the topology induced by the Wasserstein distance is the global minimizer identified in Lemma 3.1.Proof.Let ν be a stationary point of E and r = r(z) the right-continuous inverse of its distribution function as described above.We do not assume that ν is a finite measure, hence r might not have finite support.In terms of r the energy is given by The constraint 4π 3 r 3 ν(dr) = ρ corresponds to 4π 3 r 3 (z) dz = ρ.We consider now so called inner variations, that is, we consider r ε (z) := r(z + εη(z)), where Furthermore, the volume constraint Thus, in summary we find that for all ψ ∈ C ∞ 0 ((0, ∞)) and for some λ ∈ R. Let {z i }, i ∈ I ⊂ N, be the set of discontinuity points of r (which is an at most countable set, since r is decreasing).With this notation we deduce from (35) after an integration by parts that It follows that r has to be piecewise constant, say r(z) ≡ r i ∈ [z i−1 , z i ) (with the convention that z 0 = 0), and r i must be a zero of f (r) := r + γr 4 − 3λr 2 for some λ ∈ R. Furthermore, we find for all i ∈ I.A simple calculation shows that this can only be true if I = {1} and changes the configuration of stationary points.However, in our case there are no stationary states.This is due to the somewhat artificial setting, where we set particles in a bounded domain but solve the equations for the corresponding potentials in full space.From equation (20) for the evolution of particle centers we easily see that particles located at the boundary of our "cloud" of particles will be driven further outwards.
If instead we would confine particles and their corresponding potentials to a finite domain and impose periodic or Neumann boundary conditions we would obtain similar formulas for our approximate gradient flow equations and the mean-field models, but 1  |x| would be replaced by the corresponding Green's function.Then stationary states would be characterized by K ≡ const, which would imply a regular spacing of the particles.Once K is constant, the term r 3 Kν(dx dr) is also constant due to the volume constraint.The remaining discussion of stationary points is then analogous to the dilute case.
From the discussion in the mean-field theory, Section 2.4.2, we know that the particle centers move on a slower time scale d 3 than the particle radii, which change on a time scale of order R 3 .Consequently, in our simulations the movement of particles is almost negligible compared to the evolution of the radii.
Remark 3.4.(The discrete setting) Naturally, the discrete case is somewhat different from the continuous one, since the number of particles is integer and we cannot reach any arbitrary fraction of particles.As argued in Remark 3.3 above, we can concentrate on the dilute case (17).
Let us first consider the situation from the point of view of the evolution.Our initial configuration consists of N particles with total volume ρ 3 .Since during our evolution particles can only vanish, we consider stationary states with N 1 ≤ N particles.For any given N 1 there is a family of stationary states.First there is the trivial one, that is R i ≡ r for all i and r = 3ρ 4πN 1 .However, due to the discreteness of the problem there are now also stationary states with two different sizes of particles.For that let σ := N 2 N 1 ∈ (0, 1], where N 2 will be the number of particles with one radius r 1 , whereas N 1 − N 2 will be the number of particles with another radius r 2 .In view of (17) stationary states of the desired form are given by radii which are zeros of the function f (r) := 1 + γr 3 − λr for some λ ∈ R.There can be at most two different radii for which this is satisfied.We denote them by r 1 and r 2 .To say that r 1 and r 2 are zeros of the function f for some λ is equivalent to requiring that 1 r 1 + γr 2 1 = 1 r 2 + γr 2 2 holds true.In order to satisfy the volume constraint, we must have σr 3  1 + (1 − σ)r 3 2 = 3ρ 4πN 1 =: β.Using the last equation to express r 2 by r 1 we find that r 1−σ 2/3 .Since lim r→0 g(r) = +∞ and lim r→σ −1/3 g(r) = −∞ there is at least one solution.Thus we have a whole set of stationary states for given N 1 ≤ N which consist of two different sizes of particles.Due to the convexity of the energy they have higher energy than the stationary state where only one size is present and the size is given by r = 3ρ 4πN 1 .This is true for fixed N 1 .Of course, depending on the choice of parameters, it might be that a stationary state with two different sizes for one N 1 has lower energy than the trivial stationary state for a different given number of particles.
4. Simulations.For numerical simulation it is convenient to re-formulate ( 17) and ( 20)-(21) in terms of particle volumes and their velocities instead of radii.This has the advantage that the constraint of volume conservation becomes linear.For that let W i := R 3 i and ω i := d dt W i .The equations then read with λ as in (18) in the dilute case, and in case L ∼ L sc .1. Results of simulated particle systems all within the same parameter region, but different radii distribution.
To simulate the behavior of our gradient flow model with these equations we set up a simple numerical scheme consisting of an Euler predictor and a trapezoidal rule corrector step.The latter is also used as error indicator and to determine -beside the influence of vanishing particles -a suitable step size.Particles are removed when they become too small, i.e. if W i < δ W i /N where N denotes the current number of particles and δ is a user-defined tolerance.Under no circumstances we allow particle volumes to become negative.
At each time step the velocities of particle midpoints and volumes are computed via (36) and (37), respectively.While this is not much work in the first case, we have to solve a symmetric, indefinite, and dense linear system of size N + 1 in the second case.Due to its potential structure, however, the product of the system matrix with an arbitrary vector can be computed in almost linear time by means of the Fast Multipole Method [7,8].This and the approximate weak diagonal dominance of the matrix enable us to efficiently apply a Krylov subspace method to solve (37) for large systems.
We start our simulations with a number of particles distributed regularly on a mesh in some box (0, L) 3 .From the mesh size and (9) we obtain the value of the initial mean or typical radius in order that L = L sc is approximately fulfilled.We simulate the dilute and the full problem with the same radius distribution to be able to compare the results.4.1.Case I: L ≪ L sc .For any distribution of radii which we consider in our simulations, uniformly or normally distributed in some interval about the initial mean radius or special as stated below, the system reaches a stationary state.As  expected from Remark 3.4 this final state is characterized by either all particles having (almost) the same size or one of two radii with the condition stated there.The latter case, however, only appears when we already start the simulation with such a particle configuration; in general we arrive at a stationary state concentrated in one size.
However, the computed stationary state is not the global minimizer derived in Section 3, which, of course, might be infeasible due to the volume constraint and integer number of particles.In general one cannot even say that it is close to or approaches it when the number of particles grows.In fact, the computed state and its difference from the predicted minimizer of the continuous case in terms of the final radius depend on the actual distribution of the radii, as can be seen in Figure 1 and Table 1.

4.2.
Case II: L ∼ L sc .As expected, the evolution of the full problem does not reach a stationary state, i. e. a state in which all velocities are identically zero, since particles drift out of the box.But the following two observations can be made, as can be seen for instance in Figure 2. First of all, even though drifting can occur, its influence on energy minimizing is distinguishable from that of radii evolution and -as already stated in Remark 3.3 -can be neglected in case the value of L sc is not much smaller than the box size.Secondly and more important, one observes that the total energy in principle tends towards a stationary value, as does the mean particle radius.However, the final mean radius is in general not equal to that from the same initial configuration evolving under the dilute equations.Another difference is that in the full problem the particles do not tend to have equal size.
Instead, after a phase of vanishing and approaching it, all remaining particles have a radius in some interval about the final mean radius which is smaller than the initial radii interval, but whose size again seems to depend on the actual initial distribution.

4γ 1/ 3 and r 2
= 0.The point z 1 , where r jumps, is determined by the volume constraint and hence given by z 1 = 3ρ 4π r 3 * .Remark 3.3.(The inhomogeneous case) It is in principle a very interesting question, how in the case L ∼ L sc the additional inhomogeneous term in the energy,
energy diagonal term nonlocal energy off-diagonal term

Figure 2 .
Figure2.Evolution of mean radius, particle radii and energy in the full problem; compare with Figure1where the same system in the dilute case is shown.