# The Nanometric and Micrometric Scales of the Structure and Mechanics of Materials Revisited: An Introduction to the Challenge of Fully Deterministic Numerical Descriptions

Amine Ammar, Francisco Chinesta, Pierre Joyot

► To cite this version:

Amine Ammar, Francisco Chinesta, Pierre Joyot. The Nanometric and Micrometric Scales of the Structure and Mechanics of Materials Revisited: An Introduction to the Challenge of Fully Deterministic Numerical Descriptions. International Journal for Multiscale Computational Engineering, 2008, 6 (3), pp.191-213. 10.1615/IntJMultCompEng.v6.i3.20. hal-01007372

HAL Id: hal-01007372

https://hal.science/hal-01007372

Submitted on 2 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# The Nanometric and Micrometric Scales of the Structure and Mechanics of Materials Revisited: An Introduction to the Challenges of Fully Deterministic Numerical Descriptions

A. Ammar*, F. Chinesta†, P. Joyot‡

* Laboratoire de Rhéologie, 1301 Rue de la Piscine,
BP 53 Domaine Universitaire, F-38041 Grenoble cedex 9, France

amine.Ammar@ujf-grenoble.fr

† LMSP UMR 8106 CNRS-ENSAM, 151 Boulevard de l'Hôpital, F-75013 Paris, France

francisco.chinesta@paris.ensam.fr

‡ LIPSI, ESTIA technopole izarbel, F-64210 Bidart, France

p.joyot@estia.fr

## Abstract

Nano-science and nano-technology as well as the fine modelling of the structure and mechanics of materials from the nanometric to the micrometric scales use descriptions ranging from the quantum to the statistical mechanics. This paper revisits the modelling at these scales and points out the main challenges related to the numerical solution of such models that some times are discrete but involves an extremely large number of particles (as it is the case of molecular dynamics simulations or coarse-grained molecular dynamics) and other times are continuous but they are defined in highly multidimensional spaces leading to the well known curse of dimensionality issues. The curse of dimensionality concerned by some of these deterministic models will be emphasized and their numerical implications will be addressed in the second part of this work.
# 1 Introduction

## 1.1 Motivation

The fine description of the mechanics and structure of materials at the micro, nano and sub-nanometric scales introduces some specific challenges related to the impressive number of degrees of freedom required or the highly dimensional spaces in which those models are defined. Despite the fact that spectacular progresses have been accomplished in the context of computational mechanics in the last decade, the treatment of those models, as we describe in the present work, needs further developments.

The brut force approach cannot be considered as a possibility for treating this kind of models. Thus, some specialists as the Nobel Prize R.G. Laughlin, affirmed that no computer existing, or that will ever exist, can break the barriers found in quantum mechanics because it is a catastrophe of dimension [?].

