BEYOND MULTISCALE AND MULTIPHYSICS: MULTIMATHS FOR MODEL COUPLING

. The purpose of this article is to present a uniﬁed view of some multiscale models that have appeared in the past decades in computational materials science. Although very diﬀerent in nature at ﬁrst sight, since they are employed to simulate complex ﬂuids on the one hand and crystalline solids on the other hand, the models presented actually share a lot of similarities, many of those being in fact also present in most multiscale strategies. The mathematical and numerical diﬃculties that these models generate, the way in which they are utilized (in particular as numerical strategies coupling diﬀerent models in diﬀerent regions of the computational domain), the computational load they imply, are all very similar in nature.

1. Introduction. A few decades ago, multiscale modelling used to only mean modelling based on the homogenization paradigm. Since then, the field has exploded. Multiscale simulations are present in a large variety of domains. The purpose of this article is to show that, however diverse in nature they might look at first sight, multiscale simulations definitely share common features. A common denominator is the use of an increasing variety of mathematical notions and techniques.
For illustration purposes, we consider two prototypical modelling strategies, apparently very different in nature since the first one deals with fluids, the second one deals with solid materials. The former modelling strategy is the micro-macro simulation of complex fluids. In short, it consists in coupling the macroscopic equations of conservation of mass and momentum for a viscous incompressible fluid, with a kinetic description of the evolution of the microstructures within the fluid that are responsible for the non-Newtonian character. The latter modelling strategy is the discrete-to-continuum coupling in use for the simulation of solid materials. It couples an atomistic description of matter using actual positions of atoms and microscopic interaction laws between them, and a description by continuum elasticity theory in terms of gradients of deformations and energy associated with those.
Considering these two particular cases, we will show that 1. multiscale modelling arises as a response to the deficiency of single-scale modelling; its purpose is to bypass the analytical or empirical derivation of a constitutive law at the macroscopic scale, letting the numerical simulation implicitly perform the coupling between the scales; 2. multiscale modelling is computationally costly, and is a field of predilection for parallel computing; 3. multiscale modelling is actually especially useful today as a model coupling strategy: the computational domain is divided between on the one hand regions where a computationally cheap, purely macroscopic, single-scale model is sufficient to accurately reproduce the behaviour, and on the other hand hopefully small regions where the computationally expensive multiscale model has to be employed for capturing the singularities in the physical/mechanical behaviour; 4. multiscale modelling increasingly involves models that come from different areas of Physics in the broad sense; 5. but beyond all this, the true novelty of the last two decades is that multiscale modelling requires a good knowledge of the theory and the numerics of an ever increasingly wide spectrum of mathematical objects.
In a somewhat provocative (not to be taken too seriously) style, we could say that, after the times of multiscale and multiphysics, the era of multimaths has begun.
The article is organized as follows. In Section 2, we recall the most commonly used macroscopic, single-scale models for non-Newtonian fluids and elastic solid materials. We show their limitation and explain why alternative, multiscale models are useful. We next introduce some prototypical examples of such multiscale models in Section 3. We present some adapted discretization techniques in Section 4, and then discuss the challenging mathematical and numerical issues they raise in Section 5. The section illustrates the variety of mathematical notions and techniques addressed. Its purpose is to point out the questions raised by these models and briefly indicate the approaches followed to (partially) solve these questions. The article concludes with Section 6 where the current trends are emphasized. The specificity of our presentation is that, throughout the sections, we systematically, sequentially discuss each issue for the complex fluid context and for the solid material context. It is quite unusual to gather in a single article two so different contexts. We however believe this emphasizes the similarities in the scientific process, and helps to identify global achievements and challenges transverse to the field of modelling. To conclude this introduction, we emphasize that, due to space limitations, many questions and details are left aside in our presentation. More exhaustive presentations can be found in the papers [83,72] by B. Jourdain, C. Le Bris and T. Lelièvre for fluid modelling, and the articles [20,82] by X. Blanc, C. Le Bris and F. Legoll for solid modelling.
2. Deficiency of the single-scale macroscopic description. We recall in this section the basic elements of single-scale modelling of fluids, and next of solid materials. We shall in particular mention the inevitable limitations of such a single-scale strategy.
2.1. Macroscopic models for fluids. To begin with, we recall here some elements on the modelling of viscous incompressible fluids. Consider a viscous fluid with volumic mass (or density) ρ, flowing at the velocity u, and experiencing external forces f per unit mass. Denote by T the stress tensor. The equation of conservation of mass for the fluid reads ∂ρ ∂t + div (ρ u) = 0, (2.1) while the equation expressing the conservation of momentum is ∂(ρ u) ∂t + div (ρ u ⊗ u) − div T = ρf .
For such a viscous fluid, the stress tensor reads T = −p Id + τ , where p is the (hydrodynamic) pressure, and τ is the tensor of viscous stresses. The equations are supplied with appropriate initial and boundary conditions we omit. In order to close the above set of equations, a constitutive law is needed, which relates the viscous stress τ and the velocity field u, namely τ = τ (u, ρ, . . .).

(2.3)
Expression (2.3) is symbolic. A more precise formulation may involve derivatives in time, or in space, of the various fields τ , u, ρ, . . . Examples are given below. The simplest possible situation is that of Newtonian fluids (typically the case of water), for which, by definition, τ linearly depends on the velocity u. Under appropriate assumptions, it can then be shown that the relation between τ and u necessarily takes the form τ = λ (div u) Id + 2η d, where λ and η are the two Lamé coefficients, and d denotes the (linearized) rate of deformation tensor d = 1 2 (∇u + ∇u T ), with the convention (∇u) ij = ∂u i ∂x j . The derivation proceeds, and in the case of an incompressible (div u = 0) homogeneous (ρ = ρ 0 ) Newtonian viscous fluid, the equations obtained are the celebrated incompressible Navier-Stokes equations    ρ 0 ∂u ∂t + (u · ∇) u − η∆u + ∇p = ρ 0 f , div u = 0.

