Multiscaling for Classical Nanosystems: Derivation of Smoluchowski and Fokker-Planck Equations
aa r X i v : . [ c ond - m a t . m e s - h a ll ] J a n Multiscaling for Classical Nanosystems: Derivationof Smoluchowski & Fokker-Planck Equations
S. Pankavich ∗ , Z. Shreif † , P. Ortoleva † July 2, 2018
Abstract
Using multiscale analysis and methods of statistical physics, we show thata solution to the N -atom Liouville Equation can be decomposed via an expan-sion in terms of a smallness parameter ǫ , wherein the long scale time behaviordepends upon a reduced probability density that is a function of slow-evolvingorder parameters. This reduced probability density is shown to satisfy theSmoluchowski equation up to O ( ǫ ) for a given range of initial conditions. Fur-thermore, under the additional assumption that the nanoparticle momentumevolves on a slow time scale, we show that this reduced probability densitysatisfies a Fokker-Planck equation up to O ( ǫ ). This approach has applicationsto a broad range of problems in the nanosciences. Keywords: nanosystems, all-atom multiscale analysis (AMA), Gibbs Hypothesis,Smoluchowski equations, Fokker-Planck equations
Nanosystems are currently of great interest in the fundamental and applied life sci-ences. A major unresolved challenge is to develop a predictive approach to thesesystems that capture the inter-communication among the processes operating on dif-fering scales in space and time. The premise of the present work is that one canintroduce order parameters (slowly varying quantities that capture the essence oflarge-scale bionanosystem phenomena) and then, using Newton’s equations for the N -atom system, derive equations for stochastic order parameter dynamics.Examples of bionanosystems abound in nature and medicine. Viruses are supra-million atom entities with complex structural and functional characteristics, includingdramatic transitions, interactions with host cells, and self-assembly of subunits. Ri-bosomes are of size and complexity similar to viruses, and mediate an important ∗ Department of Mathematics, Indiana University, Bloomington, IN 47405, [email protected] † Center for Cell and Virus Theory; Department of Chemistry, Indiana University, Bloomington,IN 47405, [email protected] • an underlying all-atom description to evaluate the interaction of bionanosys-tems with selected molecules, membranes, or other features in their backgroundmicroenvironment • a model that does not require recalibration with each new application • an approach that builds in the detailed physical molecular laws and the predic-tive power following from them • an approach that is computationally feasible.Considering this list, we suggest that a multiscale analysis of the equations of N -atom physics will fulfill each of these requirements. Molecular Dynamics (MD)is a current state-of-the-art software package that efficiently performs simulations ofNewton’s equations for each of N atoms in a system of interest. An efficient MDcode, NAMD, has previously been used to simulate a whole virus using a 1024 CPUsupercomputer, but the process proceeds at a rate of about 1 nanosecond of simulatedtime per day. The typical timescale for a viral structural transition is on the order of amillisecond or greater. Thus, the aforementioned MD code and hardware would take3 ,
000 years or longer to attain meaningful results. As bionanosystems evolve due tothe cross-talk between processes which take place on many scales in both space andtime, a computational algorithm based on a multiscale approach seems like a naturalchoice.The use of multiscale techniques in statistical mechanics beginning with the Liou-ville equation has a long history (see [10], [6], [7], and [11], and more recently [5], [2],and [3]). In the present work, we demonstrate several new elements of the analysis.First, in our approach, the nanoparticle’s internal atomic state, as well as that of themicroenvironment, are maintained allowing for a more natural, symmetric starting2escription. Additionally, we utilize a version of the Gibbs postulated equivalenceof ensemble and long-time averages, following classical results within ergodic theory.A precise representation is obtained for the momentum factor in the normalizationconstant for the lowest-order N -atom probability density in a perturbation expansionof a solution to the Liouville equation. As a result, Fokker-Planck and Smoluchowskiequations are derived which describe the stochastic dynamics of these slow-evolvingorder parameters. These results can then be utilized in the production of an efficientsoftware module that can model nanoparticle behavior over long time scales, therebycapturing the necessary structural dynamics of a virus.In recent investigations (see [2] and [3]), the reduced probability density W wasshown to obey an unconserved equation of Fokker-Planck type up to O ( ǫ ). Thederivation of this equation was inconsistent with the mathematical framework of dif-ferential equations as the thermal average and the derivatives with respect to orderparameters ∂∂ Φ and ∂∂ Π do not commute. In the work that follows, we eliminateambiguities regarding the permutation of the thermal average and these derivatives.Additionally, the lowest order distribution was previously taken to be independentof the conjugate momentum Π. This is done in error, causing the lowest order de-pendence on the slow variable Π to be lost and propagating this throughout themultiscale analysis. In Section 3, we rigorously correct these mistakes and establishmany of the ideas of [3] on a more precise footing by showing that the correctionto the reduced probability density f W indeed satisfies a Fokker-Planck equation inconservative form. Prior to this, we show in Section 2 that if the momentum is anatomically varying quantity, rather than a slowly varying order parameter, then f W directly satisfies the Smoluchowski diffusion equation up to O ( ǫ ). In both sections,our derivations occur from the starting point of the general kinetic equation so thatthe resulting coarse-grained equations do not arrive from solubility conditions, butfrom a rational expansion of the Liouville equation. A central goal of multiscale analysis is to rigorously derive coarse-grained equationsstarting from a more fundamental, final scale theory. The Liouville equation hasbeen a common starting point. The challenge is that while the Liouville equationpreserves probability by construction, it is not guaranteed that a given truncationof the equation will be conserving. A re-examination of multiscale analysis for theLiouville equation is now carried out to identify potential difficulties of this type thatmay arise, and to set forth techniques to resolve them. In this section we resolveprobability conservation violations when the momentum is not a slow variable.Consider the Liouville equation in a multiscale framework wherein order param-eters are introduced. We consider an N -atom system consisting of a nanoparticle of3 ∗ atoms and a host medium of N − N ∗ atoms. For each atom i = 1 , ..., N , we write p i , r i ∈ R , and m i > i respectively.In addition, we use the notation Γ = { r , p , ..., r N , p N } . For each i = 1 , ..., N , definethe indicator function θ i = (cid:26) , if atom i is in the nanoparticle0 , otherwise.For the nanoparticle, we define its total mass m = N X i =1 m i θ i , the center of mass R = N X i =1 m i m r i θ i , (1)and the total momentum P = N X i =1 p i θ i . (2)To begin the multiscale analysis, we first introduce a dimensionless scaling param-eter ǫ in the mass terms by writing ǫ = b mm , (3)where b m is the mass of a typical atom. In the case that all atoms in the nanoparticlehave the same mass, m i = ˆ m for all i = 1 , ..., N ∗ , it follows that b mm = 1 N ∗ . Hence, ǫ ≈ ( N ∗ ) − . In this section, we make the following assumptions:1. The total nanoparticle momentum does not evolve slowly - P is O ( ǫ ).2. The net force on the nanoparticle is not decreased due to cancelation of atomiccontributions - f is O ( ǫ ).3. Large migration distances are not a consideration - R is O ( ǫ ).As a result, order parameters are O ( ǫ ) and need not be scaled in ǫ , even thoughNewton’s equations show that they evolve slowly as dRdt = Pm = ǫ P ˆ m = O ( ǫ ). Wenote that other scalings would be appropriate to capture different behavioral regimes.Let us assume ρ satisfies the Liouville Equation ∂ρ∂t = − N X i =1 (cid:20) p i m i · ∂∂r i + F i · ∂∂p i (cid:21) ρ ≡ L ρ (4)4here we define F i to be the force on atom i and t to be time. In addition, weassume throughout that ρ decays at infinity (a standard assumption for a probabilitydensity) so that boundary terms do not appear in the calculations from integrationby parts. Denote the collection of all atomic positions by Γ r = { r , ..., r N } . Giventhe probability density, ρ (Γ , t ), we define f W ( R, t ) = Z ∆(Γ ∗ r , R ) ρ (Γ ∗ , t ) d Γ ∗ (5)where ∆(Γ ∗ r , R ) = δ ( R − R ∗ ) , (6) R is the center of mass order parameter, and R ∗ = N X i =1 m i m r ∗ i θ i is the Γ ∗ r -dependent value of R . Then, using the dependence of f W on ρ , a solutionto the Liouville equation, we may show that f W must satisfy a conserved equation.Since ρ satisfies (4), we find ∂ f W∂t = Z ∆ ∂ρ∂t (Γ ∗ , t ) d Γ ∗ = − Z ∆ N X i =1 p ∗ i m i · ∂ρ∂r ∗ i + F i · ∂ρ∂p ∗ i ! d Γ ∗ = Z ρ N X i =1 p ∗ i m i · ∂ ∆ ∂r ∗ i + F i · ∂ ∆ ∂p ∗ i ! d Γ ∗ = − Z ρ N X i =1 p ∗ i m i · ∂R ∗ ∂r ∗ i + F i · ∂R ∗ ∂p ∗ i ! ∂ ∆ ∂R d Γ ∗ = ∂∂R (cid:18)Z ρ ∆ L R ∗ d Γ ∗ (cid:19) = − ∂∂R (cid:18)Z ρ ∆ P ∗ m d Γ ∗ (cid:19) Thus, the reduced probability density, f W , satisfies ∂ f W∂t = − ǫ ∂∂R (cid:18)Z ρ ∆ P ∗ b m d Γ ∗ (cid:19) . (7)Next, we attempt to determine ρ up to O ( ǫ ). The N-atom probability density, ρ (Γ , t ), is then assumed to be expressed as a function of an additional argument,Υ(Γ , t, · ) in such a way that when the last argument is evaluated at R , ρ is obtained,i.e ρ (Γ , t ) = Υ(Γ , t, R ). Instead of labeling this new function, we will just extend our5revious notation and refer to it as ρ (Γ , t, R ). This displays the dependence of theprobability density on multiple scales of motion. Hence, ρ depends on the all-atomdescriptive variables Γ, as well as on R defined by (1), the latter an expression of thefact that ρ has indirect dependence on the all-atom state through order parametersand thus depends on the all-atom state in several, distinct ways.We apply the Liouville operator to ρ (Γ , t, R ) and invoke the chain rule to find ∂ρ∂t = −L ρ − N X i =1 p i m i · dRdr i ∂ρ∂R . Using (1) this becomes ∂ρ∂t = −L ρ − Pm · ∂ρ∂R . (8)Here, we are writing L instead of L because these derivatives are taken at constantvalues of R . By introducing (3) into (8), the Liouville equation (4) transforms into amultiscale equation (see [2] and [3] for more details) as ∂ρ∂t = ( L + ǫ L ) ρ (9)where L = − N X i =1 (cid:20) p i m i · ∂∂r i + F i · ∂∂p i (cid:21) (10)and L = − P b m · ∂∂R . (11)Again, it must be noted that L and L , while seemingly exact in definition, differbecause the differentiation in L is performed at constant values of order parameters R . Additionally, the differentiation in L is performed at fixed values of Γ. Furtherdetails regarding the all-atom, multiscale analysis (AMA) for the Liouville equationcan be found in [2], [3], [5], and [9]. The operators (10) and (11) differ from that of thepreviously mentioned papers [2] and [3] since in that work the conjugate momentumΠ is treated as an order parameter. In the work that follows in this section, we treatthis momentum term as a micro-variable instead.Assuming the net force on the nanoparticle does not experience cancelation dueto fluctuating terms, it can be written in terms of the individual atomic forces as f = N X i =1 F i θ i . (12)We let V (Γ r ) be the N -atom potential so that ∂V∂r i = − F i for every i = 1 , ..., N . Next,we assume that ρ may be expressed as a power series in ǫ : ρ = ∞ X n =0 ρ n ǫ n (13)6 set of time variables, defined via t n = ǫ n t , is introduced to capture effects ofprocesses occurring on the various timescales. The chain rule implies ∂∂t = ∞ X n =0 ǫ n ∂∂t n . (14)We then expand (9) using (13) and (14) and separate ǫ scales. Define for n ∈ { } ∪ N ,Λ n = ∂∂t n − L n (15)where we take L n = 0 for n >
1. The expansion yields the equationsΛ ρ = 0 , and for n ∈ N , Λ ρ n = − n X i =1 Λ i ρ n − i . Assuming the statistical state of the system has quasi-equilibrium character, thelowest order distribution ρ is taken to be independent of t . Thus, to lowest orderwe find L ρ = 0 . (16)This implies ρ is a function of the conserved dynamical variables, notably the totalenergy H , as well as of R . The latter occurs because the derivatives ∂∂r i and ∂∂p i in L are to be taken at constant R , and thus L R = 0. Then, we can define H = N X i =1 (cid:18) p i m i + V (Γ r ) (cid:19) and notice that L H = 0.Using the entropy maximization principle, one arrives at the nanocanonical solu-tion to (16) from [2]: ρ = e − βH W ( R, t ) Q ≡ ˆ ρW, (17)where Q ( β, R ) = Z ∆(Γ ∗ r , R ) e − βH ∗ d Γ ∗ . (18)Here ∆ is defined as in (6) and H ∗ is the Γ ∗ -dependent value of H . For conveniencewe write t for the collection of slow time variables t = { t , t , ... } . As is standard inmultiscale theory, determination of W is delayed until higher orders in the analysis.With this, ρ is seen to factorize into the conditional probability ˆ ρ (i.e. for Γ given7 ), multiplied by the reduced probability W for the slowly evolving state of the orderparameter R . We define the thermal average of a given dynamical variable, A (Γ) by A th ≡ Z ˆ ρ ∆ A (Γ ∗ ) d Γ ∗ . (19)Now, we will assume that the nanocanonical ensemble obeys the Gibbs hypothesizedequivalence between the long-time and ensemble averages. More specifically, we utilizea classical theorem of Birkhoff ([4] can provide more detail) which states that thethermal average of a dynamical variable A (Γ) and its long-time average are equal.Using classical semigroup methods from applied partial differential equations (see[8] for more detail), one may show that the linear operator L is the infinitesimalgenerator of a strongly continuous semigroup on the function space L (Γ). Thissemigroup is then well-defined and denoted by e L t . Hence, in the analysis thatfollows, we will rely extensively on the property :lim t →∞ t Z − t e −L s Ads = A th (20)for all dynamical variables A (Γ). Thus, the long-time average or time evolution of avariable does not affect the value of its thermal average as defined in (19). The survey[4] or the classic article [1] can provide more background information and detail froman ergodic theory perspective.To O ( ǫ ) one finds Λ ρ = − Λ ρ . (21)Using the previously constructed semigroup e L t , equation (21) admits the solution ρ = e L t A − Z t e L ( t − t ′ ) Λ ρ dt ′ = e L t A − Z t e L ( t − t ′ ) (cid:20) ˆ ρ ∂W∂t + P ˆ m ˆ ρ · ∂W∂R + P ˆ m · ∂ ˆ ρ∂R W (cid:21) dt ′ The first order initial condition A is, for now, undetermined and has the dependence A (Γ , R, t ). As a consequence of the cross-level communication inherent to multiscaleanalysis, the behavior of ρ at large t provides information about the t -dependenceof W , while the analysis of (7) provides a necessary condition on A that ensures theequation determining f W is closed. Letting s = t ′ − t , one obtains ρ = e L t A − t ˆ ρ ∂W∂t + ˆ ρ (cid:20) βf th W − ∂W∂R (cid:21) · Z − t e −L s P b m ds. (22)In this equation, βf th = ∂∂R (ln Q ), so that f th is the force averaged via the nanocanon-ical ensemble. This term occurs because of the dependence of ˆ ρ on R and we will8erify the expression for ∂ ˆ ρ∂R in the Appendix. Thus, using the Gibbs Hypothesis(20), we find f th = lim t →∞ t Z − t e −L s f ds. (23)Next, we remove secular behavior from the t -dependence in ρ , thereby imposingthe additional condition that ρ remains bounded as t → ∞ . Using (22), it can beseen that if ρ grows in t , it must do so at least linearly. Hence, we may ensure that ρ does not grow in t by requiring that lim t →∞ t ρ = 0. We then divide by t , take thelimit as t → ∞ in equation (22), and use (20). Notice that ( P ) th = 0, as it involvesterms of the form Z p i exp (cid:18) p i m i (cid:19) dp i . Assuming the first order initial data is takenin the nullspace of L , that is L A = 0, use of (20) yieldslim t →∞ ρ t = − ∂W∂t . Hence, we find ∂W∂t = 0 (24)and the reduced probability density W is independent of t . Note that this propertyfollows regardless of the choice of A in the nullspace of L . Using this in (22), wefind ρ = A + ˆ ρ (cid:20) βf th W − ∂W∂R (cid:21) Z − t e −L s P b m ds (25)concluding the O ( ǫ ) analysis, although A and W are not yet determined.At this point, one would expect to conduct a O ( ǫ ) analysis of the problem anddetermine an equation for ∂W∂t . However, this is unnecessary as the correction to thereduced probability density depends only on ρ and ρ up to O ( ǫ ). Instead, definefor all n = 0 , , , ... f W n ( R, t ) = Z ∆(Γ ∗ r , R ) ρ n (Γ ∗ , t ) d Γ ∗ (26)so that, using (5) and (13), we may write f W = ∞ X n =0 ǫ n f W n . (27)Hence, we expand f W and ρ in powers of ǫ as in (13) and (27). Using (17) and (19),9he lowest order correction, f W , can be calculated as f W = Z ∆(Γ ∗ r , R ) ρ (Γ ∗ , t ) d Γ ∗ = Z ∆ ρ (Γ ∗ , t, R ∗ ) d Γ ∗ = Z ∆ ˆ ρ W ( R ∗ , t ) d Γ ∗ = W ( R, t ) . For ǫ →
0, one may see that f W → W . Hence, as the long time scales tend to zero,the correction tends to the reduced probability density. The O ( ǫ ) correction can bedetermined using (20) and (25), so that f W = Z ∆ ρ (Γ ∗ , t ) d Γ ∗ = Z ∆ (cid:20) A − ˆ ρ (cid:18) ∂W∂R − βf th W (cid:19) · Z − t dse −L s P ∗ b m (cid:21) d Γ ∗ = Z ∆ A (Γ ∗ , R ∗ , t ) d Γ ∗ − Z ∆ ˆ ρ (cid:18) ∂W∂R − βf th W (cid:19) · Z − t e −L s P ∗ b m dsd Γ ∗ (cid:12)(cid:12)(cid:12)(cid:12) R = R ∗ = Z ∆ A (Γ ∗ , R ∗ , t ) d Γ ∗ + (cid:18)Z ∆ ˆ ρ Z − t e −L s P ∗ b m dsd Γ ∗ (cid:19) · (cid:18) − ∂W∂R ( R, t ) + βf th W ( R, t ) (cid:19) = Z ∆ A (Γ ∗ , R ∗ , t ) d Γ ∗ + "Z − t (cid:18) e −L s P b m (cid:19) th ds · (cid:18) − ∂W∂R + βf th W (cid:19) = Z ∆ A (Γ ∗ , R ∗ , t ) d Γ ∗ . Now, we may write ρ (Γ , t ) in terms of its expansion up to O ( ǫ ). Using (17) and1025) in the right side of (7), we find Z ∆ P ∗ b m ρ (Γ ∗ , t ) d Γ ∗ = Z ∆ P ∗ b m ρ (Γ ∗ , t, R ∗ ) d Γ ∗ = Z ∆ P ∗ b m ( ρ + ǫρ ) d Γ ∗ = Z ∆ ˆ ρ P ∗ b m W ( R ∗ , t ) d Γ ∗ + ǫ Z ∆ P ∗ b m A (Γ ∗ , R ∗ , t ) d Γ ∗ + ǫ (cid:20)Z ∆ ˆ ρ P ∗ b m Z − t e −L s P ∗ b m ds (cid:18) βf th W ( R ∗ , t ) − ∂W∂R ( R ∗ , t ) (cid:19) d Γ ∗ (cid:21) = (cid:18) P b m (cid:19) th W ( R, t ) + ǫ Z ∆ P ∗ b m A (Γ ∗ , R ∗ , t ) d Γ ∗ + ǫ (cid:18) P ∗ b m · Z − t e −L s P ∗ b m ds (cid:19) th (cid:18) βf th W ( R, t ) − ∂W∂R ( R, t ) (cid:19) = ǫ Z d Γ ∗ ∆ P ∗ b m A (Γ ∗ , R ∗ , t ) + ǫ γ ˆ m (cid:18) βf th W ( R, t ) − ∂W∂R ( R, t ) (cid:19) where the diffusion coefficient is γ = Z − t ( P (0) · P ( s )) th ds (28)and we use the notation P ( s ) = e −L s P . Thus, (7) becomes ∂ f W∂t = − ǫ ∂∂R · (cid:18)Z ∆ P ∗ b m A (Γ ∗ , R ∗ , t ) d Γ ∗ (cid:19) + ǫ ∂∂R · (cid:20) γ ˆ m (cid:18) ∂W∂R − βf th W (cid:19)(cid:21) . Using the expressions for f W and f W , we can expand the reduced probability densityas f W = f W + ǫ f W . Then, isolating f W k terms and imposing the condition that A must stay bounded for large R , we find that this equation is closed only if A = 0.Thus, up to O ( ǫ ), the conserved equation (7) becomes the Smoluchowski equation: ∂ f W∂t = ǫ ∂∂R · (cid:20) γ ˆ m (cid:18) ∂∂R − βf th (cid:19) f W (cid:21) . (29)Hence, in the case of fast-evolving momentum, due to P and f being O ( ǫ ), the re-sulting behavior of the reduced probability density is governed by the Smoluchowskiequation (29). In the next section, we alter these assumptions on the behavior ofnanoparticle momentum and determine the corresponding changes in the structureof the equation for f W . In this section, we again use multiscale techniques to show that under similar circum-stances, the correction to the reduced probability density, f W , satisfies a Fokker-Planck11quation. In this situation, the momentum is not considered an atomistic variable,but instead as an order parameter. Hence, we define Γ , m, R , and P as before, butreformulate the problem to allow for the differing behavior of this slowly-evolvingquantity.To begin the multiscale analysis, we again introduce a dimensionless scaling pa-rameter ǫ in the mass terms by writing ǫ = b mm , (30)where b m is the mass of a typical atom. In this section, however, the assumptionson the system of interest change. We are now interested in significant migration dis-tances on the order of the nanoparticle diameter, which we take to be O ( ǫ − ), andhence scale R to be O ( ǫ − ). Additionally, under the assumption that the systemis near equilibrium, the nanoparticle kinetic energy, P m is O ( ǫ ). Using the massratio scaling, this implies that P = O ( ǫ − ), as well. Finally, we assume that the netforce on the nanoparticle is reduced due to cancelation of atomic contributions, thuscausing the momentum to evolve slowly. Hence, f is assumed to be O ( ǫ ). A moredetailed description of these assumptions can be found in [3].Under these considerations, define the scaled order parameters Φ and Π byΦ = ǫ R (31)and Π = ǫ P (32)respectively. The scaled net force f can then be written in terms of the individualatomic forces as f = ǫ − N X i =1 F i θ i , (33)and we let V (Γ r ) be the N -atom potential so that ∂V∂r i = − F i for every i = 1 , ..., N .Let us assume ρ satisfies the Liouville Equation (4) where we again consider F i to be the force on atom i and t to be time. In addition, we denote the collection ofall atomic positions by Γ r = { r , ..., r N } . Given the probability density, ρ (Γ , t ), wedefine f W (Φ , Π , t ) = Z ∆(Γ ∗ , Φ , Π) ρ (Γ ∗ , t ) d Γ ∗ . (34)where ∆(Γ ∗ , Φ , Π) = δ (Φ − Φ ∗ ) δ (Π − Π ∗ ) , (35)and the terms Φ ∗ = ǫ N X i =1 m i m r ∗ i θ i ∗ = ǫ N X i =1 p ∗ i θ i are the Γ ∗ -dependent values of Φ and Π. Then, using the dependence of f W on ρ ,a solution to the Liouville equation, we may show that f W must satisfy a conservedequation similar to that of the previous section. Since ρ satisfies (4), we find ∂ f W∂t = Z ∆ ∂ρ∂t (Γ ∗ , t ) d Γ ∗ = − Z ∆ N X i =1 p ∗ i m i · ∂ρ∂r ∗ i + F i · ∂ρ∂p ∗ i ! d Γ ∗ = Z ρ N X i =1 p ∗ i m i · ∂ ∆ ∂r ∗ i + F i · ∂ ∆ ∂p ∗ i ! d Γ ∗ = − Z ρ N X i =1 p ∗ i m i · ∂R ∗ ∂r ∗ i + F i · ∂R ∗ ∂p ∗ i ! ǫ ∂ ∆ ∂ Φ d Γ ∗ − Z ρ N X i =1 p ∗ i m i · ∂P ∗ ∂r ∗ i + F i · ∂P ∗ ∂p ∗ i ! ǫ ∂ ∆ ∂ Π d Γ ∗ = ǫ ∂∂ Φ (cid:18)Z ρ ∆ L R ∗ d Γ ∗ (cid:19) + ǫ ∂∂ Π (cid:18)Z ρ ∆ L P ∗ d Γ ∗ (cid:19) = − ǫ ∂∂ Φ (cid:18)Z ρ ∆ P ∗ m d Γ ∗ (cid:19) − ǫ ∂∂ Π Z ρ ∆ N X i =1 F i θ i d Γ ∗ ! Thus, the reduced probability density, f W , satisfies ∂ f W∂t = − ǫ ∂∂ Φ (cid:18)Z ρ ∆ Π ∗ b m d Γ ∗ (cid:19) − ǫ ∂∂ Π (cid:18)Z ρ ∆ f d Γ ∗ (cid:19) (36)which is in conservative form.Next, we conduct a multiscale analysis in order to determine ρ up to O ( ǫ ). Sim-ilar to the previous section, the N-atom probability density, ρ (Γ , t ), is assumed tobe expressed as a function of two additional arguments, Υ(Γ , t, · , · ) in such a waythat when the last two arguments are evaluated at Φ and Π, ρ is obtained, i.e ρ (Γ , t ) = Υ(Γ , t, Φ , Π). Instead of labeling this new function, we will just extendour previous notation and refer to it as ρ (Γ , t, Φ , Π).We apply the Liouville operator to ρ (Γ , t, Φ , Π) and invoke the chain rule to find ∂ρ∂t = −L ρ − N X i =1 p i m i · d Φ dr i ∂ρ∂ Φ − N X i =1 F i · d Π dp i ∂ρ∂ Π . ∂ρ∂t = −L ρ − ǫ Pm ∂ρ∂ Φ − ǫf ∂ρ∂ Π . (37)Here, we are writing L instead of L because these derivatives are taken at constantvalues of Φ and Π. By introducing (30) and (32) into (37), the Liouville equation (4)transforms into a multiscale equation as ∂ρ∂t = ( L + ǫ L ) ρ (38)where L = − N X i =1 (cid:20) p i m i · ∂∂r i + F i · ∂∂p i (cid:21) (39)and L = − Π b m · ∂∂ Φ − f · ∂∂ Π . (40)Again, it must be noted that L and L , while seemingly exact in definition, differbecause the differentiation in L is performed at constant values of order parametersΦ and Π. Additionally, the differentiation in L is performed at fixed values of Γ.Unlike the previous section, the operators (10) and (11) are now the same as thatof the previously mentioned papers [2] and [3] since the conjugate momentum Π isformulated as an order parameter.Next, we assume that ρ may be expressed as a power series in ǫ : ρ = ∞ X n =0 ρ n ǫ n (41)A set of time variables, defined via t n = ǫ n t , is introduced to capture effects ofprocesses occurring on the various timescales. The chain rule implies ∂∂t = ∞ X n =0 ǫ n ∂∂t n . (42)We then expand (38) using (41) and (42) and separate ǫ scales. Define for n ∈ { }∪ N ,Λ n = ∂∂t n − L n (43)where we take L n = 0 for n >
1. The expansion yields the equationsΛ ρ = 0 , and for n ∈ N , Λ ρ n = − n X i =1 Λ i ρ n − i . ρ is taken to be independent of t . Thus, to lowest orderwe find L ρ = 0 . (44)This implies ρ is a function of the conserved dynamical variables, notably the totalenergy H , as well as of Φ and Π. The latter occurs because the derivatives ∂∂r i and ∂∂p i in L are to be taken at constant Φ and Π, and thus L Φ = L Π = 0. As before,we define H = N X i =1 (cid:18) p i m i + V (Γ r ) (cid:19) and notice that L H = 0.Using the entropy maximization principle and proceeding as in [2], one arrives atthe nanocanonical solution to (44): ρ = e − βH W (Φ , Π , t ) Q ≡ ˆ ρW, (45)where Q ( β, Φ , Π) = Z ∆(Γ ∗ , Φ , Π) e − βH ∗ d Γ ∗ . (46)This form of the nanocanonical solution is slightly different from that of [2] since itwas stated in that article that Q is independent of Π. We find that this is not thecase and determine the exact manner in which the Π dependence can be computedin the Appendix. Here, ∆ is defined as in (35) and H ∗ is the Γ ∗ -dependent value of H . We label t as the collection of slow time variables t = { t , t , ... } . As before, wedefine the thermal average of a given dynamical variable, A (Γ) by A th = Z ∆ ˆ ρ A (Γ ∗ ) d Γ ∗ (47)and use the Gibbs Hypothesis:lim t →∞ t Z − t e −L s Ads = A th (48)for all dynamical variables A (Γ). Notice that the thermal average, and thus the state-ments (47) and (48), depend upon the new order parameter Π because of (35).To O ( ǫ ) one finds Λ ρ = − Λ ρ (49)15nd using the semigroup e L t defined in Section 2, this equation admits the solution ρ = e L t A − Z t e L ( t − t ′ ) Λ ρ dt ′ = e L t A − Z t e L ( t − t ′ ) (cid:20) ˆ ρ ∂W∂t + Πˆ m ∂ ˆ ρ∂ Φ W + Πˆ m ˆ ρ ∂W∂ Φ + f ∂ ˆ ρ∂ Π W + f ˆ ρ ∂W∂ Π (cid:21) dt ′ = e L t A − Z t e L ( t − t ′ ) (cid:20) ˆ ρ ∂W∂t − Πˆ m β ˆ ρf th W + Πˆ m ˆ ρ ∂W∂ Φ + βf ˆ ρ Πˆ m W + f ˆ ρ ∂W∂ Π (cid:21) dt ′ where we have used the results ∂ ˆ ρ∂ Φ = − β ˆ ρf th and ∂ ˆ ρ∂ Π = β ˆ ρ Πˆ m . These derivatives will be verified in the Appendix. The first order initial condition A is, for now, undetermined and has the dependence A (Γ , Φ , Π , t ). As a consequenceof the cross-level communication inherent to multiscale analysis, the behavior of ρ at large t provides information about the t -dependence of W , while the analysis of(36) provides a necessary condition on A that ensures the equation determining f W is closed. Letting s = t ′ − t , one obtains ρ = e L t A − t ˆ ρ ∂W∂t − ˆ ρβ Πˆ m W · Z − t e −L s (cid:0) f − f th (cid:1) ds − t ˆ ρ Πˆ m · ∂W∂ Φ − ˆ ρ ∂W∂ Π · Z − t e −L s f ds. (50)Next, we remove secular behavior from the t -dependence in ρ . As before, weassume that ρ remains bounded as t → ∞ . Using (50), it can be seen that if ρ grows in t , it must do so at least linearly. Hence, we may ensure that ρ does notgrow in t by requiring that lim t →∞ t ρ = 0. We then divide by t , take the limit as t → ∞ in the above equation, and use (48). Assuming the first order initial data isin the nullspace of L , this yieldslim t →∞ ρ t = − ∂W∂t − Πˆ m · ∂W∂ Φ − f th · ∂W∂ Π . Hence, we find Λ th ≡ ∂W∂t + Πˆ m · ∂W∂ Φ + f th · ∂W∂ Π = 0 (51)and W satisfies a Liouville equation in ( t , Φ , Π) space. Note that (51) follows re-gardless of the choice of A in the nullspace of L . Using this in (50), we find16 = A − ˆ ρ (cid:20) β Πˆ m W + ∂W∂ Π (cid:21) Z − t e −L s (cid:0) f − f th (cid:1) ds (52)concluding the O ( ǫ ) analysis, although A and W are not yet determined.As before, the correction to the reduced probability density depends only on ρ and ρ up to O ( ǫ ). Hence, define for every n = 0 , , , ... f W n (Φ , Π , t ) = Z ∆(Γ ∗ , Φ , Π) ρ n (Γ ∗ , t ) d Γ ∗ (53)so that, using (34) and (41), we may write f W = ∞ X n =0 ǫ n f W n . (54)In addition, we may expand f W and ρ in powers of ǫ as in (41) and (54). Using (48),the lowest order correction, f W , can be calculated as f W = Z ∆(Γ ∗ , Φ , Π) ρ (Γ ∗ , t ) d Γ ∗ = Z ∆ ρ (Γ ∗ , t, Φ ∗ , Π ∗ ) d Γ ∗ = Z ∆ ˆ ρ W (Φ ∗ , Π ∗ , t ) d Γ ∗ = W (Φ , Π , t ) . For ǫ →
0, one may see that f W → W . Hence, as the long time scales tend to zero,the correction tends to the reduced probability density. The O ( ǫ ) correction can bedetermined using (47). Notice that (48) implies (cid:0) e −L τ A (cid:1) th = A th for any finite valueof τ ∈ R . Using this and (52), we find f W = Z ∆ ρ (Γ ∗ , t ) d Γ ∗ = Z ∆ (cid:20) A − ˆ ρ (cid:18) ∂W∂ Π + β Π ∗ ˆ m W (cid:19) · Z − t dse −L s (cid:0) f − f th (cid:1)(cid:21) d Γ ∗ = Z ∆ A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ − Z ∆ ˆ ρ (cid:18) ∂W∂ Π + β Π ∗ ˆ m W (cid:19) · Z − t e −L s (cid:0) f − f th (cid:1) dsd Γ ∗ (cid:12)(cid:12)(cid:12)(cid:12) (Φ , Π)=(Φ ∗ , Π ∗ ) = Z ∆ A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ − (cid:20)Z ∆ ˆ ρ Z − t e −L s (cid:0) f − f th (cid:1) dsd Γ ∗ (cid:21) · (cid:18) ∂W∂ Π (Φ , t ) + β Πˆ m W (Φ , t ) (cid:19) = Z ∆ A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ − Z − t (cid:0) f ( s ) − f th (cid:1) th ds · (cid:18) ∂W∂ Π + β Πˆ m W (cid:19) = Z ∆ A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ f ( s ) = e −L s f . Now, we may write ρ (Γ , t ) in terms of itsexpansion up to O ( ǫ ). Using (45) and (52) in the first term on the right side of (36),we find Z ∆ Π ∗ b m ρ (Γ ∗ , t ) d Γ ∗ = Z ∆ Π ∗ b m ρ (Γ ∗ , t, Φ ∗ , Π ∗ ) d Γ ∗ = Z ∆ Π ∗ b m ( ρ + ǫρ ) d Γ ∗ = Z ∆ ˆ ρ Π ∗ b m W (Φ ∗ , Π ∗ , t ) d Γ ∗ + ǫ Z ∆ Π ∗ b m A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ − ǫ (cid:20)Z ∆ ˆ ρ Π ∗ b m (cid:18) β Π ∗ ˆ m W + ∂W∂ Π (cid:19) · Z − t dse −L s (cid:0) f − f th (cid:1) d Γ ∗ (cid:21) = Π b m W (Φ , Π , t ) + ǫ Π b m Z ∆ A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ − ǫ Π b m (cid:18) β Πˆ m W + ∂W∂ Π (cid:19) · (cid:18)Z − t ( f ( s ) − f th ) ds (cid:19) th = Π b m W (Φ , Π , t ) + ǫ Π b m Z ∆ A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ . Similarly, the second term on the right side of (7) becomes Z ∆ f ρ (Γ ∗ , t ) d Γ ∗ = Z ∆ f ρ (Γ ∗ , t, Φ ∗ , Π ∗ ) d Γ ∗ = Z ∆ f ( ρ + ǫρ ) d Γ ∗ = Z ∆ ˆ ρf W (Φ ∗ , Π ∗ , t ) d Γ ∗ + ǫ Z ∆ f A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ − ǫ (cid:20)Z ∆ ˆ ρf (cid:18) β Π ∗ ˆ m W + ∂W∂ Π (cid:19) · Z − t dse −L s (cid:0) f − f th (cid:1) d Γ ∗ (cid:21) = f th W (Φ , Π , t ) + ǫ Z ∆ f A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ − ǫ (cid:18) β Πˆ m W + ∂W∂ Π (cid:19) (cid:18)Z − t f (0) · ( f ( s ) − f th ) ds (cid:19) th = f th W (Φ , Π , t ) + ǫ Z ∆ f A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ − ǫθ (cid:18) β Πˆ m W + ∂W∂ Π (cid:19) where θ = Z − t (cid:2) ( f (0) · f ( s )) th − f th · f th (cid:3) ds. (55)Thus, (36) becomes ∂ f W∂t = − ǫ ∂∂ Φ · (cid:18) Π b m W (Φ , Π , t ) (cid:19) − ǫ ∂∂ Φ · (cid:18) Π b m Z ∆ A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ (cid:19) − ǫ ∂∂ Π · (cid:0) f th W (Φ , Π , t ) (cid:1) + ǫ ∂∂ Π · (cid:20) − Z ∆ f A (Γ ∗ , Φ ∗ , Π ∗ , t ) d Γ ∗ + θ (cid:18) β Πˆ m W + ∂W∂ Π (cid:19)(cid:21) . f W and f W , we can expand the reduced probability densityas f W = f W + ǫ f W . Then, isolating f W k terms and imposing the condition that A must stay bounded for large Φ and Π, we find that this equation is closed only if A = 0. Thus, up to O ( ǫ ), the conserved equation (36) becomes the Fokker-Planckequation: ∂ f W∂t = − ǫ ∂∂ Φ · (cid:18) Π b m f W (cid:19) − ǫ ∂∂ Π · (cid:16) f th f W (cid:17) + ǫ ∂∂ Π · (cid:20) θ (cid:18) β Πˆ m + ∂∂ Π (cid:19) f W (cid:21) . (56)Hence, under the assumption that momenta have (comparatively) large values andevolve slowly, i.e., P = O ( ǫ − ) and f = O ( ǫ ), the reduced probability density obeysa Fokker-Planck equation given by (56). We first verify the Φ derivative of ˆ ρ used in the derivation of both equations. InSection 2, the variable R is used instead of Φ (as order parameters are unscaled withrespect to ǫ ), but the statements that follow can be applied exactly to R in the samemanner as Φ. We claim ∂ ˆ ρ∂ Φ = − β ˆ ρf th . Using (17) or (45), we see that ∂ ˆ ρ∂ Φ = − e − βH Q ∂Q∂ Φ . Notice that ∂∂ Φ (cid:0) e − βH (cid:1) = 0 since Φ derivatives are to be taken at constant values ofΓ. Using (18) or (46), we see that ∂Q∂ Φ = Z ∂ ∆ ∂ Φ e − βH ∗ d Γ ∗ = − Z ∂ ∆ ∂ Φ ∗ e − βH ∗ d Γ ∗ = ǫ − Z ∆ ∂∂R ∗ (cid:0) e − βH ∗ (cid:1) d Γ ∗ Now, since this integration is performed over all values of Γ ∗ , and hence Γ ∗ is not fixed,the energy depends upon the center of mass through the potential function. Notice, wemay calculate ∂V∂R as V depends upon R implicitly in the following manner. For every i = 1 , ..., N , write the residual displacement of each atomic position as r i = σ i + Rθ i .Here the N variables { R, σ , ..., σ N − } constitute a complete set of variables as wemay write σ N in terms of each of the other σ i using the constraint N X i =1 σ i θ i = 0 , ∂V∂R = N X i =1 ∂V∂r i θ i = − ǫ f. Since the kinetic energy is independent of Γ r , we find ∂∂R ∗ (cid:0) e − βH ∗ (cid:1) = − βe − βH ∗ ∂V ∗ ∂R ∗ = βǫ f ∗ e − βH ∗ . Using this, the derivative of Q becomes ∂Q∂ Φ = ǫ − Z ∆ e − βH ∗ βǫ f ∗ d Γ ∗ = β Z ∆ e − βH ∗ f ∗ d Γ ∗ . Hence, using the modified Gibbs Hypotheses in either section (20) or (48), we find ∂ ˆ ρ∂ Φ = − e − βH Q ∂Q∂ Φ= − ˆ ρ (cid:18) Q · ∂Q∂ Φ (cid:19) = − ˆ ρβ Z Q ∆ e − βH ∗ f ∗ d Γ ∗ = − β ˆ ρ (cid:18)Z ∆ ˆ ρf ∗ d Γ ∗ (cid:19) = − β ˆ ρf th and the verification of this derivative is complete.Next, we verify the Π derivative of ˆ ρ used in Section 3. We claim ∂ ˆ ρ∂ Π = β ˆ ρ Πˆ m . First, we use (46) and rewrite Q in terms of its Π dependence using the Fouriertransform. Since the potential and kinetic energies depend exclusively on r i and p i variables respectively, we separate Q into two different integrals involving thesevariables. In addition, we separate the momentum integrals into the momenta ofparticles in the host medium, denoted Γ = { p i : θ i = 0 } and those in the nanoparticle,denoted Γ = { p i : θ i = 1 } . Using the notation K = N X i =1 p i m i (1 − θ i )20nd K = N X i =1 p i m i θ i for the kinetic energies of the host medium and nanoparticle, respectively, we write Q as Q ( β, Φ , Π) = (cid:18)Z δ (Φ − Φ ∗ ) e − βV ∗ d Γ ∗ r (cid:19) (cid:18)Z e − βK ∗ d Γ ∗ (cid:19) (cid:18)Z δ (Π − Π ∗ ) e − βK ∗ d Γ ∗ (cid:19) . (57)Keeping the Γ ∗ r and Γ ∗ integrals as they are (notice further that the Γ ∗ integral isconstant), we focus on the Γ ∗ integral. We may again split the Γ ∗ integral into oneeach in the x, y , and z directions. We will consider the integral in the x -direction,labeled I x , and state that the results we obtain will follow for the integrals in theother directions in the same manner.Now, we relabel the momenta in the nanoparticle, 1 through M , and write their x -coordinates as x through x M . Then, write the x -directional Dirac mass δ (Π x − Π ∗ x )using the Fourier transform of the function f ( k ) = 1 as δ (Π x − Π ∗ x ) = 1 √ π Z ∞−∞ exp " ik Π x − √ ǫ M X l =1 x l ! dk. Then, using the Inverse Fourier Transform on the resulting Gaussiansexp( − λk ) = 1 √ π Z ∞−∞ √ πλ exp ( − ikz ) exp (cid:18) − z λ (cid:19) dz. and (30), we find I x = 1 √ π Z Z ∞−∞ exp " ik Π x − √ ǫ M X l =1 x l ! dk · exp (cid:18) − β m x (cid:19) · · · exp (cid:18) − β m M x M (cid:19) dx · · · dx M = 1 √ π Z ∞−∞ exp ( ik Π x ) (cid:20)Z exp (cid:0) − ik √ ǫx (cid:1) exp (cid:18) − β m x (cid:19) dx (cid:21) · · ·· · · (cid:20)Z exp (cid:0) − ik √ ǫx M (cid:1) exp (cid:18) − β m M x M (cid:19) dx M (cid:21) dk = Z ∞−∞ e ik Π x r m β e − ǫm k β · · · r m M β e − ǫm M k β dk = r m πβ · · · r m M πβ Z ∞−∞ exp ( ik Π x ) exp (cid:20) − ǫk β m (cid:21) dk = C exp (cid:18) − β m Π x (cid:19) C = r m πβ · · · r m M πβ · r πβ ˆ m . We extend this in the y and z directions and multiply to find Z δ (Π − Π ∗ ) e − βK ∗ d Γ ∗ = ( C ) exp (cid:18) − β m Π (cid:19) . Thus, we can write Q ( β, Φ , Π) = Q (Φ)( C ) C exp (cid:18) − β m Π (cid:19) where Q (Φ) = (cid:18)Z δ (Φ − Φ ∗ ) e − βV ∗ d Γ ∗ r (cid:19) and C = (cid:18)Z e − βK ∗ d Γ ∗ (cid:19) . Finally, ˆ ρ can be expressed in the formˆ ρ = e − βH Q (Φ)( C ) C exp (cid:18) β m Π (cid:19) . (58)Taking a Π derivative in (58), which must be done at fixed values of Γ, we find ∂ ˆ ρ∂ Π = β ˆ ρ Πˆ m and the verification of this derivative is complete. References [1] Hopf, E., Proof of Gibbs’ Hypothesis on the Tendency Toward Statistical Equi-librium.
Proc. Natl. Acad. Sci. U. S. A. (4):333-340 (1932).[2] Miao, Y. and Ortoleva, P., All-atom multiscaling and new ensembles for dynam-ical nanoparticles. J. Chem. Phys. , J. Chem. Phys. , Revista Brasileira de Ensino de Fisica , (2):189-201 (2007).[5] Ortoleva, P., Nanoparticle Dynamics: A Multiscale Analysis of the LiouvilleEquation. J Phys. Chem , J. Chem. Phys.
528 (1998).[7] Peters, M.H.,
J. Stat. Phys.
557 (1999).[8] Renardy, M. and Rogers, R., An Introduction to Partial Differential Equations.(2000)[9] Shreif, Z. and Ortoleva, P., Curvilinear All-Atom Multiscale (CAM) Theory ofMacromolecular Dynamics,
J. Stat. Physics , to appear (2007).[10] Shea, J.E. and Oppenheim, I.,
J. Phys. Chem
Physica A247