We can understand the catastrophe of dimension by assuming a model defined in a hyper-cube in a space of dimension $D$, $\Omega =]-L, L[^D$. Now, if we define a grid to discretize the model, as it is usually performed in the vast majority of numerical methods (finite differences, finite elements, finite volumes, spectral methods etc.), consisting of $N$ nodes on each direction, the total number of nodes will be $N^D$. If we assume that for example $N \approx 10$ (an extremely coarse description) and $D \approx 80$ (much lower than the usual dimensions required in quantum or statistical mechanics), the number of nodes involved in the discrete model reaches the astronomical value of $10^{80}$ that represents the presumed number of elementary particles in the universe! We shall come back to the analysis of these systems later.

Thus, progresses on this field need further developments on the physical modelling as well as the introduction of new ideas and methods in the context of computational physics. In this work we are exploring different modelling scales, starting from the finest one, the quantum mechanics, to derive molecular dynamics, Brownian dynamics and finally kinetic theory models. The main particularities of such models and the main numerical difficulties associated with their solution will be emphasized.

This work does not pretend to introduce advances in computational physics, but only revisit the main ideas of models defined in the finest scales. Most part of the models are well established and well known in the context of both the computational physics and the computational mechanics community, and consequently all of them will be summarized in the first section that represents a broad introduction.

In the second part of this work we address the curse of dimensionality problematic that some of
these models involve, pointing out the advantages that separated representations could represent for the numerical solution of such models. Finally, some examples will be considered to illustrate the potentiality of such approaches, but also to emphasize the remaining challenges, others than their highly multidimensional character.

## 1.2 From quantum mechanics to statistical mechanics: a walk on the frontier of the simulable world

### 1.2.1 The finest description: the quantum approach

The quantum state of a given electronic distribution could be determined by solving the Schrödinger equation. This equation has been for longtime considered as one of the finest descriptions of the world. However, before focusing on the challenges of its numerical solution, we would like to recall that this equation is not relativistic and then it fails when it is applied to describe heavy atoms. Moreover, the Pauli's principle constraint was introduced in the Schrödinger formalism in an "ad hoc" way and constitutes, as we illustrate in the last section of this work, the main difficulty in its solution.

Some simplificative hypotheses are usually introduced, as for example the Born-Oppenheimer that states that the nuclei can be in first approximation assumed as classical point-like particles, that the state of electrons only depends on the nuclei positions and that the electronic ground state corresponds to the one that minimizes the electronic energy for a nuclear configuration. This equation defines a multidimensionnal problem whose dimension increases linearly with the number of the electrons in the system.

Thus, the knowledge of a quantum system reduces to the determination of the wavefunction $\Psi(x_1, x_2, \dots, x_N, t; X_1, \dots, X_M)$ (that establishes that the electronic wavefunction depends parametrically on the nuclei positions) whose evolution is governed by the Schrödinger equation:

$$i\hbar \frac{\partial \Psi}{\partial t} = -\frac{\hbar^2}{2m_e} \sum_{e=1}^{e=N} \nabla_e^2 \Psi + \sum_{e=1}^{e=N-1} \sum_{e'=e+1}^{e'=N} V_{ee'} \Psi + \sum_{e=1}^{e=N} \sum_{n=1}^{n=M} V_{en} \Psi \quad (1)$$

where $N$ is the number of electrons and $M$ the number of nuclei, the last ones assumed located and fixed at positions $X_j$. Each electron is defined in the whole physical space $x_j \in \mathbb{R}^3$, $i = \sqrt{-1}$, $\hbar$ represents the Planck's constant divided by $2\pi$ and $m_e$ is the electron mass.

The differential operator $\nabla_e^2$ is defined in the conformation space of each particle, i.e.: $\nabla_e^2 \equiv \partial^2/\partial x_e^2 + \partial^2/\partial y_e^2 + \partial^2/\partial z_e^2$. The Coulomb's potentials accounting for the electron-electron and electron-electron potentials.
nuclei interactions write:

$$V_{ee'} = \frac{(q_e)^2}{\|\mathbf{x}_e - \mathbf{x}_{e'}\|} \qquad (2)$$

$$V_{en} = - \frac{q_n q_e}{\|\mathbf{x}_e - \mathbf{X}_n\|} \qquad (3)$$

The electron charge is represented by $q_e$ and the nuclei charge by $q_n = |q_e| \times Z$ (where Z is the atomic number).

The time independent Schrödinger equation (from which one could determine the ground state, perform quantum static computations or accomplishing separated representations of the time-dependent solution) writes:

$$-\frac{\hbar^2}{2m_e} \sum_{e=1}^{e=N} \nabla_e^2 \Psi + \sum_{e=1}^{e=N-1} \sum_{e'=e+1}^{e'=N} V_{ee'} \Psi + \sum_{e=1}^{e=N} \sum_{n=1}^{n=M} V_{en} \Psi = E \Psi \qquad (4)$$

where the ground state corresponds to the eigenfunction $\Psi_0$ associated with the most negative eigenvalue $E_0$.

Several techniques have been proposed for solving this equation. Some of them lie in the direct solution of the (time-independent or time-dependent) Schrödinger equation. Due to the curse of dimensionality its solution was only possible for very reduced populations of electrons.

Other solution strategies are based on the Hartree-Fock (HF) approach and its derived approaches (post-Hartree-Fock methods). The main assumption of this approach lies in the approximation of the joint electronic wavefunction (related to the N electrons) as a product of N 3D-functions (the molecular orbitals) verifying the antisymmetry restriction derived from the Pauli's principle. Thus, the original HF approach consists of writing the joint wavefunction from a single Slater's determinant. The Schödinger equation allows computing the N molecular orbitals after solving the resulting strongly non-linear problem. This technique has been extensively used in quantum chemistry to analyze the structure and behavior of molecules involving a moderate number of electrons. Of course, the HF assumption represents sometimes a too crude approximation which invalidate the derived results.

To circumvent this crude approximation different multi-determinant approaches have been proposed. Interested readers can refer the excellent overview of Cancès et al. [?] as well as the different chapters of the handbook on computational chemistry [?]. The simplest possibility consists in writing the solution as a linear combination of some Slater determinants built by combining n molecular orbitals, with $n > N$. These molecular orbitals are assumed known (e.g. the orbitals related to the hydrogen atom) and the weights are searched to minimize the electronic energy. When the molecu-
lar orbitals are built from the Hartree-Fock solution (by employing the ground state and some excited eigenfunctions) the technique is known as Configuration Interaction method (CI). A more sophisticated technique consists in writing this many-determinants approximation of the solution by using a number of molecular orbitals $n$ (with $n > N$) assumed unknown. Thus, the minimization of the electronic energy leads to compute simultaneously the molecular orbitals as well as the associated coefficients of this many-determinants expansion. Obviously, each one of these unknown molecular orbitals are expressed in an appropriate functional basis (e.g. gaussian functions, ...). This strategy is known as Multi-Configuration Self-Consistent Field (MCSCF).

All the just mentioned strategies (and others like the coupled cluster or the Moller-Plesset perturbation methods) belong to the family of the wavefunction based methods. In any case all these methods can be only used to solve quantum systems composed of a moderate number of electrons. As we confirm later the main difficulty is not in the dimensionality of the space, but in the use of the Slater determinants (needed to account for the Pauli's principle) whose complexity scales on the factorial of the number of electrons, i.e. in $N!$.

The second family of approximation methods, widely used in quantum systems composed of hundreds, thousands and even millions of electrons, are based on the density functional theory (DFT). These models, more than looking for the expression of the wavefunction (with the associated multi-dimensional issue) look for the electronic distribution $\rho(\mathbf{x})$ itself. The main difficulties of this approach are related to the expressions of both the kinetic energy of electrons and the inter-electronic repulsion energy. The second term is usually modelled from the electrostatic self-interaction energy of a charge distribution $\rho(\mathbf{x})$. On the other hand the kinetic energy term is also evaluated in an approximate manner (from the electronic distribution itself in the Thomas-Fermi and related orbital-free DFT models or from a system of $N$ non-interacting electrons –Kohn-Sham models–). Obviously, due to the just referred approximations introduced in the kinetic and inter-electronic interaction energies, a correction term is needed, the so-called exchange-correlation-residual-kinetic energy. However, no exact expression of this correction term exists and then different approximate expressions have been proposed and used. Thus, the validity and accuracy of the computed results will depend on the accuracy of the exchange-correlation term that must be fitted for each system.

The models related to the Thomas-Fermi, less accurate in the practice because the too phenomenological expression of the kinetic energy coming from the reference system of an uniform non-interacting electron gas, allows to consider large multi-electronic systems. In a recent work, Gavini et al. [?] per-
formed multi-million atom simulations by employing the Thomas-Fermi-Weizsacker family of orbital-free kinetic energy functionals. On the other hand, the Kohn-Sham based models are a priori more accurate, but they need the computation of the *N* eigenfunctions related to the *N* lowest eigenvalues of a non-physical atom composed of *N* non-interacting electrons.

Transient solutions are very common in the context of quantum gas dynamics (physics of plasma) but are more infrequent in material science when the structure and properties of molecules or crystals are concerned. For this reason, in what follows, we are focusing on the solution of the time-independent Schrödinger equation which leads to the solution of the associated multidimensional eigenproblem, whose eigenfunction related to the most negative eigenvalue constitutes the ground state of the system.

Quantum chemistry calculations performed in the Born-Oppenheimer setting consist either (i) in solving the geometry optimization problem, that is, to compute the equilibrium molecular configuration (nuclei distribution) that minimizes the energy of the system, finding the most stable molecular configuration that determines numerous properties like for instance infrared spectrum or elastic constants; or (ii) in performing an "ab initio" molecular dynamics simulation, that is, to simulate the time evolution of the molecular structure according to the Newton law of classical mechanics. Molecular dynamics simulations allow to compute various transport properties (thermal conductivity, viscosity, ...) as well as some others non-equilibrium properties.

For more details on the mathematical aspects of these models the interested reader can refer to [?], [?] and [?] and the references therein.

### 1.2.2 From "ab initio" to molecular dynamics

Depending on the choice of the method, on the accuracy required, and on the computer facility available, the *ab initio* methods allow today for the simulations of systems up to ten, one hundred or some million atoms. In time dependent simulations, they are only convenient for small-time simulations, say not more than a picosecond. However, some times larger systems are concerned, and for this purpose one must focus on faster approaches, obviously less accurate. Two possibilities exist: the semi-empirical and the empirical approaches. The semi-empirical approaches speed up the *ab initio* methods by profiting of the information coming from experiments or previous simulations. Empirical methods go on by considering explicitly only the nuclei, by introducing "empirical" potentials leading to the forces acting on the nuclei. Thus, in the stationary setting only the stable configuration is searched, and for this a geometrical optimization (to computed the nuclei equilibrium distribution) is
addressed leading to the so-called molecular mechanics. The transient setting results in the classical molecular dynamics but now the computation is speed up of many orders of magnitude with respect to the molecular dynamics where the potentials are computed at the *ab initio* level.

Thus, if we assume a population of *M* nuclei (of mass *m*<sub>n</sub>) and a two-body potential (many-body potentials are also available), now the Newton's law writes for a generic nuclei *n*:

$$m_n \frac{d^2 \mathbf{X}_n}{dt^2} = \sum_{k=1, k \neq n} \mathbf{F}_k^n, \quad \forall n \in [1, \dots, M] \qquad (5)$$

where $F_k^n$ denotes the force acting on nucleus *n* originated by the presence of nucleus *k*. Obviously these forces can be computed from the gradient of the assumed inter-particles potentials.

Accurate algorithms for integrating these equations exist. The simplectic Verlet's scheme is one of the most used. Molecular dynamics simulations are confronted, despite its conceptual simplicity, with diverse difficulties of different nature:

* The first and most important comes, as previously indicated, from the impossibility of using an "exact" interaction potential derived from quantum mechanics. This situation is particularly delicate when we are dealing with some irregular nuclei distributions as the ones encountered in the neighborhood of defaults in crystals (dislocations, crack tips, etc.), interfaces between different materials or in zones where different kinds of nuclei coexist.

* The second one comes from the units involved in this kind of simulations: the nuclei displacements are in the nanometric scale, the energies are of the order of the electron-volts, the time steps are of the order of $10^{-15}$s. Thus, because of the limits in the computers precision, a change of units is required, which can be easily performed.

* In molecular dynamics the behavior of atoms and molecules is described in the framework of classical mechanics. Thus, the particles energy variations are continuous. The applicability of MD depends on the validity of this fundamental hypothesis. When we consider crystals at low temperature the quantum effects (implying discontinuous energy variations) are preponderant, and in consequence the matter properties at these temperatures cannot be determined by MD simulations. The use of MD is restricted to temperatures higher than the Debye's temperature (for example the Debye's temperature for the Fe $\alpha \approx 460K$). This analysis is in contrast to the vast majority of MD simulations carried out nowadays. In fact, higher is the temperature (kinetic energy) and higher results the velocity of particles, requiring shorter time steps in order
to ensure the stability of the integration scheme. For this reason, nowadays most of the MD
simulations in solid mechanics are carried out at zero degrees Kelvin or at very low temperatures
but, as just pointed out, at these temperatures the validity of the computed MD solutions are
polluted by the non negligible quantum effects.

* The prescription of boundary conditions is another delicate task. If the analysis is restricted to systems with free boundary conditions, then the MD simulation can be carried out without any particular treatment. On the other case we must consider a system large enough to ensure that in the analyzed region the impact of the free surfaces can be neglected. Other possibility lies in the prescription of periodic boundary conditions, where an atom leaving the system for example through the right boundary is re-injected in the domain through the left boundary. Moreover, the particles located in the neighborhood of a boundary are influenced by the ones located in the neighborhood of the opposite boundary. The imposition of other boundary conditions is more delicate from both the numerical and the conceptual points of view. For example, what is the meaning of prescribing a displacement on a boundary? Each situation requires a deep analysis in order to define the best (the most physical) way to prescribe the boundary conditions.

* There are other difficulties related to the transient analysis. We consider a thermal system in equilibrium, i.e. a system in which the velocities follows the Maxwell-Boltzmann distribution. Now, we proceed to heat the system. One possibility lies in increasing suddenly the kinetic energy of each particle. Obviously, even if the resulting velocities define a Maxwell-Boltzmann's distribution, the system remains off equilibrium because the partition between kinetic and potential energies is not the appropriate one. For this reason we must proceed to relax the system that evolves from this initial (non-physical) state to the equilibrium one. Other (more physical) possibility lies in the incorporation of a large enough ambient region around the analyzed system, whose particles are initially in equilibrium at the highest temperature. Now, both regions (the system and the ambient) interact, and the system initiates its heating process that reaches its equilibrium some time latter. The final state of both evolutions is the same, but the time required to reach it depends on the technique used to induce the heating. The first transient is purely numerical whereas the second one is more physical allowing the identification of some transport coefficients (for example the thermal conductivity).

* Finally the CPU time continues to be the main limitation of MD simulations. The strongest
handicap is related to the necessity of considering at each time step and for each particle the
influence of all the others particles. Thus, the integration method seems to scale with the square of
the number of particles. Even if some computing time savings can be introduced in the neighbors
search, the extremely small time steps and the extremely large number of particles required to
describe real scenarios, limit considerably the range of applicability of this kind of simulations,
that has been accepted to be nowadays -in 2008- of the order of a cubic micrometer, even when
the systems are considered at very low, and then non-physical, temperatures (close to zero degrees
Kelvin). We can notice that, despite the impressive advances in the computational availabilities,
the high performance computing and the use of massive parallel computing platforms, the state of
the art does not allow the treatment of macroscopic systems encountered in practical applications
of physics, chemistry and engineering.

The above mentioned difficulties to perform fully molecular dynamics (MD) simulations motivated
the proposal of hybrid techniques that apply MD in the regions where the unknown field varies in a
non-uniform way (molecular dynamics model) and a standard finite element approximation in those
regions where the unknown field variation can be considered as uniform (continuous model). The
main questions concerned by these bridging strategies concern: (i) the kinematics representations in
both models; (ii) the transfer conditions on the MD and continuous models interface and (iii) the
macroscopic constitutive equation employed in the continuous model.

Different alternatives exist, and the construction of such bridges is nowadays one of the most active
topics in computational mechanics. The spurious reflection of the high frequency parts of the waves is
one of the main issues. We would like only to mention three "families" of bridging techniques, giving
some key references in which the interested reader could find other extremely valuable references: (i)
the quasi-continuum method proposed by Tadmor and Ortiz [?]; (ii) the multi-scale method proposed
by Wagner and Liu [?] and (iii) the methods based on the "Arlequin" approach [?] like the one proposed
by Belytschko in [?].

1.2.3 Coarse grained modelling: Brownian dynamics

Some times one is interested in analyzing the behavior of a system composed by a series of microscopic
entities (particles assumed with a null extension) dispersed into another fluid (the solvent). The
kinematic of such particles depends, of course, on their interactions with the solvent particles. A
real molecular dynamics simulation is definitively forbidden (the nowadays MD simulation feasibilities
rarely exceeds the number of particles contained within a cube of one micron of side).

One possibility to reduce the size of the discrete models lies in considering only the particles of interest. The other particles (the ones that constitute the solvent) are not considered explicitly and only their averaged effects are retained in the modelling.

Thus, the motion equation of a particle whose position is described by $\mathbf{x}_i$, is governed by the Langevin's equation:

$$m \frac{d^2 \mathbf{x}_i}{dt^2} = \xi \left( \frac{d\mathbf{x}_i}{dt} - \mathbf{v}_f(\mathbf{x}_i) \right) + \mathbf{F}_i^{\text{ext}}(t); \quad \forall i \qquad (6)$$

where $m$ denotes the particle mass, $\mathbf{x}_i$ the position of particle $i$, $\xi$ the friction coefficient, $\mathbf{v}_f(\mathbf{x}_i)$ the fluid velocity at position $\mathbf{x}_i$ and $\mathbf{F}_i^{\text{ext}}(t)$ all the other forces acting on the particle $i$ (coming from a external potential or from the solvent particles bombardment). We can notice that even if this model doesn't incorporate explicitly the solvent particles population, their effects are taken into account from the drift term $\xi \left( \frac{d\mathbf{x}_i}{dt} - \mathbf{v}_f(\mathbf{x}_i) \right)$ as well as by the impact forces.

In the last expression the drift term is quite standard, however the external forces acting on each particle deserve some additional comment. In what follows and without any detriment of generality we assume that there is not other forces than the one coming from the solvent particles bombardment and that the solvent is macroscopically at rest, i.e. $\mathbf{v}_f = \mathbf{0}$. The random nature of the interaction force is modelled from a statistical distribution function that becomes fully defined as soon as its mean value and its standard deviation are fixed. In that concerns the mean, one expect a null value if the microscopic dynamics is isotropic. Concerning the standard deviation one must proceed within the statistical mechanics framework. In what follows we summarize the main ideas for the derivation of the standard deviation expression.

If we define $B_{\Delta t}$

$$B_{\Delta t} = \int_{0}^{\Delta t} \frac{F^{\text{ext}}(t)}{m} dt = \sum_{k=1}^{k=p} \frac{F^{\text{ext}}(t_k)}{m} \delta t \qquad (7)$$

where we assume that in $\Delta t$ the particle is subjected to $p$ impacts from the solvent particles, with $p \gg 1$. These impacts are modelled by a constant force $F^{\text{ext}}(t_k)$ that applies for a time $\delta t$ ($\Delta t = p\delta t$). By invoking the central limit theorem we conclude that $B_{\Delta t}$ follows a gaussian distribution $N(0, q\Delta t)$ because the number of impacts scales linearly with $\Delta t$.

Now, to compute $q$, one could integrate the Langevin's equation to obtain the equilibrium velocity
distribution W:

$$W\left(\frac{dx}{dt}, t \to \infty\right) = \sqrt{\frac{\xi}{m\pi q}} e^{-\frac{\xi \left(\frac{dx}{dt}\right)^2}{m\pi q}} \quad (8)$$

that must coincide with the Maxwell-Boltzmann one (canonical ensemble), from which we deduce the expression of q:

$$q = \frac{2\xi k_b T}{m^2} \qquad (9)$$

where $k_b$ is the Boltzmann constant and T the absolute temperature.

Thus, the Langevin's equation is fully defined, by writing:

$$B_{\Delta t} = \int_{0}^{\Delta t} \frac{F^{ext}(t)}{m} dt = N \left( 0, 2 \frac{\xi k_b T}{m^2} \Delta t \right) \qquad (10)$$

expression that runs also in presence of other forces like the ones coming from a gradient of a potential or when $v_f \neq 0$.

Thus, one could track the movement of each particle $x_i$ by considering a standard drift term and a random force whose distribution is perfectly known. This stochastic approach has been traditionally also used for solving deterministic models as the advection-diffusion one, because one could compute some moments of the resulting distribution by tracking a moderate population of particles instead to discretize the deterministic counterpart of the advection-diffusion model by using one of the standard mesh-based discretization techniques that could involve in the case of 3D models an excessive number of degrees of freedom.

To illustrate this procedure we consider the simplest form of the advection-diffusion equation:

$$\frac{\partial C}{\partial t} + \mathbf{v} \cdot \nabla C = D \Delta C \qquad (11)$$

where $C = C(\mathbf{x}, t)$ is the concentration field and $D$ is the so called diffusion coefficient.

Obviously, this simple parabolic equation could be solved by using any standard technique (finite differences, finite elements, spectral methods, finite volumes, the method of particles, ...), but in what follows we are solving it using a stochastic approach. For this purpose we assume the initial condition represented by N-Dirac masses:
$$C(\mathbf{x}, t=0) \approx C^0(\mathbf{x}) = \sum_{i=1}^{i=N} c_i \delta(\mathbf{x} - \mathbf{x}_i(t=0)) \quad (12)$$

Here, we are not discussing the choice of the coefficients $c_i$ and the locations $\mathbf{x}_i$. Several possibilities exist to perform a choice trying to represent, as precisely as possible, the initial concentration distribution. Most of them proceed by regularizing the Dirac distribution and then enforcing the minimization of $\|C(\mathbf{x}, t=0) - C^0(\mathbf{x})\|$. From now on, $c_i$ and $\mathbf{x}_i(t=0)$ are assumed known. Moreover, as the considered advection-diffusion equation does not contain source terms, the weights $c_i$ remain unchanged during the evolution of the pseudo-particles positions $\mathbf{x}_i(t)$.

Now, the simplest explicit integration algorithm proceeds by updating the particle position considering both the deterministic and the random contributions:

$$\mathbf{x}_i^{n+1} = \mathbf{x}_i^n + \mathbf{v}(\mathbf{x}_i^n)\Delta t + \mathcal{N}(0, 2D\Delta t)\mathbf{u}, \quad \forall i \quad (13)$$

where $\mathbf{x}_i^{n+1} \equiv \mathbf{x}_i(t = (n + 1) \times \Delta t)$ and $\mathbf{u}$ is a unit random vector.

Now the distribution moments can be easily computed, and the concentration field could be reconstructed by employing some appropriate smoothing.

It is usual to find this kind of advection-diffusion equations in many branches of science and engineering. In particular they are encountered when one models macromolecular materials within the statistical mechanics framework, as we describe later. Despite its intrinsic simplicity, some times they arise defined in highly dimensional spaces including the physical and the conformation coordinates. To avoid the curse of dimensionality drawback characteristic of mesh-based techniques, different authors proposed the use of stochastic techniques exploiting the equivalence between the so-called Fokker-Planck equation and the Ito's stochastic equation [?], similar to the just described equivalence between models (??) and (??). Thus, if one is only interested in computing some moments of the resulting distribution function a moderate population of particles is enough to describe accurately the evolution of such moments. The size of the pseudo-particles population that must be considered for computing accurately the different moments of the solution scales linearly with the dimension of the space, however, if one want to reconstruct the distribution itself an impressive number of particles is required which don't scale anymore linearly with the dimension of the space.

In general the technique just presented, is conceptually very simple and then easy to implement in a computer or in a parallel computing platform. Explicit integration schemes are usually employed,
needing for a careful choice and control of the time step.

It is nowadays widely recognized that the main drawbacks of stochastic techniques are: (i) the control of the statistical noise that makes difficult the use of the stochastic approach to perform inverse parameter identification or optimization, because the poor accuracy in the sensitivity analysis; (ii) the difficulty to reconstruct the model solution itself even in moderate multidimensional models; and (iii) the necessity to solve always the transient model, even if one is only interested in the steady state.

Moreover, in complex flows simulation using for example a finite element solver for the flow kinematics computation, one must ensure that all the elements contain a number of particles, at least enough to allow computing the virial stress (the usual micro-macro bridge). Different possibilities exist, but all of them have a non negligible impact on the solution. Thus, if new particles are added (and probably others removed) the size of the model is changing, tracking procedures are time consuming, and the initialization of the just introduced particles induces a noticeable numerical diffusion.

One could think that the aforementioned difficulties could be circumvented by employing a Lagrangian description of the flow combined with a Lagrangian description of the microstructure evolution (stochastic approach), however usual mesh-based strategies fails because the high distortion of the meshes during the flow. A first tentative of coupling a meshless Lagrangian description of the flow kinematics (using the natural element method widely described in [?][?]) and a Lagrangian microstructure description have been recently performed in [?]. This technique could be extended for coupling a Lagrangian flow description with a stochastic description of the microstructure evolution. In any case the difficulties just mentioned will persist.

To reduce the computational cost of numerical simulations different model reduction techniques have been recently proposed [?]. However, the coupling of such techniques (based on the proper orthogonal decomposition -also known as Karhunen-Loève decomposition-) with a Lagrangian description of the microstructure evolution remains nowadays an open problem. However, sometimes the stochastic model can be written in an Eulerian form (Brownian Configuration Fields) and in that case, as illustrated in [?], the model reduction runs opening some interesting perspectives.

### 1.2.4 Coming back to continuous descriptions: kinetic theory models

The next level of description concerns the statistical mechanics where the knowledge of individuals is substituted by a kind of averaged knowledge described by a probability distribution function that
depends on some conformational coordinates depending on the considered model. In what follows we are addressing some models involving more or less large conformational spaces.

**The Vlasov-Poisson-Boltzmann equation.** Firstly, we consider the dynamics of N electrically charged particles of masse m. When N becomes too large, direct molecular dynamics simulations result prohibitive from the computing time viewpoint. Thus, more than describing the system from the position and velocities of all the particles, one could introduce the function $f(t, \mathbf{x}, \mathbf{v})$ given the number of particles that at time t, are located within an elementary volume $d\mathbf{x} = (dx, dy, dz)^T$ placed at position $\mathbf{x}$ and having velocities within the volume defined by $d\mathbf{v} = (du, dv, dw)^T$ around $\mathbf{v}$. Now, the density balance writes:

$$ \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla_x f + \mathbf{a}(\mathbf{x}, t) \cdot \nabla_v f = S(t, \mathbf{x}, \mathbf{v}) \quad (14) $$

where $\nabla_x$ and $\nabla_v$ represent the gradient operator in the physical and velocity spaces respectively. We assumed that the acceleration $\mathbf{a} = \frac{d\mathbf{v}}{dt}$ does not depend on the velocity (by this reason it is not affected by the velocity-gradient $\nabla_v$). The source term $S(t, \mathbf{x}, \mathbf{v})$ represents the so called collision term and can be derived from an appropriate physical analysis.

We don't need any physics to model the velocity field $\mathbf{v}$ because now the velocity field is a real coordinate, like the spatial ones. On the contrary, we need to define the acceleration field $\mathbf{a}(\mathbf{x}, t)$. For this purpose we consider the Newton's law $\mathbf{a} = \frac{\mathbf{F}}{m}$, and compute the force acting on the particles by taking into account the nature of the system, that in the case here addressed consists of a population of charged particles interacting by means of the Coulomb's potential. The electrostatic potential $U(\mathbf{x}, t)$ can be computed by solving the associated Poisson's problem:

$$ \Delta U(\mathbf{x}, t) = -4\pi k Q(\mathbf{x}, t) \quad (15) $$

where $k$ is the Coulomb's law constant and $Q(\mathbf{x}, t)$ is the electrical charge inside the volume $d\mathbf{x}$ around $\mathbf{x}$ that can be computed from:

$$ Q(\mathbf{x}, t) = q \int_{\mathbb{R}^3} f(\mathbf{x}, \mathbf{v}, t) d\mathbf{v} \quad (16) $$

with $q$ the particle charge. The force acting on a particle can be finally computed by using:
$$F(\mathbf{x}, t) = -q \nabla_x U(\mathbf{x}, t) \qquad (17)$$

Remarks:

1. When the particles are not charged the steady state solution of Eq. (17) leads to the Maxwell-Boltzmann distribution when the appropriate choice of the collision term representing the particles collisions is made.

2. When the interaction potential does not comply the Coulomb's law, model (17)-(18) cannot be employed. In this case we must compute the force by using:

$$F(\mathbf{x}, t) = \int_{\mathbb{R}^3} F(\mathbf{x}, t; \mathbf{x}') \tilde{f}(\mathbf{x}', t) d\mathbf{x}' \qquad (18)$$

where

$$\tilde{f}(\mathbf{x}, t) = \int_{\mathbb{R}^3} f(\mathbf{x}, \mathbf{v}, t) d\mathbf{v} \qquad (19)$$

and $F(\mathbf{x}, t; \mathbf{x}')$ represents the force at position $\mathbf{x}$ originated by a particle located at position $\mathbf{x}'$. In any case this analysis fails when (i) the inter-particle potential leads to a non-definite integral (19); and (ii) in the case of dense systems where the movements of particles is correlated.

3. When the acceleration does not depend on the velocity (as was assumed in Eq. (19)) and the collision terms vanishes, an equivalence between the conservation equation and the Liouville's theorem in the phase space can be established.

4. In general the establishment of collision terms is quite difficult. To circumvent this difficulty and assuming that the equilibrium distribution is known $f_{eq}(\mathbf{x}, \mathbf{v}, t)$, one could approximate the collision term by

$$S(t, \mathbf{x}, \mathbf{v}) = \frac{f_{eq}(\mathbf{x}, \mathbf{v}, t) - f(\mathbf{x}, \mathbf{v}, t)}{\tau} \qquad (20)$$

where $\tau$ represents a relaxation time. This approximation leads to the so called BFK models.

5. Equation (19), also known as the Vlasov-Poisson-Boltzmann equation, is widely used to model quantum plasmas. One could imagine the extension of this formalism (within the BBGKY hierarchy) to a variety of physical models: colloids, ferrofluids, coarse grained molecular dynamics,
crystallization, demixing, etc. The associated kinetic theory descriptions constitute the Vlasov-Fokker-Planck models.

6. The kinetic theory formalism allows transforming a discrete model into its continuous counterpart. However, in general, the continuous descriptions involve highly multidimensional spaces and their evolutions are governed by hyperbolic non-linear partial differential equations. To solve these kind of models appropriate stabilized solvers, able to proceed in highly multidimensional spaces, are needed.

In the previous paragraphs we introduced some ideas related to kinetic theory models of systems composed of particles. We are introducing in the next paragraphs some models describing the microscopic modelling of polymeric liquids.

**Kinetic theory description of polymeric liquids.** For the sake of simplicity we are focusing in polymer solutions (the entangled systems related to the polymer melts can be also described in the kinetic theory framework [?][?]). We are addressing the Bead-Spring-Chain (BSC) model of polymer solutions. The BSC chain consists of $S + 1$ beads connected by $S$ springs. The bead serves as an interaction point with the solvent and the spring contains the local stiffness information depending on local stretching (see [?][?] for more details). From now on we are also assuming a fully developed homogeneous flow. Thus, the microstructural state does not depend on the space coordinates.

The dynamics of the chain is governed by viscous drag, Brownian and connector forces. If we denote by $\dot{\mathbf{r}}_k$ the velocity of bead $k$ and by $\dot{\mathbf{q}}_k$ the velocity of the spring connector $\mathbf{q}_k$, we have

$$ \dot{\mathbf{q}}_k = \dot{\mathbf{r}}_{k+1} - \dot{\mathbf{r}}_k \quad \forall k = 1, \dots, S \qquad (21) $$

The dynamics of each bead can be written as:

$$ \underbrace{-\zeta(\dot{\mathbf{r}}_k - \mathbf{v}_0 - \nabla\mathbf{v} \cdot \mathbf{r}_k)}_{\text{Viscous drag}} - \underbrace{k_b T \frac{\partial \ln(\psi)}{\mathbf{r}_k}}_{\text{Brownian effects}} + \underbrace{\mathbf{F}_k^c - \mathbf{F}_{k-1}^c}_{\text{Interaction Forces}} = 0, \qquad (22) $$

where $\zeta$ is the drag coefficient, $\mathbf{v}$ is the velocity field, $\mathbf{v}_0$ is an average velocity, $k_b$ is the Boltzmann constant, $T$ is the absolute temperature and $\psi$ is the distribution function $\psi(\mathbf{r}_1, \dots, \mathbf{r}_{S+1}, t)$. From
equations (??) and (??) we obtain:

$$ \dot{\mathbf{q}}_k = \nabla \mathbf{v} \cdot \mathbf{q}_k - \frac{1}{\zeta} \sum_{l=1}^{S} \mathbf{A}_{kl} \cdot \left( k_b T \frac{\partial \ln(\psi)}{\partial \mathbf{q}_l} + \mathbf{F}_l^c \right), \quad (23) $$

where $\mathbf{A}_{kl}$ are the components of the Rouse matrix (see [?]) for more details).

In the Rouse model the connector force $\mathbf{F}^c$ is a linear function of the connector vector, but we could use Finitely Extensible Non-Linear Elastic (FENE) springs, with a dimensionless connector force given by:

$$ \mathbf{F}^c(\mathbf{q}_k) = \frac{1}{1 - \mathbf{q}_k^2/b} \mathbf{q}_k, \quad (24) $$

where $\sqrt{b}$ is the maximum dimensionless length of each spring connector of the chain.

Introducing Eq. (??) in the equation governing the evolution of the distribution function, and considering homogeneous flows,

$$ \frac{\partial \psi(\mathbf{q}_1, \dots, \mathbf{q}_S, t)}{\partial t} = - \sum_{k=1}^{S} \left( \frac{\partial}{\partial \mathbf{q}_k} (\dot{\mathbf{q}}_k \psi(\mathbf{q}_1, \dots, \mathbf{q}_S, t)) \right), \quad (25) $$

we obtain

$$ \frac{\partial \psi}{\partial t} = - \sum_{k=1}^{S} \left( \frac{\partial}{\partial \mathbf{q}_k} \left( (\nabla \mathbf{v} \cdot \mathbf{q}_k - \frac{1}{\zeta} \sum_{l=1}^{S} \mathbf{A}_{kl} \cdot \mathbf{F}_l^c) \psi \right) \right) + \frac{k_b T}{\zeta} \sum_{k=1}^{S} \sum_{l=1}^{S} \mathbf{A}_{kl} \frac{\partial^2 \psi}{\partial \mathbf{q}_k \partial \mathbf{q}_l}. \quad (26) $$

The micro-macro bridging is performed by computing the virial stress, that within the rheology community is known as the Kramer's rule. The main difficulty in using this description is the highly multidimensional problem defined by Eq. (??) that needs for specific advanced solvers as the ones that we proposed in some of our former works and that we revisit in the next section.

**Kinetic theory description of rod-like suspensions in complex flows.** In the case of a dilute suspension of rod-like particles (short fibers, nanofibers, functionalized carbon nanotubes or even rod-like molecules), the configuration distribution function (also known as orientation distribution function) gives the probability of finding the particle in a given direction. Obviously, this function depends on the physical coordinates (space and time) as well as on the configuration coordinates, that taking into account the rigid character of the particles, are defined on the surface of the unit sphere. Thus, we can write $\psi(\mathbf{x}, t, \mathbf{p})$, where $\mathbf{x}$ defines the position of the rod center of mass, $t$ the time and $\mathbf{p}$ the unit vector
defining the rod orientation. The evolution of the distribution function is given by the Fokker-Planck equation

$$ \frac{d\psi}{dt} = -\frac{\partial}{\partial \mathbf{p}}(\psi \dot{\mathbf{p}}) + \frac{\partial}{\partial \mathbf{p}}\left(D_r \frac{\partial \psi}{\partial \mathbf{p}}\right) \qquad (27) $$

where $d/dt$ represents the material derivative, $D_r$ is a diffusion coefficient and $\dot{\mathbf{p}}$ is the particle rotation velocity. The orientation distribution function must verify the normality condition:

$$ \int_{S(0,1)} \psi(\mathbf{p}) d\mathbf{p} = 1 \qquad (28) $$

where $S(0, 1)$ denotes the surface of the unit sphere.

For ellipsoidal particles and when the suspension is dilute enough, the rotation velocity can be obtained from the Jeffery's equation

$$ \dot{\mathbf{p}} = \mathbf{\Omega} \cdot \mathbf{p} + k \mathbf{D} \cdot \mathbf{p} - k(\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \mathbf{p} \qquad (29) $$

where $\mathbf{\Omega}$ and $\mathbf{D}$ are the vorticity and the strain rate tensors respectively, associated with the fluid flow undisturbed by the presence of the suspended particles, and $k$ is a scalar which depends on the particle aspect ratio $\lambda$ (ratio of its length and diameter)

$$ k = \frac{\lambda^2 - 1}{\lambda^2 + 1} \qquad (30) $$

that for rod-like particles can be assumed $k \approx 1$. In complex flows simulations the solution of the Fokker-Planck equation involves some numerical difficulties related to: (i) its multidimensional character, i.e. $\psi(\mathbf{x}, t, \mathbf{p}) : \Omega \subset \mathbb{R}^3 \times \mathbb{R}^+ \times S(0, 1) \to \mathbb{R}^+$; (ii) the geometrical complexity of the physical domain $\Omega$; (iii) its purely advective character in the physical space; and (iv) the advection effects in the conformation space. Both last behaviors need appropriate numerical stabilizations (e.g. upwinding) of discrete models.

## 2 Advanced solvers for multi-dimensional models

The solution of models like the one just addressed in Eq. (??) needs new advanced numerical strategies, because the standard ones suffer the curse of dimensionality. Some strategies have been recently proposed for solving models defined in multi-dimensional spaces. The sparse grid techniques are one
of the most popular [?], but as we described in the first section, separated representations were also used for solving the models encountered in quantum mechanics (Hartree-Fock and post-Hartree-Fock techniques).

The sparse grids [?] are restricted (as argued in [?]) for treating models involving up to twenty dimensions. On the other hand, separated representations like the one considered in the multi-configuration-self-consistent-fields (MCSCF) [?] [?], fix the number of products (anti-symmetrized) and also assumes that those products are constructed by combining a certain number of unknown one-dimensional functions.

In the next section we are describing another alternative and efficient technique recently proposed and based on the use of separated representations within a variational framework.

## 2.1 Circumventing curse of dimensionality by using separated representations

Recently we proposed a new strategy able to solve highly multi-dimensional models circumventing the curse of dimensionality. This technique was successfully applied for treating some multi-dimensional models encountered in the kinetic theory description of complex fluids [?] [?]. It allows defining the optimal number of products containing the optimal one-dimensional functions. In what follows we are summarizing the main ideas that the aforementioned technique involves. For the sake of simplicity, we are illustrating the solution procedure by solving the Poisson problem defined in a space of dimension $D$:

$$ \Delta T = -f(x_1, x_2, \dots, x_D), \qquad (31) $$

where $T$ is a scalar function of $x_1, x_2, \dots, x_D$. Problem (31) is defined in the domain $\Omega =] -L, +L]^D$ with vanishing boundary conditions.

The problem solution can be written in the form:

$$ T(x_1, x_2, \dots, x_D) = \sum_{j=1}^{\infty} \alpha_j \prod_{k=1}^{D} F_{kj}(x_k), \qquad (32) $$

where $F_{kj}$ is the $j^{th}$ approximation function, with unit norm and that depends only on the $k^{th}$ coordinate.
It is well known that the solution of numerous problems can be accurately approximated using a finite (sometimes very reduced) number of approximation functions, i.e.:

$$ T(x_1, x_2, \dots, x_D) \approx \sum_{j=1}^{Q} \alpha_j \prod_{k=1}^{D} F_{kj}(x_k). \qquad (33) $$

The previous expression implies the same number of approximation functions in each dimension, but each one of these functions could be expressed in a discrete form using different number of parameters (nodes of the one-dimensional meshes).

Now, an appropriate numerical procedure is needed for computing the coefficients $\alpha_j$ as well as the $Q \times D$ one-dimensional approximation functions. The proposed numerical scheme consists of an iteration procedure that solves at each iteration $n$ the following three steps:

### Step 1: Projection of the solution in a discrete basis

If we assume the functions $F_{kj}, \forall j \in [1, \dots, n]; \forall k \in [1, \dots, D]$ known (verifying the boundary conditions), the coefficients $\alpha_j$ can be computed by introducing the approximation of $T$ into the Galerkin variational formulation associated with Eq. (33):

$$ \int_{\Omega} \nabla T^* \cdot \nabla T d\Omega = \int_{\Omega} T^* f \ d\Omega. \qquad (34) $$

Introducing the approximation of $T$ and $T^*$:

$$ T(x_1, x_2, \dots, x_D) = \sum_{j=1}^{n} \alpha_j \prod_{k=1}^{D} F_{kj}(x_k), \qquad (35) $$

$$ T^*(x_1, x_2, \dots, x_D) = \sum_{j=1}^{n} \alpha_j^* \prod_{k=1}^{D} F_{kj}(x_k), \qquad (36) $$

into Eq. (33), we have

$$ \int_{\Omega} \nabla \left( \sum_{j=1}^{n} \alpha_j^* \prod_{k=1}^{D} F_{kj}(x_k) \right) \cdot \nabla \left( \sum_{j=1}^{n} \alpha_j \prod_{k=1}^{D} F_{kj}(x_k) \right) d\Omega = $$
$$= \int_{\Omega} \left( \sum_{j=1}^{n} \alpha_{j}^{*} \prod_{k=1}^{D} F_{kj}(x_{k}) \right) f \, d\Omega. \qquad (37)$$

Now, we assume that $f(x_1, \dots, x_D)$ can be written in the form

$$f(x_1, \dots, x_D) = \sum_{h=1}^{m} \prod_{k=1}^{D} f_{kh}(x_k). \qquad (38)$$

Eq. (37) involves integrals of products of $D$ functions each one defined in a different dimension. Let $\prod_{k=1}^{D} g_k(x_k)$ be one of these functions to be integrated. The integral over $\Omega$ can be performed by integrating each function in its definition interval and then multiplying the $D$ computed integrals according to:

$$\int_{\Omega} \prod_{k=1}^{D} g_k(x_k) \, d\Omega = \prod_{k=1}^{D} \int_{-L}^{L} g_k(x_k) dx_k, \qquad (39)$$

which makes possible the numerical integration in highly dimensional spaces.

Now, due to the arbitrariness of the coefficients $\alpha_j^*$, Eq. (37) allows to compute the approximation coefficients $\alpha_j$, solving the resulting linear system of size $n \times n$. This problem is linear and moreover rarely exceeds the order of tens of degrees of freedom in standard elliptic models. Thus, even if the resulting coefficient matrix is densely populated, the time required for its solution is negligible with respect to the one required for performing the approximation basis enrichment (step 3).

### Step 2: Checking convergence

From the solution of $T$ at iteration $n$ given by Eq. (37) we compute the residual $Re$ related to Eq. (37):

$$Re = \frac{\sqrt{\int_{\Omega} (\Delta T + f(x_1, \dots, x_D))^2}}{\|T\|}. \qquad (40)$$

The integral in Eq. (37) can be written as the product of one-dimensional integrals by performing a separated representation of the square of the residual.

If $Re < \epsilon$ (epsilon is a small enough parameter) the iteration process stops, yielding the solution $T(x_1, \dots, x_D)$ given by Eq. (37). Otherwise, the iteration procedure continues.

### Step 3: Enrichment of the approximation basis
From the coefficients $\alpha_j$ just computed the approximation basis can be enriched by adding the new function $\prod_{k=1}^{D} F_{k(n+1)}(x_k)$. For this purpose we solve the Galerkin variational formulation related to Eq. (??):

$$ \int_{\Omega} \nabla T^* \cdot \nabla T d\Omega = \int_{\Omega} T^* f \, d\Omega, \qquad (41) $$

using the approximation of $T$ given by

$$ T(x_1, x_2, \dots, x_D) = \sum_{j=1}^{n} \alpha_j \prod_{k=1}^{D} F_{kj}(x_k) + \prod_{k=1}^{D} R_k(x_k). \qquad (42) $$

and the associated weighting functions

$$ \begin{aligned} T^*(x_1, x_2, \dots, x_D) &= \\ &= R_1^*(x_1) R_2(x_2) \cdots R_D(x_D) + \cdots + R_1(x_1) R_2(x_2) \cdots R_D^*(x_D) \end{aligned} \qquad (43) $$

leading to a non-linear variational problem, whose solution are the $D$ functions $R_k(x_k)$, $k = 1, \dots, D$. After convergence of the non-linear solver, functions $F_{k(n+1)}(x_k)$ are finally obtained by normalizing the functions $R_1, R_2, \dots, R_D$.

The discretization needs an appropriate discrete approximation of functions $R_k(x_k)$ and $F_{kj}(x_k)$. Each one of these functions is approximated using a 1D finite element description. If we assume than $p_k$ nodes are used to construct the interpolation of functions $R_k(x_k)$ and $F_{kj}(x_k)$ in the interval $[-L, L]$, then the size of the resulting discrete non-linear problem is $\sum_{k=1}^{k=D} p_k$. The price to pay for avoiding a whole mesh in the multidimensional domain is the solution of a non-linear problem whose size scales linearly with the dimension of the space. However, even in high dimension the size of the non-linear problems remains moderate and no particular difficulties have been found in its solution.

Concerning the computation time, even when the non-linear solver converges quickly, this step consumes the main part of the global computing time. Different non-linear solvers have been analyzed: fixed-point, Newton or one based on an alternating directions scheme [?] [?] [?].

Some preliminary convergence analysis, one of them presented in the next section, reveal that as expected, when one uses piecewise linear one-dimensional finite elements to approximate all the one-dimensional functions involved by the separated representation, the convergence rate (in the $L^2$-norm)
is second order. At the same time one could modify the one-dimensional approximations to enhance
the approximation consistency and then the convergence rates by using spectral approximations or
use wavelet approximation basis to take advantage of its multi-resolution properties, specially useful
to define adaptive multi-scale strategies. Other possibility to enhance the convergence rate could be
the use of moving least squares -MLS- based approximations (widely used in the context of meshless
methods [?] [?]). We are summarizing the main ideas related to this approach in the next section.

## 2.2 Consistency enhancement by using MLS-based one-dimensional approximations

A first possibility to enhance the approximation consistency lies in the use of higher degree piecewise
polynomials (quadratic or cubic). A potentially better strategy could consist in using one-dimensional
approximation bases involving orthogonal polynomials, as usually encountered in spectral techniques
making use of Chebyshev, Fourier or Legendre polynomials. However, the use of these kind of bases
induce the apparition of fully populated matrices. Moreover, sometimes, there is some "a priori"
knowledge about the solution, that one would like to introduce in the approximation basis. One
possibility to perform this enrichment lies in the use of the partition of unity paradigm [?], like in
the extended finite element method (X-FEM) [?], however this strategy needs for the introduction of
additional degrees of freedom related to the enrichment functions. The use of partition of unity for
enhancing the approximation of the electronic density in the context of the density functional theory
of quantum mechanics was employed in [?].

A nice compromise between all the previous techniques lies in the definition of $C^\infty$ polynomial approximation, of any degree and with compact support, within the context of the reproducing kernel particle method -RKPM- [?]. Moreover, this framework allows us introducing any enrichment function in the reproduction conditions without the introduction of additional degrees of freedom. In what follows we are describing the construction of such enriched approximation, focusing in a simple one-dimensional model.

Let $\Omega$ be a 1D domain where the problem is defined. The points within this domain will be noted by $x$ or $s$.

### 2.2.1 Reproduction conditions

The approximation $u^h(x)$ of $u(x)$ is built from the convolution integral
$$u^h(x) = \int_{\Omega} w(x - s, h)u(s)ds \qquad (44)$$

where $w(x - s, h)$ is the kernel function and $h$ a parameter defining the size of the approximation support.

The main idea in the enriched RKPM method [?] is to enforce the reproduction of a generic function that we can write in the form of a polynomial plus another function noted by $u^e(x)$:

$$u^h(x) = a_0 + a_1x + \dots + a_nx^n + a_{n+1}u^e(x) \qquad (45)$$

In the following paragraphs we analyze the required properties of the kernel function $w(x - s, h)$ for reproducing a function expressed by Eq. (3).

From Eq. (3), the reproduction of a constant function $a_0$ writes

$$\int_{\Omega} w(x - s, h)a_0ds = a_0 \qquad (46)$$

which implies

$$\int_{\Omega} w(x - s, h)ds = 1 \qquad (47)$$

which constitutes the partition of unity.

Now, the required condition to reproduce a linear function $u^a(x) = a_0 + a_1x$ is

$$\int_{\Omega} w(x - s, h)(a_0 + a_1s)ds = a_0 + a_1x \qquad (48)$$

By using the partition of unity (3), Eq. (3) can be rewritten as

$$\left\{ \begin{array}{l} \int_{\Omega} w(x - s, h) ds = 1 \\ \int_{\Omega} w(x - s, h)s ds = x \end{array} \right. \qquad (49)$$

which implies the linear consistency of the approximation. Repeating this reasoning, we can write the $n$-order consistency as
$$ \left\{ \begin{array}{l} \int_{\Omega} w(x-s, h) \, ds = 1 \\ \int_{\Omega} w(x-s, h) s \, ds = x \\ \vdots \\ \int_{\Omega} w(x-s, h) s^n \, ds = x^n \end{array} \right. \qquad (50) $$

and consequently, the reproduction of the function given by (50) implies

$$ \begin{aligned} & \int_{\Omega} w(x-s, h)(a_0 + a_1 s + \dots + a_n s^n + a_{n+1} u^e(s)) \, ds = \\ & = a_0 + a_1 x + \dots + a_n x^n + a_{n+1} u^e(x) \end{aligned} \qquad (51) $$

from which it results

$$ \left\{ \begin{array}{l} \int_{\Omega} w(x-s, h) \, ds = 1 \\ \int_{\Omega} w(x-s, h) s \, ds = x \\ \vdots \\ \int_{\Omega} w(x-s, h) s^n \, ds = x^n \\ \int_{\Omega} w(x-s, h) u^e(s) \, ds = u^e(x) \end{array} \right. \qquad (52) $$

In the original procedure proposed by Liu et al. [?] only n-order consistency was imposed, but it can not be directly used to enforce the reproduction condition associated with the non-polynomial function $u^e(x)$.

### 2.2.2 The moment matrix

We denote by $u^r(x)$ the approximation function verifying the conditions (50). Usually a cubic spline is considered as kernel function, and consequently the conditions given by Eq. (50) are not satisfied. Liu et al. [?] proposed the introduction of a correction function $C(x, x-s)$ for satisfying the reproduction conditions. In our case we consider the more general form $C(x, s, x-s)$ whose pertinence was discussed in some of our former works [?] [?]. Thus $u^r(x)$ will be expressed by

$$ u^r(x) = \int_{\Omega} C(x, s, x-s) w(x-s, h) u(s) \, ds \qquad (53) $$
where $C(x, s, x-s)$ is assumed to have the following form

$$C(x, s, x-s) = \mathbf{H}^T(x, s, x-s)\mathbf{b}(x) \qquad (54)$$

where $\mathbf{H}^T(x, s, x-s)$ represents the vector containing the functions considered in the approximation basis, and $\mathbf{b}(x)$ is a vector containing unknown functions that will be determined for satisfying the reproduction conditions. Thus, Eq. (54) can be rewritten as

$$\left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \mathbf{H}^T(x, s, x-s) \mathbf{b}(x) w(x-s, h) \, ds = 1 \\[1.5ex]
\displaystyle \int_{\Omega} \mathbf{H}^T(x, s, x-s) \mathbf{b}(x) w(x-s, h) s \, ds = x \\[1.5ex]
\vdots \\[1.5ex]
\displaystyle \int_{\Omega} \mathbf{H}^T(x, s, x-s) \mathbf{b}(x) w(x-s, h) s^n \, ds = x^n \\[1.5ex]
\displaystyle \int_{\Omega} \mathbf{H}^T(x, s, x-s) \mathbf{b}(x) w(x-s, h) u^e(s) \, ds = u^e(x)
\end{array}
\right.
\qquad (55)$$

In fact, the reproduction conditions must be enforced in a discrete form. For this purpose we consider $N$ points (also refereed as nodes) which allow to compute the discrete form of Eq. (54), i.e.

$$\left\{
\begin{array}{l}
\displaystyle \sum_{i=1}^{N} \mathbf{H}^T(x, x_i, x-x_i) \mathbf{b}(x) w(x-x_i, h) \Delta x_i = 1 \\[1.5ex]
\displaystyle \sum_{i=1}^{N} \mathbf{H}^T(x, x_i, x-x_i) \mathbf{b}(x) w(x-x_i, h) x_i \Delta x_i = x \\[1.5ex]
\vdots \\[1.5ex]
\displaystyle \sum_{i=1}^{N} \mathbf{H}^T(x, x_i, x-x_i) \mathbf{b}(x) w(x-x_i, h) x_i^n \Delta x_i = x^n \\[1.5ex]
\displaystyle \sum_{i=1}^{N} \mathbf{H}^T(x, x_i, x-x_i) \mathbf{b}(x) w(x-x_i, h) u^e(x_i) \Delta x_i = u^e(x)
\end{array}
\right.
\qquad (56)$$

that in a matrix form results

$$\left[ \sum_{i=1}^{N} \mathbf{R}(x_i) \mathbf{H}^T(x, x_i, x-x_i) w(x-x_i, h) \Delta x_i \right] \mathbf{b}(x) = \mathbf{R}(x) \qquad (57)$$

where $\mathbf{R}(x)$ is the reproduction vector

$$\mathbf{R}^T(x) = [1, x, \dots, x^n, u^e(x)] \qquad (58)$$

Eq. (54) allows computing the vector $\mathbf{b}(x)$,

$$\mathbf{b}(x) = \mathbf{M}(x)^{-1} \mathbf{R}(x) \qquad (59)$$
where the moment matrix $M(x)$ is defined by

$$M(x) = \sum_{i=1}^{N} R(x_i)H^T(x, x_i, x - x_i)w(x - x_i, h)\Delta x_i \quad (60)$$

This moment matrix differs from the usual moment matrix proposed in [?], and in our case it becomes non symmetric.

### 2.2.3 Discrete form of the approximation function

The discrete form $u^r(x)$ of $u^h(x)$ derives from Eqs. (??), (??) and (??)

$$ \begin{aligned} u^r(x) &\approx \sum_{i=1}^{N} H^T(x, x_i, x - x_i) M(x)^{-1} R(x) w(x - x_i, h) u(x_i) \Delta x_i = \\ &= \sum_{i=1}^{N} \psi_i(x) u_i \end{aligned} \quad (61) $$

where $\psi_i$ is the enriched RKP approximation shape function

$$\psi_i(x) = H^T(x, x_i, x - x_i) M(x)^{-1} R(x) w(x - x_i, h) \Delta x_i \quad (62)$$

As in the classical RKPM we consider $\Delta x_i = 1$. Different quadrature rules exist and they have been tested without a significant incidence on the reproducing condition accuracy.

In principle this procedure allows: (i) defining $C^\infty$ approximations of any order of consistency; (ii) the associated shape functions have a compact support avoiding fully populated linear systems; and (iii) the approximation can be enriched by considering any knowledge about the solution through the use of the enrichment function $u^e(x)$. The analysis of the just described approximation is being performed and preliminary results will be shortly reported.

## 3 Numerical examples

In this section we are illustrating the application of separated representations to the solution of different models defined in highly multidimensional spaces suffering of the so-called curse of dimensionality.

In [?] [?] we described the application of separated representations to the solution of highly multidimensional models arising from kinetic theory modelling of molecular systems. It was noticed that
the separated representation avoided the curse of dimensionality. In this section we explore the applicability of this technique for addressing more complex systems.

Firstly we address a problem whose exact solution is known in order to quantify the capabilities of the discretization technique. Then, the solution of a rod-like suspension flow (whose model was introduced in section ??) will be considered. Finally, we come back to the numerical treatment of the Schrödinger equation (introduced in section ??), where the main challenge related to its solution will be pointed out.

## 3.1 A multidimensional Poisson problem

We consider the problem

$$ \Delta T(\mathbf{x}) = -f(\mathbf{x}); \mathbf{x} \in \Omega =]0, \pi^{[D} \quad (63) $$

where $T(\mathbf{x}) : \mathbb{R}^D \to \mathbb{R}$. The components of $\mathbf{x}$ will be denoted by $(x_1, x_2, \dots, x_D)$. In what follows we set $D = 100$, i.e. the resulting model is defined in a space of hundred dimensions.

If the source term $f(\mathbf{x})$ in Eq. (??) takes the form:

$$ f(x_1, \dots, x_D) = \sum_{k=1}^{k=D} \left\{ (1+k)(\sin(x_k))^{k-1} ((\sin(x_k))^2 - k(\cos(x_k))^2) \prod_{l=1, l \neq k}^{l=D} (\sin(x_l))^{l+1} \right\} \quad (64) $$

then, the exact solution of Eq. (??) writes:

$$ T^{\text{exact}}(x_1, \dots, x_D) = \prod_{k=1}^{k=D} (\sin(x_k))^{k+1} \quad (65) $$

The numerical solution was searched by assuming the separated representation

$$ T^h(x_1, x_2, \dots, x_D) \approx \sum_{j=1}^{k=n} \alpha_j \prod_{k=1}^{k=D} F_{kj}(x_k), \quad (66) $$

The one-dimensional functions $F_{kj}$ where approximated using 1D linear finite elements. Different nodal spacings were considered: $h = \pi/5$, $h = \pi/10$, $h = \pi/20$, $h = \pi/40$, $h = \pi/80$ and $h = \pi/160$. Thus, when the finest mesh was considered, the number of nodes in each direction was $p_k = 161$, that for $D = 100$ would have implied, in the context of mesh-based techniques, the use of much more than $10^{200}$ nodes.
<figure><img src="image_1.png" /><figcaption>Figure 1: Convergence analysis of the separated representation technique for a multidimensional Poisson problem</figcaption></figure>

In all the cases the solution was reached after a single iteration, i.e. $n = 1$, and associated error was computed according to

$$E = \|T^{\text{exact}}(\mathbf{x}) - T^{h}(\mathbf{x})\|_{2} \qquad (67)$$

Figure ?? illustrates the evolution of the error (E) as a function of the nodal spacing (h) proving the expected second order rate of convergence.

## 3.2 Kinetic theory description of rod-like aggregating suspensions

In this section we come back to the suspension model summarized in section ?? . For increasing the model complexity we are assuming that the rods can flocculate creating large aggregates that due to the shear induced by the flow, are continuously broken. Thus, aggregation and disaggregation mechanisms
coexist and two populations can be identified: the one related to free rods (pendant population) and the one associated with the aggregated rods (active population). Figure ?? depicts both populations and the flow induced aggregation/disaggregation.

<figure><img src="image_2.png" /><figcaption>Figure 2: Flow induced aggregation/disaggregation</figcaption></figure>

The kinetic theory description of such systems contains two coupled Fokker-Planck equations given the orientation distribution of rods belonging to each one of these populations, $\Psi(\mathbf{x}, t, \mathbf{p})$ and $\Phi(\mathbf{x}, t, \mathbf{p})$ for the active and pendant respectively:

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial \mathbf{p}}(\dot{\mathbf{p}}\Psi) + D_{r1}\frac{\partial^2\Psi}{\partial\mathbf{p}^2} - V_d\Psi + V_c\Phi \quad (68) $$

$$ \frac{d\Phi}{dt} = -\frac{\partial}{\partial \mathbf{p}}(\dot{\mathbf{p}}\Phi) + D_{r2}\frac{\partial^2\Phi}{\partial\mathbf{p}^2} + V_d\Psi - V_c\Phi \quad (69) $$

where $V_d$ and $V_c$ represent the velocity of destruction and construction of the active population respectively.

The normality condition writes
$$ \int_{S(0,1)} (\Psi(\mathbf{x}, t, \mathbf{p}) + \Phi(\mathbf{x}, t, \mathbf{p})) \, d\mathbf{p} = 1, \quad \forall \mathbf{x}, \forall t \qquad (70) $$

We consider that the flow takes place in a converging channel. The steady state flow kinematics (assumed undisturbed by the presence of the suspended particles) was computed by solving the Stokes equations.

As we are interested in computing the steady state solution, and because the advection character of the Fokker-Planck equations in the spatial coordinates (x), we decided to integrate both coupled Fokker-Planck equations along some particular flow streamlines. The separated representation of both orientation distribution functions writes:

$$ \begin{pmatrix} \Psi_{st}(s, \mathbf{p}) \\ \Phi_{st}(s, \mathbf{p}) \end{pmatrix} = \sum_{j=1}^{n} \begin{pmatrix} \alpha_j^{st} E_j^{st}(\mathbf{p}) F_j^{st}(s) \\ \beta_j^{st} G_j^{st}(\mathbf{p}) H_j^{st}(s) \end{pmatrix} \qquad (71) $$

where s denotes the curvilinear coordinate defining the streamlines, and the index "st" refers to the particular streamline along which the integration is performed.

Figure ?? depicts the resulting orientation distributions of both populations at some points on some flow streamlines. In this figure the orientation distribution is directly depicted on the unit surface, and the color indicate the intensity of the orientation in each direction.

## 3.3 On the separated representation of the Schrödinger equation

Finally we are considering a challenging problem, the one associated with the solution of the Schrödinger equation introduced in section ??.

In this case to circumvent the curse of dimensionality that this model involves we could try to perform, as in the previous examples, a separated representation of the wavefunction Ψ.

In what follows we consider a system composed of N electrons and a single nucleus with atomic number Z = N. In this case one could assume the time-independent wavefunction approximation:

$$ \Psi(\mathbf{x}_1, \dots, \mathbf{x}_N; \mathbf{X}) \approx \sum_{j=1}^{j=n} \alpha_j \phi_j^1(\mathbf{x}_1) \phi_j^2(\mathbf{x}_2) \dots \phi_j^N(\mathbf{x}_N) \qquad (72) $$

where X denotes the nucleus position.

However, the Pauli's exclusion principle that applies for fermions (electrons belong to the fermions family) implies the antisymmetry of the wavefunction. To ensure the antisymmetry one could proceed
<figure><img src="image_3.png" /><figcaption>Figure 3: Orientation distribution of active and pendant populations in a contraction flow</figcaption></figure>
by introducing the Slater's determinants:

$$D_j = \begin{vmatrix} \phi_j^1(\mathbf{x}_1) & \phi_j^2(\mathbf{x}_1) & \cdots & \phi_j^N(\mathbf{x}_1) \\ \phi_j^1(\mathbf{x}_2) & \phi_j^2(\mathbf{x}_2) & \cdots & \phi_j^N(\mathbf{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_j^1(\mathbf{x}_N) & \phi_j^2(\mathbf{x}_N) & \cdots & \phi_j^N(\mathbf{x}_N) \end{vmatrix} \quad (73)$$

that introduced in the separated representation

$$\Psi(\mathbf{x}_1, \cdots, \mathbf{x}_N; \mathbf{X}) \approx \sum_{j=1}^{j=n} \alpha_j D_j \quad (74)$$

ensures the antisymmetry of the resulting wavefunction.

However, the last representation hides two difficulties. The first one related to the multidimensional character of the model (we must recall that the wavefunction is defined in a space of dimension $3 \times N$) and the second one related to the complexity in evaluating the Slater determinants that scales in $N!$. In order to isolate each difficulty, we are considering an "imaginary world" in which the electrons are not subjected to the Pauli's exclusion principle. In this case the separated representation (??) works and it could be used to solve the Schrödinger equation.

For the sake of simplicity, in what follows, we assume that each electron "lives" in a 1D dimensional space, i.e. $\Psi(x_1, \cdots, x_N; X)$, where $x_k \in \mathbb{R}, \forall k$.

Now, we consider different systems containing the more and more electrons: $N = 1, 2, 3, 4, 5, 10, 20, 50$. In the last case the Schrödinger equation is defined in a space of dimension $D = 50$.

For each system, the time independent Schrödinger equation (??) is solved using the separated representation (??), and the ground state eigenvalue $\Psi_0$ and its associated energy $E_0$ computed. Now, the ground state electronic distribution $\rho^0(x)$ could be obtained by calculating the ground state distribution of each electron $e$, $\rho_e^0$:

$$\rho_e^0(x) = \int_{\mathbb{R}^{N-1}} |\Psi_0(x_1, \cdots, x_{e-1}, x_{e+1}, \cdots x_N)|^2 dx_1 \cdots dx_{e-1} dx_{e+1} \cdots dx_N \quad (75)$$

and adding all the electrons contributions

$$\rho^0(x) = \sum_{e=1}^{e=N} \rho_e^0 \quad (76)$$

Figure ?? depicts the ground-state electronic distribution for $N = 1, N = 2, N = 3, N = 5$
<figure><img src="image_4.png" /><figcaption>Figure 4: Ground state electronic density for systems composed of 1, 2, 3, 5 and 10 electrons -the Pauli's principle has not been applied-</figcaption></figure>

and $N = 10$. Obviously the electronic distribution is concentrated around the nucleus position and increases as the number of electrons increases. In that figure, to ensure a good quality in the curves resolution we do not include the electronic distributions of the systems composed of 20 and 50 electrons.

It is easy to verify that because the Pauli's exclusion principle does not apply, $\rho^0(N) = N\rho^0(1)$, where $\rho^0(N)$ and $\rho^0(1)$ denote the ground state densities of systems composed of $N$ and 1 electrons respectively. Thus, we can conclude that all the electrons are occupying the same orbital, a kind of $s$-orbital (that in 3D systems would be spherical).

Figure ?? depicts for these "imaginary" systems the evolution of the ground state energy (the most negative eigenvalue related to the solution of Eq. (??)) as a function of the number of electrons involved in the system.

From this analysis we conclude that separated representation allow to circumvent the curse of dimensionality related to quantum mechanics systems, making possible the treatment of multidimensional models. However, the previous analysis concerns "non real" systems, because electrons are fermions, and for fermions the Pauli's principle applies.
<figure><img src="image_5.png" /><figcaption>Figure 5: Evolution of the ground state energies with the number of electrons involved in the system -the Pauli's principle has not been applied-</figcaption></figure>
Even if separated representations could be applied for solving bosonic quantum systems (the Pauli's principle does not apply for bosons) in material science such systems have a limited interest. For this reason we must come back to the "real world" of electronic systems.

In the present case, the separated representation to be considered is the one making use of the Slater determinants (??). All the practical details concerning the discretization of the Schrödinger equation by using a wavefunction separated representation within the Slater's determinants formalism can be found in [?]. In this case, as the complexity scales with the factorial of the number of electrons involved in the quantum system, at present, in our knowledge, only small systems can be solved (containing up to ten electrons).

Figure ?? depicts the ground state electronic distribution of systems composed of 1 until 5 electrons and a single nucleus with atomic number $Z = 3$. It also shows the differences between each two successive configurations. We can notice that the first two electrons are located in a $s$-type orbital (both electrons have different spin which makes it possible). For more than two electrons we appreciate the appearance of $p$-type orbitals. These $p$-type orbitals can be identified by subtracting the electronic density functions related to both kinds of populations (up to two electrons and more than two electrons).

Figure ?? compares the evolution of the ground state energy with the number of electrons $N$ involved in the quantum system composed on a single nucleus with $Z = N$, when the Pauli's principle is or not taken into account. For comparison purposes, when the Pauli's exclusion principle was activated, the spin coordinate was removed, because like this each electron is enforced to occupy a different orbital (as was the case when the Pauli's principle was not taken into account). We can notice that its consideration increases the value of the ground state energy.

From the previous analysis we can affirm that the remaining unsolved difficulty related to the solution of the Schrödinger equation lies in the antisymmetry restriction that the Pauli's principle addresses for fermions, and that as previously argued, its complexity scales as $N!$, which becomes extremely large even for systems composed of a moderate number of electrons.

# 4 Conclusions

In this paper we revise different physical descriptions at different scales involved in material science. At the lowest scale, quantum mechanics involves high-dimensional models whose solution, when fermions are considered, must be anti-symmetric. We have pointed out that in fact, more than the multidimensional character of quantum mechanics models, the real challenge is the antisymmetry con-
<figure><img src="image_6.png" /><figcaption>Figure 6: Ground state electronic density for systems composed of 1, 2, 3, 4 and 5 electrons -the Pauli's principle applies-</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 7: Evolution of the ground state energy when the Pauli's principle is activated</figcaption></figure>

constraint.

Coarser modelling involves molecular dynamics or Brownian dynamics but despite its conceptual simplicity its use is not exempted of computational and conceptual difficulties. The next level is the one that corresponds to statistical mechanics descriptions leading to continuous models described by non-linear and coupled multidimensional partial differential equations. Some of these models have been successfully solved in some of our former works by using the adaptive separated representation revisited in the present work. This technique seems a promising alternative for addressing highly-multidimensional models, but as argued in the previous section, its use in quantum mechanics systems composed of fermions is limited, at present, by the antisymmetry restriction.

## References

[1] Y. Achdou, O. Pironneau, Computational methods for option pricing, SIAM Front. Appl. Math., 2005.
[2] A. Ammar, D. Ryckelynck, F. Chinesta, R. Keunings, On the reduction of kinetic theory models related to finitely extensible dumbbells, 134, 2006, 136-147.

[3] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, J. Non-Newtonian Fluid Mech., 139, 2006, pp. 153-176.

[4] A. Ammar, B. Mokdad, F. Chinesta et R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space-time separated representations, J. Non-Newtonian Fluid Mech., 144, 2007, pp. 98-121.

[5] A. Ammar, F. Chinesta, Circumventing curse of dimensionality in the solution of highly multidimensional models encountered in quantum mechanics using meshfree separated representations, Lecture Notes on Computational Science and Engineering, Springer, In press.

[6] T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, International Journal for Numerical Methods in Engineering, 37, 1994, pp. 229-256.

[7] H. Ben Dhia, Multiscale mechanical problems: the Arlequin method, C. R. Acad. Sci., Paris 326, Ser-II b, 1998, pp. 899-904.

[8] G. Beylkin, M. Mohlenkamp, Algorithms for numerical analysis in high dimensions, SIAM J. Sci. Commun., 26/6, 2005, pp. 2133-2159.

[9] B.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager, Dynamics of polymeric liquids, in: Kinetic Theory, vol. 2, John Wiley & Sons, 1987.

[10] C. Le Bris, P. L. Lions, From atoms to crystals: a mathematical journey, Bulletin of the American Mathematical Society, 42/3, 2005, pp. 291-363.

[11] C. Le Bris, Computational chemistry from the perspective of numerical analysis, Acta Numerica, 2005, pp. 363-444.

[12] C. Le Bris, Mathematical and numerical analysis for molecular simulation: accomplishments and challenges, Proceedings of the International Congress of Mathematicians, Madris, Spain, 2006, pp. 1507-1522.
[13] H.J. Bungartz, M. Griebel, Sparse grids, Acta Numer., 13, 2004, pp. 1-123.

[14] E. Cancès, M. Defranceschi, W. Kutzelnigg, C. Le Bris, Y. Maday, Computational Quantum Chemistry: a primer, Handbook of Numerical Analysis, Vol. X, Elsevier, 2003, pp. 3-270.

[15] J.S. Chen, W. Hu, M. Puso, Orbital HP-clouds for solving Schrödinger equation in quantum mechanics, Comput. Methods Appl. Mech. Engrg., 196/37-40, 2007, pp. 3693-3705.

[16] F. Chinesta, A. Ammar, A. Falco, M. Laso, On the reduction of stochastic kinetic theory models of complex fluids, Modelling and Simulation in Material Science and Engineering, 15, 2007, pp. 639-652.

[17] V. Gavini, K. Bhattacharya, M. Ortiz, Quasi-continuum orbital-free density-functional theory: A route to multi-million atom non-periodic DFT calculation, Journal of the Mechanics and Physics of Solids, 55/4, 2007, pp. 697-718.

[18] D. Gonzalez, E. Cueto, F. Chinesta, M. Doblare, An updated Lagrangian strategy for free-surface fluid dynamics, Journal of Computational Physics, 223, 2007, pp. 127-150.

[19] Handbook of Numerical Analysis, Vol. X: Computational Chemistry, C. Le Bris editor, Elsevier, 2003.

[20] R.G. Larson, The structure and rheology of complex fluids, Oxford University Press, New York, 1999.

[21] R.B. Laughlin, The theory of everything, Proceeding of the U.S.A. National Academy of Sciences, 2000.

[22] W.K. Liu, S. Jun, Y.F. Zhang, Reproducing kernel particle methods, Int. J. Numer. Methods Fluids, 21, 1995, pp. 1081-1106.

[23] M.A. Martinez, E. Cueto, M. Doblare, F. Chinesta, Natural element meshless simulation of injection processes involving short fiber suspensions, Journal of Non-Newtonian Fluid Mechanics, 115, 2003, pp. 51-78.

[24] M.A. Martinez, E. Cueto, I. Alfaro, M. Doblare, F. Chinesta, Updated Lagrangian free surface flow simulations with the natural neighbour Galerkin methods, International Journal for Numerical Methods in Engineering, 60/13, 2004, pp. 2105-2129.
[25] J.M. Melenk, I. Babuska, The partition of unity finite element method: basic theory and applications, Computer Methods in Applied Mechanics and Engineering, 4, 1996, pp. 289-314.

[26] N. Möes, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering 46, 1999, pp. 131-150.

[27] B. Nayroles, G. Touzot, P. Villon, Generalizing the finite element method: diffuse approximation and diffuse elements, Computational Mechanics, 10, 1992, pp. 307-318.

[28] H.C. Öttinger, M. Laso, Smart polymers in finite element calculation, In P. Moldanaers and R. Keunings (Editors) Theoretical and Applied Rheology, Vol. 1, Proceedings of the XIth International Congress on Rheology (Elsevier), 1992, pp. 286-288.

[29] V.B. Shenoy, R. Millera, E.B. Tadmor, D. Rodney, R. Phillips, M. Ortiz, An adaptive finite element approach to atomic-scale mechanics – the quasicontinuum method, Journal of the Mechanics and Physics of Solids, 36, 1999, pp. 500-531.

[30] J. Trunzler, P. Joyot, F. Chinesta, Discontinuous derivative enrichment in RKPM, Lectures Notes on Computational Science and Engineering, Springer, 43, 2005, pp. 93-108.

[31] J. Trunzler, P. Joyot, F. Chinesta, P. Reuter, Enriched reproducing kernel particle approximation for simulating problems involving moving interfaces, Lectures Notes on Computational Science and Engineering, Springer, 57, 2007, pp. 149-164.

[32] G.J. Wagner, W.K. Liu, Coupling of atomistic and continuum simulations using a bridging scale decomposition, Journal of Computational Physics, 190, 2003, pp. 249-274.

[33] S.P. Xiao, T. Belytschko, A bridging domain method for coupling continua with molecular dynamics, Comput. Methods Appl. Mech. Engrg., 193, 2004, pp. 1645-1669.