(2.4)
When the fluid is not Newtonian, that is, in many practically relevant situations (most of the fluids of the everyday life are indeed not Newtonian: blood, wine, oil, . . . ), a constitutive equation (2.3) needs to be identified. It is of course specific to the fluid under consideration. Equations (2.1)- (2.2) together with the constitutive equation (2.3) (and possibly equations for the energy and the temperature) then form the model for the fluid. For non-Newtonian fluids, many constitutive laws, and thus many purely macroscopic models, exist. All are based upon considerations of continuum mechanics. It is usual to decompose the stress tensor τ as the sum τ = τ n + τ p , where τ n denotes the Newtonian contribution and τ p denotes the part of the stress (called non-Newtonian or extra stress) that cannot be modelled in the Newtonian manner. The bottom line is then to write an equation, in the vein of (2.3), ruling the evolution of the non-Newtonian contribution τ p and/or encoding a relation between the latter and other quantities characterizing the fluid dynamics, such as the deformation tensor d, or ∇u itself. One famous example is the Oldroyd B model , which, in non-dimensional form, writes: The Reynolds number Re > 0, the Weissenberg number (ratio of the characteristic relaxation time of the microstructures in the fluid to the characteristic time of the fluid) We > 0 and ǫ ∈ (0, 1) are the non-dimensional parameters of the model. The Oldroyd B model is not capable of reproducing many experimentally observed behaviors. Refined macroscopic models for viscoelastic fluids have thus been derived, allowing for a better agreement between simulation and experiments: the Giesekus model, the Phan-Thien Tanner model, the FENE-P model, etc. Each model correspond to a particular constitutive law. Overall, they yield better results than the Oldroyd B model, and satisfactorily agree with several prototypical experiments on simple flows. Of course, in terms of scientific computing, solving the three-field problem (2.5) is much more difficult and computationally demanding than the 'simple' Newtonian problem (2.4). However, the major scientific difficulty is neither a mathematical one nor a computational one. The major difficulty is to derive a constitutive equation (2.3). It requires a deep qualitative and quantitative understanding of the physical properties of the fluid under consideration. And there are many fluids, and many experimental situations. For many non-Newtonian fluids, complex in nature, reaching such an understanding is therefore a challenge. Moreover, even if such an equation is approximately known, evaluating the impact of its possible flaws on the final outcome of the simulation is not an easy issue. It can only be completed a posteriori, comparing the results to actual experimental observations, when the latter exist, and they do not always exist. The difficulty is all the more prominent that non-Newtonian fluids are very diverse in nature.
All this, in its own rights, motivates the need for alternative strategies, based on an explicit microscopic modelling of the fluid. This gives rise to the so-called micromacro models (see Figure 1). The lack of information at the macroscopic level is then circumvented by a multiscale strategy consisting in searching for the information at a finer level (where reliable models do exist, based on general conservation equations, posed e.g. on the microstructures of the fluids). The latter information is then inserted in the equations of conservation at the macroscopic level.
We shall now see that similar considerations fully apply to our second example, although coming from an entirely different context: multiscale modelling for solid materials. 2.2. Continuum elasticity theory. We now briefly recall the standard mechanical description of a solid material subjected to forces. For simplicity and brevity, the setting is static. We refer e.g. to the monograph [36] of Ph. G. Ciarlet for a complete mathematically oriented presentation including time-dependent situations. We denote by D the reference domain that the material occupies at rest, by ϕ the deformation it is subjected to, i.e. the map from D to R 3 that gives the current position of the material. We also denote by u(x) = ϕ(x) − x the displacement, and by F = ∇ϕ : D −→ M 3 the gradient of deformation, where M 3 denotes the space of square matrices of size 3 × 3. The general equations that describe the equilibrium of our sample material, when subjected to body forces f and boundary forces g, read − div T = f in D, (2.6) with the boundary condition T · n = g on ∂D. Here, T denotes the stress tensor (more precisely the first Piola-Kirchhoff stress tensor), and n is the unitary outward normal vector on ∂D. Equation (2.6) plays, in our present discussion, the role of equations (2.1)-(2.2) from the previous section.
In order to close the equation (2.6), a relation is needed between the stress tensor T and the kinematic description of the material, provided by the fields ϕ, u or F . In contrast to equations (2.6) which are general, the relation linking T to, say, ϕ, depends on the material considered. In such a relation is indeed encoded the physical and mechanical nature of the material. Formally, such a closure relation reads T = T (x, ϕ(x), ...), (2.7) and is called, similarly to (2.3), a constitutive relation, or a law of behavior. Again, equation (2.7) is symbolic: derivatives of ϕ may also be inserted, as well as other points than x (or, in a time-dependent setting, times anterior to the time t at which the stress tensor is evaluated). The relation may be a differential equation, a partial differential equation, an integral equation, etc.
A usual assumption in deriving such a law is the elasticity of the material: T = T (x, F (x)). Hyperelasticity is often additionally assumed: the material is assumed to dissipate no energy during a cyclic deformation. From this is inferred the existence of a density W of mechanical energy such that where A is the set of all admissible deformations ϕ. We omit several technical (and actually also fundamental) issues in this expository presentation, see for example the works by J. M. Ball [9] and Ph. G. Ciarlet [36]. The most famous example of density W is provided by linearized elasticity. Then, W is the simple quadratic form This model has the advantage that the Euler-Lagrange equations of the variational problem (2.9), with W given by (2.10), are linear. Note however that this density W is not invariant under rigid rotations. For a given material, the derivation of the constitutive law (2.7), or equivalently in the hyperelasticity setting, that of the density W of mechanical energy, is a central theoretical question, and a challenging practical issue. Then, the resolution of (2.6), or the minimization of (2.9), is the purpose of the numerical simulations performed. Despite all these efforts, and all the expertise accumulated, the strategy for the derivation of a constitutive law for a solid material has the exact same flaws as that previously mentioned for fluids. Therefore, similarly to the context of fluid mechanics above, it is useful to develop an alternative strategy: multiscale modelling.
In conclusion to this section, we emphasize that bypassing the macroscopic constitutive law is one of the major purpose of multiscale modelling in materials science. In principle, the constitutive law is expected to encode the physical nature of the material, so that the phenomena taking place at all the scales finer than the macroscopic ones are accounted for. A multiscale model is also aimed at encoding such a relation, but, in contrast, without translating it into an explicit mathematical relation. The physical nature of the material is inserted via an explicit microscopic description, which is in turn coupled with the usual macroscopic description of the material. Loosely speaking, the derivation of the constitutive law is implicitly performed by the simulation itself. The amount of physical intuition needed is expected to be smaller, the computational effort will compensate for it. Likewise, it is hoped that less modeling assumptions will be needed: ideally, a universal microscopic model is inserted in the universal macroscopic description. In doing so, the sources for inaccuracies are easier to identify, and the assessment of the quality of the result is simpler.
3. Multiscale modelling for fluids and solids. The variety of situations in which multiscale modelling is now employed is so large that, for our exposition, we have to pick specific examples of materials. For fluid modelling, we choose infinitely dilute solutions of polymers, and for solid modelling, crystalline solids. We mention that, for instance for fluids, the modelling of liquid crystals (see the reference books by M. Doi and S.F. Edwards [49] or H.C.Öttinger [103], and the articles by H. Zhang and P.W. Zhang [121], P. Constantin, I. Kevrekidis, and E.S. Titi [39] or G. Forest, Q. Wang, and R. Zhou [58] for some mathematical and numerical studies), or suspensions (see the paper [63] by P. Hébraud and F. Lequeux for an interesting model and the series of works [30,31,32,15] Figure 2. Each polymeric chain, itself possibly consisting of thousands of atoms, is admittedly well modeled using a coarse-grained description. The simplest possible such description consists in a single end-to-end vector, called a dumbbell, that models the total length and overall direction of the chain, see Figure 3. in each grain, atoms are located on (or close to) a perfect lattice, whose orientation is different from one grain to the other. Each grain can be described by a continuum model. In constrast, the interface between two grains is in general too small and too heterogeneously deformed for a continuum model to hold.
Our purpose is now to explain how the derivation of an explicit constitutive law can be bypassed in this context. The success of the enterprise owes to the existence of a well established kinetic theory for solutions of polymeric chains, see the monographs [17,18]  Let us now denote ψ(t, x, r) the probability density for the end-to-end vector r modeling the polymer chains at macropoint x and time t. The variation of ψ follows from three different phenomena: 1. a hydrodynamic force: the dumbbell is elongated or shortened because of the interaction with the fluid; its two ends do not necessarily see the same X Figure 3. The dumbbell model: the end-to-end vector X is the vector connecting the two beads that model the entire polymeric chain. It is supposed to accurately model the typical behaviour of the entire chain.
fluid velocity, the slight difference in velocities (basically ∇u(t, x) r) results in a force elongating the dumbbell ζ∇u(t, x) r, where ζ denotes a friction coefficient; 2. an entropic force F issued from the coarse-graining procedure and which is reminiscent of the actual, much more complex, geometry of the entire polymeric chain; 3. a Brownian force, modelling the permanent collisions of the polymeric chain with solvent molecules, which (randomly) modifies its evolution.
The equation of conservation of momentum reads as the following evolution equation on ψ: x, r). (3.1) The three terms of the righthand side of (3.1) respectively correspond to the three phenomena listed above. A crucial point is that, in this right-hand side, all differential operators acting on ψ are related to the variable r of the configuration space, not of the ambient physical space. In contrast, the gradient of the left-hand side is the usual transport term in the physical space u · ∇ x . In the absence of such a transport term (this will indeed be the case for extremely simple geometries, such as that of a Couette flow presented below), (3.1) is simply a family of Fokker-Planck equations posed in variables (t, r) and parameterized in x. These equations are coupled only through the macroscopic field u. When the transport term is present, (3.1) is a genuine partial differential equation in all variables (t, x, r). It is intuitively clear that the latter case is much more difficult, computationally and mathematically. Once ψ is obtained, its contribution to the total stress, and, further, its impact on the macroscopic flow, need to be formalized. Elementary considerations of continuum mechanics (see Figure 4) show that the contribution to the stress is given by the so-called Kramers formula, where n p denotes the total number of polymeric chains per unit volume.
The complete system of equations combines the equation of conservation of momentum at the macroscopic level, the incompressibility constraint, the Kramers formula, and the Fokker-Planck equation for the distribution of the end-to-end vector. In non-dimensional form, it reads: ∆ r ψ(t, x, r). An alternative description of the evolution of polymeric chains is provided by the stochastic viewpoint. This viewpoint is actually extremely useful in practice, because it allows to circumvent the difficulties related to the high-dimensionality of the Fokker Planck equation (3.1). It is indeed to be borne in mind that when the geometric description of the chain is richer than a coarse, dumbbell model, then the end-to-end vector is replaced by a possibly highly multidimensional vector, and the dimensionality of the Fokker-Planck equation grows correspondingly.
As an alternative to (3.1), we may model the mesoscopic part of the system by the set of stochastic differential equations: (3.4) where X t (x) denotes the stochastic process modeling the conformation of the polymeric chain at x at time t. The stress is then given by where n p is the concentration of polymers. The coupled system is thus: In both the above systems (3.3) and (3.6), the force F needs to be made specific.
In full generality, it is assumed of the form F (X) = π ′ ( X ) X X , for a given potential π. The simplest potential π is the quadratic potential π Hook (l) = H l 2 2 . A peculiarity of this choice is that the multiscale model is then equivalent to the Oldroyd B model (2.5). Other choices of forces allow to capture important physical behaviours and yield models that are genuinely multiscale: the FENE force corresponding to the potential π FENE (l) = − bkT 2 ln 1 − l 2 bkT /H is the most commonly used type of such forces. In general, it is believed that a multiscale model is more accurate than a purely macroscopic model, and that the description of the microstructure need not be sophisticated to give excellent results, capturing the right qualitative physics being the only important issue (see the FENE force in contrast to the Hookean force).

