friends. They are Xia Qing

Wen Jing

Zhang Bo

Zheng Shuo

Yan Haile

Liu Xiaorui

Fan Jiangkun

Wang Peng

Li Fan

Tian Haitao

Huang Zhicheng

Zhao Yajun

Xu Shun

Huang Qun

Zhang Yuqi

Zhang Lin They

Joanna Akrouche

Sobhe Hamiye

Olfat Deeb

Sara Al Jawhari

Nadège Meyer They

First of all, I would like to thank Prof.

Chapter 1 Introduction

This thesis investigates mechanical properties such as elasticity, viscoelasticity, and glass transition of model systems by means of molecular simulations.

Elasticity is one of the fundamental properties of materials. It shows the deformation capability of a substance under the external force field. In macroscopic scale and in the linear response regime, Hooke's law is used to characterize the elasticity of materials, by introducing the elastic constants. With the development of simulation techniques, we are able to link microscopic aspects of materials with these macroscopic constants. The fluctuation formalism method is one of them, in which the elastic constants are related to the fluctuations of the stress tensor components. This approach was first introduced by Squire et al. [START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF] to calculate the isothermal elastic constants for argon crystal at finite temperatures by Monte Carlo (MC)simulations. The formalism was also used in Molecular Dynamics (MD) simulations by Barrat et al. [START_REF] Barrat | Elastic response of a simple amorphous binary alloy near the glass transition[END_REF] to calculate elastic response of a binary alloy near the glass transition. The remarkable advantage of this approach is that the unstrained systems are used, thus the elastic constants can be obtained via equilibrium molecular simulations. In reference [START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF], the constants are expressed by the sum of three contributions: the kinetic term, the Born term and the stress-fluctuation term. The kinetic term is the ideal gas term. It is zero for T = 0. The authors suggested that the fluctuation term arises only with non zero temperatures. Thus at zero temperature, only the Born term would subsist, leading to the well known formulas by Born and Huang [START_REF] Born | Dynamical theory of crystal lattices[END_REF] for the elasticity of solids at T = 0. The issue was later discussed by Lutsko [START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF], who pointed out that even at T = 0, the fluctuation term can be non zero if the solid is complex, for example 1 having more than one atom per primitive unit cell. The Born term accounts for the average stress resulting from an affine uniform microscopic strain field of the system in response to a macroscopic strain. The fluctuation term, at T = 0, corresponds hence to non-affine displacements of the particles in order to relax the internal stress. At T > 0, it represents indeed the thermal fluctuations of the stress tensor elements. This prediction was verified for model glassy systems at T = 0, where the fluctuation term was found as important at half the magnitude of the Born term [START_REF] Wittmer | Vibrations of amorphous, nanometric structures: When does continuum theory apply? EPL[END_REF]. Another interest of the method is that calculations can be done for P = 0 [START_REF] Jarić | Density-functional theory of elastic moduli: Hardsphere and Lennard-Jones crystals[END_REF], making it an attractive tool to study the elasticity of materials for T = 0, and P = 0. In crystals, there are usually more than two independent elastic constants [START_REF] Ashcroft | Solid state physics[END_REF]. For isotropic bodies, there are only two of them. The convenient set can be the so-called elastic moduli [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF], consisting of the bulk modulus K and the shear modulus G. Their macroscopic definitions are well known: K = -(1/V )(∂P/∂V ) T and G = δτ /δγ, where δτ is the average shear stress and δγ a small shear strain. Their "microscopic" expressions are well defined too within the stress fluctuation formalism and shall be explained in detail in the thesis.

Although this formalism was initially proposed for solids, the formulas for computing the shear modulus G and the bulk modulus K are found to be valid for the liquid state, where one obtains consistently G = 0 and K being the same expression as the one obtained by Rowlinson in the context of an isotropic liquid state [START_REF] Rowlinson | Liquids and Liquid Mixtures[END_REF]. This observation is important as one studies the liquid-solid phase transition, as we do in this thesis.

We are thus able to use the stress fluctuation formalism on both phases and through the transition. Beyond characterizing G and K through a liquid-crystal transition, we are also interested in low temperature crystals. We shall verify indeed that for simple crystals, the non-affine contribution to the elasticity is zero, whereas for more complex crystals, such as the AB13 system [START_REF] Xu | A density functional study of superlattice formation in binary hard-sphere mixtures[END_REF], this contribution is on the other hand non-zero, as suggested by Lutsko. All materials display some viscoelastic response. The mathematical formulation of viscoelasticity theory makes it possible to predict the material response to any load history. Viscoelasticity is beneficial in a number of fields of materials science, metallurgy, and solid state physics as it is causally connected to various microphysical processes and can be used as an experimental probe of those processes [START_REF] Lakes | Viscoelastic materials[END_REF]. The causal links between viscoelasticity and microstructure are exploited in the use of viscoelastic tests as an inspection tool as well as in the design of materials [START_REF] Lakes | Viscoelastic materials[END_REF]. In complex fluids, viscoelasticity is also an important property [START_REF] Larson | The structure and rheology of complex fluids[END_REF][START_REF] Colby | Polymer physics[END_REF]. It can reveal different time scales involved in the stress relaxation process. It is thus a precious probe of dynamics of complex fluids. The shear relaxation modulus G(t) gives the measure of the viscoelasticity of a system. This function can be measured in rheological experiments. More usually, the dynamic moduli G ′ (ω) and G ′′ (ω) (with ω the frequency) are measured, thanks to oscillatory stresses. As G ′ and G ′′ are related to G(t) in a simple way via sine and cosine Fourier transformations, knowledge of one set leads to the other set.

As we know, computer simulations can bring a lot for the understanding of materials properties. It is important that we can compute these functions by using statistical physics tools. In liquids, it is well known [START_REF] Hansen | Theory of simple liquids[END_REF] that G(t) is given by the shear stress autocorrelation function (SACF), which can be computed by equilibrium molecular dynamics simulations. For example, recently Heyes et al. [START_REF] Heyes | Scaling of Lennard-Jones liquid elastic moduli, viscoelasticity and other properties along fluid-solid coexistence[END_REF] investigated the elastic moduli, viscoelasticity properties of the Lennard-Jones fluid along the fluid-solid coexistence line by equilibrium Molecular Dynamics simulation. For solids, however, generally one cannot obtain G(t) directly from the SACF [START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF]. Recently, a formula [START_REF] Wittmer | Simple average expression for shearstress relaxation modulus[END_REF] has been proposed and tested for various model systems. It has been proven to be robust and general for systems as different as permanent elastic bodies [START_REF] Wittmer | Simple average expression for shearstress relaxation modulus[END_REF], or selfassembled transient networks [START_REF] Wittmer | Shear-stress fluctuations in self-assembled transient elastic networks[END_REF]. In this thesis, we shall apply this formula to our model systems. Our aim is to compare the dynamic moduli in various states: liquids, crystals and glasses, in order to reveal some general trends, related to different states.

The third topic of the thesis is the glass transition. Generally speaking, the glass transition is a transition by which a liquid changes into an amorphous state without crystallization. The most ordinary way of glass transition is quenching, namely, supercooling a liquid. By this method, lots of glassy substances can be obtained. The remarkable phenomenon upon supercooling is that the viscosity increases dramatically from the liquid side to glass. The glass is similar to the relevant liquid in structure. In terms of mechanical properties, it is solid-like. The glass transition has attracted much attention. Many models and theoretical approaches are proposed to explain it, such as Adam-Gibbs model, Kob-Andersen model, random first order transition theory, and mode coupling theory. In opposition to the freezing transition, the glass transition is continuous in specific volume change and the transition temperature T g is cooling-rate dependent. The Kob-Andersen (KA) model has been investigated extensively since it was first proposed to test mode-coulping theory [START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture. ii. intermediate scattering function and dynamic susceptibility[END_REF]. The KA model is a binary mixture of Lennard-Jones (LJ) particles. Another popular glass-former is the polydisperse LJ systems [START_REF] Wittmer | Vibrations of amorphous, nanometric structures: When does continuum theory apply? EPL[END_REF]. With these LJ glasses, different aspects are studied, such as static yield stress [START_REF] Varnik | A study of the static yield stress in a binary Lennard-Jones glass[END_REF], aging effects [START_REF] Kob | Aging effects in a Lennard-Jones glass[END_REF], fluctuation dissipation ratio [START_REF] Barrat | Fluctuation-dissipation ratio in an aging Lennard-Jones glass[END_REF], single particle jumps [START_REF] Vollmayr-Lee | Single particle jumps in a binary Lennard-Jones system below the glass transition[END_REF], local elasticity map and plasticity [START_REF] Tsamados | Local elasticity map and plasticity in a model Lennard-Jones glass[END_REF]. The examined assumption impact factors of glass transition contain static length scale [START_REF] Biroli | Comparison of static length scales characterizing the glass transition[END_REF], point-to-set length scales [START_REF] Li | Growing point-to-set length scales in Lennard-Jones glass-forming liquids[END_REF], heterogeneity [START_REF] Sillescu | Heterogeneity at the glass transition: a review[END_REF]. As we can see, the elasticity is not enough investigated and we focus our investigations on this topic. As is well known, the shear modulus is zero in liquids and non zero in solids (crystal or glass). This quantity can thus be a signature of the glass transition, as an order parameter. This is for example the criterion used in [START_REF] Klix | Glass elasticity from particle trajectories[END_REF]. In this case, one can question about the nature of the change in G at T g , whether it is a continuous change, or a jump. Since there is still debate on this question, we carry out investigations on KA2d model system with focus on this topic.

In order to examine the above mentioned properties using the stress fluctuation formalism, we choose a few simple model systems in this thesis. The glass transition is especially studied with the KA2d system. The companion system, the monodisperse LJ2d system, which crystallizes, is also studied for the sake of comparison. Two crystals are investigated for their elasticity and viscoelasticity, the first one is the monodisperse LJ3d system which forms a fcc crystal. The second is the AB13 system which forms a cubic complex superlattice. These studies are done by using Molecular dynamics simulations provided by the LAMMPS code [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF]. To a minor extent, monodisperse and polydisperse hard sphere fcc crystals are also studied by constant pressure Monte Carlo simulations [START_REF] Allen | Computer simulation of liquids[END_REF], in order to test the computation of the bulk modulus by volume fluctuations.

The thesis is organized as follows: in chapter 2, we present the methodology used in this thesis. After an overview of classical simulation methods, Monte Carlo and Molecular Dynamics, we discuss the different statistical ensembles that are relevant to our studies. They are microcanonical (NVE), canonical (NVT) and isothermal and isobaric (NPT) ensembles. A large portion of this chapter is devoted to the presentation of the stress fluctuation formalism for the calculation of the elastic constants. Chapter 3 contains our results on static properties of the model systems we studied. Three topics have been considered. The first is the calculation of bulk moduli of monodisperse and polydisperse hard sphere crystals by volume fluctuations, the second is a study of a glass former in two dimensions: 80-20 Kob-Andersen model, and the last topic concerns shear stress fluctuations in simple and complex crystals. Chapter 4 investigates timedependent properties of our systems. Several questions have been examined. The first is the sampling-time dependence of the stress fluctuations, then the shear-stress autocorrelation function and the shear-stress mean-square displacement are considered.

The last part concerns the viscoelastic properties of our model systems, in liquid, glass and crystal states. Conclusions are gathered in chapter 5.

Chapter 2 Methodology

In this chapter we shall introduce the theoretical foundation of our work. First we shall describe briefly the simulation methods, Monte Carlo (MC) and Molecular dynamics (MD) simulations. A brief discussion on the relevant statistical ensembles is given. Then we shall present the method of computing the elastic properties of a model material using the stress fluctuation formalism.

Monte Carlo methods

Monte Carlo methods used in numerical analysis are mainly aimed at calculating multivariable integrals. In molecular simulations, it is used to simulate statistical systems and produce a representative ensemble of configurations to reach thermodynamical quantities without the necessity (or the possibility) of computing them analytically [START_REF] Walter | An introduction to monte carlo methods[END_REF].

Monte Carlo integration

Monte Carlo technique can be a method of integration. Take a simple integral for example

F = x 2 x 1 dxf (x).
(2.1)

This integral can be calculated approximately by

F ′ = x 1 -x 2 N N 1 f (x i ) (2.2)
where the points x i are evenly distributed in the range (x 1 , x 2 ). As long as we select enough points, that is, the number N is big enough, the F ′ will be sufficiently close to the exact F . However, the points can be also selected randomly and in this case the integral can be written as follows

F = x 2 x 1 dx f (x) ρ(x) ρ(x) (2.3)
with ρ(x) an arbitrary probability density function (ρ(x) = 0). Consider conducting lots of trials tr, and every time a random number ξ tr is chosen from the distribution ρ(x) in the interval (x 1 , x 2 ). We obtain that

F = f (ξ tr ) ρ(ξ tr ) trials (2.4)
where the brackets show an average over all trials. The probability density function ρ(x) could be as simple as a uniform distribution (sample mean integration),

ρ(x) = 1 x 2 -x 1 , x 1 ≤ x ≤ x 2 , (2.5)
following this, the integral can be approximated by

F ≈ x 2 -x 1 tr max trmax tr=1 f (ξ tr ), (2.6)
where tr max is the total trial number. We notice that for a one-dimensional integration, this method is not as efficient as many others (Simpson, Gauss, etc). But its use is more indicated for multi-dimensional integrations.

Importance sampling

In Monte Carlo algorithms, lots of trials may give a very small contribution to the Monte Carlo average if we use an uniform sample mean integration method (above).

Importance sampling techniques can solve this problem, by choosing random numbers from a non-uniform distribution [START_REF] Allen | Computer simulation of liquids[END_REF]. They concentrate the samplings in the regions of interest which make important contributions to the integral and hence enhance the overall efficiency. This technique is essential to molecular simulations. For example, in the canonical ensemble (NV T ), if we want to calculate the average of a physical quantity A, we write

A N V T = dΓρ N V T (Γ)A (Γ), (2.7)
where Γ is a point of the phase space, and ρ N V T (Γ) is the canonical distribution.

Comparing with eq.(2.1), we have f = ρ N V T A . In general, the integrand f will be important where ρ N V T is important, and neglectable where ρ N V T is neglectable. By repeating the previous MC integration method we obtain the integral

A N V T = A ρ N V T /ρ trials . (2.8)
if we choose ρ = ρ N V T , we perform an efficient estimate of the integral. Then

A N V T = A trials (2.9)
In order to generate a series of random states at the end of the simulations where every state has appeared with the appropriate probability, a Markov chain of states of the system is set up. This is done using the Metropolis algorithm below.

The Metropolis method

The Metropolis algorithm is used in molecular simulations. It generates a sequence of states according to the equilibrium probability distribution of a given system. The following gives a presentation according to Chandler's book [START_REF] Chandler | Introduction to modern statistical mechanics[END_REF]. We only consider the change of particle's positions. Suppose a state m of the system, and ρ m the equilibrium probability, consistent with the Boltzmann distribution ρ m ∝ exp (-βV m ), (β = 1/k B T , and V m is the potential energy of the system at state m). If T mn is the transition probability (per unit time) that the system transit from state m to n, then the kinetic equation for ρ m writes ρm = n (ρ n T nm -ρ m T mn ) .

(2.10)

At equilibrium, we have ρm = 0 (the equilibrium state is a stationary solution of the kinetic equation). From eq.(2.10), we see that this is the case when

ρ n T nm = ρ m T mn . (2.11)
Eq.(2.11) is referred to as the detailed balance condition. More precisely, we can write

T mn T nm = ρ n ρ m = exp (-β∆V mn ) (2.12)
where ∆V mn = V n -V m . In practice, the Metropolis Monte Carlo (MMC) scheme proposes

T mn =      α nm if ρ n ≥ ρ m , α nm ρn ρm if ρ n < ρ m .
(2.13)

where α mn = α nm being the elements of a symmetric stochastic matrix. Here is an MC step: suppose the system at state m, a trial move is made, by displacing one atom randomly, tempting to bring the system to state n. If ∆V mn ≤ 0, then the move is accepted. If ∆V mn > 0, then exp (-β∆V mn ) is compared to a random number ξ picked up uniformly in [START_REF]nos méthodes de calcul des constantes élastiques pour étudier les transitions de phase solide-solide sous haute pression, ou en fonction de la température. Des potentiels plus généraux, tels que des potentiels à trois corps ou du type EAM devraient être examinée, afin d'étudier les propriétés élastiques de matériaux plus réalistes, surtout en fonction de la température et de la pression. Les fonctions viscoélastiques peuvent être calculées pour de nombreux systèmes de la matière molle, pour des comparaisons avec les résultats expérimentaux[END_REF][START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF]. The move is only accepted when exp(-β∆V mn ) > ξ. In other words, in the Metropolis algorithm, the acceptance probability for a move from state m to state n is

P (m → n) = min (1, exp (-β (V n -V m ))) (2.14)
It is important that the Markov chain generated in such a fashion is ergodic. That means from any state, any other state should be reached after a finite number of trial moves. An MC simulation contains a very large number of MC steps, typically more than a million.

MC simulation of hard spheres

The hard sphere (HS) model is a simple model to represent dense fluids. The HS potential is

βu(r) =      0 if r ≥ σ ∞ if r < σ (2.15)
where σ is the HS diameter. This singular potential leads to a very simple Metropolis MC algorithm. Consider the system be in state m. The state n consists of displacing one particle i a small distance:

r (m) i → r (n) i .
Then, substituting eq.(2.15) in eq.(2.14), the acceptance probability for a move from state m to state n is [START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF] i.e. it is an overlap test. The scheme is easily generalized to HS mixtures, where each species has a diameter σ µ . If the particle i belongs to species µ , its position is r iµ .

P (m → n) =      1 if | r (n) i -r (m) j | ≥ σ for all j = i 0 if | r (n) i -r (m) j | < σ for one j = i, (2.
The MC acceptance probability is then

P (m → n)      = 1 if | r (n) iµ -r (m) jν | ≥ σ µν for all jµ = iν = 0 if | r (n) iµ -r (m) jν | < σ µν for one jµ = iν, (2.17)
where σ µν = (σ µ + σ ν )/2. This is the scheme we used in one of our studies in the next chapter.

To end this section, we point out that the MC method is a powerful tool to study equilibrium properties of a system. On the other hand, it does not follow the Newtonian laws for the particles. That's why the molecular dynamics method (following section) is used when dynamic properties of a system are investigated.

Molecular dynamics

Molecular Dynamics simulation is a technique that computes the equilibrium and transport properties of a classical many-body system. It generates the microscopic information of the investigated system such as particles positions, velocities and so on. This information can be converted to the macroscopic properties of the system by statistical mechanics [START_REF] Allen | Computer simulation of liquids[END_REF]. The motion of the constituent particles conforms to the laws of classical mechanics. In the following, we display our equations in a 3d system.

Equations of motion

Consider a system of N particles which interact with the pair potential u. The forces imposed on particles derive from the potential energy (the system thus conserves the total energy). All the particles move in all directions of the system. The Lagrangian equation can be used to describe the motion of the particles. The equation (for the

variable q k (t)) is d dt (∂L /∂ qk ) -(∂L /∂q k ) = 0, k = 1, • • • 3N (2.18)
where q k are the set of coordinates of all the particles (k = 1, ...3N) and qk are the time derivatives of q k (the momentum of the particles). The Lagrangian function

L ({ qk , q k }) is L = K -V (2.19)
with kinetic energy term K = (1/2) 3N k=1 q2 k /2m and potential energy term V . The total potential is

V = (1/2) N i=1 j =i u(r ij ) = N -1 i=1 j>i u(r ij ), (2.20)
where we have r i = (x i , y i , z i ), x 1 = q 1 , y 1 = q 2 , z 1 = q 3 , etc. The force on the particle i is

F i = -∇ r i V . (2.21)
The force that atom j exerts on atom i is

f ij = -∇ r i u(r ij ).
And for F i we have

thus F i = j =i f ij = -j =i ∇ r i u(r ij ). According to Newton's third law, f ji = -f ij ,
meaning that each pair force needs only to be calculated once. If we take u to be the Lennard-Jones potential, then the pair force is

- → f ij,LJ = 48ǫ σ 2 σ r ij 14 - 1 2 σ r ij 8 - → r ij (2.22)
The equation of motion for the particle i is

m ¨ i r = F i = j =i f ij (2.23)
where the sum is over all N particles except i itself and m is the mass of the particle.

These equations must usually be integrated numerically.

Finite difference methods

In molecular simulations, usually, the equations of motion can only be solved numerically because the force acting on particle i depends not only on its position r i but also on the positions of other particles r j (j = i), and in a non-linear manner (see for example eq.(2.22) for the Lennard-Jones pair force). An integrator is used to calculate the trajectory step by step with a finite time increment δt. All the particles in a system evolve with time. What we want to know is the trajectory of the system during a period of time in the phase space. There are a lot of time integration algorithms to generate the MD trajectory. The most popular ones are the Gear predictor-corrector algorithm and the Verlet algorithm. We first show the Gear predictor-corrector algorithm as in [START_REF] Allen | Computer simulation of liquids[END_REF]. It proceeds in two steps. The first one is the prediction step. Suppose that at time t, the positions and velocities of the particles are known. Thus the accelerations a i (t) = F i /m are known too. Suppose that we have also (a numerical estimate of) the third time derivative of r i at t, b i (t) and the fourth derivative c i (t). For a 5-value predictor-corrector algorithm, the predicted positions and velocities, etc, at time t + δt are given by a simple Taylor expansion

r p i (t + δt) = r i (t) + v i (t)δt + a i (t)δt 2 /2 + b i (t)δt 3 /3! + c i (t)δt 4 /4! v p i (t + δt) = v i (t) + a i (t)δt + b i (t)δt 2 /2 + c i (t)δt 3 /3! a p i (t + δt) = a i (t) + b i (t)δt + c i (t)δt 2 /2 b p i (t + δt) = b i (t) + c i (t)δt (2.24)
The second step is the correction step. We calculate the accelerations a i (t + δt) from the new (predicted) positions { r p i (t + δt))}, then the difference ∆ a i = a i (t + δt) -a i (t). Using this, the corrected quantities are

r c i (t + δt) = r p i (t + δt) + G 0 ∆ a i δt 2 /2 v c i (t + δt) = v p i (t + δt) + G 1 ∆ a i δt/2 b c i (t + δt ) = b p i (t + δt) + G 3 ∆ a i (3/δt) c c i (t + δt) = c p i (t + δt) + G 4 ∆ a i (12/δt 2 ) (2.25)
where, according to Gear, G 0 = 19/120, G 1 = 3/4, G 3 = 1/2 and G 4 = 1/12. The corrected quantities allow to cumulate the statistics of the system, and constitute the starting point of the next MD step.

A more direct and widely used method is the Verlet algorithm. We make sole use of the equation of the motion, no need for the extra numerical coefficients. Its velocity version (velocity Verlet) proceeds also by two steps. It has been shown to be very stable, and is one of the most popular MD schemes. First, from the positions, velocities , and the forces (thus the accelerations) of the particles at time t, we compute the positions at time t + δt and the velocities at time t + δt/2,

r i (t + δt) = r i (t) + v i (t)δt + a i (t)δt 2 /2 v i (t + δt/2) = v i (t) + a i (t)δt/2, (2.26)
then, using the new positions { r i (t + δt))}, the forces (and the accelerations) at time t + δt are computed, and the velocities are advanced to t + δt:

v i (t + δt) = v i (t + δt/2) + a i (t + δt)δt/2. (2.27)
The eqs.(2.26-2.27) constitute one step of the velocity Verlet algorithm. A comparison of the two methods has been made in [START_REF] Allen | Computer simulation of liquids[END_REF]. It appears that if the forces only depend on the positions (not on the velocities) as in our systems, the velocity Verlet method can be preferred, since for relatively large δt, it conserves better the total energy.

Furthermore, it is simple to implement.

Whatever the integrator, our aim is to accumulate the statistic information of the system, and extract macroscopic properties. The results will depend on the ensemble used. The following gives some basic information about statistical ensembles that are relevant to our work.

Statistical ensembles

Statistical physics is used to study the thermodynamics of model materials. From the microscopic movement of particles composing the systems, it provides us with the relevant macroscopic properties . In the macroscopic world, the state of the system is characterized by some macroscopic quantities such as temperature, volume and pressure. In microscopic world, the classic statistical physics views systems as a mechanic system containing a very large number of particles whose movements conform to the Newtonian laws.

In statistical physics, ensembles are regarded as a large set of individual systems which have the same properties and structures in given macroscopic conditions. In our simulations, several ensembles are often encountered, including the microcanonical ensemble, the canonical ensemble and the isothermal-isobaric ensemble.

Microcanonical ensemble

For the microcanonical ensemble, the number of the particles in the system, the volume and the total energy of the system are fixed to specific values (it is the NVE ensemble). This is the natural ensemble for a standard MD simulation presented in the section 2.2. The probability density for this ensemble is proportional to

δ (H (Γ) -E) (2.28)
in this expression, Γ represents the phase space and H (Γ) is the Hamiltonian. The function describes the probability density of obtaining specific states of an N-particle system in a space of volume V which have the desired energy E. The partition function is as follows

Q N V E = Γ δ (H (Γ) -E) (2.29)
its expression for a classical N-particles continuous system is

Q N V E = 1 N! 1 h 3N drdpδ (H (r, p) -E) (2.30)
with h being the Planck constant. The correspondant thermodynamic potential is (minus) entropy:

-S = -k B ln Q N V E . (2.31)
If, in the MD simulations, we were able to compute Q N V E , then according to eq. (2.31) we can get access to the useful thermodynamics of our systems [START_REF] Callen | Thermodynamics and an Introduction to Thermostatistics[END_REF]. For example, for the temperature, we have 1/T = ∂S/∂E| V and for the pressure, we have P/T = ∂S/∂V | E . However, in an MD simulation, we do not compute the partition function

Q N V E
, because we cannot possibly cover sufficiently the phase space in our runs. So the entropy S is not computed. Nevertheless, many thermodynamic quantities can be obtained by simply averaging their microscopic expressions during the run. For example, the average kinetic energy is given by

E K = K M D , (2.32)
leading to, according to the equipartition principle [START_REF] Chandler | Introduction to modern statistical mechanics[END_REF], the temperature

T = 2 3k B (E K /N) (2.33)
Eq.(2.33) is the standard way of calculating the temperature of our systems.

In the NVE ensemble, E is fixed during the MD run. This means that when one starts a simulation, one has the initial configuration, consisting of the initial positions and initial velocities of the particles. If we set the initial velocities according to a target temperature T 0 by using the Maxwell-Boltzmann distribution, the initial positions usually do not correspond to the equilibrium situation. They do fix the total energy, giving E = E 0 . During the MD run of equilibration, the system evolves to the equilibrium state (maximum entropy) corresponding to given (E, ρ), during which there is exchange between the kinetic energy and potential energy. Consequently, the equilibrium temperature will not be the target temperature T 0 , but a different value.