3.2.
Atomistic-to-continuum modelling for crystals. When a model for the density of mechanical energy (2.8) for a given solid material is unknown, or is inapplicable because all the assumptions of regularity of stress and strain for continuum mechanics to apply are not satisfied, one has to return to the microscopic scale, and use appropriate atomistic models. The task is then to deduce from the atomistic scale an appropriate density of mechanical energy, that is, a specific form for a constitutive law (2.7).
To this end, we proceed as follows: we give ourselves a microscopic description of the sample, next we fix the deformation ϕ that is imposed to the material at the macroscopic level and directly apply it to the atomic sites at the microscopic level, and we search for the macroscopic limit of the energy obtained.
Somewhat similarly to what has been performed above for fluids using the kinetic description and Kramers formula to obtain the stress, we expect to derive an explicit form for the macroscopic mechanical energy of the sample deformed by ϕ, thus a link between the density of mechanical energy W , and the energy model used at the microscopic level.
For simplicity, we consider an elementary two-body interaction potential at the atomic scale, and place the atoms on a particularly simple geometrical arrangement. At the microscopic scale, the material occupies the domain D, and we assume its substructure is the truncation D∩L of a periodic lattice L. Every pair of sites (x i , x j ) in D∩L interacts with a pair potential V (x i −x j ). We assume the periodic lattice has a cubic unit cell, and that the atomic site stands at the center of this cell. The length of the cell is denoted by ε, and vanishes in the macroscopic limit. The potential V , taken radially symmetric, is assumed to be smooth, and have compact support.
No body force, nor boundary force, is applied to the material. The deformation we apply to the sample is assumed to be smooth. All these assumptions aim at avoiding unnecessary technicalities. Let us now apply the deformation ϕ (see Figure 5). The microscopic energy formally reads: In fact the above sum is truncated over a finite domain, containing N = ε −1 sites per dimension. It is normalized by the total number of particles (since the energy is an extensive quantity). Note also that we need to rescale all distances by a factor 1 N = ε, so that the equilibrium length (i.e., say, the length r that minimizes the function r −→ V (r)) is also of order ε for consistency. The purpose is now to identify the limit of (3.7) when N −→ ∞. This means that we both let the truncated lattice go to the whole infinite lattice at the microscopic scale, and let the lattice size vanishes, so as to pass to the macroscopic limit.
In this simple case, the analysis is obvious. As ϕ is smooth, we write its Taylor expansion at the first order and neglect the terms of higher order: We insert this in the potential V , and in turn search for the first order. Arguing formally, omitting some technicalities in particular related to boundary terms that can be easily handled, and using the periodicity of the lattice, we see that the limit formally reads Denoting by we observe that (3.9) is a Riemann sum in Ψ, and converges as N −→ +∞ to 1 where |D| is the volume of the domain D. At the macroscopic level, the (nondimensional) density of mechanical energy therefore reads We observe that this energy is indeed the energy of the original periodic lattice deformed by the linear map ∇ϕ(x) at the macroscopic point x. A simple example is the case of the quadratic potential V ∝ |r − r eq | 2 where r eq denotes some equilibrium interatomic distance. It is easily seen, at least formally and in dimension one (and if particles stay ordered), that the above derivation then yields the linearized elasticity model (2.10) as a limit in (3.12). For the simple case treated here (periodic lattice, pair-potential interaction), this derivation of a macroscopic density of energy has been known for long. The work [23]  It is illustrative to mention a slight extension of the simple case presented above, because it establishes an enlightening connection between the elements of complex fluids modelling introduced in the previous section and the present setting. We consider, instead of a periodic lattice of atomic sites as the microscopic model for the crystalline structure, a random perturbation of this arrangement. More explicitly, we assume that, before rescaling, the atomic sites now stand at points x i (ω) = i + X i (ω), where i is the three-dimensional integer index that also denotes an arbitrary point of Z 3 , and where X i (ω) is a set of ergodic stationary random variables. Then, arguing as above, but this time additionally using the ergodic theorem, the density of mechanical energy obtained in the macroscopic limit reads instead of (3.12). We refer to the original paper [23] by X. Blanc, C. Le Bris and P.-L. Lions for the details of the derivation. The similarity of expression (3.14) with the Kramers formula (3.5) is evident. The stress (or equivalently its primitive, the density of mechanical energy) is inferred from an average taken on microscopic configurations. Another, more formal similarity is that the reader now understands how the random, microscopic description is genuinely coupled with the deterministic, macroscopic description. The mathematics behind the models necessarily embraces both aspects. The above derivation of the continuum mechanics energy (3.11) from the atomistic energy (3.7) is based on a smoothness assumption on the deformation ϕ. In many situations of interest (for instance when dislocations appear, as in the nanoindentation simulation shown in Figure 6), such an assumption does not hold in the whole domain, and one cannot use a model based on (3.11). Using an atomistic model in the whole domain is not possible either, due to its prohibitive computational cost. We now describe a coupling method, whose motivation is based on the observation that, actually, the deformation that we are after is not smooth in only a small part of the solid. So, a natural idea is to try to take advantage of both models, the continuum mechanics one and the atomistic one, and to couple them, in a domain decomposition spirit. The description below is a toy-version of the QuasiContinuum Method (QCM), as presented in its initial version [116,115] by E.B. Tadmor, M. Ortiz and R. Phillips. The method has been next amended, and we will describe a more mature formulation in Section 4.2. The microscopic energy of a given deformation ϕ is (3.7), that we recast as denotes the energy of atom i. We now split the computational domain D into two non-overlapping subdomains, where D reg is a domain where the deformation is expected to be smooth (see Figures 7 and 8). Consequently, in D reg , we can approximate the atomistic energy by a continuum mechanics expression, in the spirit of (3.11)-(3.12). In D sing , we cannot make this approximation, and we keep the original atomistic model. For a given deformation ϕ, we hence write where Ψ is defined by (3.10). The expression (3.17) is next approximated by Recall that Ψ(x) = W (∇ϕ(x)) is the energy of an atom in an infinite system deformed by the linear map ∇ϕ(x). The domain D reg is hence often called the local zone, since only the knowledge of the deformation ϕ in a neighborhood of x is needed to compute the energy of an atom located at x. The situation is completely different in D sing . We make the standard assumption that, beyond a cut-off radius r c , the potential V vanishes. Hence only the atoms inside the ball of center x and of radius r c interact with an atom located at x ∈ D. In view of (3 .18)   In practice, a finite element method is employed in D reg to approximate the deformation ϕ. The first term of (3.18) is computed by a numerical quadrature, involving a few quadrature points in each finite element. In the case of a piecewise linear approximation, the degrees of freedom are the values of ϕ on the vertices of the mesh in D reg , and the positions ϕ (x j /N ) of the atoms belonging to D sing (that is, those such that An equilibrium configuration ϕ eq is defined as the global minimizer of an approximation of (3.18), along the above lines (note that some methods alternatively focus on local minimizers, rather than global ones, or on critical points of the energy). This configuration only makes sense if it is smooth in D reg . Indeed, in that region, we replaced the original atomistic energy by a continuum mechanics energy, which is possible only under some regularity assumptions on the deformation. If, according to a certain criterion which is part of the multiscale method, the deformation ϕ eq is considered not smooth in a subdomain D irreg of D reg , then the partition (3.16) is updated accordingly, and a new partition is considered: the domain D irreg , which was described at the continuum scale, is now described at the atomistic scale. A coupled energy of type (3.18), based on this new partition, is considered. This yields the following iterative procedure. Start with a given partition of D, and then iterate over the steps: • on the basis of the current partition, define the coupled energy (3.18); • solve the variational problem associated with that energy (we denote ϕ eq its global minimizer); • update the partition on the basis of ϕ eq , along the lines of the above discussion. We now assess the computational gain. Evaluating the energy of atom i is more expensive when using the nonlocal model (3.15) than when using the local model (3.10). Indeed, in the local formulation (3.10), it is easy to compute the position of all the atoms in the ball of cutoff radius around atom i. On the other hand, computing the energy of atom i according to the nonlocal expression (3.15) requires to know all the positions ϕ (x j /N ) of the atoms j in the cutoff radius ball centered around atom i. This is a computationally demanding task, since one first has to determine in which finite element each of these atoms j is. Besides, in practice, finite elements of the subdomain D reg contain a very large number of atoms n e . As a consequence, computing the contribution of the local zone D reg to the total energy is much cheaper with the coarse-grained model than with the reference model. Indeed, in the coarse-grained model, this contribution, which is exactly the first term of (3.18), can be computed by evaluating the energy of an extremely small number of atoms (those at the quadrature points, where Ψ needs to be evaluated). In constrast, in the reference model, one has to compute the energy of a large number n e of atoms. So the computational gain comes from a two-fold argument: in D reg , the energy of much fewer atoms needs to be evaluated, and each of these evaluations is cheaper, since a local model is used rather than a nonlocal one.

Polymeric fluids.
We now turn to the discretization of the systems (3.3) or (3.6). Most of the numerical methods employed for the simulation of non Newtonian fluids are based upon a finite element discretization in space and Euler schemes in time, using a semi-explicit scheme: at each timestep, the velocity is first solved for a given stress, and then the stress is updated, for the corresponding velocity.
Three types of difficulties typically arise: 1. An inf-sup condition must be satisfied by the spaces respectively used for the discrete velocity, the pressure and the stress (in order for the discretization to be stable for ǫ close to 1). 2. The advection terms need to be appropriately discretized, in the equation of conservation of momentum, in the equation on τ p in (2.5), in the equation on ψ in (3.3), on in the SDE in (3.6). 3. The nonlinear terms require, as always, special attention. On the one hand, some nonlinear terms stem from the coupling: ∇uτ p + τ p (∇u) T in (2.5), ∇uX t in (3.6) or div r (∇u r ψ(t, x, r)) in (3.3). On the other hand, for rheological models more complicated than the Oldroyd-B or Hookean dumbbell models, some nonlinear terms come from the model itself (see the entropic force F (X t ) in (3.6) for FENE models for example).
Besides, for both micro-macro models and purely macroscopic models, one central difficulty of the simulation of viscoelastic fluids is the so-called High Weissenberg Number Problem (HWNP), see for instance R. Keunings [74,75]. It is indeed observed that numerical simulations do not converge when We is too large and that the problem gets all the more delicate as the mesh is refined. Because the reader might be more familiar with a purely deterministic formulation, we begin by outlining the numerical approach for problems of type (3.3). The discretization of the Fokker-Planck equation in (3.3) is typically performed using spectral methods (see A. Lozinski [92], J.K.C. Suen, Y.L. Joo and R.C. Armstrong [114] or D.J. Knezevic and E. Süli [77]). It is not easy to find a suitable variational formulation of the Fokker-Planck equation, and an appropriate discretization that satisfies the natural constraints on the probability density ψ (namely non negativity, and normalization). We refer to C. Chauvière and A. Lozinski [34,93] for suitable discretizations in the FENE case. A major difficulty is the possible high dimensionality of the Fokker-Planck equation. In the context of polymeric fluid flow simulation, when the polymer chain is modelled by a chain of N + 1 beads linked by N springs, the Fokker-Planck equation is a parabolic equation posed on a 3Ndimensional domain. Some numerical methods have been developed to discretize such high dimensional problems. The idea is to use an appropriate Galerkin basis, whose dimension remains limited when dimension grows. We refer to P. Delaunay, A. Lozinski and R.G. Owens [41], T. von Petersdorff and C. Schwab [119], H.-J. Bungartz and M. Griebel [29] for the sparse-tensor product approach, to L. Machiels, Y. Maday, and A.T. Patera [95] and D.J. Knezevic and A.T. Patera [78] for the reduced basis approach and to A. Ammar, B. Mokdad, F. Chinesta and R. Keunings [3,4] and to C. Le Bris, T. Lelièvre, and Y. Maday [84] for a method coupling a sparse-tensor product discretization with greedy algorithms used in approximation theory.
The numerical approach for the coupled system (3.6) involving the stochastic differential equations formulation is different. A Monte Carlo method is employed to discretize the expectation: at each macroscopic point x (i.e. at each node of the mesh once the problem is discretized), many replicas (or realizations) (X k,K t ) 1≤k≤K of the stochastic process X t are simulated, driven by independent Brownian motions (W k t ) k≥1 , and the stress tensor is obtained as an empirical mean over these processes: In this context, the discretization method coupling a finite element method and a Monte Carlo technique is called CONNFFESSIT for Calculation Of Non-Newtonian Flow: Finite Elements and Stochastic SImulation Technique (see M. Laso and H.C.Öttinger [79], and Figure 9). We describe below the implementation of this method in a simple geometry.

Numerical analysis of SDEs
Numerical analysis in fluid mechanics  To give a flavour of the numerical difficulties involved when solving coupled problems of the type (3.3) or (3.6), we momentarily consider the simple situation of a start-up Couette flow (see Figure 10). The fluid flows between two parallel planes. Such a model is typically obtained considering a flow in a rheometer, between two cylinders, and taking the limit of large radii for both the inner and the outer cylinders (see Figure 10). At initial time, the fluid is at rest. The lower plane (y = 0, modelling the inner cylinder of the rheometer) is then shifted with a velocity V (t), which, for simplicity, will be set to a constant value V (sinusoidal velocities may also be applied). On the other hand, the upper plane (y = L, modelling the outer cylinder of the rheometer) is kept fixed. Such a setting is called a start-up flow, and because it is confined between two parallel planes, a Couette flow.
We denote by x and y the horizontal and vertical axes, respectively. The flow is assumed invariant in the direction perpendicular to (x, y).  Figure 10. Schematic representation of a rheometer. On an infinitesimal angular portion, seen from above, the flow is a simple shear flow (Couette flow) confined between two planes with velocity profile (u(t, y), 0, 0).
After appropriate assumptions (the flow is laminar, the velocity writes u = u(t, y) e x , the entropic force is taken Hookean, etc), the coupled system (3.3) simplifies into                      Re ∂u ∂t (t, y) = (1 − ǫ) ∂ 2 u ∂y 2 (t, y) + ∂τ ∂y (t, y), τ (t, y) = ǫ We R 2 P Q ψ(t, y, P, Q) dP dQ, ∂ψ ∂t (t, y, P, Q) = − ∂ ∂P ∂u ∂y (t, y)Q − 1 2We P ψ(t, y, P, Q) Q ψ(t, y, P, Q) + 1 2We where P and Q are the two components of the end-to-end vector r, along the x and y axes respectively. In the above system, τ (t, y) denotes the xy entry of the tensor τ p . Actually, the pressure field, and the other entries of the stress tensor may be then deduced, independently.
We emphasize at this stage the tremendous simplifications that the Couette model allows for. Owing to the simple geometric setting and the fact that the flow is assumed laminar, the divergence-free constraint (3.3) is fulfilled by construction of the velocity field and can be eliminated from the system. In addition, both transport terms (u · ∇)u and (u · ∇)ψ vanish, again because of evident geometrical considerations. This explains the extremely simple form of the equation of conservation of momentum in this context, which indeed reduces to a simple onedimensional heat equation. This set of simplifications is specific to the Couette flow. Substantial difficulties arise otherwise.
The macroscopic equation is discretized with finite elements: P1 finite elements for the velocity and P0 finite elements for the stress tensor.
If the mesoscopic scale is modelled by the Fokker Planck equation, the latter may be discretized using finite difference in time (taking explicit the transport term and implicit the parabolic terms) and a spectral method for the space variable. More precisely, we first introduce the equilibrium solution for the last line of (4.1) when , rewrite the Fokker-Planck equation using ϕ = ψ ψ ∞ as the primary unknown function, and next proceed with the discretization. Because of the specific form of ψ ∞ , the most appropriate Galerkin basis consists of (tensor products of) Hermite polynomials H i , which indeed satisfy 1 √ 2π R H i (P )H j (P ) exp (−P 2 /2) dP = δ ij . Note that the use of such a spectral basis allows to circumvent the practical difficulty related to the fact that the equation is posed on the whole space.
Alternatively to the last two lines of (4.1), the mesoscopic scale may be modelled using the stochastic differential equations where V t and W t are two mutually independent one-dimensional Brownian motions, and the stress is given by τ (t, y) = ǫ We R 2 P Q ψ(t, y, P, Q) dP dQ = ǫ We E(P (t, y)Q(t)).

(4.3)
The system is discretized using a forward Euler scheme in time and a standard Monte-Carlo method for replacing the above expectation value by an empirical mean: for 1 ≤ k ≤ K (number of realizations of the random variables), where V n k and W n k are independent normal random variables, and This discretization is the CONNFFESSIT approach mentioned above, implemented in a simple case (see Figure 11).
A crucial remark (which not only applies to the Couette flow, but also to more general situations) is the following. Since the stress (τ h ) n+1 i is an empirical mean (4.6), it is thus also a random variable. It follows that the discretized macroscopic velocity itself is a random variable. In contrast, at the continuous level, K → ∞ in (4.6), and the stress and the velocity are both deterministic quantities (since the expectation value is a deterministic quantity). Consequently, computing the velocity or the stress using the stochastic approach implies performing a collection of simulations, and averaging out the results. This immediately brings variance Figure 11. The CONNFFESSIT method in a shear flow.
issues into the picture. Appropriate variance reduction methods can be applied.
We refer to C. Le Bris and T. Lelièvre [83] (and more specifically to B. Jourdain, C. Le Bris, and T. Lelièvre [67]) for more details. The Matlab codes corresponding to the above description are available at http://hal.inria.fr/inria-00165171 The numerical approach to treat other types of entropic forces (like the FENE force) and other geometries of flows follows the same line, but the extension can be quite involved, especially for the Fokker-Planck approach.
It is evident from the above discussion that the work load implied by such a coupling may be overwhelming in many practical, especially three-dimensional, situations. Thus the micro-macro simulations are, to date, limited in applications. They serve as a backroom strategy to validate or derive appropriate constitutive laws. They may also be employed on a limited portion of the computational domain (typically in a layer close to the boundaries where non-Newtonian effects are expectedly important; see for example A. Ern and T. Lelièvre [52] for results in this direction).

4.2.
Solids. The method analyzed in Section 3.2 is a model example for more advanced methods, such as the QuasiContinuum Method (QCM). In its initial version [115,116] by E.B. Tadmor, M. Ortiz and R. Phillips, the approach first considers the continuum scale, with a standard continuum mechanics model, discretized by a finite element method. The multiscale feature of the method appears when the elastic energy of an element is computed. Depending on some criteria, some elements are declared to be too heterogeneously strained for a macroscopic description to be valid. They are considered henceforth as a set of discrete particles. The associated energy is then computed according to an underlying atomistic model. Otherwise, for an element that is smoothly deformed, the standard original continuum mechanics model is used to compute the energy.
In the second, later, version of the method [112] by V.B. Shenoy, R. Miller, E.B. Tadmor, D. Rodney, R. Phillips and M. Ortiz, that we describe below, the opposite viewpoint is adopted. The starting point is a multibody atomistic energy, sum of the energies E i (ϕ) of each individual atom i when the current configuration of the atomistic system is defined by ϕ. We define the equilibrium configuration as the solution to inf {E micro (ϕ); ϕ ∈ A} , (4.8) with A = ϕ ∈ R dN , ϕ satisfies some boundary conditions , where d is the space dimension. In practice, the system under consideration is composed of an extremely large number N of atoms. Hence, the evaluation of (4.7), for a given ϕ, is already a challenging task. Furthermore, the variational problem (4.8) is set in a high-dimensional space.
To drastically diminish the number N of degrees of freedom, N r atoms are selected, with N r ≪ N . They are called the representative atoms, abbreviated repatoms. Let i α , 1 ≤ α ≤ N r , denote their indices. Their current positions ϕ iα Nr α=1 are the only remaining degrees of freedom of the reduced system. The positions of the N − N r non-representative atoms are obtained by interpolation. The idea of interpolation is related to the Cauchy-Born rule; see the work [60] by G. Friesecke and F. Theil for a seminal paper on the validity of this assumption for a spring lattice system. More precisely, a mesh is built upon the repatoms in the reference configuration. Let ϕ i 0 be the reference position of atom i, and S α (x) be the piecewise affine function associated with the node α (for simplicity, we henceforth consider a P1 finite element method). In a one-dimensional setting, we thus have The position of any atom i in the deformed configuration is obtained from the current positions of the repatoms using the interpolation formula Otherwise stated, a Galerkin approximation inf {E micro (ϕ); ϕ ∈ A and satisfies (4.9)} (4.10) of (4.8) is performed. We now turn to the practical evaluation of the energy (4.7), and explain how to handle the large number of terms it involves. Assume that the energy of atom i reads for some interaction potential V with some cutoff radius r cut , and assume that the reference configuration is a periodic lattice. Consider an atom i that only interacts with the atoms j of the same finite element. The key point is to observe that its energy E i (ϕ) actually does not depend on i. Indeed, when atoms i and j belong to the same finite element, on which S α is an affine function, we infer from (4.9) that where g ℓ α is the gradient of S α in the finite element ℓ, which is a constant vector. Hence, ϕ j − ϕ i is a function of ϕ j 0 − ϕ i 0 , parameterized by the finite element ℓ and the current positions ϕ iα Nr α=1 of the repatoms: Inserting the above relation in (4.11), we observe that E i (ϕ) is a sum of terms that only depend on ϕ j 0 − ϕ i 0 , that we write The reference configuration being a perfect lattice, we obtain that the above sum actually does not depend on i. Hence, all atoms i in a finite element ℓ that only interact with atoms j in the same finite element share the same energy, which we may thus denote by E ℓ , and which only depends on the current positions of the repatoms (actually, only on those related to the vertices of the finite element): E i (ϕ) = E ℓ ϕ iα Nr α=1 . By choice, the exact same expression E ℓ ϕ iα Nr α=1 is taken as an approximation of the energy of atoms that interact with atoms belonging to different finite elements. We hence approximate the energy (4.7) using where n ℓ is the number of atoms included in the finite element ℓ. Note that there are much fewer terms in the sum (4.12) than in the sum (4.7). We next approximate (4.10) by the problem inf E micro (ϕ); ϕ ∈ A and satisfies (4.9) , (4.13) which can be solved in practice, since it is posed on a space of moderate dimension (there are dN r ≪ dN degrees of freedom), and it involves energies E micro (ϕ) that can be computed practically. So, in its second version, the QuasiContinuum method somewhat consists in an efficient quadrature rule to compute (4.7). This second formulation leads to similar equations as the first version, presented in Section 3.2, in the case of a discretization using Lagrangian P1 finite elements. Indeed, in that case, the degrees of freedom in both the first and second versions are the current positions of the atoms at the mesh vertices. In addition, using the fact that ∇ϕ is constant over each triangle T ℓ , the first term of (3.18) reads where ∪ ℓ T ℓ is the triangulation of D reg . The volume ratio |T ℓ |/|D| is equal to n ℓ /N , where n ℓ is the number of atoms included in the finite element T ℓ , and N is the total number of atoms in the system. Besides, by definition, W (∇ϕ |T ℓ ) is the energy of an atom in a lattice deformed by the linear map ∇ϕ |T ℓ , which is exactly equal to the energy E ℓ ϕ iα Nr α=1 . Hence, the first term of (3.18) reads Up to a multiplicative constant, we recover the expression (4.12) in the subdomain D reg . The QuasiContinuum method has been applied in a number of practical situations (see e.g. the works [99,111,117] [6] by M. Anitescu, D. Negrut, P. Zapol and A. El-Azab for an application of a similar idea to another context, and the paper [94] by M. Luskin and C. Ortner for some discussion of the consistency and accuracy of several quadrature rules that have been proposed for the practical evaluation of (4.7), along the above lines.

5.
Mathematical and numerical analysis. The questions raised by the above multiscale models in terms of mathematical analysis and numerical analysis are diverse in nature, because the mathematical objects involved are themselves diverse.
Before proceeding with our overview of these questions, let us mention that some knowledge on the mathematical nature of the (single-scale) macroscopic models is a prerequisite. The situation cannot be expected to be simpler than, say, models like the Oldroyd B model (2.5) for complex fluids, or the nonlinear elasticity models (2.9) for solid materials. There is a two-fold reason for this. Firstly, the multiscale models are indeed, for some simple particular cases, equivalent to some single-scale models. We have mentioned the Hookean model equivalent to the Oldroyd B model, or the nearest neighbour harmonic interaction potential that leads to linearized elasticity (2.10). Secondly, multiscale systems at least involve macroscopic equations, either because they couple general equations of conservation (like the first two lines of (3.3)) with equations at finer scales, or because they use in practice two different sets of equations in different regions of the computational domain, and one of such set of equations is often purely macroscopic. Without even speaking of the coupling issues, considering separately each region, or each equation, is necessary.
As a reference, and because, as we have seen, this is implicitly or explicitly embedded into the analysis of the multiscale problems, we each time begin our discussion below with briefly recalling the mathematical and numerical challenges for the single-scale equations.

Fluids.
5.1.1. Mathematical analysis. Systems like (2.5) modelling purely macroscopically non-Newtonian fluids typically include the Navier-Stokes equations, with the additional term div τ p in the right-hand side. The equation on τ p is essentially a transport equation and, formally, τ p has at best the regularity of ∇u (this formal observation is important for the choice of appropriate functional spaces for the mathematical setting, and of the discretization spaces for numerical methods). The term div τ p in the right-hand side of the momentum equation is not likely to bring more regularity on u. It is thus expected that the study of these coupled systems contains at least the well-known difficulties of the Navier-Stokes equations. Recall that for the three-dimensional Navier-Stokes equations, it is known that globalin-time weak solutions exist but the regularity, and thus the uniqueness, of such solutions for appropriate data is only known locally in time.
Besides the difficulties already contained in the Navier-Stokes equations (which essentially originate from the Navier term u·∇u), the coupling with the equation on τ p raises additional problems. First, these equations contain a transport term u · ∇τ p without any diffusion term in the space variable. They are hyperbolic in nature. The regularity on the velocity u is typically not sufficient to rigorously treat this transport term by a characteristic method. Moreover, these equations involve a nonlinear multiplicative term ∇u τ p . Finally, for the most sophisticated models, the equations defining τ p generally contain additional non-linearities. These difficulties of course limit the state-of-the-art mathematical well-posedness analysis to mainly local-in-time existence and uniqueness results. They also have many implications on the numerical methods (in terms of choice of the discretization spaces, stability of the numerical schemes, ...).
Example of results for such macroscopic models that may be found in the literature are contributions by M. Renardy in [107] (local-in-time existence for an abstract model that covers the specific Oldroyd-B case), by C. Guillopé and J.C. Saut [61,62] (existence results for less regular solutions for non-zero viscosity of the solvent), by E. Fernandez-Cara, F. Guillen and R.R. Ortega [54] and references therein (local-intime well-posedness in Sobolev spaces), by F.-H. Lin, C. Liu and P.W. Zhang [87] (local-in-time existence and uniqueness results and global-in-time existence and uniqueness results for small data), by P.-L. Lions and N. Masmoudi [90] (the only global-in-time existence result we are aware of, for a specific Oldroyd-like model with the corotational convective derivative on the stress tensor rather than the upper convected derivative). For reviews on the mathematical analysis of macroscopic models citing many more contributions, we refer to M. Renardy [109], or E. Fernandez-Cara, F. Guillen and R.R. Ortega [54].
We now turn to the multiscale models. The analysis of micro-macro models for polymeric fluids has begun with an early work by M. Renardy [108], where a micromacro model in its Fokker-Planck formulation is studied. There is now a growing literature on such models, presumably because they are prototypical of a broad class of multiscale models, where some parameters inserted in the macroscopic equations are computed using models at finer scales.
The difficulties present for the purely macroscopic models discussed above are also present mutatis mutandis in the multiscale models. They are related to the transport terms (in addition to u · ∇u, we now have u · ∇X t and u · ∇ψ), the nonlinear terms either coming from the coupling between variables (in addition to the variables (u, p) and τ p , we now have ∇uX t and div r (∇u r ψ)), or inherently contained in the equations defining τ p (due to the non-linear entropic force F ).
Several existence and uniqueness results for the coupled system involving the Fokker-Planck equation have been obtained. For some local existence and uniqueness results, we refer to M. Renardy [108], T. Li, H. Zhang and P.W. Zhang [85], H. Zhang and P.W. Zhang [122], N. Masmoudi [96]. The latter author has also shown a global in time existence result for initial data close to equilibrium (see also F.-H. Lin, C. Liu and P.W. Zhang [88] for a prior, more restricted result). Moreover, N. Masmoudi has also recently obtained global existence of weak solutions to the FENE model in [97].
Global existence results have also been obtained for some closely related problems: • Existence results for a regularized version: In a series of papers [12,13,14], J.W. Barrett and E. Süli obtain global existence results using space regularized versions. For example, the velocity u in the Fokker Planck equation is replaced by a smoothed velocity, and the same smoothing operator is used on the stress tensor τ p in the right-hand side of the momentum equations. Alternatively, a small diffusion term (with respect to the space variable) is introduced in the Fokker-Planck equation. See also L. Zhang, H. Zhang and P.W. Zhang [123]. . More precisely, in [91], a global-in-time existence result of weak solutions is obtained in dimension 2 and 3, while in [96], it is proved that in dimension 2, there exists a unique global-in-time strong solution. A related recent result by F.-H. Lin, P. Zhang and Z. Zhang is [89].
We would like also to mention the related works [37,38,40] by P. Constantin, C. Fefferman, N. Masmoudi and E.S. Titi, on existence results for coupled Navier-Stokes Fokker-Planck micro-macro models.
We now turn to system (3.6), which is quite difficult to study in full generality. Two simplifications of this general setting are usually considered, at least for preliminary arguments: homogeneous flows and shear flows. In homogeneous flows, ∇u does not depend on the space variable, and therefore X t (and thus τ p ) does not depend on the space variable either. It greatly simplifies the mathematical situation. For illustration purposes, we prefer to consider here a shear flow, studied in particular in M. Laso and H.C.Öttinger [79], J.C. Bonvin and M. Picasso [27], C. Guillopé and J.C. Saut [62], W. E, T. Li and P.W. Zhang [50]. The equations have already been introduced in (4.1) and (4.2). For a general force (instead of an Hookean force), the system reads: where (P t (y), Q t (y)) are the two components of the stochastic process X t (y), and (F P , F Q ) the two components of the force F . A standard formal manipulation on this system, involving some elements of stochastic calculus, yields the following a priori estimate: where Π is the potential of the force F = ∇Π. Notice that the right-hand side of this equality is typically a positive term (recall that, in practice, the potential Π is convex). The situation is different from the usual a priori estimates for, say, the Navier-Stokes equations where the right-hand side is zero. Here, some energy is brought to the system by the finer scales. On finite time intervals, this is however not a difficulty for the mathematical analysis. A typical result is, in the case of Hookean dumbbells in a shear flow, the globalin-time existence and uniqueness result proved for the first time in B. Jourdain, C. Le Bris and T. Lelièvre [70]. The solution (u, X t ) on the interval [0, T ] satisfies the estimate This setting (Hookean dumbbell in a shear flow) is actually extremely specific. A global-in-time existence and uniqueness result is obtained since the coupling term ∇uX t of the original problem simplifies to ∂u ∂y Q t in (5.1), where Q t is known independently of (u, P t ). In other words, this coupling term is, serendipitously, no more nonlinear.
For the FENE force, two new difficulties have to be addressed: first, the SDE contains an explosive drift term and second, even in a shear flow, the coupling term ∇uX t is genuinely nonlinear. Advanced questions of stochastic analysis arise, the model being really multimaths. We refer to B. Jourdain and T. Lelièvre [69] and to B. Jourdain, C. Le Bris and T. Lelièvre [71] for the first studies in this setting and in particular the proof of a local-in-time existence and uniqueness result. For a similar result in a more general setting (3-dimensional flow) and forces with polynomial growth, we refer to W. E, T. Li and P.W. Zhang [51]. The authors prove a local-in-time existence and uniqueness result in high Sobolev spaces. We also refer to A. Bonito, Ph. Clément and M. Picasso [26] for existence results for Hookean dumbbells, neglecting the advection terms.
More generally, it is important to notice that, when the velocity field is not sufficiently regular, and similarly to the situation seen for the Fokker-Planck setting above, it is difficult to give a sense to the transport term in the SDE (which is actually a Stochastic Partial Differential Equation). We refer to C. Le Bris and P.-L. Lions [80,81] for works in this direction.
On the other hand, the presence of an energy source, the right-hand side of (5.2), affects the analysis of the long-time behaviour, like questions related to return to equilibrium. For such questions, the appropriate notion to introduce is that of free energy rather than energy. Assume zero Dirichlet boundary conditions on the velocity u. The expected stationary state (equilibrium) is u(∞, x) = 0, ψ(∞, x, X) = ψ eq (X) ∝ exp(−Π(X)).
The free energy sum of the kinetic energy plus the relative entropy with respect to the equilibrium ψ eq , can be shown to satisfy: 3) Comparing with (5.2), we observe that the introduction of the entropy allows to eliminate the positive right-hand side. Standard techniques of kinetic theory (like Logarithmic Sobolev inequalities, see [5]) allow then to conclude that, under appropriate conditions, the fluid returns to equilibrium after perturbations. We refer to B. Jourdain, C. Le Bris, T. Lelièvre and F. Otto [68], or A. Arnold, J.A. Carrillo and C. Manzini [7]. The multimaths character of the setting is evident.
We would like to also mention that these estimates on the micro-macro system can be used as a guideline to derive new estimates on related macro-macro models (see D. Hu and T. Lelièvre [64]) and also to derive new approximation schemes (see S. Boyaval, T. Lelièvre and C. Mangoubi [28]). This is an interesting (perhaps general) byproduct of mathematical studies of multiscale systems to actually contribute to better understand and approximate the associated purely macroscopic models.

Numerical analysis.
Of course, the difficulties raised by the discretization of the models are, as always, reminiscent of the difficulties of the mathematical analysis. Here again, as mentioned above, the treatment of the multiscale problem necessarily requires a good knowledge of the treatment of the purely macroscopic model. An overview of the numerical difficulties encountered when simulating purely macroscopic models for non Newtonian fluids may be found in R. Keunings [74], F.P.T. Baaijens [8], and R. Owens and T. Phillips [105]. As for multiscale problems, we will only, for brevity, address here the stochastic formulation of the equations, that is, system (3.6). A typical result of numerical analysis, proved in B. Jourdain, C. Le Bris and T. Lelièvre [70] and W. E, T. Li and P.W. Zhang [50], deals with Hookean dumbbells in a shear flow. The error estimate for the discretization approach reads: The main difficulties for the proof (see Figures 9, 10 and 11 for the context) originate from the following facts: • The velocity u n h is a random variable. The energy estimate at the continuous level cannot be directly translated into an energy estimate at the discrete level (which in turn would yield the stability of the scheme).
• The end-to-end vectors (P n i,k , Q n k ) 1≤k≤K are coupled random variables (even though the driving Brownian motions (V n k , W n k ) 1≤k≤K are independent). • The stability of the numerical scheme requires an almost sure bound on Q n k . For an extension of these results to a more general geometry and discretization by a finite difference scheme, we refer to T. Li and P.W. Zhang [86]. A convergence result in space and time may be found in A. Bonito, Ph. Clément and M. Picasso [25]. Many other studies now exist in the literature.

Solids.
As announced at the beginning of Section 5, before studying multiscale models of solids, we first briefly review nonlinear elasticity models. A very nice review by J.M. Ball on this type of mathematical problems may be found in [9]. The (important) difficulties present in these models are also present in multiscale models. The coupling with discrete models brings additional problems which we study in the sequel.
We return to the study of problem (2.9): it defines the deformation ϕ of a solid subject to the body force f and the surface force g. The set A may be for instance of the form: where M is a given matrix, and p > 1 should be related to W (see (5.5) below). The appropriate mathematical notion involved in the study of (2.9)-(5.4) is quasiconvexity. We refer to the papers [11] [11]. Assume that W is quasiconvex and bounded from below and that f ∈ L 1 (D). Assume in addition that W satisfies the growth condition: and that A is defined by (5.4) with the same p as in (5.5). Then, problem (2.9) has a minimizer ϕ ∈ W 1,p (D). This result however does not apply to most physically relevant situations. Indeed, since matter cannot interpenetrate itself, a natural assumption is that W (M ) −→ +∞ as det(M ) −→ 0 + . This clearly contradicts (5.5).
In the same spirit, an expression of the form (3.12) implies some invariances that prevent W from being quasiconvex, unless it is a function of the determinant of ∇ϕ(x). We refer to I. Fonseca [56] for the proof, and to M. Chipot and D. Kinderlehrer [35] and I. Fonseca [57] for related results. A good simple model to understand the underlying difficulties is the so-called one-well problem: is the group of rotations in R 3 . Obviously, W is nonnegative and vanishes only on SO(3). This functional satisfies the frame invariance. According to D. Kinderlehrer [76] (see also S. Müller [101]), we have, if f = g = 0: • If ϕ is a minimizer of (2.9), with A defined by (5.4), then ϕ(x) = M x almost everywhere in D. • Any minimizing sequence (ϕ j ) j∈N of (2.9) satisfies ∇ϕ j −→ M in measure.
The second assertion is delicate. The gradient of a minimizing sequence is easily shown to converge almost everywhere to a matrix of the form RM , where R ∈ SO(3) may depend on x ∈ D. The point is to show that this rotation R is actually constant. The result is closely related to fundamental results called rigidity lemmas, contained in [65,66] by F. John, [59] by G. Friesecke, R.D. James and S. Müller, and [110] by Y.G. Reshetnyak. More involved is the case when W in (2.9) has two minima (up to rotational invariance). More precisely, assume that W ≥ 0, and that there exists two matrices A = B such that (3)) . (5.6) We then have the following property (see J.M. Ball and R.D. James [10]): if ϕ ∈ W 1,∞ (D) is such that for all x ∈ D, ∇ϕ(x) ∈ {A, B}, then: where a ∈ R 3 and n ∈ R 3 , then there exists a Lipschitz function h such that h ′ ∈ {0, 1} and ϕ(x) = Ax+ah(x·n)+b, for some b ∈ R 3 . Hence, in the first case, i.e. if A and B are not rank-one connected, then ∇ϕ must be constant, and the behavior is similar to the one-well problem case. On the contrary, in the second case, i.e. if A and B are rank-one connected, minimizing sequences may behave differently. Indeed, if the boundary condition M satisfies M = λA + (1 − λ)B, for some λ ∈ (0, 1), then it is possible to construct a minimizing sequence (ϕ j ) j≥0 of (2.9) such that ∇ϕ j = A and B alternately on strips of width λ/j and (1 − λ)/j, respectively. A suitable tool to study such oscillating minimizing sequences is gradient Young measures. We refer to S. Müller [101], P. Pedregal [106] and L.C. Young [120].
Let us now turn to the problem of minimizing the coupled energy (3.18). We perform the analysis in a one-dimensional setting, as in X. Blanc, C. Le Bris and F. Legoll [20,19]. The generalization to higher dimensions is, to date, unclear. As we will see below, serious difficulties already arise in the one-dimensional case. We assume that the interaction potential V is the Lennard-Jones potential V (r) = 1 r 12 − 2 r 6 , and we only account for nearest neighbor interactions.
We assume that D = (0, 1), that the equilibrium configuration corresponds to a constant distance between the atoms: that D = D reg ∪ D sing is a nontrivial partition of D, and, for simplicity, that We denote by B N = ⌊N b⌋ ∈ N the largest integer that is lower than N b. The coupled energy (3.18) reads 8) and the problem we consider is inf ϕ∈A E c (ϕ), (5.9) where the set A reads A = ϕ; ϕ |Dreg ∈ W 1,1 (D reg ), ϕ |Dsing is the discrete set of variables ϕ i N , In principle, D reg and D sing should be chosen such that the expected singularities of the deformation all belong to D sing . However, as established by X. Blanc, C. Le Bris and F. Legoll [20], if N is large enough and if ϕ(1) is set to a > 1, the minimizers of (5.9) are of the form: (5.12) with I ⊂ N and i∈Iṽ i = a − 1 (see Figure 12). Here, H denotes the Heaviside function: Figure 12. A typical minimizer of (5.9), when N is large enough and ϕ(1) is set to a > 1.
Equations (5.11)-(5.12) indicate that the material breaks if a > 1, and that the fracture systematically occurs in the regular zone D reg , that is in the zone where ϕ was expected to be smooth. Loosely speaking (see the papers [20,82,21] by X. Blanc, C. Le Bris, F. Legoll and P.-L. Lions for more details), one easily constructs test-functions ϕ 1 and ϕ 2 such that ϕ 1 has a fracture in D sing and ϕ 2 has a fracture in D reg , and E(ϕ 1 ) = E(ϕ 2 ) + 1 N + o 1 N > E(ϕ 2 ). Note however that the difference of energy is of order 1 N , which is tiny, since N is the typical number of atoms in the system under consideration. In addition, the above-mentioned argument assumes an exact evaluation of the energies. In practice, the energy of the continuous model is evaluated using a discretization method (say a Finite Element Method for simplicity), which has a finite mesh size H. Hence, As explained by X. Blanc, C. Le Bris, F. Legoll and P.-L. Lions [21,82], it is easily seen that this approximation introduces an error term of order 1/N H ≫ 1/N , which turns out to counterbalance the above-mentioned effect. It follows that, although the theory establishes the fracture always occurs in the continuous zone D reg , the discretized approximation of the continuous energy shows a different behavior. A numerical artifact somehow saves the situation. However, this indicates that one should be careful with the method: it might break down if the mesh is too much refined. Let us also point out that in many situations, global minimizers are not physically relevant. In such a case, the above phenomenon does not take place. In any event, the above discussion shows that, when coupling two different types of mathematical models (here, discrete and continuous), new difficulties arise. Another example of such difficulties is the existence of ghost forces, to which we now turn.
The notion has been first introduced and discussed by V.B. Shenoy, R. Miller, E.B. Tadmor, D. Rodney, R. Phillips and M. Ortiz in [112], see also [42,44] by M. Dobson and M. Luskin, [82] by F. Legoll and [113] by T. Shimokawa, J. Mortensen, J. Shiotz and K. Jacobsen. The phenomenon is not present in the oversimplified case of one dimensional nearest neighbor interactions. The simplest case in which it occurs is a one-dimensional system with second nearest neighbor interactions and no external force. The atomistic energy reads where V 1 and V 2 are the nearest neighbor and the second nearest neighbor interaction potentials, respectively. The energy (5.13) is a sum of pair-wise interactions. As a consequence, when the configuration is homogeneously deformed, e.g. ϕ i = a i + b for any i, the force on any atom that does not interact with the boundaries (namely, any atom i with 2 ≤ i ≤ N − 2) vanishes. We now consider the partition D = D reg ∪ D sing , with D sing defined by (5.7), and the associated coupled energy which is the generalization of (5.8) to the case of second nearest neighbor interactions. Note that atom B N − 1 has two discrete neighbors on its left, but only one on its right. A simple computation shows that the derivative of E NNN c , evaluated on a homogeneously deformed configuration, does not vanish for this atom. Some spurious force, the so-called ghost force, appears at atom B N − 1, due to the coupling of a discrete model with a continuous one.
As mentioned by F. Legoll in [82], one can either work with the coupled energy (5.14), hoping that the effect of ghost forces is not too important, or correct these forces by various means (see the works by M. Dobson, M. Luskin and C. Ortner [44,43,45,47] and [42], respectively). The former option is the energy based formulation, whereas the latter option is the force based formulation (see M. Dobson, M. Luskin and C. Ortner [42,46,48], R. Miller, E.B. Tadmor [98] and V.B. Shenoy, R. Miller, E.B. Tadmor, R. Phillips and M. Ortiz [111]). The latter consists in approximating the atomistic forces (rather than the energy) by a coupled formulation. In that setting, by construction, the forces cancel for any homogeneously deformed configuration (up to forces at the boundaries of the computational domain, that were already present in the reference model). It turns out that the corresponding set of equations is not the Euler-Lagrange equation of a variational problem (there is no underlying energy minimization problem), although the whole strategy is supposed to approximate a problem originally posed as a variational problem (see e.g. (4.8)). Specific techniques are to be used to address this case, both for the standpoint of mathematical analysis and numerical algorithms.
The problem of ghost forces is similar to the problem of defining boundary or interface conditions for multiphysics models. Such conditions need to be appropriately defined in order to avoid spurious effects (like ghost forces here). Note also that the above discussion only considered static problems. Similar problems are encountered in the dynamics. This is reminiscent of well-known problems in, for instance, wave propagations, where dedicated methods such as perfectly matching layers (see J.-P. Bérenger [16], or F.L. Teixeira and W.C. Chew [118] for instance) and/or transparent boundary conditions are needed (see Y. Achdou, C. Sabot, and N. Tchou [1], U.E. Aladl, A.S. Deakin, and H. Rasmussen [2], J.B. Keller and D. Givoli [73], or D.P. Nicholls and N. Nigam [102]).
The above discussion certainly demonstrates that appropriately setting problem (5.9) is delicate. Although the above study is only carried out in a one-dimensional geometry, there is no reason to think that similar difficulties cannot be expected in more general situations. The simple mathematical observation performed above, along with the absence of a complete mathematical analysis of the situation (in particular regarding the systematic study for choosing D reg ), indicate that all results obtained with atomistic-to-continuum coupling methods must be taken with great care. The definite practical success of numerical approaches involving hybrid problems of the type (5.9) should motivate further mathematical efforts. The state of the mathematical understanding is certainly lagging behind the success of the numerical simulations. 6. Trends for multiscale modeling. Polymeric fluids modelling and crystalline solids modelling, which have been jointly examined here, examplify a general trend for multiscale modelling: the inclusion of an increasing number of mathematical objects often different in nature, in order to account for all the phenomena at play at the various scales.
Schematically, traditional single-scale modelling is based on a description of the form Dτ p Dt = G(τ p , u). (6. 2) The structure of system (6.2) is a common denominator to many multiscale models in many contexts. A global macroscopic equation is coupled with a local (microscopic) equation, via an averaging formula. For instance, the reader familiar with homogenization theory for materials formally recognizes in (6.2) the homogenized equation, the value of the homogenized tensor, and the corrector equation, respectively. On the numerical front, it is also a structure shared with multiscale algorithmic approaches: a global coarse solver coupled to a local fine one using an averaging process (think of the Godunov scheme for solving the Riemann problem in computational fluid dynamics).
Already the consideration of systems of the type (6.2) requires several mathematical notions different in nature. But, lately, an additional ingredient appeared. A particular feature of the most recent multiscale approaches is the consideration of random quantities. Randomness may be caused by the need to describe in a concise manner many possible behaviours at a given scale. This is the case with the statistical description seen above of the assembly of thousands of polymeric chains floating in the fluid considered. It may be also related to the wish to get closer to real materials, as is the case in the modelling of solids. Periodicity for solid is often an idealization. In real life, solid materials consist of different grains, are polycrystals, contain defects, etc. Randomness is a way to appropriately account for these flaws in the material.
Understanding mathematically this new category of multiscale systems, and improving the corresponding numerical approaches, therefore require an appropriate background, at the intersection of many mathematical domains.