This is a well-known difficulty of the standard MD simulation. It can be overcome by rescaling the velocities regularly during the run (Andersen's method), in order to obtain an equilibrium configuration corresponding to the target temperature, meaning that we run one kind of NVT MD for the equilibrium stage. Afterwards, NVE runs can be carried out to study the system's properties.

As the stress fluctuation formalism which we use is most directly case in the canonical (NVT) ensemble, we give in the following some basis of this ensemble.

Canonical ensemble

For the canonical ensemble, the number of particles in the system and the temperature of the system are specified. The probability density is proportional to

exp (-H (Γ) /k B T ) , (2.34)
the partition function is

Q N V T = Γ exp (-H (Γ) /k B T ) , (2.35)
the expression for an N-particle continuous system is

Q N V T = 1 N! 1 h 3N d N rd N p exp (-H (Γ) /k B T ) (2.36)
The corresponding thermodynamic potential is the Helmholtz free energy

F F = -k B T ln Q N V T . (2.37)
It is also convenient to define the configuration partition function, or the configuration integral:

Z N V T = d N r exp(-V /k B T ). (2.38)
Z N V T is the relevant partition function, when we want to focus on the effects of particle interactions. As we mentioned earlier, the standard MC simulations (see section 2.1) are performed in the NVT ensemble. There is, of course, no question of calculating the partition function Q N V T , for the same reason as our not calculating Q N V E . However, we can compute plenty of physical quantities as MC averages, for example the energy, and its fluctuations. The latter is of course related to the constant-volume heat capacity of the system, as we shall see later.

Although the standard MD simulations are done in the NVE ensemble, nowadays it is common practice to perform them in the NVT ensemble. This can be done in several ways. One of them is the above-mentioned Andersen's method. It however does not correspond to any deterministic dynamics, as opposed to the extended system Nose-Hoover method. In the latter, we consider the system in contact with a thermal

reservoir. An extra degree of freedom is introduced, which represents the reservoir.

Energy is allowed to flow between the system and the reservoir. The introduced degree of freedom is expressed as s and the conjugate momentum is p s . An extra potential energy is related to s as follows

V s = (f + 1)k B T ln s (2.39)
Here f is the number of degrees of freedom and T is the specified temperature. The relevant kinetic energy is

K s = 1 2 Q ṡ2 = p 2 s /2Q (2.40)
Here Q is the thermal inertia parameter. It controls the rate of temperature fluctuations. The Lagrangian of the system is

L s = K + K s -V -V s (2.41)
The equations of motion for the system can be derived

r = f ms 2 - 2 ṡṙ s (2.42)
and

Qs = i m ṙ2 i s - 1 s (f + 1)k B T (2.43)
The Hamiltonian H s for the extended system is conserved

H s = K + K s + V + V s (2.44)
and the density function for the extended system is microcanonical

ρ N V Es (r, p, s, p s ) = δ(H s -E s ) drdpdsdp s δ(H s -E s ) (2.45)
This NVT scheme is used in our simulations. The NVT ensemble is important for our work. We shall show later how the stress fluctuation formalism is derived within this ensemble. Before this, we present another useful ensemble to our work.

Isothermal-isobaric ensemble

The isothermal-isobaric ensemble (NPT ensemble), is convenient for many studies because it fixes the temperature and the pressure of the system. This is similar to many experimental approaches. For this ensemble, the probability density is proportional to

exp (-(H + P V ) /k B T ) (2.46)
The partition function is

Q N P T = Γ V exp (-(H + P V ) /k B T ) = V exp (-P V /k B T ) Q N V T (2.47)
The expression for an N-particle system is

Q N P T = 1 N! 1 h 3N 1 V 0 dV drdp exp (-(H + P V ) /k B T ) (2.48)
The appropriate thermodynamic potential is the Gibbs free energy

G G = -k B T ln Q N P T (2.49)
Again, we do not actually compute G in our simulations. But we do calculate the averages and the fluctuations of physical quantities.

In order to perform MD simulations in fixed pressure P , we use extended system methods. Andersen proposed coupling the system to an external variable V , which is the volume of the simulation box. The coupling imitates the function of a piston on a real system. The kinetic energy for the piston is

K V = 1 2 Q V 2 (2.50)
Here Q is the mass of this piston. The extra potential energy for the system is

V V = P V (2.

51)

Here P is the specified pressure. The kinetic energy and potential energy related to the particles of the system are

V = V V 1 3 s (2.52) K = 1 2 m i v 2 i = 1 2 mV 2 3 i ṡ2 i (2.53) Here r = V 1 3 s and v = V 1 3
ṡ. The Lagrangian of the system is

L V = K + K V -V -V V (2.54)
The equations of motion for the system can be derived

s = f mV 1 3 - 2 ṡ V 3V V = P -P Q (2.55)
Here f is the force and P is the pressure. They are calculated using unscaled coordinates and momenta. For the system, the Hamiltonian H V is conserved.

H V = K + K V + V + V V . (2.56)
Eq.(2.55) corresponds actually to the dynamics of a constant NPH ensemble (H = E + P V is the enthalpy). In order to carry out NPT runs, this scheme is coupled to the Nose-Hoover constant-temperature scheme shown earlier. This is the NPT-MD we used, within LAMMPS [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF].

To perform isothermal-isobaric (NPT) MC simulations, the scheme of section 2.1 is also modified. We recall [START_REF] Allen | Computer simulation of liquids[END_REF] that the configuration integral in this case is,

Z N P T = dV d N r exp (-β(V + P V )) , (2.57)
with V the volume variable. By performing the scaling s = V -1/3 r, we can rewrite Z N P T as

Z N P T = dV V N d N s exp (-β(V + P V )) = dV d N s exp (-β(V + P V ) + N ln V ) .
(2.58)

The new MC scheme then generates states consistent with the probability

ρ N P T ∝ exp (-β(V + P V ) + N ln V ) , (2.59)
by setting H = V +P V -k B T N ln V and by consdering the ensemble {s 1 , • • • , s N , V } as the (3N + 1) variables to which to apply random changes and the Metropolis algorithm

P (m → n) = min (1, exp (-β (H n -H m ))) (2.60)
where

H m = V m + P V m -k B T N ln V m and

Simple averages and fluctuations

Our simulations lead us to calculate various physical quantities. Some are simple averages, other are the fluctuations of relevant quantities which are related to the response functions of our systems. As explained in [START_REF] Allen | Computer simulation of liquids[END_REF], for the same thermodynamic state, simple averages are independent of the ensemble used, whereas the fluctuations depend on the ensemble. For example, if we take a system defined by its density ρ = N/V and its temperature T in the NVT ensemble, at equilibrium, the energy is

given by

E = H N V T . (2.61)
Equivalently, we can run NVE simulations, with the same E, and we shall obtain the (average)equilibrium temperature by eq.(2.33), which is the same as the one of our NVT ensemble. This idea also applies to the calculation of the pressure P. Running (equilibrium) NVT simulations gives the pressure P . If we run NPT simulations, then we shall obtain average volume V . Again, if we have P = P , then V = V . More precisely, here are some common simple averages

E = H E K = K = (1/2) i m i v 2 i E P = V = N -1 i=1 j>i u(r ij ) T = (2/3)E K /Nk B (2.62)
for total energy, kinetic energy, potential energy and temperature respectively.

The pressure P is an important quantity in our studies. Its expression is well-known for pairwise potentials. To show this, we start from the definition P = -(∂F/∂V ) T .

As F = -k B T ln Q N V T , we have

P = k B T Q N V T ∂Q N V T ∂V T (2.63)
From eq.(2.36) and eq.(2.38)we have

Q N V T = A(T )Z N V T = A(T ) d N r exp(-V /k B T ). (2.64)
where A(T ) = (2πmk B T /h 3 ) N /N! is independent of the volume V. By the scaling method r = V 1/3 s, we can rewrite Z N V T as

Z N V T = V N d N s exp(-V /k B T ).
(2.65)

As V = N -1 i=1 j>i u(r ij ) = N -1 i=1 j>i u(V 1/3 s ij ), we have (∂V /∂V ) T = N -1 i=1 j>i u ′ (r ij )s ij (1/3)V -2/3 = (1/3V ) N -1 i=1 j>i u ′ (r ij )r ij (2.66)
Using eq.(2.63) to eq.(2.66), we obtain finally,

P = ρk B T -(1/3V ) N -1 i=1 j>i u ′ (r ij )r ij . (2.67)
This is the famous virial expression for the pressure. The first term is the ideal gas term, P id , the second one is the "excess" term P ex , due to the interactions between the particles.

Another common simple average we compute is the radial distribution function g(r), which characterizes the local structure in our systems. The function ρg(r) represents the probability density of finding a particle at r, given that another is at the origin.

Furthermore, in an isotropic system, g(r) is related to the thermodynamic quantities of the system, E P (isotropic)/N = 2πρ drr 2 g(r)u(r) (2.68)

P (isotropic) = ρk B T -(2π/3)ρ 2 drg(r)u ′ (r)r 2 . (2.69)
This function will be studied in detail in the next chapter. Now we turn to the fluctuations. The fluctuations depend sensibly on the ensemble.

Take for example the energy E. In an NVE ensemble, E does not fluctuate, i.e.

∆E ≡ (H -E) 2 = 0. On the other hand, we have [START_REF] Allen | Computer simulation of liquids[END_REF] ∆E

N V T = k B T 2 C V > 0, (2.70)
where C V = (∂E/∂T ) V is the specific heat capacity of the system. So the fluctuations of E in the NVT ensemble is a way of calculating the constant-volume specific heat capacity. (Functions like C V are often called "response functions"). Another example, more related to our work, is the fluctuations of the volume V . Obviously, in an NVT ensemble, ∆V = 0, whereas in an NPT ensemble, we have [START_REF] Allen | Computer simulation of liquids[END_REF] ∆V

N P T ≡ (V -V ) 2 = V k B T κ T > 0, (2.71)
where κ T = -V -1 (∂V /∂P ) T is the isothermal compressibility of the system. Using the bulk modulus K = 1/κ T , eq.(2.71) can be rewritten as

K = V k B T ∆V N P T . (2.72)
We see that eq.( 2.72) provides a means of computing the bulk modulus. Thus this relation is very useful to our work, as shown in the next chapter.

Calculation methods of elastic properties

Here we present the basic formulas of the calculation of elastic properties by the stress fluctuation formalism. We are in general situations of P = 0, and T = 0. The simulations are equilibrium simulations, i.e. our systems are never actually deformed or submitted to shear stress. On the other hand, we are restricted to pairwise potentials.

Elasticity of solids under pressure

The formalism is most conveniently derived in the canonical ensemble.

We follow the notations in [START_REF] Wallace | Thermoelasticity of stressed materials and comparison of various elastic constants[END_REF]. Given a system (N, V, T ), let X be the initial configuration, x the final configuration. The displacement gradient u αβ is

x α -X α = u α ( X); u αβ = ∂u α /∂X β (2.73)
The conventional strain tensor is

ǫ αβ = 1 2 (u αβ + u βα ) (2.74) The distance R ij = | R i -R j | changes to r ij according to r 2 ij = R ij (1 + ǫ T )(1 + ǫ)R ij = R ij (1 + 2η)R ij (2.75)
where

η = 1 2 ǫ + ǫ T + ǫ T ǫ (2.76)
is the Lagrangian strain. The free energy per unit of (the undeformed) volume f = F/V 0 , can be expanded in powers of η:

f ( x, T ) = f ( X, η αβ , T ) = f ( X, 0, T ) + C αβ η αβ + 1 2 C αβχκ η αβ η χκ +• • • (2.77)
where C αβχκ are the elastic constants.

The stress σ αβ is related to the 1st derivative of f by

σ αβ ( x) = det(h) -1 h αχ ∂f ∂η χκ ( x) h κβ (2.78)
where h = I + ǫ. And the stress-strain relation is:

σ αβ = B αβχκ ǫ χκ (2.79)
where B αβχκ are the Birch coefficients [START_REF] Wallace | Thermoelasticity of stressed materials and comparison of various elastic constants[END_REF]. They are the elastic constants implied in the Hooke's law eq.(2.79).

In case σ αβ ( X) = -P δ αβ (isotropic initial stress), we have

B αβχκ = C αβχκ + P (δ αβ δ χκ -γαχ γβκ -δ ακ δ βχ ) (2.80)
Obviously, if P = 0, we have B = C. This is often assumed in many textbooks on the elasticity of the solids.

By the Voigt notation, we write 1 = xx, 2 = yy, 3 = zz, 4 = yz (or zy), 5 = zx (or xz), 6 = xy (or yx). Thus, the strain tensor reads,

σ =     σ 1 σ 6 σ 5 σ 6 σ 2 σ 4 σ 5 σ 4 σ 3     (2.81)
And the stress tensor reads, by convention

ǫ =     ǫ 1 ǫ 6 /2 ǫ 5 /2 ǫ 6 /2 ǫ 2 ǫ 4 /2 ǫ 5 /2 ǫ 4 /2 ǫ 3     (2.82)
Special case one: 3-dimensional Cubic lattice

For a cubic lattice, the B tensor is written as

B =              B 11 B 12 B 12 0 0 0 B 12 B 11 B 12 0 0 0 B 12 B 12 B 11 0 0 0 0 0 0 B 44 0 0 0 0 0 0 B 44 0 0 0 0 0 0 B 44              , (2.83)
where B 44 = G is the shear modulus. The bulk modulus is given by K = (B 11 +2B 12 )/3.

The matrix inverse to B (the compliance tensor S) is thus (for cubic symmetry):

S = B -1 =              b 1 b 2 b 2 0 0 0 b 2 b 1 b 2 0 0 0 b 2 b 2 b 1 0 0 0 0 0 0 G -1 0 0 0 0 0 0 G -1 0 0 0 0 0 0 G -1              (2.84)
where

b 1 = B 11 + B 12 (B 11 -B 12 )(B 11 + 2B 12 ) b 2 = - B 12 (B 11 -B 12 )(B 11 + 2B 12 )
.

(2.85)

Thus the Young modulus is

E = 1 b 1 = (B
B 11 = d.K d.K + (d -2)E d 2 .K -E B 12 = d.K d.K -E d 2 .K -E B 66 = G.
(2.93)

Isotropic d-dimensional systems

For an isotropic material, the Lamé coefficients λ and µ are usually introduced.

They are related to the C constants by In terms of λ and µ, the bulk modulus K = -V ∂P ∂V | T is given by

C αβχκ = λδ αβ δ χκ + µ(δ αχ δ βκ + δ ακ δ βχ ). ( 2
K = 1 d (B 11 + (d -1)B 12 ) = λ + 2µ/d + (d -2)P/d, (2.96)
the shear modulus G is

G = µ -P. (2.97)
Inversely, we have λ = K -2G/d -P , µ = G + P . Thus, from eq.( 2.95), we can express the relevant Bs in terms of K and G:

B 11 = K + 2(d -1)G/d B 12 = K -2G/d B 66 = G (2.98)
Of course, here we have (B 11 -B 12 )/2 = G = B 66 . In this case, the Young modulus can be obtained from eq.(2.93) and eq.(2.98), leading to

E = 2d 2 K • G d(d -1)K + 2G = 2(µ -P ) d • λ + 2µ + (d -2)P (d -1)λ + 2µ + (d -3)P , (2.99)
and the Poisson ratio, from eq.(2.92), gives

ν = d • K -2G d(d -1)K + 2G = λ + P (d -1)λ + 2µ + (d -3)P . ( 2

Stress fluctuation formalism for crystalline solids

Within the stress fluctuation formalism [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF], we have, supposing the Hamiltonian

H = Σ i p 2 i /(2m) + Σ i<j u(r ij )
, the first term being the kinetic part, the second being the potential energy part, (where Σ i<j is a short-handed notation for Σ i Σ j>i ), and using eqs.(2.75) -(2.80), the elastic constants C expressed as the sum of three terms:

C αβχκ = C K αβχκ + C B αβχκ -C F αβχκ (2.101)
where C K is the kinetic part

C K αβχκ = 2k B T ρ(δ ακ δ βχ + δ αχ δ βκ ) (2.102)
with the density ρ. The Born part is well-known:

C B αβχκ = 1 V i<j u ′′ (r ij ) - u ′ (r ij ) r ij r α ij r β ij r χ ij r κ ij r 2 ij (2.103)
And the fluctuation part is

C F αβχκ = V k B T [ σαβ σχκ -σαβ σχκ ] (2.104)
with the stress tensor element σαβ given by, for pairwise potentials:

σαβ = 1 V i<j u ′ (r ij )r α ij r β ij r ij - i p α i p β i m (2.105)
We note that in eq.(2.104), the fluctuation term is defined as the opposite of that in [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF].

This choice is made for numerical convenience (see next chapter). We can split C F into "kinetic" and "configurational" parts:

C F K αβχκ = 1 V k B T i j p α i p β i p χ j p κ j /m 2 i p α i p β i j p χ j p κ j /m 2 = ρk B T (δ ακ δ βχ + δ αχ δ βκ ) (2.106)
and

C F C αβχκ = 1 V k B T i<j k<l u ′ (r ij )u ′ (r kl )r α ij r β ij r χ kl r κ kl r ij r kl + i<j u ′ (r ij )r α ij r β ij r ij i<j u ′ (r ij )r χ ij r κ ij r ij (2.107)
As we can see, the 1st term in formula (2.107) implies 3-particle and 4-particle distribution functions. We have now

C = C B -C F C + C KK with C KK = C K -C F K , i.e. C KK αβχκ = ρk B T (δ ακ δ βχ + δ αχ δ βκ ) (2.108)
To split completely the "ideal gas" term and the "excess" term (due to the interactions), we can rewrite C as

C αβχκ = C B αβχκ -C F C αβχκ + C KK αβχκ (2.109) C id αβχκ = C KK αβχκ = ρk B T (δ ακ δ βχ + δ αχ δ βκ ) (2.110) C ex αβχκ = C B αβχκ -C F C αβχκ . (2.111)
As we shall see later, in most of our (dense) systems, the contribution of the ideal term (proportional to T ) is negligible. The Born term is usually more important than the fluctuation term, with the latter not at all neglectable, except for simple crystals at

T = 0.
As we always compute the bulk modulus K and the shear modulus G in our systems, we give them specific expressions in the following.

Recalling that K = -V (∂P/∂V ) T = [C 11 + (d -1)C 12 + (d -2)P ] /d, we can show
that after some algebra, the bulk modulus K can be expressed as

K = P + χ V -δP 2 ex ( V k B T ) (2.112)
here the "hypervirial function" [START_REF] Allen | Computer simulation of liquids[END_REF], also referred to as the "Born-Lamé coefficient" [START_REF] Wittmer | Compressibility and pressure correlations in isotropic solids and fluids[END_REF] is

χ = 1 d 2 i<j r ij d (r ij u ′ (r ij )) dr ij (2.113)
and we have P = P id + P ex with P id = ρk B T , and

P ex = P ex = - 1 d • V i<j r ij u ′ (r ij ) (2.114)
is the "virial equation". Thus δP ex = P ex -P ex . We see that eq.( 2.112) is an isotropic expression, as expected, for the bulk modulus.

For the shear modulus G, we have

G = C 66 -P = C B 66 -C F C 66 + C KK 66 -P = C B 66 -C F C 66 -P ex (2.115)
In general, this is not an isotropic property, depending on the symmetry of the crystal.

We shall examine it in more details later.

Coming back to K, and following our recent works, it is convenient to write

K = P id + η A,ex -η F,ex (2.116)
where η A,ex is the elastic bulk modulus corresponding to the excess affine elasticity of the system. More precisely, we have

η A,ex = P ex + χ V (2.117) and η F,ex = δP 2 ex V k B T (2.118)
For pairwise potentials and T > 0, obviously, χ /V and P ex (thus the affine parts) can be expressed from the radial distribution function (RDF) g(r) and the pair potential

u(r), i.e. χ V = ρ 2 2d 2 dr S d g(r) r d (ru ′ (r)) dr (2.119)
where S d is the surface of a d-sphere (S 3 = 4πr 2 , S 2 = 2πr). and

P ex = - ρ 2 2d dr S d g(r) r u ′ (r) (2.120)
In a crystal at T = 0, eq.( 2.113) can be simply written as

χ V = ρ 2d 2 k s k r k d (r k u ′ (r k )) dr k (2.121)
where r k are the coordination shells, and s k are the coordination numbers for each shell.

The summation terminates when r k ≥ r c , the cutoff distance of the pair potential u(r).

For example, the triangular lattice with lattice constant a, we have r k = (a, √ 3a, 2a, ...)

and s k = (6, 6, 6, ...). The lattice constant, related to the density by ρ = 2/ √ 3a 2 , is of course determined by the fixed pressure P , where for T = 0, P = P ex , from eq. (2.114)

is

P ex = - ρ 2d k s k r k u ′ (r k ) (2.122)
In this way, we can compute η A,ex at T = 0. However, η F,ex cannot be computed using (2.118). It will be estimated from extrapolation of low temperature results.

Now some details about G. The computation of G in a cubic crystal follows eq.( 2.115), with the symmetry C 66 = C 55 = C 44 (numerically, it is possible to average these three elements). Next we turn to the case of a 2d triangular lattice. In such a symmetry, we also have a second formula for the shear modulus, denoted by G 2 ,

G 2 = (C 11 -C 12 )/2 -P = (C B 11 -C B 12 -C F C 11 + C F C 12 )/2 -P ex (2.123)
By denoting G 1 the result given by eq.( 2.115), we can write the averaged G as

G 2d = (G 1 + G 2 )/2.
This averaging is applied to the study of low temperature monodisperse 2d LJ system, forming a triangular lattice [START_REF] Abraham | The phases of two-dimensional matter, their transitions, and solidstate stability: a perspective via computer simulation of simple atomic systems[END_REF].

Elastic moduli of the glassy state

The glassy phase is an isotropic solid. Its elastic moduli K and G fully characterize its elasticity. For K, the equation is given by eq.( 2.116), i.e.

K = P id + η A,ex -η F,ex .
As for G, starting from eq.( 2.115), we can express it similarly as

G = µ A,ex -µ F,ex (2.124) with µ A,ex = C B 66 ang -P ex = d d + 2 χ V -P ex = d d + 2 (η A,ex -2P ex ) , (2.125)
with C B 66 ang being the angular average of C B 66 in the d-space, and

µ F,ex = C F C 66 = V k B T δσ 2 xy,ex . (2.126)
We see from eq.( 2.125) that η A,ex and µ A,ex are related. This is not the case for the non-affine part of K and G, η F,ex and µ F,ex . These quantities, involving three-and four-particles correlations, must be computed from the stress fluctuations. As the case of η F,ex , µ F,ex given by eq.( 2.126) can only be computed for T > 0. So its value at T = 0 will be an extrapolation of the results for T > 0. We shall see that, depending on the complexity of the crystal, it is not necessarily zero.

Our glass formers are often binary mixtures. We adopt here eqs.(2.119) and (2.120)

to the case of a binary mixture:

χ V = 2 a=1 2 b=1 ρ a ρ b 2d 2 drS d (r)g ab (r)r d (ru ′ ab (r)) dr (2.127)
where a and b run over the species, and

P ex = - 2 a=1 2 b=1 ρ a ρ b 2d drS d (r)g ab (r)ru ′ ab (r) (2.128)
As for the equation of σαβ , eq.( 2.105), its extension to mixtures is straightforward, since the sums run over each particle of the system.

Elastic moduli of the liquid state

In the liquid state, there is no well-defined displacement field for the particles. What show that P id = ρk B T , and

we know is K = -(1/V )(∂P/
P ex = k B T dL d-1 Z ∂Z ∂L T (2.130)
with Z the canonical configuration partition function, defined by eq. (2.38). This gives the virial expression of P (see eq. 2.67). As the bulk modulus K = -V (∂P/∂V ) T = ρk B T -V (∂P ex /∂V ) T , we need to differentiate once more the second member of eq.( 2.130) with respect to V (thus L). Using again the "trick" eq.( 2.129) and after some algebra, eq.( 2.112) can be proven. We note that the proof can be extended to a rectangular box with fixed aspect ratios.

In a liquid, obviously, G = 0. The fluctuation expression (2.124) is indeed consistent with this fact. As Zwanzig showed [START_REF] Zwanzig | High-frequency elastic moduli of simple fluids[END_REF], we have

δσ 2 xy,ex = σ2 xy,ex -σxy,ex 2 = 1 V 2 i x i F y i k x k F y k = - k B T V 2 i ∂ x i k x k F y k /∂y i = - k B T V 2 i x i ∂ k l>k x kl F y kl /∂y i = k B T V 2 i j>i x 2 ij ∂ 2 u(r ij )/∂y 2 ij = k B T V 2 i j>i x 2 ij u ′ r ij + (u ′′ - u ′ r ij ) y 2 ij r 2 ij .
(2.131)

Note that in eq.( 2.131), we assumed σxy,ex = 0. This is true in a static liquid, and is well verified in simulations. After spherical average (in the d-space), we obtain

δσ 2 xy,ex = k B T dV 2 i j>i r ij u ′ + u ′′ - u ′ r ij r 2 ij d + 2 = k B T d(d + 2)V 2 i j>i (d + 1)r ij u ′ + u ′′ r 2 ij = d k B T (d + 2)V -P ex + χ V (2.132)
Putting eq.( 2.132) in eq.(2.124). One shows indeed G = 0 for liquids. We note that using the radial distribution function g(r), we can also write, for a one-component system

βV δσ 2 xy = ρ 2 2d(d + 2) dr S d g(r) (d + 1)ru ′ (r) + r 2 u ′′ (r) (2.133)
and for mixtures

βV δσ 2 xy = a b ρ a ρ b 2d(d + 2) dr S d g ab (r) (d + 1)ru ′ ab (r) + r 2 u ′′ ab (r) (2.134)
where a and b run over the species. Thus, in liquids, the fluctuation term does not involve higher order (3 or 4) distribution functions.

Additional theoretical aspects

Before closing the section, we point out that the stress fluctuation formulas in subsections 2.5.2 -2.5.4 are obtained (and only valid) within the NVT ensemble, more precisely, the NV γT ensemble, with γ the shear deformation fixed and equal to zero. This is a consequence of the fluctuations being dependent on the ensemble used (see Section 2.4). As the Lebowitz-Percus-Verlet transformation [START_REF] Percus | Ensemble dependence of fluctuations with application to machine computations[END_REF] relates the fluctuations in different (conjugated) ensembles, the elastic moduli can be computed using other ensembles than the one chosen here. This issue has been illustrated in some recent works [START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF][START_REF] Wittmer | Communication: Pressure fluctuations in isotropic solids and fluids[END_REF].

It is also interesting to point out that the stress fluctuation equations for the elastic moduli, both for the liquid and solid states, can be deduced from general thermostatistical considerations, without introducing a local displacement field, as in ref. [START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF]. This has been shown in detail in [START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF] for an isotropic system, where an affine canonical transformation of the positions and the momenta of the particles was made, and the shear modulus obtained via the second derivative of the free energy with respect to the shear deformation, yielding the same equation as eq.(2.124).

Some technical issues 2.6.1 Periodic boundary conditions and minimum images method

In a molecular simulation, the number of particles is always limited. If we only take one simulation box, there will be important surface effects. In order to calculate macroscopic bulk properties of the model system studied, periodic boundary conditions are to be applied [START_REF] Allen | Computer simulation of liquids[END_REF]. In this way, a small number of particles is extended to an infinite system. This can remove surface effects. As the Fig. 2.1 shows, the shaded box represents the original simulation box, while the surrounding boxes are exact copies in every detail. Whenever an atom leaves the simulation cell, it is replaced by another with exactly the same velocity, entering from the opposite cell face so that the number of atoms in the cell is conserved. In simulations with relatively short range interactions, we only need the nearest neighbours around the simulation box. We calculate the interactions between a particle i in the box and all the other particles j within the range of the potential, by applying the minimum image convention [START_REF] Allen | Computer simulation of liquids[END_REF]. For a given range of the potential, with a cut-off distance, the simulation box side must thus be larger than twice the cut-off distance.

Truncation of the pair potential and related corrections

In simulations, it is necessary to truncate the pair potential u(r) at some finite distance r = r c (the cut-off distance), such that u(r) = 0 for r > r c . If there is discontinuity of u(r) at the cut-off distance, then there must be an "impulsive correction"

to the virial expression of the pressure [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]. This can be avoided by shifting the pair potential, imposing u sh (r) = u(r) -u(r c ) for r ≤ r c and u sh (r) = 0 for r > r c . This shift will not alter the forces, so will not affect the particles trajectories. In this thesis, we use indeed shifted pair potentials. However, this shift will not generally prevent the first derivative of the potential, u ′ (r), to be discontinuous at the cut-off. As the Born terms of the elastic constants involve second derivatives of the pair potential, it is necessary to make another "impulsive correction" to these calculations. The issue has been discussed in [START_REF] Xu | Impulsive correction to the elastic moduli obtained using the stress-fluctuation formalism in systems with truncated pair potential[END_REF]. Relevant expressions will be given in Chapter 3.

To conclude this chapter, we point out that we have laid out here the theoretical foundation for our applications in the next chapter (Ch.3). They are all concerned with static properties. The time-dependent properties are all gathered in Chapter 4.

Chapter 3

Static elastic properties

In this chapter, we gather static properties computed using the formalism of the previous chapter. We have carried out three investigations. The first is the computation of the bulk moduli of hard sphere crystals by volume fluctuations. The second is a study of a two-dimensional glass-former: the Kob-Andersen 2d model. In the third part we investigate stress fluctuations in model crystals, one is a simple LJ fcc crystal, the other is an AB13 superlattice, formed by a binary mixture of repulsive particles, in the proportion 1:13 of large and small spheres.

Bulk moduli of hard sphere crystals by volume fluctuations

The bulk modulus is an important property of materials. From its definition,

K = -V ∂P ∂V | T ,
it is possible to compute it by performing a finite (but small) pressure change ∆P , and measure the change in volume ∆V , or vice versa. However, this procedure can imply large numerical errors, especially in the regime where the pressure is high, and where it is highly non-linear with respect to the volume change. Within constant pressure simulations, it is possible to evaluate K from the volume fluctuations during the simulation [START_REF] Allen | Computer simulation of liquids[END_REF], via the relation

K = k B T V δV 2 ,
where k B is the Boltzmann constant, V is the average volume at the given thermodynamic state point (NP T ) and δV 2 = (V -V ) 2 is the variance of the volume (Notice that ... indicates the thermostatistical averages). This idea is shown to give accurate estimates of K, when compared with other more elaborate stress-fluctuation approaches [START_REF] Wittmer | Compressibility and pressure correlations in isotropic solids and fluids[END_REF]. But the cases studied there involve situations with quite weak pressure. Here we apply this method to a hard sphere solid, both monodisperse and polydisperse, up to high pressure, to show first its feasibility and second its predictive power in general physically relevant situations.

The hard sphere crystal phase has been well studied since long time [START_REF] Hoover | Use of computer experiments to locate the melting transition and calculate the entropy in the solid phase[END_REF]. The analytical fit of the equation of state (EOS) βP (ρ), (β = 1/k B T ), proposed by Young and Alder [START_REF] Young | Studies in molecular dynamics. xvii. phase diagrams forsteppotentials in two and three dimensions[END_REF] has proven to accurately cover the whole solid region, and is robust enough when compared to more recent simulations and EOS based on much larger systems [START_REF] Speedy | Pressure and entropy of hard-sphere crystals[END_REF]. Once the EOS is known, the bulk modulus can be calculated by direct derivative of P with respect to the density ρ. There exists another route to compute the elastic constants, via the stress-fluctuation formalism [START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF]. This formalism has been successfully implemented to model glass forming colloidal systems [START_REF] Wittmer | Communication: Pressure fluctuations in isotropic solids and fluids[END_REF]. The special case of hard spheres has been treated by Farago and Kantor [START_REF] Farago | Fluctuation formalism for elastic constants in hardspheres-and-tethers systems[END_REF] for the computation of its elastic constants within the stress-fluctuation formalism. However, the implementation of the method is quite elaborate. Furthermore, the formalism is limited to classical pairwise potentials. This narrows its potential use for some real materials. If only the bulk modulus is needed, then other more direct and simpler simulations can be preferred. In the present work, we show that the widely used constant pressure (NP T ) simulations can indeed produce correctly the bulk modulus K. The method is mostly interesting in cases where an accurate fit of the EOS does not exist, and where the pressure varies quite non-linearly with the density, rendering a finite difference method inaccurate.

Monodisperse hard sphere solid

The case of monodisperse hard sphere solid is well known. It will be used here to check the relevance of the method. We perform constant pressure Monte Carlo simulations [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] of hard spheres solid, consisting of an fcc crystal. The simulation box contains N = 864 hard particles of diameter σ. The reduced pressure P * = βP σ 3 varies from about 12 to 50. As P is proportional to T (at given ρ) for hard spheres, the physical controlling parameter is indeed βP , instead of P and T separately. Although the number of particles is not as large as in some recent simulations [START_REF] Speedy | Pressure and entropy of hard-sphere crystals[END_REF], it is significant enough for the purpose of the present work (cf. analysis in [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]). Starting from an fcc crystal, with an initial volume V 0 , the system is equilibrated during 5 × 10 5 MC cycles under the constant pressure P * 0 . Statistics are further gathered during typically 2 × 10 6 MC cycles. Information such as the average volume V , the volume variance δV 2 , and the volume distribution (V-histogram) are the main results we obtain. This allows us to compare our simulation results with the EOS of ref. [START_REF] Young | Studies in molecular dynamics. xvii. phase diagrams forsteppotentials in two and three dimensions[END_REF]. For the sake of completeness, we display this EOS in the following form.

βP ρ = 3 n=-1 a n y n , (3.1)
where y = ρ cp /ρ -1 with the close-packing density ρ cp = √ 2σ -3 (fcc structure), and

a -1 = 3, a 0 = 2.566, a 1 = 0.55, a 2 = -1.
19 and a 3 = 5.95. The bulk modulus K is thus given by the equation

βK ρ = 3 n=-2 b n y n . (3.2)
The coefficients b n are given in Table 3.1.1. Following our simulation, the first direct result is the average density ρ = N/ V . In Fig. 3.1, we plot ρ * = ρσ 3 vs P * and compare our results with eq.(3.1). We see that indeed, the agreement is excellent. The (reduced) bulk modulus from simulations, K * = βKσ 3 , is plotted vs P * in Fig. 3.2. The simulation results are compared with eq.(3.2). Again, excellent agreement is observed. However, unlike the average density ρ, the bulk modulus K-results bear some non-negligible measurement uncertainty, which is represented by error-bars in Fig. 3.2. They are estimated by standard deviations of K from several independent runs. In Fig. 3.3, the volume histogram is shown for a relatively high pressure P * = 35.5. Indeed, the curve follows a Gaussian distribution, as it should, concerning equilibrium fluctuations [START_REF] Hill | An introduction to statistical thermodynamics[END_REF]. We stress that rather long runs are necessary in order to have Gaussian-like volume histogram, meaning a correct sampling of the volume fluctuations.

Polydisperse system

Having shown the feasibility of the method, we try to make some predictions with it. To do this, we take a polydisperse hard sphere system, containing three components, with the composition 1/3 for each of them, and diameters σ 1 = (1 -ǫ)σ 0 , σ 2 = σ 0 and

= ρ ∂P ∂ρ . (b -2 = -a -1 , b 3 = -2a 3 , and b n = (1 -n)a n -(1 + n)a n+1 for n = -1 to 2). n -2 -1 0 1 2 3 a n 3
σ 3 = (1 + ǫ)σ 0 .
The mean diameter (the first moment of the size-distribution) σ = σ 0 , the polydispersity index I = σ2 /σ 2 -1 = (2/3)ǫ 2 ( σ2 is the second moment of the sizedistribution). This system has been studied previously by one of us, in the framework of the crystallisation [START_REF] Huang | Polydisperse hard sphere mixtures: equations of state and the melting transition[END_REF]. Here we do not consider the fluid phase, because for this phase, the well-known EOS by Mansoori et al. provides a quite accurate description of the polydisperse fluid, at least up to the crystallisation densities. The crystallisation is out of our scope too, because of the fractionation phenomenon accompanying the freezing of polydisperse systems [START_REF] Huang | Polydisperse hard sphere mixtures: equations of state and the melting transition[END_REF]. We thus limit ourselves to the crystal phase, in the same range of the reduced pressure as for the monodisperse system, i.e. 12 ≤ P * ≤ 50 (P * = βP σ3 ). For low to moderate (reduced) pressure, i.e. 1.0 ≤ ρσ 3 ≤ 1.2, an EOS has been proposed by Bartlett [49], based on the hypothesis that the thermodynamics of a polydisperse systems can be mapped to that of a binary mixture (through the identification of the first three moments of the size distribution), and taking advantage of the extensive simulations on binary hard sphere solid make by Kranendonk and

Frenkel [START_REF] Kranendonk | Computer simulation of solid-liquid coexistence in binary hard sphere mixtures[END_REF][START_REF] Kranendonk | Thermodynamic properties of binary hard sphere mixtures[END_REF]. This EOS (eq.( 34) of ref. [START_REF] Bartlett | Thermodynamic properties of polydisperse hard spheres[END_REF]) leads naturally to a approximate analytical expression for the bulk modulus K, which we shall test with our simulation results. The EOS of ref. [START_REF] Bartlett | Thermodynamic properties of polydisperse hard spheres[END_REF] contains 20 double-precision fitting coefficients. In the case of a symmetric size-distribution, the number of coefficients is reduced to 9. This number is still large, compared to the monodisperse case, showing the difficulty of fitting the EOS of polydisperse high density systems. In the following, we display the EOS of ref. [START_REF] Bartlett | Thermodynamic properties of polydisperse hard spheres[END_REF], in the case of symmetric size-distribution (relevant to our system).

βP p σ3 = βP m σ3 + i,j C ij0 (ρσ 3 ) i (10 4 I) j (3.3)
whereP p stands for the pressure of the polydisperse system, P m that of the monodisperse system (given by eq.(3.1) ) as a function of the dimensionless density ρσ 3 , the coefficients C ijk are displayed in Table 3 of ref. [START_REF] Bartlett | Thermodynamic properties of polydisperse hard spheres[END_REF], and I is the above-defined polydispersity index. From eq.(3.3), the bulk modulus K is readily expressed to be

βK p σ3 = βK m σ3 + i,j i C ij0 (ρσ 3 ) i (10 4 I) j (3.4)
where the K p and K m stand respectively for the bulk modulus of the polydisperse and monodisperse systems. K m is given by eq.(3.2). Our simulated system consists always of N = 864 particles, with 288 particles for each species. In this work, a weak polydispersity is taken, i.e. ǫ = 0.02. This will allow us to investigate how ρ and K are sensitive to the polydispersity, at given P * . In Fig. 3.1, we plot ρ * = ρσ 3 vs P * for the polydisperse system. Two points of observation can be made. First, ρ is indeed sensitive to the polydispersity. At a given pressure, the density is sensibly lower in a polydisperse system. This is consistent with the commonly known fact that the crystal is destabilized by polydispersity (see e.g. ref [START_REF] Kranendonk | Computer simulation of solid-liquid coexistence in binary hard sphere mixtures[END_REF]). The second observation is that the EOS of ref. [START_REF] Bartlett | Thermodynamic properties of polydisperse hard spheres[END_REF] coincide with our simulation results in the range of its initial domain of fit. And not surprisingly, it deviates from simulation results outside this domain. So extrapolation of this EOS is not correct. In Fig. 3.2, we have plotted K * = βK σ3 vs P * for the polydisperse system. Here again, we see that the analytical expression eq.(3.4) cannot be extrapolated beyond P * σ 3 ≈ 20. On the other hand, at given P * , the value of K * do not differ much with the monodisperse case. This result can be interpreted by two opposite effects of the polydispersity. Actually, in a (slightly) polydisperse solid, at given P * , V is larger, but δV 2 too. These two effects somehow compensate each other in K. In Fig. for the polydisperse case. The histograms are scaled such that their maximum value is one.

Topical summary

In summary, we tested the feasibility of directly computing the bulk modulus from NP T simulations, by taking the hard sphere solids. The method can cover a very large range of the solid phase, including high pressure solid, where finite difference approach is not appropriate. In the polydisperse case, we show first that the EOS is sensitive to the polydispersity, and secondly, the EOS of Bartlett cannot be extrapolated beyond its initial fitting range (of low to moderate pressures). By computing K up to P * = 50, we demonstrate that these simulations can indeed provide correct K values for high pressure polydisperse solids. Finally, we point out two opposite effects of the polydispersity, rendering K much less sensitive to the polydispersity than the average density does (at a given pressure). This method can find many applications for various systems and situations, for example in the case of structural phase transition in metallic systems, or in the case of self-organisation in complex liquids. The classical 80-20 Kob-Andersen model (KA) [START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture. ii. intermediate scattering function and dynamic susceptibility[END_REF], a binary Lennard-Jones (LJ) bead mixture with a fraction f = 0.8 of large spheres, is an important numerical reference model for the understanding of the glass transition in three dimensions [START_REF] Rowlinson | Liquids and Liquid Mixtures[END_REF][START_REF] Larson | The structure and rheology of complex fluids[END_REF][START_REF] Hansen | Theory of Simple Liquids 3rd Ed[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory: A Mode-Coupling Theory[END_REF][START_REF] Stillinger | Glass transition thermodynamics and kinetics[END_REF]. Its recently proposed two-dimensional (2d) version, called below "KA2d model", has been much less investigated [START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF], despite the high experimental relevance of 2d glass-forming systems [START_REF] Klix | Glass elasticity from particle trajectories[END_REF][START_REF] Cicuta | Shearing or compressing a soft glass in 2d: time-concentration superposition[END_REF][START_REF] Hoffmann | Partial clustering in binary two-dimensional colloidal suspensions[END_REF][START_REF] Ebert | The experimental realization of a two-dimensional colloidal model system[END_REF][START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF][START_REF] Nagamanasa | Direct measurements of growing amorphous order and non-monotonic dynamic correlations in a colloidal glass-former[END_REF]. This may be due to the disappointing finding reported in [START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF] that the standard KA fraction of large spheres did not show a glass transition, but rather a crystal-like low-temperature phase. Various different fractions f have thus subsequently been studied [START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF][START_REF] Berthier | The role of attractive forces in viscous liquids[END_REF][START_REF] Sengupta | Adam-gibbs relation for glass-forming liquids in two, three, and four dimensions[END_REF][START_REF] Marcotte | Nonequilibrium static growing length scales in supercooled liquids on approaching the glass transition[END_REF][START_REF] Hocky | A small subset of normal modes mimics the properties of dynamical heterogeneity in a model supercooled liquid[END_REF][START_REF] Puosi | Time-dependent elastic response to a local shear transformation in amorphous solids[END_REF][START_REF] Dailidonis | Mechanical properties and plasticity of a model glass loaded under stress control[END_REF][START_REF] Rottler | Predicting plasticity with soft vibrational modes: from dislocations to glasses[END_REF][START_REF] Cubuk | Identifying structural flow defects in disordered solids using machine-learning methods[END_REF][START_REF] Flenner | Fundamental differences between glassy dynamics in two and three dimensions[END_REF] with the notable exception of [START_REF] Sengupta | Adam-gibbs relation for glass-forming liquids in two, three, and four dimensions[END_REF] where the classical 80-20 KA model has also briefly been considered in two dimensions. In the present work we readdress some of the results reported in [START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF]. Here are general goals of present work. One central point we want to make is merely that most computational studies [START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF][START_REF] Berthier | The role of attractive forces in viscous liquids[END_REF][START_REF] Sengupta | Adam-gibbs relation for glass-forming liquids in two, three, and four dimensions[END_REF][START_REF] Marcotte | Nonequilibrium static growing length scales in supercooled liquids on approaching the glass transition[END_REF][START_REF] Hocky | A small subset of normal modes mimics the properties of dynamical heterogeneity in a model supercooled liquid[END_REF][START_REF] Puosi | Time-dependent elastic response to a local shear transformation in amorphous solids[END_REF][START_REF] Rottler | Predicting plasticity with soft vibrational modes: from dislocations to glasses[END_REF][START_REF] Cubuk | Identifying structural flow defects in disordered solids using machine-learning methods[END_REF][START_REF] Flenner | Fundamental differences between glassy dynamics in two and three dimensions[END_REF] compare configurations prepared and sampled at an imposed constant volume V and this in addition corresponding to an unrealistically large number density ρ ≈ 1.2(in units of the large sphere's diameter) [START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF]. Obviously, from the experimental and application point of view one should rather control the pressure P (or more generally the stress tensor) while the systems are quenched into the amorphous state. Surprisingly, there exist at present only a few computational studies [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Dailidonis | Mechanical properties and plasticity of a model glass loaded under stress control[END_REF] where the 2d glass transition is investigated allowing the volume to fluctuate at an imposed constant pressure. Following our recent study of the three-dimensional (3d) KA model (KA3d) [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF], the aim of the current work is thus to characterize the standard KA2d model at a moderate pressure corresponding to a much smaller density ρ ≈ 1.0 in the low-temperature limit. We show that under these conditions the KA2d model is in fact a reasonably good glass-former, just as its 3d counterpart,and no indications of a crystallization or other forms of long-range orientational correlations [START_REF] Flenner | Fundamental differences between glassy dynamics in two and three dimensions[END_REF] have been found. As a companion system we investigate in parallel a monodisperse Lennard-Jones (mdLJ) system under the same external constraints. As one expects [START_REF] Abraham | The phases of two-dimensional matter, their transitions, and solidstate stability: a perspective via computer simulation of simple atomic systems[END_REF], this model forms a triangular lattice below a freezing temperature T f . We compare thus various properties of the KA2d model glass with this reference. There exists a cuspsingularity of shear modulus. The present study focuses on simple static and quasistatic thermodynamic properties such as the number density ρ(T ), the compression modulus K(T) or the shear modulus G(T) as a function of the temperature T on both sides of the glass transition temperature T g . The equilibrium shear modulus G is obviously an important order parameter characterizing in general the transition from the liquid/sol (G = 0) to the solid/gel state (G > 0) where the particle permutation symmetry of the liquid state is lost for the time window probed [START_REF] Alexander | Amorphous solids: their structure, lattice dynamics and elasticity[END_REF][START_REF] Ulrich | Elasticity of highly cross-linked random networks[END_REF]. (Please note that all reported shear moduli are quasi-static or transient in this sense.) As in other related numerical studies [START_REF] Barrat | Elastic response of a simple amorphous binary alloy near the glass transition[END_REF][START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF][START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF][START_REF] Xu | Impulsive correction to the elastic moduli obtained using the stress-fluctuation formalism in systems with truncated pair potential[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF][START_REF] Barrat | Microscopic elasticity of complex systems[END_REF] we shall determine the only two relevant elastic moduli K(T) and G(T) using the stress-fluctuation formalism [START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF][START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF][START_REF] Rowlinson | Liquids and Liquid Mixtures[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]. In contradiction to the additive jump discontinuity predicted by the mode-coupling theory [START_REF] Klix | Glass elasticity from particle trajectories[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory: A Mode-Coupling Theory[END_REF][START_REF] Szamel | Emergence of long-range correlations and rigidity at the dynamic glass transition[END_REF][START_REF] Ozawa | Jamming transition and inherent structures of hard spheres and disks[END_REF], we show that the shear modulus G(T) of the KA2d model reveals a continuous cusp-singularity

G(T ) ≈ G g (1 -T /T g ) α for T /T g < 1 (3.5)
with α ≈ 0.6. This result is in qualitative agreement with related numerical [START_REF] Barrat | Elastic response of a simple amorphous binary alloy near the glass transition[END_REF][START_REF] Wittmer | Communication: Pressure fluctuations in isotropic solids and fluids[END_REF] and theoretical [START_REF] Zaccone | Disorder-assisted melting and the glass transition in amorphous solids[END_REF] work, albeit in conflict to some recent experimental work [START_REF] Klix | Glass elasticity from particle trajectories[END_REF][START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF].

We begin by presenting in subsection 3.2.2 the numerical model systems considered and summarize the relevant stress-fluctuation relations for the determination of K(T )

and G(T ). We describe then in subsection 3.2.3 our simulation results and summarize finally this work in subsection 3.2.4.

Algorithmic Details

Model Hamiltonians. For comparison we have studied two soft bead models. As a reference we compute systems of monodisperse Lennard-Jones beads (mdLJ) [START_REF] Allen | Computer simulation of liquids[END_REF], i.e.

we take advantage of a generic model for "simple liquids" [START_REF] Hansen | Theory of Simple Liquids 3rd Ed[END_REF] with a perfectly known 2d phase diagram [START_REF] Abraham | The phases of two-dimensional matter, their transitions, and solidstate stability: a perspective via computer simulation of simple atomic systems[END_REF]. Our system consists of N = 1250 particles interacting through a shifted and truncated LJ potential. The central model of the present study is the 2d version (KA2d) [START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF] of the standard 3d KA model [START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture. ii. intermediate scattering function and dynamic susceptibility[END_REF]. It consists of a binary mixture of LJ particles. We have sampled these model Hamiltonians by means of molecular dynamics (MD) simulations [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] using the LAMMPS code [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] and taking advantage of the standard Nosé-Hoover thermostat and barostat provided by this algorithm.

The temperature coupling constant is set to 10 and the isobaric coupling constant to 25 (with both values in simulation units). A rectangular box with L y = √ 3L x was used in order to be compatible with the triangular crystal phase of the mdLJ model.

For both models the average normal pressure P is kept at a constant P = 2 for all temperatures considered as in [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF]. For the mdLJ model we have carried out both cooling and heating cycles. For the KA2d model the systems have been quenched

with a constant cooling rate of 10 -5 starting from the liquid limit at T = 1. This corresponds to the main rate discussed in [START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF]. We remind that areas of hexagonally crystallized A particles in a matrix of amorphous AB material were observed for exactly this rate and the standard fraction f = 0.8 of large spheres. Table 3.2.1 summarizes some properties of the KA2d model discussed in the next section. We use stressfluctuation relations for elastic moduli. The compression modulus K(T ) and the shear modulus G(T ) indicated in table 3.2.1 have been obtained using the relevant stressfluctuation relations for simulations at constant volume V and constant shear strain γ = 0 [START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF][START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF]. While the systems are quenched imposing a constant mean pressure, i.e. the volume fluctuates while we cool, the elastic moduli are subsequently determined for each T in the NVT-ensemble by switching off the barostat after some tempering [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF]. The compression modulus is given by the Rowlinson stress fluctuation formula [START_REF] Wittmer | Compressibility and pressure correlations in isotropic solids and fluids[END_REF][START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]]

K = η A,ex -η F,ex + P id (3.6)
The first "affine" contribution η A,ex is a sum of moments of first and second derivatives of the pair potentials used. Characterizing the average energy change under an affine dilatation strain it is related to the well-known "hypervirial" [START_REF] Allen | Computer simulation of liquids[END_REF]. See [START_REF] Wittmer | Compressibility and pressure correlations in isotropic solids and fluids[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF] for details.

The second contribution η F,ex stands for the reduced fluctuation η F,ex ≡ βV δ P 2 ex of the instantaneous excess pressure Pex (with β being the inverse temperature) and the third contribution P id to the mean ideal pressure. As in related studies [START_REF] Barrat | Elastic response of a simple amorphous binary alloy near the glass transition[END_REF][START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF][START_REF] Tsamados | Local elasticity map and plasticity in a model Lennard-Jones glass[END_REF][START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF][START_REF] Xu | Impulsive correction to the elastic moduli obtained using the stress-fluctuation formalism in systems with truncated pair potential[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF] the shear modulus G may be obtained most readily using

G = µ A,ex -µ F,ex = µ A -µ F (3.7)
with µ A,ex being the excess contribution to the affine shear elasticity µ A = µ A,id + µ A,ex , a simple average of moments of first and second derivatives of the pair potential characterizing the mean energy under an affine pure shear strain. The second contribution µ F,ex ≡ βV δτ 2 ex stands for the excess contribution to the total shear stress fluctuation µ F = µ F,id + µ F,ex with τex being the excess contribution to the instantaneous shear stress. Since for a shear strain at constant volume the ideal free energy contribution does not change, the explicit kinetic energy contributions must be irrelevant for G. (An ideal gas can not elastically support a finite shear stress.) As one thus expects, the kinetic contributions µ A,id = µ F,id = P id cancel and can be dropped when G is determined using the first equation of eq.(3.7). We note finally that since second derivatives are relevant for η A,ex and µ A,ex , impulsive corrections must be taken into account for the truncated and shifted potentials considered in the present work as stressed in [START_REF] Xu | Impulsive correction to the elastic moduli obtained using the stress-fluctuation formalism in systems with truncated pair potential[END_REF].

Otherwise the shear modulus G, determined using eq.(3.7), does not vanish in the liquid limit as it must. G are of order one in these units, all other properties are known to higher precision.

Numerical results

Specific

Boltzmann's constant is set to unity and all properties are given in LJ units [START_REF] Allen | Computer simulation of liquids[END_REF]. 3.2.1) and the specific volume v(T ) = 1/ρ(T ) per particle must decrease. This is shown in Fig. 3.4 for both models. A jump singularity is observed for v(T ) at the freezing temperature T f ≈ 0.62 of the mdLJ model. This jump is of course a consequence of the first order phase transition from the liquid phase to the crystalline state [START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF][START_REF] Callen | Thermodynamics and an Introduction to Thermostatistics[END_REF][START_REF] Abraham | The phases of two-dimensional matter, their transitions, and solidstate stability: a perspective via computer simulation of simple atomic systems[END_REF]. At variance with this, the KA2d model changes more gradually with linear slopes fitting reasonably both the low-and the high-T limits. The observed (more or less) sudden change of the tangent slopes at the glass transition may be used to operationally define T g [START_REF] Larson | The structure and rheology of complex fluids[END_REF] by matching the indicated two lines as done in various recent numerical studies [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF]. Using this calorimetric criterion we obtain

T ρ L x R P id η A,ex K µ A,ex G 0.
T g ≈ 0.
35 for the KA2d model. This is slightly smaller than the corresponding value

T g ≈ 0.
41 for the KA3d model at P = 1 [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF]. We remind that for the poly-disperse purely repulsive LJ system (pdLJ) considered in [START_REF] Wittmer | Vibrations of amorphous, nanometric structures: When does continuum theory apply? EPL[END_REF][START_REF] Xu | Impulsive correction to the elastic moduli obtained using the stress-fluctuation formalism in systems with truncated pair potential[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF][START_REF] Tanguy | Continuum limit of amorphous elastic bodies: A finite-size study of low-frequency harmonic vibrations[END_REF] the same criterion yields a glass transition temperature of only T g ≈ 0.26 at the same pressure P = 2.

In qualitative agreement with [START_REF] Berthier | The role of attractive forces in viscous liquids[END_REF] this suggests that, having a much larger T g , the KA2d model should be a more promising numerical model for investigations of the glass transition in two dimensions.

Radial distribution functions. One of the conclusions of [START_REF] Brüning | Glass transitions in one-, two-, three-, and four-dimensional binary Lennard-Jones systems[END_REF] is that the KA2d model with a fraction f = 0.8 of large spheres should crystallize instead of forming an amorphous glass at low temperatures. In order to clarify this issue for our constant pressure systems, we compare in the main panel of Fig. This criterion suggests a glass transition temperature T g ≈ 0.32 which is consistent, albeit slightly smaller than the one obtained above. The excess normal pressure fluctuation contribution η F,ex (T) shows a striking peak at the freezing temperature T f ≈ 0.62 of the mdLJ model and vanishes then rigorously for T → 0. Since in this limit the ideal pressure contribution P id vanishes also, the compression modulus K(T = 0) is completely determined from the affine contribution η A,ex . This is expected for sim- and the fact that the shear modulus of a liquid must vanish. µ F,ex (T ) shows a maximum at, respectively, T f ≈ 0.62 for the liquid-solid and T g ≈ 0.32 for the liquid-glass transition. While µ F,ex (T ) vanishes rigorously in the low-T limit for the mdLJ model due to the simple triangular lattice adopted [START_REF] Barrat | Microscopic elasticity of complex systems[END_REF], µ F,ex (T ) is seen to level off for the KA2d model (horizontal arrow). In agreement with the discussion of the compression modulus K above, this confirms that an affine shear does neither lead for an amorphous body to a mechanically stable configuration and a finite energy per volume µ F,ex (T = 0) ≈ µ A,ex (T = 0)/2 of the strained ground state system can be relaxed by non-affine displacements as discussed elsewhere [START_REF] Wittmer | Vibrations of amorphous, nanometric structures: When does continuum theory apply? EPL[END_REF][START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF][START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF][START_REF] Xu | Impulsive correction to the elastic moduli obtained using the stress-fluctuation formalism in systems with truncated pair potential[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF][START_REF] Barrat | Microscopic elasticity of complex systems[END_REF][START_REF] Tanguy | Continuum limit of amorphous elastic bodies: A finite-size study of low-frequency harmonic vibrations[END_REF]. As one and G g ≈ 16.5, T g ≈ 0.32 and α ≈ 0.6 seems to fit our data. This is better seen from the zoom on the glass transition region given in the inset of Fig. 3.8 (b). We remind that a slightly lower exponent α ≈ 1/2 was fitted for the KA3d model in our previous work [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF] as may be seen from the rescaled data (diamonds) included in the inset.

Whether the weak difference of the exponents reflects a fundamental effect due to the different spatial dimensions, as suggested by the recent work of Flenner and Szamel [START_REF] Flenner | Fundamental differences between glassy dynamics in two and three dimensions[END_REF], or whether it should merely be attributed to the small system sizes used, can currently not be answered. Larger systems and, more importantly, longer trajectories are warranted to clarify this question.

Sampling time dependence of stress fluctuations. As already emphasized elsewhere [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF], the convergence with sampling time of the stress fluctuation relations for the elastic moduli eqs. (3.6) and (3.7) may be slow even for permanent elastic networks and one thus has to check whether sufficiently large trajectories have been used to determine reliable long-time estimates. The reason for this is that, while the simple means η A,ex and µ A,ex converge essentially immediately, the stress fluctuation contributions η F,ex and µ F,ex become sampling-time dependent properties ηF,ex (∆t) and μF,ex (∆t) if sampled over trajectories of finite lengths ∆t. Please note that this has apriori nothing to do with an insufficient equilibration or ageing of the systems, but stems from the fact that the stress-fluctuations simply need time to explore the available phase space. Our results for this topic will be shown in the next chapter, with other time-dependent properties. We note here briefly that indeed, a strong dependence has been observed, especially around the glass transition.

Topical conclusions

To summarize, we have investigated by means of molecular dynamics simulations the glass transition of the two-dimensional version of the Kob-Andersen model with a Discussions. We note here that all the presented plateau values G(T), especially for the KA2d model [START_REF] Sausset | Do solids flow?[END_REF], correspond strictly speaking to intermediate (albeit perhaps rather large) shoulders of µ A,ex -μF,ex (∆t). Note that G(T ) is only finite below T g ≈ 0.3 since on the time scale ∆t available numerically the particle permutation symmetry of the liquid state is lost [START_REF] Alexander | Amorphous solids: their structure, lattice dynamics and elasticity[END_REF]. As discussed in [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF], if trajectories of increasing lengths are compared, G(T ) becomes monotonously smaller, T g is shifted to smaller values and the glass transition becomes sharper. The phenomenological exponent α(∆t) characterizing the cusp-singularity eq.(3.5) should thus decrease somewhat with the sampling time ∆t. It is thus possible that the jump discontinuity (α → 0) predicted by mode coupling theory may be recovered by appropriately extrapolating α(∆t) for larger sampling times ∆t.

It is of course necessary to investigate in future work larger systems. Preliminary tests with systems comprising N = 10000 particles suggest, however, that all the reported properties are system size independent, especially no indication for a crystallization has been observed. Even more importantly, as already stated above, future studies should focus on sampling time effects. This should also allow to verify whether the amorphous low-temperature state (Fig. 3.5) remains stable and no additional static long-range correlations (of perhaps Kosterlitz-Thouless-type as suggested in [START_REF] Flenner | Fundamental differences between glassy dynamics in two and three dimensions[END_REF]) may become relevant. In view of recent experimental trajectory measurements on 2d glass-forming colloidal systems [START_REF] Klix | Glass elasticity from particle trajectories[END_REF][START_REF] Ebert | The experimental realization of a two-dimensional colloidal model system[END_REF][START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF][START_REF] Nagamanasa | Direct measurements of growing amorphous order and non-monotonic dynamic correlations in a colloidal glass-former[END_REF] it is also important to analyse the Fourier transform of the displacement fields and other wavevector dependent static and dynamical properties [START_REF] Flenner | Fundamental differences between glassy dynamics in two and three dimensions[END_REF].

Shear stress fluctuations in simple and complex crystals

In this section we calculate the elastic constants of two model crystals, using the stress fluctuation formalism. The two model are the Lennard-Jones system, and a binary mixture of repulsive particles in 1:13 portion, the AB13 system. Our focus will be mainly to point out the importance of the fluctuation terms in crystals, since this term is usually neglected in classical calculations of the elasticity of crystals [START_REF] Ashcroft | Solid state physics[END_REF].

Lennard-Jones model system

The Lennard-Jones model is the most popular one in molecular simulations. At low temperatures, it forms a fcc crystal [START_REF] Hansen | Phase transitions of the Lennard-Jones system[END_REF]. It is the one used by Squire et al. for the first calculation of the elastic constants in solids by the stress fluctuation formalism [START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF]. We revisit here the model by considering it as a simple crystal.

Simulation details

In our molecular dynamics simulation, we use NPT and NVT ensembles. The Lennard-Jones potential between a pair of particles at a distance r is

u LJ (r) = 4ǫ σ r 12 - σ r 6 . (3.8)
In practice, a shifted version is used:

u LJsh (r) =      4ǫ σ r 12 -σ r 6 -u 0 , r < r c 0, r ≥ r c (3.9)
where r c is the cutoff distance of the potential, and u 0 = u LJ (r c ). All units are reduced by σ, ǫ and m (m being the mass of a particle). We chose r c = 2.5 in our simulations. For the heating or cooling or equilibrium processes, NPT ensemble is used while for the statistical calculations NVT ensemble is used. At the beginning, the configuration of our system is obtained at T = 0.4. Then we have done two types of work including cooling and heating. The processes are as follows : First, cooling or heating the system. We cool or heat the system from the initial configuration (T=0.4)

to the next temperature (for example, T = 0.3 or T=0.5) using NPT ensemble and the cooling or heating rate is 10 -5 . The MD integration time step is 0.005 in reduced units. Second, we equilibrate the system by 2 × 10 6 timesteps using NPT ensemble. In order to keep the pressure at the fixed value (here 1.0), we calculate the average volume during the equilibrium process and change the volume of the last configuration of the equilibrium process to this average volume before running again 1 ×10 6 timesteps. The MD integration time step is still 0.005 in reduced units. Third, we collect the statistics by running 2 × 10 6 timesteps using NVT ensemble. The MD integration time step is 0.0025 in reduced units.

On the one hand, the cooling process of the configurations as the Table 3.3.1 shows.

At every temperature point, the same three steps are repeated. On the other hand, we heat the configuration at T = 0.4 to the configuration at T = 0.5 and the heating rate is 10 -5 , then we reproduce the same three steps as above. The same run is done for other temperatures as the Table 3.3.1 shows. By checking the specific volume (see below), we see that the melting temperature T m is located at 0.8 < T m < 0.85. Thus more temperature points are examined in this range.

Numerical results

Specific volume. We first examine the evolution of the specific volume v = V /N as a function of the temperature. From the graph Fig. 3.9, we can see that there is a Table 3. [START_REF] Xu | Impulsive correction to the elastic moduli obtained using the stress-fluctuation formalism in systems with truncated pair potential[END_REF]. The impulsive correction formula is as follows :

C B αβχκ = C Bcut αβχκ + ∆C B αβχκ , (3.10)
where C Bcut αβχκ are results without impulsive corrections, and ∆C B αβχκ are the impulsive corrections. For a cubic crystal, they can be written, in an "isotropic" approximation, as

∆C B 44 = ∆µ B , ∆C B 11 = 3∆µ B , ∆C B 12 = ∆µ B (3.11)
where

∆µ B = - 8π 3/2 5Γ(3/2) (ρσ 3 ) 2 ǫ σ 3 f (r c )g(r c ) (3.12) with f (r c ) = [1 -(r 0 /r c ) 6
]/r 3 with r 0 = 2 1/6 at which the potential u LJ is minimum.

The Γ(3/2) is a gamma function and equal to √ π/2. f (r c ) is the force strength between particles at distance r c , the cutoff distance. g(r c ) is the radial pair distribution function at cutoff distance.

Our results for the elastic constants (in the solid range) are gathered in Fig. 3.11, where we have plotted the elastic constants C 11 , C 12 and C 44 as functions of the temperature, as well as their fluctuation contributions. Whereas the constants themselves increase when the temperature is lowered, the opposite trend is observed for their fluctuation contributions. We can see that all the fluctuation terms go to 0 when the temperature decreases to 0. This result is obtained by making a linear extrapolation of our results at low temperatures. When T = 0, according to the formulas given in the previous section, the Born terms of the elastic constants are equal to the value of the elastic constants. The results from the linear extrapolation for the Born terms at low temperatures and the direct calculation at T = 0 are almost the same. Our numerical results are also displayed in Table 3.3.2. We have compared our results with some previous (few) calculations, for example in [START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF]. The agreement is excellent.

The shear modulus G can be computed by G = µ -P = C 44 -P . From Fig. 3.12, we can see that the shear modulus of the system decreases linearly when temperature increases from 0.0 to 0.83, then it suddenly becomes zero at T = 0.84, which corresponds Table 3.3.2 Elastic constants and the relevant fluctuation terms of elastic constants at different temperatures for LJ system (LJ3d) in solid state (fcc crystal). We recall that P = 1 throughout, and the melting temperature is T m = 0.84.

Fluctuation terms

Elastic constants The pressure is P = 1. Also plotted is the second "shear" modulus G ′ defined by This means that the LJ fcc crystal is quite anisotropic. The bulk modulus is calculated by K = (C 11 + 2C 12 + P )/3. We have plotted in Fig. 3.13 its evolution with the temperature, both at the solid and the liquid states. We can see that the bulk modulus for the system decreases almost linearly when the temperature increases. But the slope are quite different in each state, being strong in the solid, and quite weak in the liquid.

T C11 f C12 f C44 f C11 C12 C44
G ′ = (C 11 -C 12 )/2 -P .
At T = 0, we have K ≈ 80, whereas at T = 0.8 (solid), K ≈ 20. In the liquid side, T > 0.84, K does not exceed 10. At the transition, K is not well defined, so its numerical value in the plot may be accidental. Our numerical results are shown in system using the stress fluctuation formalism. These results will serve as references when we tackle a more complex crystal, the AB13 superlattice.

AB13 model system

The superlattice formed by a binary mixture of hard spheres with a diameter ratio of σ B /σ A = 0.58, with the composition of one A-particle for 13 B-particles is a quite complex crystal [START_REF] Xu | A density functional study of superlattice formation in binary hard-sphere mixtures[END_REF]. The A particles form a simple cubic lattice. Inside each cube there is a centered icosahedra of 13 B particles. Because of a rotational freedom, each unit cell contains 8 cubes, that is 112 particles. This structure was found by experiments, theory and simulations of colloids, at high packing fraction. We use here this crystal as an example of complex crystals. In Figs.

Simulation details

Again, in our molecular dynamics simulations, we use NPT and NVT ensembles.

The potential between a pair of particles i and j, of species α and β, is

u αβ (r iα,jβ ) =      4ǫ αβ σ αβ r iα,jβ 12 - σ αβ r iα,jβ 6 + ǫ αβ 0 < r iα,jβ < r 0αβ 0, r iα,jβ > r 0αβ (3.13)
where r 0αβ = 2 1/6 σ αβ is the position of the minimum of the LJ pair potential. Thus the pair potential is purely repulsive. Both the potential and its first derivative are zero We follow the same procedure as the case of LJ system. Here the starting temperature is T = 1 (solid). It is cooled down to T = 0, and heated up to T = 3.5. Between 0 and 3.5, several temperatures have been examined (see Table 3.3.4). The pressure was always fixed at the same value. So the heating and cooling are done using the NPT ensemble (rate 10 -5 ). The equilibration is done using NPT as well. The statistics are gathered using the NVT ensemble.

Numerical results

We present the results in the same order as the LJ system.

First, we look at the specific volume. From the graph Fig. 3.17, we can see that there is a jump point along the curve v(T ), which implies the phase transition from crystal to liquid for AB13 model system. The melting temperature is about T m = 2.1. Fig. 3.17 Specific volume (v = V /N) as a function of temperature for AB13 system. The pressure is fixed at P = 114. All units refer to the large spheres. Next, we calculate the RDF for a few temperatures. Fig. 3.18 shows the RDFs at very low temperature, T = 0.1. We see that the curves of different pairs have pronounced peaks, and they subsist at long distances. This means that the system is well structured (low temperature crystal).

In graph Fig. 3.19, the radial distribution functions for two higher temperatures are displayed . When the temperature is 2.0, the number of peaks reduces and the peaks become wider and lower than the T = 0.1 case, which reflects some disorder in the system. But still, the RDFs are crystal-like. So we are close to the melting. When the temperature is 2.1, the radial distribution functions of A particle and B particle become about 1 at long range, which is the feature of a liquid. As a result, the system undergoes the melting phase transition at this temperature. This is consistent with the specific volume result. Now we characterize the elastic constants for the solid state. As before, we have three independent constants, the system being a cubic superlattice. From the graph Fig. 3.20, we can see that all the elastic constants for AB13 model decrease linearly when the temperature increases. Their values are much higher than the LJ system, because of the very high pressure here. On the other hand, we have, as previously, C 12 The shear modulus is of course examined. We recall that G = C 44 -P , and the alternative shear modulus is G ′ = (C 11 -C 12 )/2 -P . From Fig. 3.22, we see that the shear modulus for AB13 system decreases when temperature increases from zero to 2.0. Then it becomes zero when temperature is above 2.1, which means that the melting transition happens at this temperature. This is again consistent with the specific volume behaviour. If we compare G and G ′ , we see that they are very close, with G slightly larger. This shows that the AB13 crystal is much more isotropic than with µ A being the affine modulus, µ F the fluctuation of the stress. We have, by definition

µ F = βV (τ -τ ) 2 = βV τ 2 -βV τ 2 = μF -µ * (4.2)
where μF = βV τ 2 , and µ * = βV τ 2 = βV τ 2 . Now we show that µ F depends on the sampling time used. In order to explain this and simplify the notations, we define a physical observable a, with instantaneous values â. If a is measured over a time window ∆t, meaning n records at interval δt, then we can write the time-average of a

as ā = 1 n n i=1 âi (4.3)
where âi = â(t i ) with t i = (i -1)δt for i ranging from 1 to n. If we have m trajectories of length n, then the expectation value of a is the ensemble average given by

E(a) = ā = 1 m m j=1 āj (4.4)
where āj denotes the time-average resulting from the number-j trajectory. For many observables, the ensemble and the time averages commute. They are called "simple averages". Examples are the pressure P , the average shear stress τ , the affine shear constant µ A , and the mean value of τ 2 denoted by μF . On the other hand, the fluctuations of the observables do depend on ∆t. They are thus not simple averages. To show this, we imagine a long simulation run, consisting on N times δt. This trajectory can be divided into m sub-trajectories each of length ∆t. The expectation value for the fluctuations of a is

F l(a) = â2 -ā2 = 1 m m j=1 (â 2 j -ā2 j ). (4.5)
Obviously, F l(a) will depend on ∆t (because the second term in eq. (4.5) depends on ∆t). It is actually a monotonically increasing function of ∆t, starting from zero when ∆t = 0. In practice, when we compute µ F , we have â = √ βV τ . Thus

µ F (∆t) = βV τ 2 -τ 2 (4.6)
is an increasing function of ∆t (see Fig. Its statistical value can be estimated when µ F (∆t) approaches its plateau-value asymptotically. The real relaxation time of the system can be defined by t * such that the plateau-value as close as 1%. Fig. 4.2 shows our results for LJ2d and KA2d systems. We can observe that the time to reach the plateau depends on the system, and on the thermodynamic state. For the LJ2d system, t * is quite short (about 100).

For the KA2d system, a binary mixture, it is certainly longer (more than 500 for the temperatures investigated).

Our results for 3d systems are gathered in Figs. (4.3-4.4). In these systems, t * appears to be quite short, either for the crystal or the liquids states. For AB13, it is of the order 10, for LJ3d system, about 100. Here the mixture reaches the plateau faster than the one-component system. This may be explained by the important sizedisparity between the spheres, so that the small one have more freedom to move and help relax the system faster.

From eq. (4.1) and eq. (4.6), we can write, for the shear modulus

G(∆t) = µ A -µ F (∆t) = µ A,ex -µ F,ex (∆t), (4.7)
as a function of the sampling time ∆t. around T g , than for the LJ2d model. One thus has to check whether sufficiently large trajectories have been used to determine reliable long-time estimates. Please note that this has apriori nothing to do with an insufficient equilibration or ageing of the systems, but stems from the fact that the stress-fluctuations simply need time to explore the available phase space. The elastic moduli G(T ) and K(T ) presented previously thus stem from reasonable plateau values as one expects.

Shear-stress autocorrelation function

We study the shear-stress auto-correlation function (SACF) of the shear stress, defined by, as in [START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF],

C(t) = β V δτ (t)δτ (0) . (4.8)
This function is important. It allows a better theoretical understanding of the dynamics in liquids [START_REF] Balucani | Dynamics of the liquid state[END_REF], it leads to measurable quantities, such as the viscosity [START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF], and to the shear modulus [START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF].

Two comments on C(t) can be made. The first is the way we compute this cor- relation function. It is important to remind that the time-translational invariance is always assumed. This means that we have

C(t) = β V δτ (t + t 0 )δτ (t 0 ) , (4.9)
where the horizontal bar denotes the gliding average over the time-origins t 0 taken within the trajectory. The second comment remarks that C(t) is a correlation of the stress fluctuations. Thus it depends on the sampling time ∆t. We can designate it as C(t, ∆t). What we expect is that if ∆t << t * (t * is a typical relaxation time), then C(t) is not correctly evaluated. In practice, we divide a long trajectory into m sub-trajectories of ∆t, (j = 1, ...m). Thus the function C(t) obtained from the jth sub-trajectory is

C j (t) = β V δτ (t + t 0 )δτ (t 0 ) j = β V τ (t + t 0 )τ (t 0 ) j -µ * j (4.10)
where µ * = βV τ 2 is obviously ∆t-dependent, since this time-average is taken within the time-interval ∆t. The final result C(t) is the ensemble average of C j (t):

C(t) = 1 m m j=1 C j (t) = 1 m m j=1 Cj (t) -µ * j (4.11)
where C(t) = β V τ (t + t 0 )τ (t 0 ) is the unshifted shear-stress SACF. At given ∆t, we can commute the time and ensemble averages for C(t), it is thus a simple average.

Details for C(t). We denote f (t) = √ βV τ (t), the correlation time t k = kδt, the correlation function F (k). The time-average in j-th sub-trajectory is

F (k) j = [f (t 1,j )f (t 1+k,j ) + ... + f (t n-k,j )f (t n,j )]/(n -k) (0 ≤ k ≤ n -1) (4.12)
and the final correlation function is, by ensemble average,

F (k) = j F (k) j /m = j i f (t i,j )f (t i+k,j )/[(n -k) * m], (4.13)
as the sums commute. C is indeed a simple average. On the other hand, its precision can depend on ∆t. To show this, we suppose that we start with ∆t 0 , corresponding to n 0 data points, coming from one of the m 0 sub-trajectories of a long run of N records.

Thus, we have

F 0 (k) = j i f (t i,j )f (t i+k,j )/[(n 0 -k) * m 0 ]. (4.14)
If we set now m = m 0 * n s segments, each having n = n 0 * n s data points, then the new estimate of F (k) is

F s (k) = j i f (t i,j )f (t i+k,j ) (n -k) * m = j i f (t i,j )f (t i+k,j ) [(n 0 -k/n s ) * m 0 ] . (4.15)
Obviously, except for k = 0, we have F s (k) different of F 0 (k). But the difference can be numerically neglectable, if we ensure that for each k, there are at least 10 2 data for the statistics, between the time and ensemble averages.

On the other hand, µ * is not a simple average. We have, for a given ∆t,

µ * = j µ * j m = j i f (t i,j )/n 2 /m, (4.16)
showing that the sum over i (n terms) and the one over j (m terms) are not commutable. Furthermore, µ * depends strongly on ∆t, being a maximum for ∆t = 1.δt, and decreasing to its plateau value for n >> 1.

From the above discussion, it is clear that the SACF, for a finite ∆t, should be noted C(t, ∆t). However, we can drop the second argument if ∆t is large enough (∆t >> t * ), after the function has converged.

Fig. 4.6 -4.9 show results on this topic. We see that indeed, they all depend on ∆t. But the convergence can be reached after a reasonable sampling time, of the order ∆t ≈ 100, for the state points explored. Another observation is that clearly, the crystals (LJ2d, LJ3d and AB13) show strong non monotonic SACF at short times, whereas in the liquid, the function is a monotonic decreasing function of time. In the low temperature glass (KA2d), the function is also monotonic, although less smooth as in the liquid. Another general feature is that the function goes to zero at long times, for all systems and all states. This means that it is impossible to predict the shear modulus of a solid from this limit, as some authors believe [START_REF] Klix | Glass elasticity from particle trajectories[END_REF]. The sequence µ F (∆t) can be generated by C(t), as shown in [START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF]. We have (for more details, see the next section),

µ F (∆t) µ F = 1 -(2/∆t 2 ) 1 µ F ∆t 0 dt(∆t -t)C(t) ≡ 1 -y(∆t) (4.17)
We can thus define a time θ such that µ F (θ) is very close to µ F , for example at 95%, or y(θ) = 0.05. In Fig. 4.10, the SACF C(t) and the function y(∆t) of the monodisperse LJ2d system are presented, at the liquid and solid states, respectively. We see that for both states, the C(t) curves decay to 0 quite quickly, in about one time unit. The "real" relaxation time of the SACF, as suggested in ref [START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF], could be defined by the time for which y(∆t) reaches 0.05 (the blue line in the inset), corresponding here to about 5 time units for the states explored. The functions, however, are much longer lasting in the glass (the KA2d model), near T g , as shown in Fig. 4.11. A good estimate of the SACF relaxation time can again be obtained by the method suggested in ref [START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF].

Using this definition, we obtained results displayed in Table 4.1, where we see indeed that the relaxation in the glass (near T g ) is about 500, compared to (about) 5 in the crystal.

The shear-stress mean-square displacement

In [START_REF] Wittmer | Simple average expression for shearstress relaxation modulus[END_REF], we introduced the "mean square displacement" (MSD) of the shear-stress, defined by

h(t) = βV [τ (t) -τ (0)] 2 /2. (4.18)
It is straightforward to show that

h(t) = C(0) -C(t) = C(0) -C(t) (4.19)
As C(0) = μF and C(0) = µ F , we can invert eq. (4.19) to write

C(t) = μF -h(t); C(t) = µ F -h(t). (4.20)
Now we show that h(t) is a simple average. We have, for a given ∆t, where the large bar denotes the gliding average over t 0 within the a ∆t package. In detail, for the j-th sub-trajectory, we have

h(t, ∆t) = βV [τ (t + t 0 ) -τ (t 0 )] 2 /2,
h(k) j = [((f (t 1+k,j ) -f (t 1,j )) 2 + ... + (f (t n,j ) -f (t n-k,j )) 2 ]/(n -k) (0 ≤ k ≤ n -1), (4.22) thus h(k) = m j=1 h(k) j m , (4.23)
showing that the time-average and the ensemble average commute. On the other hand, similarly to the C function, the numerical precision may depend on ∆t chosen. But this can be controlled ensuring enough statistics for each k value, as mentioned previously. While h(t) reaches the asymptote rapidly at t ≈ t A , the sampling-time dependent stress fluctuation µ F (∆t) reaches its plateau value µ F much more slowly, at a time scale ∆t ≈ t * . Remarkably, we can understand it by the relation between these two time-dependent properties. Starting from µ F , we have

µ F (∆t) = (1/n) i (a i -ā) 2 = a 2 -ā2 = (1/2n 2 ) i l (a i -a l ) 2 (4.24)
where ... stands for the ensemble-average over the trajectories (of length n). The last line can be written as

µ F (∆t) = (1/n 2 ) n i=l+1 n-1 l=1 (a i -a l ) 2 = (1/n 2 ) n-1 l=1 n-l k=1 (a l+k -a l ) 2 = (1/n 2 ) n-1 k=1 n-k l=1 (a l+k -a l ) 2 = (1/n 2 ) n-1 k=0 (n -k)g(k) (4.25)
where "mean-square advance" g(k) = (a l+k -a l ) 2 is independent of l (time translational invariance). Noting g(k) = 2h(k) (the MSD of the shear stress) and supposing δt << 1, we can write eq. (4.25) into an integral:

µ F (∆t) = (2/

The shear relaxation modulus

The shear relaxation modulus, G(t), is a central rheological property of materials.

It is defined by the linear response macroscopic relation

G(t) = ∂τ (t; γ) ∂γ (4.30)
where γ is a small step strain imposed at t = 0, and τ (t) is the averaged stress increment of the system. G(t) gives the measure of the viscoelasticity of a system. It is important that we can compute this function by using statistical physics tools. Previously, we showed that the static shear modulus G can be computed using the stress fluctuations formalism, by Monte Carlo (MC) or molecular dynamics (MD) simulations. This quantity is denoted by G eq here, to distinguish it from G(t). Obviously, lim t→∞ G(t) = G eq . For an elastic body, G eq > 0, and for a fluid, G eq = 0. Given that

τ (t) = τ (t)
, what is the relation between G(t) and the SACF C(t), or C(t). Many authors assumed G(t) = C(t) = C(t). We shall demonstrate that this is only true for the fluids, where we have τ = 0 (at equilibrium), and G eq = 0. In more general case, we start from

τ (t) = t -∞ dsG(t -s)dγ(s)/ds = G(t -s)γ(s)| t -∞ - t -∞ dsdG(t -s)/dsγ(s) (4.31)
4.6 Dynamic moduli 4.6.1 Definitions and simple models

In a simulation, the function G(t) can be collected from the dynamics of the system. This quantity can also be measured in rheological experiments. However, more frequently, people can use oscillatory stresses to measure the dynamic modulus G * (see refs. [START_REF] Larson | The structure and rheology of complex fluids[END_REF][START_REF] Colby | Polymer physics[END_REF]). This quantity is frequency dependent G * = G * (ω) with ω the angular frequency. More precisely, if we impose an oscillatory shear strain (γ 0 sin(ωt)), given, in complex representation, by

γ(t) = γ 0 exp(iωt), (4.36)
in the linear response regime the average shear stress will be

τ (t) = G * (ω)γ(t) ≡ (G ′ (ω) + iG ′′ (ω))γ 0 exp(iωt), (4.37)
The real part of G * , G ′ (ω), is the storage modulus, corresponding to the part of the response of the system in phase with γ(t), i.e. the elastic response; the imaginary part G * , G ′′ (ω), is the loss modulus, corresponding to the part that is dephased π/2 w.r.t. γ(t), i.e. the viscous response.

The physical meaning of G ′′ can be shown by computing the average power that the stresses should provide to maintain the periodic oscillation of the system. We write, per unit volume, this average as

Power T = (1/T ) T 0 τ (t) * γ(t)dt, (4.38)
where T = 1/2πω is the period of the sollicitation γ(t). Using γ(t) = γ 0 sin(ωt) and

τ (t) = Im(G * .γ 0 exp(iωt)) = G ′ sin(ωt) + G ′′ cos(ωt), we obtain Power T = G ′′ .γ 0 .ω/2 ∝ G ′′ .
As this is the power lost by the system, G ′′ is indeed the "loss modulus".

In order to relate G * with G(t), let's consider first a viscoelastic liquid, with G eq = 0, and compute G * l (l for "liquid"). We start from the Boltzmann superposition principle [START_REF] Colby | Polymer physics[END_REF], which gives

τ (t) = t -∞ G l (t -t ′ ) γ(t ′ )dt ′ . (4.39)
From eq. (4.36), we obtain γ(t) = iωγ 0 exp(iωt). After a change of the variable in eq. (4.39), by posing s = t -t ′ , we obtain

τ (t) = iωγ 0 exp(iωt) ∞ 0 G l (s) exp(-iωs)ds = γ(t)ω ∞ 0 G l (s) (sin(ωs) + i cos(ωs)) ds (4.40)
By identification, we obtain

G ′ l (ω) = ω ∞ 0 G l (t) sin(ωt)dt G ′′ l (ω) = ω ∞ 0 G l (t) cos(ωt)dt. (4.41)
In general situations where G eq = 0, eq. ( 4.41) is generalized to

G ′ (ω) = G eq + ω ∞ 0 (G(t) -G eq ) sin(ωt)dt G ′′ (ω) = ω ∞ 0 (G(t) -G eq ) cos(ωt)dt. (4.42)
This shows that G ′ corresponds to the Fourier-sine transform of G(t), and G ′′ to its Fourier-cosine transform. Their low frequency limits are well known:

lim ω→0 G ′ (ω) = G eq lim ω→0 G ′′ (ω)/ω = ∞ 0 (G(t) -G eq ) dt = ∞ 0 C(t)dt = η (4.43)
where we used eq. (4.33), and the Green-Kubo relation to obtain viscosity coefficient η (see [START_REF] Hansen | Theory of simple liquids[END_REF]) . For a liquid, η is the usual viscosity. For a solid, η represents the slope of

G ′′ for ω → 0.
We can also draw some conclusions on the high frequency limit. Supposing, for the stress ACF C(t), lim t→∞ C(t) = 0, and using C(0) = µ F , from eq. (4.42) we can deduce (by integration by parts),

lim ω→∞ G ′ (ω) = µ A = G(0) lim ω→∞ G ′′ (ω) = 0. (4.44)
This tells us that the system, either in liquid or solid states, possesses a high frequency non-zero elasticity, given by the affine modulus µ A . On the other hand, G ′′ vanishes at high frequencies.

Now we can compare these general trends with the ones obtained from two simple classical rheological model, the Maxwell model and the Kelvin-Voigt (K-V) model [START_REF] Colby | Polymer physics[END_REF].

They consist on the association of two elementary rheological elements, the elastic element, a spring of constant G, and a viscous element of viscosity η. The Maxwell model is a serial association and the K-V model a parallel association (see Fig. To determine G(t), we can make a "relaxation" experiment: we set

γ(t) =      0 if t < 0, γ 0 if t ≥ 0, (4.46)
and calculate σ(t) for t > 0 according to eq. ( 4. [START_REF] Speedy | Pressure and entropy of hard-sphere crystals[END_REF], given that at t = 0 + , it is the elastic element which can react instantly: σ(0

+ ) = Gγ E (0 + ) = Gγ(0 + ) = Gγ 0 . The solution is G(t) = Ge -t/τ M (t ≥ 0) (4.47)
with τ M = η/G the Maxwell relaxation time. This result corresponds to a viscoelastic fluid: G eq = lim t→∞ = 0, and G(0) = G.

The dynamic moduli can be obtained in two ways: a) by solving directly eq. (4.45), using γ(t) = γ 0 e iωt and σ(t) = G * γ(t); b) by eq. (4.42). Both give

G ′ (ω) = G (ωτ M ) 2 1 + (ωτ M ) 2 ; G ′′ (ω) = G ωτ M 1 + (ωτ M ) 2 .
(4.48)

We see that

G ′ (0) = G ′′ (0) = 0, lim ω→0 G ′′ (ω)/ω = Gτ M = η, lim ω→∞ G ′ (ω) = G, lim ω→∞ G ′′ (ω) = 0.
For ω < τ M -1 , G ′ < G ′′ (the system is more viscous at low frequencies), for ω > τ M -1 , G ′ > G ′′ (the system is more elastic at high frequencies),

at ω = τ M -1
, there is a crossover, G ′ = G ′′ , thus τ M -1 is the characteristic frequency of the system. For the K-V model, we have γ = γ E = γ V , σ = σ E + σ V , thus the rheological equation is:

Gγ + η γ = σ (4.49)
To determine G(t), we can again make the "relaxation" experiment, by solving eq. (4.49)

for t > 0, where γ = 0. The solution is

G(t) = G; (t > 0) (4.50)
This result corresponds to a solid: G eq = lim t→∞ = G. For this model, G(0) is unphysical, because it is singular (looking at eq. (4.49) ).

The dynamic moduli can be obtained by solving directly eq. (4.49), using γ(t) = γ 0 e iωt and σ(t) = G * γ(t). This gives

G ′ (ω) = G; G ′′ (ω) = Gωτ M , (4.51)
Obviously, G ′′ is unphysical at high frequencies, because it tends to infinity. By using eq. ( 4.42), we can complete the singular part of G(t), to obtain

G(t) = Gτ M δ(t) + G = ηδ(t) + G. (4.52)
From eq. (4.52), we have the stress ACF, C(t = 0) → ∞, in contradiction to C(0) = µ F (finite) for any real system. Another measurable mechanical quantity is the compliance J * (see for example ref. [START_REF] Colby | Polymer physics[END_REF]), defined by γ(t) = J * (ω)σ(t), where σ(t) = σ 0 e iωt is imposed. Clearly J * = 1/G * . It is customary to write J = J ′ -iJ ′′ . For the Maxwell model, we have

J ′ (ω) = 1 G ; J ′′ (ω) = 1 Gωτ M . (4.53)
Again, we see that for ω < τ M -1 , J ′ < J ′′ (the viscosity dominates at low frequencies),

for ω > τ M -1 , J ′ > J ′′ (the elasticity dominates at high frequencies). As for the compliance, we obtain, for the K-V model, And the K-V model does not apply to our solids either, because, clearly we have two limit values for the storage modulus (G eq and µ A ) instead of only one in the K-V model (see Fig. 4.22). Furthermore, our G ′′ are very different from G ′′ by the K-V model, which is unphysical for high frequencies. Nevertheless, we will (later) compare the compliance of our solids with the K-V model, to see whether there is some agreement for certain frequency zones.

J ′ (ω) = 1 G(1 + (ωτ M ) 2 ) ; J ′′ (ω) = ωτ M G(1 + (ωτ M ) 2 ) . ( 4
Another point about the dynamic moduli of our solids is that the loss modulus G ′′ is not at all negligible for some frequencies. Its maximum can even be as important as G eq .

This shows the potential importance of this quantity (representing the dissipation), even for a solid. It highlights certainly the interest of our method, allowing accurate computation of this quantity. Furthermore, we can see that G ′′ (ω) displays more complex features for KA2d and AB13 solids than the LJ2d and LJ3d solids. This can certainly be traced back to the former systems being binary mixtures, rather than monodisperse systems, as the LJ2d and LJ3d systems. (circles) and 2 × 10 -4 (squares). This allows to cover both the low frequency and the high frequency zones of the dynamic moduli. Chapter 5

Conclusions

Summary

In this thesis, we have studied elasticity, viscoelasticity, and glass transition of model systems by computer simulations, using equilibrium stress and volume fluctuations.

The model systems considered are as follow: monodisperse hard sphere fcc crystal, polydisperse hard sphere fcc crystal, a glass-former consisting of the two dimensional Kob-Andersen model (KA2d system), its companion system, a monodisperse LJ2d system. Then the monodisperse LJ3d system and a binary repulsive mixture AB13 system, forming respectively simple fcc and complex superlattice crystals, are studied.

From these systems we obtained some new and interesting results.

First, as a simple test of the relation linking the bulk modulus and the volume fluctuations under constant pressure, we computed the bulk modulus of (slightly) polydisperse hard sphere crystals at high pressure by Monte Carlo NP T simulations [START_REF] Li | Testing a simple method for computing directly the bulk modulus by NPT simulation: The case of polydisperse hard sphere solids[END_REF]. We showed that the equation of state proposed by Bartlett [49] is not valid for high pressures, whereas our results give the effects of the polydispersity on the average density and the bulk modulus. Furthermore, the volume distribution displays, as expected, a Gaussian histogram, indicating good statistics of our simulations. This approach is a relatively simple way to calculate the bulk modulus of general systems. It is thus potentially interesting for various purposes (testing liquid theories, studying systems near a phase transition, or under high pressure, etc).

Second, we investigated the two dimensional 80-20 Kob-Andersen model [START_REF] Li | Glass transition of two-dimensional 80-20 Kob-Andersen model at constant pressure[END_REF] under constant moderate pressure with focus on the elastic moduli calculation. Our results

show that the KA2d model is a reasonably good glass-former. We examined three ways of determining the glass transition temperature T g , and found that they give consistent results. The shear modulus as a function of the temperature G(T ) has been characterized. The result has been found to depend strongly on sampling time, with very slow convergence near the glass transition. We observed a continuous cusp-type variation of G(T ) across the glass transition. This finding is in qualitative agreement with recent numerical work using similar glass-forming colloidal systems in two and three dimensions [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF], but in contradiction with predictions from the mode-coupling theory.

Along with the KA2d model, we also investigated a simple companion model provided by monodisperse 2dLJ beads, which crystallizes at low temperature. In contrast with the KA2d model, we observed a discontinuous jump of G(T ) at the freezing transition.

Third, we investigated the shear stress fluctuations in simple and complex crystals including Lennard-Jones (LJ3d) model systems and a binary AB13 model systems.

For each system, the pressure is fixed. The evolution of the systems with temperature is examined. We were able to locate the melting transition temperature by the discontinuous jump of the elastic constants. This jump is consistent with the jump in specific volume at the transition. In the low temperature region, the fluctuation terms of elastic constants for these two model systems are scrutinized. While for the simple fcc crystal, these terms go to zero for T = 0, they are notably non-zero for the more complex AB13 crystal, contributing to about 10% for the shear modulus. This is less than in a glass, but still significant. The degree of anisotropy of the two crystals can be examined by looking at the two ways of computing the shear modulus, using G = C 66 -P , or G 2 = (C 11 -C 12 )/2 -P , by recalling that they are equal for an isotropic system. We observe that the fcc crystal is quite anisotropic, with G ≈ 3G 2 .

On the other hand, the AB13 superlattice seems to be much more isotropic, since we have in this case G ≈ G 2 , with G only slightly above G 2 at very low temperatures.

Fourth, we examined the time dependent properties of our four model systems such as the sampling time dependence of the stress fluctuations, the shear-stress autocorrelation function C(t), the shear stress mean square displacement h(t), the shear-stress relaxation modulus G(t), the relationship between these functions, and their behaviour with the thermodynamic conditions. The dynamic moduli, G ′ (ω) and G ′′ (ω), and the compliance moduli J ′ (ω) and J ′′ (ω), are also computed. From the sampling time ∆t dependence of the stress fluctuations, we can see that the sampling time more than 100 in reduced units is sufficient for the LJ2d system, the AB13 system and the LJ3d system, while it is more than 1000 for KA2d system at least, when the temperature is close to T g . Clearly, this glassy system (near T g ) evolves much more slowly than crystal systems. As for the shear stress autocorrelation function C(t; ∆t), we have confirmed, as expected, the sampling time dependence of this function. We note that in most cases, ∆t ≈ 100 gives convergent C(t). For the shear stress mean square displacement h(t), the AB13 system show oscillation at low temperature during the relaxation time.

The function is seen to be non-monotonic in the crystal, at low temperatures. This feature has not been observed in the permanent elastic bodies [START_REF] Wittmer | Simple average expression for shearstress relaxation modulus[END_REF]. The long-time plateau value of h(t) increases with the increment of temperature at solid state while the trend is the opposite at liquid state. This is not surprising since it corresponds to µ F . We note that h(t) is smoother than C(t), being a simple average. In general, it indeed is a better function for the calculation of G(t) via G(t) = µ A -h(t), as pointed out in [START_REF] Wittmer | Simple average expression for shearstress relaxation modulus[END_REF]. Here, for some state points, we have computed shear relaxation modulus G(t). The function displays a rapid decrease from G(0) = µ A (the affine elasticity), and it reaches its equilibrium modulus G eq = µ A -µ F at long times. From G(t), we have computed the dynamic moduli, and compared our results with the Maxwell and the Kelvin-Voigt rheological models. The four systems studied have qualitatively the same behaviour as the Maxwell rheological model when they are liquid (with some small differences discussed in Chapter 4). At solid state, the AB13 system and the KA2d system show more complicated features than the other two systems. We speculate that the potentials for the AB13 system and the KA2d system are more complex and they have two kinds of particles so that they have different response to the external fields. Here, for the storage modulus G ′ , clearly the Maxwell model does not apply since we have two distinct limit G ′ values, i.e. G ′ 0 = G eq for ω = 0 and

G ′ ∞ = µ A for ω → ∞, with G ′ ∞ > G ′ 0 .
Because of these two values, the Kelvin-Voigt model does not apply either. Interestingly, we have again a peak of G ′ , at some intermediate frequency, where G ′ max > µ A (resonance). As for the loss modulus G ′′ , while it is zero for low and high frequencies, as expected, its value is by no means neglectable for some intermedi-ate frequencies, with its maximum even comparable in magnitude to G eq . This shows that even in these simple model systems, the energy dissipation in the solid state can be important for some frequencies. The chapter ends with a brief discussion on the compliance moduli J ′ and J ′′ of our models in solid state.

Outlook

Many future extensions of our work are possible:

Concerning the elasticity, we can use the elastic constants to probe solid-solid phase transitions under high pressure, or with temperature. More general potentials, such as three-body, or EAM potentials should be examined. The viscoelastic functions, G(t), G ′ and G ′′ can be calculated for many soft matter systems and, by comparison with experimental results, allow a better modelling of these systems. For the glass transition, our study of the KA2d system must be carried further. For example, the trajectory analysis can be implemented. More independent configurations should be used, in order to reduce statistical errors in our moduli, especially near the glass transition. In parallel, it would be interesting to obtain similar expressions as G(t) = µ A -h(t) for hard sphere systems, since these systems are objects of many theoretical and experimental investigations (hard sphere colloids).

Résumé de la thèse en français Cette particularité n'a pas été observée dans les corps élastiques permanents [START_REF] Wittmer | Simple average expression for shearstress relaxation modulus[END_REF]. La valeur du plateau à long terme de h(t) augmente avec l'augmentation de la température à l'état solide tandis que la tendance est à l'opposé à l'état liquide. Cela n'est pas surprenant puisqu'il correspond à la fluctuation des contraintes de cisaillement µ F . Ensuite, pour certains états thermodynamiques, nous avons calculé le module de relaxation de cisaillement G(t). Comme indiqué dans [START_REF] Wittmer | Simple average expression for shearstress relaxation modulus[END_REF], cette fonction peut se calculer via G(t) = µ A -h(t), où µ A est le module de l'élasticité affine. La fonction affiche une

  .94) Thus B αβχκ = (λ + P )δ αβ δ χκ + (µ -P )(δ αχ δ βκ + δ ακ δ βχ ) (2.95) yielding B 11 = λ+2µ-P , B 12 = λ+P , and B 66 = µ-P = G. These are d-independent relations.

  ∂V ) T > 0 and G = 0. It seems difficult to transpose previous results to liquids, since they are based on explicit displacement fields. For the compression modulus, Rowlinson [9] demonstrated eq.(2.112) for K via a discussion of the fluctuations of P ex . The system was put in a cubic box of length L . One writes for the particle i r iα = s iα L (2.129) where α = x, y, z (d = 3), (x, y for d = 2), the reduced coordinates 0 ≤ s iα ≤ 1. As the pressure P = -(∂F/∂V ) T = P id + P ex , with V = L d , F the free energy, one can

Fig. 2 .

 2 Fig. 2.1 A two-dimensional simulation box (shaded) containing 5 particles with its nearest periodic images.

Fig. 3 . 1

 31 Fig.3.1 The reduced average density ρ * vs the reduced pressure P * for monodisperse and slightly polydisperse solids. In the monodisperse system, ρ * = ρσ 3 and P * = βP σ 3 .In the polydisperse case, ρ * = ρσ 3 and P * = βP σ3 . The continuous line is the analytical EOS of the monodisperse system given by eq.(3.1). The dots represent NP T simulation results of the monodisperse system. The dashed line is the analytical EOS of the polydisperse system given by eq.(3.3), for P * ≤ 20. The dotted line is its extrapolation for P * ≥ 20. The triangles represent NP T simulation results of the polydisperse system.

Fig. 3 . 2

 32 Fig.3.2 The reduced bulk modulus K * vs the reduced pressure P * for monodisperse and slightly polydisperse solids. In the monodisperse system, K * = βKσ 3 and P * = βP σ 3 .In the polydisperse case, K * = βK σ3 and P * = βP σ3 . The continuous line is given by eq.(3.2). The circles (with error-bars) represent NP T simulation results of the monodisperse system. The dashed line is the analytical expression for the polydisperse system given by eq.(3.3), in the range 12 ≤ P * ≤ 20. The dotted line is its extrapolation for P * ≥ 20. The triangles represent NP T simulation results of the polydisperse system (error-bars, similar to the monodisperse case, are not shown).

Fig. 3 .

 3 Fig. 3.3 Volume-histogram at the reduced pressure P * = 35.5 for monodisperse and slightly polydisperse solids. The dots represent NP T simulation results of the monodisperse system. The triangles represent those of the polydisperse system. The continuous and dashed lines represent fits by Gaussian distribution for respectively the monodisperse and polydisperse cases. V * = V /σ 3 for the monodipserse case, and V * = V /σ 3

3. 2 A

 2 glass former in two dimensions: 80-20 Kob-Andersen model 3.2.1 Background

  volume. As explained above, we impose a normal pressure P = 2 for all Table 3.2.1 Several properties of the KA2d model at pressure P = 2 as a function of temperature T: number densityρ = N/V = 1/v, linear box size √ 3L x = L y ,Abraham parameter R = g min /g max[START_REF] Abraham | An isothermal-isobaric computer simulation of the supercooledliquid/glass transition region: Is the short-range order in the amorphous solid fcc?[END_REF], ideal pressure contribution P id = T ρ to the imposed normal pressure P = P id + P ex , hypervirial contribution η A,ex to the compression modulus, compression modulus K obtained using the Rowlinson formula eq.(3.6)[START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Wittmer | Compressibility and pressure correlations in isotropic solids and fluids[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature, and sampling time[END_REF], excess contribution µ A,ex to the affine shear elasticity µ A = µ A,id + µ A,ex and shear modulus G obtained using the stress-fluctuation formula eq.(3.7). The last five columns have the dimension energy per volume. The error bars for K and

Fig. 3 . 4

 34 Fig. 3.4 Rescaled specific volume v(T ) vs. temperature T for both models considered. The freezing temperature T f ≈ 0.62 for the mdLJ model and the glass transition temperature T g ≈ 0.35 for the KA2d model are indicated by vertical lines.

3 . 5

 35 the radial pair-distribution function (RDF) of the mdLJ crystal at T = 0.3 (spheres) with the (total) RDF g(r) of the KA2d system for T = 0.1(solid line). Whereas the former system becomes clearly crystal-like with long-range spatial correlations (at least if boundary conditions with L y = √ 3L x are used), we find for the KA2d model that g(r) ≈ 1 for r > 4 as is typical for an amorphous substance. The amorphous, non-crystalline behaviour is also confirmed using the corresponding structure factor S(q) in reciprocal space (not shown) and more readily by inspection of snapshots such as the one for T = 0.2 shown in the inset of Fig.3.5. We use R-parameter criterion. An additional empirical criterion for the liquid-to-solid/glass transition has been proposed by Abraham[START_REF] Abraham | An isothermal-isobaric computer simulation of the supercooledliquid/glass transition region: Is the short-range order in the amorphous solid fcc?[END_REF]. It uses the temperature dependence of the so-called R-parameter defined by R ≡ g min /g max with g max being the main peak of the radial pair-distribution function (RDF) and g min its subsequent minimum. As shown in Fig.3.6, R(T) shows for the mdLJ model a jump discontinuity at T f ≈ 0.62 consistent with Fig.3.4. Although the criterion was suggested originally only for a one-component system, we extend it here to our binary mixture in two ways by using either the total RDF g(r) sampled over all beads or the RDF g 11 (r) considering only the large spheres. The R-parameter obtained from the latter distribution is shown by the diamonds. The glass transition temperature may be obtained by matching the low-and the high-temperature linear slopes indicated in the figure. We get T g ≈ 0.35 as above from the specific volume v(T ).Compression modulus. We turn now to the two elastic moduli, the compression modulus K(T ) and the shear modulus G(T ), which characterize completely the linear elastic response of both models even in the triangular phase of the mdLJ model. The temperature dependence of the compression modulus K (circles) and its contributions η A,ex (squares) and η F,ex (triangles) is presented in Fig.3.7. The affine hypervirial contribution η A,ex (T ) is seen in both cases to increase more or less linearly with decreasing temperature. As seen in panel (b) for the KA2d model, one may determine again a transition temperature from the intercept of the low-and the high-T tangent slopes.

Fig. 3 . 5 Fig. 3 . 6 R

 3536 Fig. 3.5 Main panel: Radial pair-distribution function g(r) vs. distance r for two lowtemperature states of both models. While long-range correlations are visible for the mdLJ model (small circles), no long-range order is seen for the KA2d model where g(r) → 1. Inset: Snapshot for T = 0.2 ≪ T g confirming the amorphous structure of the KA2d model. Open/filled circles correspond to the large/small beads.

Fig. 3 . 7

 37 Fig. 3.7 Compression modulus K(T) and its contributions η A,ex (T ) and η F,ex (T ) according to the Rowlinson stress-fluctuation formal eq.(3.6) : (a) Data for the mdLJ model. K(T ) vanishes at the phase transition at T f ≈ 0.62. The stress fluctuation contribution η F,ex (T ) vanishes for T → 0. (b) Data for the KA2d model showing a similar, albeit slightly smaller transition temperature T g ≈ 0.32 as above. η F,ex (T ) does not vanish for T → 0.

Fig. 3 . 8

 38 Fig. 3.8 Affine shear elasticity µ A,ex (squares), shear stress fluctuation µ F,ex (triangles) and shear modulus G = µ A,ex -µ F,ex (circles) vs. temperature T: (a) For mdLJ systems G(T ) reveals a striking jump discontinuity at T f ≈ 0.62. (b) The continuous increase of G(T) for KA2d systems below T g ≈ 0.32 shows a cusp-singularity with α ≈ 0.6.µ F,ex (triangles) does not vanish for T → 0. Inset: Zoom on the glass transition region plotting G(T ) = G g vs. x = T /T g for the KA2d model (T g = 0.32,G g = 16.5) and the KA3d model (T g = 0.41, G g = 24.5). The lines indicate equation eq.(3.5) with α = 0.6 for the KA2d model and α = 0.5 for the KA3d model. 55

3 . 1 2 Fig. 3 . 9 0 Fig. 3 .

 3123903 Fig.3.9 The specific volume as a function of the temperature for LJ model system.The pressure is fixed at P = 1.

83 Fig. 3 .

 833 Fig. 3.11 The elastic constants and the fluctuation terms of the elastic constants for LJ model system as a function of the temperature (in the solid region). The upper group of three (curves with hollow shapes) represent the total elastic constants, in decreasing order C 11 , C 44 and C 12 . The lower group of three (curves with solid shapes) represent the fluctuation parts, in decreasing order C 11f , C 12f and C 44f . The pressure is P = 1.

0 Fig. 3 .

 03 Fig.3.12 The shear modulus G as a function of the temperature for LJ model system.

0 Fig. 3 .

 03 Fig.3.[START_REF] Colby | Polymer physics[END_REF] The bulk modulus K as a function of the temperature for LJ model system.The pressure is P = 1.

  3.14 and 3.15, the superlattice structure at T = 0 is shown.

Fig. 3 .

 3 Fig. 3.14 Basic units of the primitive cell of the AB13 superlattice.

Fig. 3 .

 3 Fig. 3.15 The primitive cell of the AB13 superlattice. (a) Front view. (b) Side view.

Fig. 3 .

 3 Fig. 3.16 Potentials for AB13 model system. Species 1 stands for A particles, 2 for B particles. Distances are measured by units of σ A , potentials by ǫ AA .

Fig. 3 .

 3 Fig. 3.18 Radial distribution functions of AB13 system at temperature T = 0.1.

Fig. 3 .

 3 Fig. 3.19 Radial distribution functions of AB13 system at two different temperatures (2.0, 2.1).

Fig. 3 . 0 Fig. 3 .

 303 Fig. 3.20 Elastic constants vs. temperature for AB13 crystal.

5 Fig. 3 .

 53 Fig. 3.23 Bulk modulus vs. temperature for AB13 model

  4.1).

8 Fig. 4 . 1

 841 Fig. 4.1 The stress fluctuations µ F (∆t) as a function of the sampling time ∆t (LJ2d system).

Fig. 4 .Fig. 4 . 2

 442 Fig. 4.2 The stress fluctuations µ F (∆t) as a function of the sampling time ∆t (LJ2d and KA2d systems).

2 Fig. 4 . 3 0 Fig. 4 . 4

 243044 Fig. 4.3 The stress fluctuations µ F (∆t) as a function of the sampling time ∆t (LJ3d system).

8 |Fig. 4 . 5

 845 Fig. 4.5 Sampling-time-dependence of the shear modulus G F (∆t) = µ A,ex -μF,ex (∆t) as a function of the sampling time ∆t for LJ2d and KA2d models and different temperatures as indicated. In all cases we reach a horizontal plateau G(T ) for sufficiently large ∆t as indicated by the dashed horizontal lines. (Here μF,ex (∆t) is defined by eq. (4.6), using only the excess shear stress, and µ A,ex = µ A -ρk B T ).

1 C 1 C 8 Fig. 4 . 6

 11846 Fig. 4.6 The shear-stress autocorrelation function (SACF) of the LJ2d system as a function of the time, for several sampling times ∆t. The left panel corresponds to a low temperature crystal, T = 0.2, the right one to a liquid, T = 0.8.

1 C 8 Fig. 4 . 7

 1847 Fig. 4.7 The SACF of the KA2d system as a function of the time, for several sampling times ∆t. The left panel corresponds to a low temperature glass, T = 0.1, the right one to a liquid, T = 0.8.

1 C 1 C 8 Fig. 4 . 8

 11848 Fig. 4.8 The SACF of the LJ3d system as a function of the time, for several sampling times ∆t. The left panel corresponds to a low temperature crystal, T = 0.2, the right one to a liquid, T = 0.8.

1 C 2 Fig. 4 . 9

 1249 Fig. 4.9 The SACF of the AB13 system as a function of the time, for several sampling times ∆t. The left panel corresponds to the AB13 superlattice, T = 2.0, the right one to a liquid, T = 2.2.

Table 4 . 1 Fig. 4 .

 414 Fig. 4.10 The SACF C(t) and the function y(∆t) (inset) for the monodisperse LJ2d system. The red (dashed) curve corresponds to a crystal (T = 0.6), while the black (continuous) curve represents a liquid state (T = 0.7). The blue dash-dotted line in the inset indicates the level y = 0.05.

Fig. 4 .

 4 Fig.4.11 The SACF C(t) and the function y(∆t) (inset) for the glass-forming KA2d system. The red (dashed) curve corresponds to a glass (T = 0.3), near T g = 0.32, while the black (continuous) curve represents a liquid state (T = 0.4). For the liquid, C(0) ≈ 34 (the upper black dot), while for the glass, C(0) ≈ 33 (the lower black dot).The blue dash-dotted line in the inset indicates the level y = 0.05.

Figs. 4 .

 4 Figs.[START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF].12 -4.[START_REF] Heyes | Scaling of Lennard-Jones liquid elastic moduli, viscoelasticity and other properties along fluid-solid coexistence[END_REF] show examples of h(t) function. We see that the function starts at 0, and saturates at long times to its plateau value. For all curves but the third one, the relaxation time t A is about 1. The evolution of h(t) is much slower for the KA2d liquid and glass close to the glass transition temperature, and is oscillatory for the low temperature AB13 crystal. The first observation can be understood by longer stress

Fig. 4 .

 4 Fig. 4.12 The shear stress mean square displacement as a function of the time for the LJ2d model systems. The first time point is t 1 = 0.005).

Fig. 4 .

 4 Fig. 4.13 The shear stress mean square displacement as a function of the time for the KA2d model systems. Here, the first time point is t 1 = 0.05 (except for T = 0.28, where t 1 = 0.1).

Fig. 4 .

 4 Fig. 4.14 The shear stress mean square displacement as a function of the time for the LJ3d model systems.

Fig. 4 .

 4 Fig. 4.15 The shear stress mean square displacement as a function of the time for the AB13 model systems.

4

 4

  .20).

Fig. 4 .

 4 Fig. 4.20 Rheological models. σ is the shear stress, γ is the shear strain.

Fig. 4 .

 4 Fig. 4.21 shows G(t) and G * of the Maxwell model.

Fig. 4 .

 4 Fig. 4.21 Rheological moduli of the Maxwell model, τ M = η/G.

Fig. 4 .

 4 Fig. 4.22 shows G(t) and G * of the K-V model.

  G(t) of K-V model.

Fig. 4 .

 4 Fig. 4.22 Rheological moduli of the Kelvin-Voigt model, τ M = η/G.

Fig. 4 .

 4 Fig. 4.23a shows J * of the Maxwell model.

. 54 )Fig. 4 .

 544 Fig. 4.23b shows J * of the K-V model. One can see that the J * moduli display better (than G * ) the viscoelastic character of the solid represented by this model.

0 G′′ ≈ 7 . 5 G′′ ≈ 0 G′ ≈ 13 . 6 Fig. 4 .

 07501364 Fig. 4.24 Dynamic moduli of LJ2d model in the triangular crystal structure T = 0.6.Two recording time intervals for the stress tensor are used to generate G(t), i.e. 5×10-3

Fig. 4 .

 4 Fig. 4.25 Dynamic moduli of LJ2d model system in liquid state T = 0.8. The fit uses τ M = 0.10, G = µ A and eq. (4.48).

0 G′ ≈ 14 . 6 Fig. 4 .

 01464 Fig. 4.26 Dynamic moduli of KA2d model system as a low temperature glass T = 0.1.

8 G′ ≈ 29 G′′ ≈ 0 G′′ ≈ 14 G′, G′′ ≈ 0 G′ ≈ 26

 829014026 model G′′ of Maxwell model G′ and G′′ for KA2d model systems N = 1024, P = 2.0, T = 0.

Fig. 4 .

 4 Fig. 4.27 Dynamic moduli of KA2d model system in liquid state T = 0.8. The fit uses τ M = 0.091, G = µ A and eq. (4.48).

8 G′ ≈ 41 . 6 G′ ≈ 37 . 6 G′′ ≈ 0 G′′ ≈ 10 . 0 G′ ≈ 19 . 3 G′′ ≈ 0 Fig. 4 .

 8416376010019304 Fig. 4.28 Dynamic moduli of LJ3d model system in crystal state T = 0.8.

9 G′ ≈ 23 . 4 G′ ≈ 20 . 6 G′′ ≈ 0 G′′ ≈ 12 . 3 G′, G′′ ≈ 0

 923420601230 model G′′ of Maxwell model G′ and G′′ for LJ3d model systems N = 1372, P = 1.0, T = 0.

Fig. 4 .

 4 Fig. 4.29 Dynamic moduli of LJ3d model system in liquid state T = 0.9. The fit uses τ M = 0.088, G = µ A and eq. (4.48).

0 Fig. 4 .≈ 0 Fig. 4 .

 0404 Fig. 4.30 Dynamic moduli of AB13 model system in crystal state T = 2.0.

0. 1

 1 Cadre généralL'élasticité et la viscoélasticité sont des propriétés essentielles de beaucoup de matériaux, en état solide ou fluides complexes. Il est important de les connaître pour de nombreuses applications pratiques. Ces grandeurs sont sensibles aux conditions thermodynamiques, notamment la température et la pression. En outre, leur brusque changement est souvent la signature d'une transition de phase. Ces aspects montrent qu'il est souhaitable de pouvoir étudier ces propriétés par simulation sur ordinateur, afin de modéliser les matériaux, et de prédire leur comportement dans des conditions potentiellement difficiles d'accès par expérience. Dans cette thèse, l'élasticité, viscoélasticité et la transition vitreuse de liquides vitrifiants et de solides modèles sont étudiés par des simulations moléculaires à l'équilibre en utilisant le formalisme de fluctuation de contraintes. Ce formalisme permet des études de l'élasticité et de la viscoélasticité sans déformer le système, mais en utilisant les fluctuations des contraintes à l'équilibre. Bien que proposé à la fin des années 1960[START_REF] Squire | Isothermal elastic constants for argon. Theory and Monte Carlo calculations[END_REF], il n'a réellement attiré l'attention des chercheurs du domaine qu'à la fin des années 1980[START_REF] Barrat | Elastic response of a simple amorphous binary alloy near the glass transition[END_REF][START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF]. Et certaines avancées théoriques sont survenues assez récemment (dans les années 2000)[START_REF] Wittmer | Vibrations of amorphous, nanometric structures: When does continuum theory apply? EPL[END_REF][START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Wittmer | Shear-stress fluctuations in self-assembled transient elastic networks[END_REF]. Notre objectif est d'appliquer ce formalisme à des systèmes modèles représentatifs de systèmes colloïdaux ou atomiques pour répondre à des questions d'ordre fondamental, surtout liées à des transitions de phase. Les systèmes modèles étudiés sont composés de systèmes de Lennard-Jones à deux et trois dimensions, un mélange binaire de Lennard-Jones (LJ) des particules constituant la version 2d du modèle bien connu de Kob-Anderson (modèle KA2d), un mélange binaire AB13 de particules répulsives et un mélange ternaire de sphères dures. Nous étudions principalement l'effet de la température T sur les propriétés mécaniques 123 des systèmes, en particulier les transitions de phase telles que la cristallisation ou la transition vitreuse. Les résultats sont présentés en deux parties, à savoir les propriétés statiques et les propriétés dépendantes du temps. La plupart des résultats sont obtenus avec le code LAMMPS, de simulation en dynamique moléculaire à léquilibre. Les propriétés statiques que nous avons étudiées sont, pour les systèmes isotropes (verres et liquides), le module de cisaillement G et le module de compression K, et pour les cristaux, les constantes élastiques. Ce sont des grandeurs sensibles à une transition de phase (solidification). Ainsi, nous examinons leur évolution en fonction de la température, et leur comportement à une transition de phase. Ainsi, nous avons montré une nette différence de comportement du module de cisaillement lors d'une transition vitreuse et lors d'une cristallisation. Alors que dans le 2ème cas, il y a une nette discontinuité de G, dans le 1er, l'évolution semble être continue. Un autre aspect intéressant est la contribution de la partie non-affine à l'élasticité, qui est calculable avec notre formalisme de fluctuation des contraintes (via le terme correspondant à ces fluctuations). Ce que nous voulons examiner est d'une part l'évolution de cette contribution en fonction de la température, et d'autre part son importance à température nulle. Nos résultats montrent que pour le cristal simple (du type cfc), ce terme est nul à T = 0, comme attendu, par contre, il est non nul pour un super-réseau du type AB13, contribuant jusqu'à 10% de la valeur du module élastique. Pour un verre, naturellement, ce terme est encore plus important (environ la moitié du module). Dans la partie consacrée aux propriétés dépendantes du temps, nous avons étudié plusieurs aspects et fonctions. Tout d'abord, nous avons exploré la dépendance du temps d'échantillonnage de nos résultats statiques, et montré que cette influence est forte pour les fluctuations de contrainte. Deuxièmement, la dynamique d'équilibre est étudiée par le biais de l'autocorrélation des contraintes de cisaillement et déplacement quadratique moyen de la contrainte, ainsi que la relation entre ces deux fonctions temporelles. L'objet principal, cependant, est le module de relaxation de la contrainte de cisaillement G(t) et les modules dynamique associée G ′ (ω) et G ′′ (ω), avec ω la fréquence. Ces fonctions caractérisent la viscoélasticité de nos systèmes. En utilisant une formule proposée récemment (Wittmer et al, Phys. Rev. E 2016), nous avons pu déterminer ces fonctions avec précision et étudier leur comportement pour nos systèmes 124 en états liquide, cristallin et vitreux. 0.2 Thèmes étudiés et résultats Tout d'abord, comme un simple test de la relation reliant le module de compression aux fluctuations de volume sous pression constante, nous avons calculé le module de compression des cristaux de sphères dures polydisperses (légèrement) à haute pression par des simulations Monte Carlo dans l'ensemble NP T [82] (voir Figure 1). Nous avons montré que l'équation d'état proposée par Bartlett [49] n'est pas valable pour des pressions élevées. En plus, nos résultats donnent les effets de la polydispersité sur la densité moyenne et le module de compression. De plus, la distribution du volume affiche, comme prévu, un histogramme gaussien, indiquant de bonnes statistiques de nos simulations. Cette approche est un moyen relativement simple de calculer le module de compression des systèmes généraux. Il est donc potentiellement intéressant à diverses fins (test des théories des liquides, étude des systèmes près d'une transition de phase, ou sous haute pression, etc.).

  Densité versus la pression.

  Module de compression.

Fig. 1

 1 Fig. 1 Densité et module de compression en fonction de la pression pour un cristal de sphères dures monodisperse (ronds) et polydisperse (triangles) respectivement, obtenus par des simulations Monte-Carlo à pression constante (ensemble NPT). Ces résultats sont comparés aux équations d'état de Speedy[45] (continu) pour le cas monodisperse et de Bartlett[49] (pointillés) pour le cas polydisperse.

6 α = 0 . 5 (

 605 b) Module de cisaillement vs T .

Fig. 2

 2 Fig. 2 Structure et le module de cisaillement G pour le modèle KA2d.

  , P = 114, T=0.05~2.0 (a) Constantes élastiques du solide AB13.

  , P = 114, T = 0.05~2.0 (b) Termes de fluctuation du solide AB13.

Fig. 3

 3 Fig. 3 Constantes élastique et termes de fluctuation du solide AB13 (super-réseau) en fonction de la température.

  , P = 114, T = 0.05~3.5 (b) Module G du solide AB13.

Fig. 4

 4 Fig. 4 Module de cisaillement G en fonction de la température des solides LJ3d (cfc) et AB13 (super-réseau) respectivement. Deux modes de calcul sont représentés (voir texte).

  diminution rapide de G(0) = µ A , et il atteint son module d'équilibre G eq = µ A -µ F à au temps long. A partir de G(t), nous avons calculé les modules dynamiques, et avons comparé nos résultats aux modèles rhéologiques simples de Maxwell et de Kelvin-128 Voigt. Les quatre systèmes étudiés ont qualitativement le même comportement que le modèle de Maxwell quand ils sont liquides (avec quelques petites différences discutées au chapitre 4). A l'état solide, le système AB13 et le système KA2d présentent des caractéristiques plus compliquées que les deux autres systèmes, qui sont monodisperses. Nous supposons cela peut s'expliquer par le caractère binaire de KA2d et AB13. Pour nos solides, concernant le module de conservation G ′ , il est clair que le modèle Maxwell ne s'applique pas puisque nous avons deux valeurs limites de G ′ non nulles (à basse et à haute fréquences). En raison de ces deux valeurs, le modèle de Kelvin-Voigt ne s'applique pas non plus. Fait intéressant, nous avons de nouveau un pic de G ′ , à une fréquence intermédiaire, où G ′ max > µ A (résonance). Quant au module de perte G ′′ , alors qu'il est nul pour les fréquences basses et hautes, comme prévu, sa valeur n'est nullement négligeable pour certaines fréquences intermédiaires, son maximum est par-fois même comparable en ordre de grandeur à G eq . Cela montre que même dans ces systèmes de modèles simples, la dissipation d'énergie à l'état solide peut être importante pour certaines fréquences. Le chapitre se termine par une brève discussion sur les modules de complaisance de nos systèmes à l'état solide (voir Figure5).

Fig. 5

 5 Fig. 5 Les modules dynamiques du modèle KA2d en état liquide et solide respectivement. G ′ (ω) est le module de conservation, G ′′ (ω) est le module de perte.

  They are B 11 = B 22 and B 12 = B 21 . For the shear modulus B 66 , we have

	Two-dimensional triangular lattice
	For a triangular lattice (2d), B is a 3×3 tensor. Only two independent B constants
	exist. B 66 = (B 11 -B 12 )/2. The other elements are zero. We see that the situation (concerning the elasticity) is similar to the case of a two-dimensional isotropic solid (glass).
	.100)

Table 3 .

 3 1.1 Coefficients entering in Eqs.(3.1) and (3.2). {a n } are taken from the work of Young et al. [44], {b n } are computed using K

  01 1.068 23.52 0.03 0.01 88.2 83.5 42.1 15.5

	0.10 1.057 23.66 0.04 0.11 85.3 74.6 40.7 11.2
	0.20 1.043 23.81 0.05 0.21 82.1 64.9 39.2 8.5
	0.25 1.035 23.90 0.06 0.26 78.9 54.9 37.7 5.7
	0.30 1.024 24.02 0.07 0.31 78.3 45.7 37.5 2.2
	0.32 1.020 24.07 0.07 0.33 75.1 40.2 35.9 0.2
	0.34 1.015 24.14 0.08 0.35 73.7 37.1 35.2	0
	0.36 1.011 24.18 0.08 0.36 74.4 37.7 35.5	0
	1.00 0.819 26.87 0.29 0.82 42.7 9.2 20.1	0

systems considered. If the systems are cooled down, the number density ρ(T ) must thus increase (Table

Table 3 .

 3 3.3. To summarize, we characterized the elastic properties of the LJ model

Table 3 .

 3 3.3 Shear modulus G, bulk modulus K and specific volume v for LJ3d model system, as a function of the temperature. The "second" shear modulus, G ′ , is defined by G ′ = (C 11 -C 12 )/2 -P . Results for T = 0 are extrapolations of T > 0 data. The pressure is P = 1 throughout, and the melting temperature is T m = 0.84.

	T	G	G ′	K	v
	0.00	60.45	22.73	77.61	0.92
	0.01	59.94	22.54	76.85	0.92
	0.05	57.92	21.73	73.71	0.92
	0.10	55.01	20.49	69.46	0.93
	0.20	49.19	17.97	60.55	0.94
	0.30	44.01	15.88	53.56	0.96
	0.40	39.13	13.99	46.14	0.98
	0.50	34.55	12.17	41.39	0.99
	0.60	29.83	10.24	32.98	1.02
	0.70	24.78	8.04	26.15	1.04
	0.80	19.31	5.55	20.22	1.08
	0.81	18.63	5.26	17.22	1.08
	0.82	17.89	4.97	16.39	1.09
	0.83	17.23	4.61	15.80	1.09
	0.84	-0.15	-0.15	13.71	1.22
	0.85	0.07	-0.01	10.43	1.24
	0.86	0.03	-0.05	11.54	1.25
	0.88	0.01	0.03	10.31	1.26
	0.90	0.04	-0.01	10.45	1.27
	1.00	-0.01	0.07	8.59	1.33
	1.10	0.01	0.11	6.48	1.40
	1.20	0.04	-0.03	5.29	1.47

Table 3 . 3

 33

	Process	Temperature
	Cooling 1.0 0.8 0.6 0.4 0.2 0.1 0.05
	Heating 1.0 1.2 1.4 1.6 1.8 2.0 2.02
	2.02 2.04 2.06 2.08 2.1 2.2 2.3
	2.3 2.

.4 AB13 system. Temperature points studied by cooling or heating a reference equilibrium system at T = 1.0 (crystal). The pressure is P = 114 throughout.

Table 3 .

 3 3.5 Elastic constants and fluctuation terms of elastic constants at different temperatures for AB13 model system in solid state (AB13 superlattice). The pressure is P = 114. The melting temperature is T m = 2.1.

		Fluctuation terms	Elastic constants
	T	C11 f	C12 f	C44 f	C11	C12	C44
	0.0	37.41	-12.56	49.88	1603.74 600.84	538.41
	0.05	53.15	-7.54	54.64	1579.67 592.86	531.20
	0.1	68.89	2.52	59.40	1555.60 584.88	523.99
	0.2	97.37	9.86	66.43	1513.31 565.63	511.11
	0.4	158.40	41.34	80.88	1422.28 518.30	482.82
	0.6	223.59	78.55	95.14	1325.24 464.71	454.16
	0.8	280.68	106.39	110.23	1236.04 421.10	425.21
	1.0	333.95	130.81	125.26	1155.05 382.71	398.09
	1.2	410.07	178.18	138.03	1045.65 319.79	371.59
	1.4	455.57	198.38	150.77	969.55	284.93	345.94
	1.6	499.56	213.93	163.65	899.79	257.09	322.47
	1.8	534.82	219.87	176.45	835.99	237.70	297.87
	2.0	572.70	227.52	188.31	767.15	216.13	273.62

Table 3 .

 3 3.6 Shear modulus G, bulk modulus K, system volume V and box length L for AB13 model system for 0 ≤ T ≤ 3.5. We recall that N = 896, P = 114, T m = 2.1, and G ′ = (C 11 -C 12 )/2 -P . Results at T = 0 are extrapolations of the data for T > 0.

	T	G	G ′	K	V	L
	0.00	424.45	387.50	973.12	173.06	5.58
	0.05	417.24	379.45	959.78	173.45	5.58
	0.10	410.03	371.40	946.44	173.84	5.58
	0.20	397.15	359.88	919.51	174.65	5.59
	0.40	368.97	338.14	857.58	176.38	5.61
	0.60	340.42	316.52	789.47	178.24	5.63
	0.80	311.47	293.73	730.66	180.25	5.65
	1.00	283.99	272.07	678.19	182.32	5.67
	1.20	257.70	249.04	599.71	184.68	5.69
	1.40	232.18	228.55	551.06	187.18	5.72
	1.60	208.48	207.36	509.32	189.75	5.75
	1.80	184.06	185.33	475.07	192.69	5.78
	2.00	160.75	164.85	430.36	196.05	5.81
	2.02	158.47	161.69	424.43	196.15	5.81
	2.04	155.63	159.48	439.01	196.49	5.81
	2.06	152.25	157.07	407.50	196.79	5.82
	2.08	149.88	154.17	377.94	197.23	5.82
	2.10	1.30	-0.18	455.14	206.91	5.91
	2.20	-0.92	-0.67	436.22	212.41	5.97
	2.30	-1.50	-0.50	427.94	213.76	5.98
	2.40	-0.83	-1.07	427.98	215.16	5.99
	2.50	-2.82	-1.20	428.09	216.44	6.00
	3.00	-1.82	-0.94	418.19	222.99	6.06
	3.50	-0.18	-0.36	388.49	229.14	6.12

  ∆t 2 ) , knowing t A of the function h(t). For this, we define y(∆t) = 1 -µ F (∆t)/µ F , where µ F is the plateau value of µ F (∆t). ≈ 100t A . So we have now an order of magnitude of t * . As C(t) = C(0) -h(t).For a given ∆t this implies C(t, ∆t) = µ F (∆t) -h(t), i.e.,

	i.e. t C(t, ∆t) = (2/∆t 2 )	0	∆t	dt(∆t -t)h(t) -h(t).	(4.29)
	Thus we have a converged C(t) only if ∆t >> t * . This fact is important, especially
	when we want to compute the shear stress relaxation modulus G(t), as we shall see in
	the next section.
				∆t
				0	dt(∆t -t)h(t)	(4.26)
	eq. (4.26) can give us a fair idea of t According to
	eq. (4.26), we have
				∆t
	y(∆t) = 1 -(2/∆t 2 )	0	dt(∆t -t)h(t)/µ F .	(4.27)
	Assuming t * is such that y(t * ) = 1%, and t * >> t A , from eq. (4.27), we obtain,
	0.01 ≈ t A /t *	(4.28)

* *

Acknowledgements

Thank China Scholarship Council for its very important support for my study.

Computer resources from the SIMPA centre of the University of Lorraine (Metz) are also gratefully acknowledged. Thank Wuhan University of Technology for its recommendation and I am grateful to Prof. Liu Hanxing and Prof.

The last elastic quantity shown is the bulk modulus K. From the graph Fig. 3. 23, we can see that the bulk modulus for the AB13 system decreases when the temperature increases. There exists a singularity when temperature is about 2.1 (transition temperature). Again, the slope in the solid is much larger than the one in the liquid.

At T = 0, this modulus is very important, K ≈ 1000, about ten times the pressure.

Near the transition, it is about 400, roughly four times the pressure. We also notice very little variation of K just below and just above T m . This should be related to the more isotropic feature of the solid (than the fcc crystal). Numerical results for these quantities are gathered in Table 3.3.6.

To summarize, our comparison of simple and complex crystals shows a difference in behaviour of the fluctuation terms at low temperatures. We also detected a more isotropic character of the AB13 superlattice. To our knowledge, the present results on Chapter 4 Time-dependent properties

Introduction

In this chapter we investigate time-dependent properties. Some of them are numerical issues, such as the sampling time dependence of the physical quantities we compute, especially the stress fluctuations. Others are the intrinsic dynamical properties of the systems investigated, for example the shear-stress autocorrelation function C(t), the shear-stress mean square displacement h(t), the shear-stress relaxation modulus G(t), the relationship between them, and their behaviour with the thermodynamic conditions. Another interesting property is the dynamic moduli, which can be obtained by sine and cosine transforms of G(t). These are important mechanical properties of materials that are accessible to the experiments. We shall show that they can also be computed quite accurately and efficiently by our equilibrium MD simulations. Most of the theoretical formalism is based on ref. [START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF][START_REF] Wittmer | Simple average expression for shearstress relaxation modulus[END_REF][START_REF] Wittmer | Shear-stress fluctuations in self-assembled transient elastic networks[END_REF][START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF]. Our investigations are focused on four systems: monodisperse two-dimensional Lennard-Jones (LJ2d) system, KA2d system, monodisperse three-dimensional Lennard-Jones (LJ3d) system and the AB13 system.

Sampling-time dependence of the stress fluctuations

In the stress fluctuation formalism, the static shear modulus G is calculated according to

where we have used integration by parts. Here we suppose γ(t) a step function (γ(t) = γ for t ≥ 0 and γ(t) = 0 for t < 0),

where, in the first line, we define χ(t) ≡ -G ′ (t) = G ′ (-t) is the response function (or the "after-effect function", and in the second line, G eq appears naturally as an integration constant. As the response function is related to

(see [START_REF] Hansen | Theory of simple liquids[END_REF]), we obtain

or

Furthermore, from eq. ( 4.33), we see that G(t) can be related to the stress MSD in a very simple way:

Eq. (4.35) shows that G(t) can be evaluated from a simple averages µ A and h(t). As h(t) relaxes much faster than µ F (∆t), eq. (4.35), G(t) relaxes actually much faster than C(t) itself, reaching G eq faster than G F (∆t) = µ A -µ F (∆t). It also shows clearly that G(0) = µ A , the affine shear elasticity, which is always positive, either the system is solid, or even liquid.

The Figs. 4.16 -4.19 show the function G(t) as obtained from eq. (4.35) for our model systems. The common feature is that all the curves shown (except one) reach their limit values quite rapidly, at a time scale of 1. The exception concerns the low temperature AB13 crystal, which oscillates for some time before being stabilized. This is of course for the same reason as its h(t) curve.

To summarize, eq. (4.35) not only gives the function G(t) in a quite accurate way (simple averages), it also provides the equilibrium shear modulus without the necessity of long sampling times ∆t. So it is an advantageous way of studying the elasticity and the viscoelasticity by simulations. However, some tests have to be done using eq. (4.34)

with converged C(t). The resulting G(t) seems quite comparable. From the simulation data G ′ (ω) and G ′′ (ω), the expressions of J ′ and J ′′ are

(4.55)

These relations will be useful for us to convert our G * moduli to J * ones. Their limits are obvious: J ′ (0) = 1/G eq , J ′ (∞) = 1/µ A , and J ′′ (0) = J ′′ (∞) = 0. But J ′′ is not zero for all ω (see the next sub-section). As J ′ (0) must be finite, we shall only discuss this quantity when our system is in solid state.

Simulation results

Turning to our systems, it is an easy task to calculate G * once G(t) is known. It would be interesting to investigate how our G * compare with those predicted by the simple rheological models. For the somewhat delocalization high frequency limit, we use Filon's method [START_REF] Allen | Computer simulation of liquids[END_REF] to perform the sine and cosine transformations. Our results are gathered in Figs.

J′ and J′′ for LJ2d model systems N = 1250, P = 1.0, T = 0.6 Fig. 4.32 Compliance moduli of the LJ2d model system at T = 0.6 (crystal). Two sets of fit are performed using eq. (4.54). Both take the time τ M = 0.22. For J ′ 1 and J ′′ 1 , G = G eq is used. For J ′ 2 and J ′′ 2 , G = µ A is used, and J ′ 2 is shifted by 1/µ A . Now it's interesting to plot the compliance of our systems in the solid state, in order to compare them with the predictions of the K-V model. They are displayed in Figs. 4.32 -4.35. Again, the time parameter τ M can be extracted from the simulation J ′′ curve, at its maximum, i.e. τ M = 1/ω max , where ω max is the frequency for the maximum of J ′′ . This give τ M = 0.22, 0.33, 0.14 and 0.063 for LJ2d, KA2d, LJ3d and AB13 systems respectively. As for G to be used in eq. (4.54), we have the choice of G eq , or µ A , with G eq be the natural choice. The fit using G eq is satisfactory for the mixtures (KA2d and AB13), for J ′′ overall and J ′ in the range ω < ω max . Whereas for the one-component systems, J ′′ is too large. If we use µ A then J ′′ is well fitted, but J ′ must be shifted by 1/µ A . These general trends could depend on temperature. Further investigations are needed.

Topical conclusions

To summarize, we presented in this chapter general considerations on time dependent properties. We first examined the sampling-time dependence of the stress fluctuations. Then we studied the shear stress autocorrelation function C(t), the shear stress mean square displacement h(t), the shear stress relaxation modulus G(t), and the relations between them. The key relation G(t) = µ A -h(t) is used to calculate the shear stress relaxation modulus for our model systems at liquid and solid states.

The resulting curves show satisfactory accuracy, and allow for the computation of the dynamic moduli. This quantity, of high experimental relevance, is shown to be now available by equilibrium MD simulations.