A Pineau

B Chimier

S X Hu

G Duchateau

Modeling the electron collision frequency during solid-to-plasma transition of polystyrene ablator for direct-drive inertial connement fusion applications

Within the inertial confinement fusion (ICF) framework, the solid-to-plasma transition of the ablator arouses increasingly interest, in particular due to the laser-imprint issue. Phase evolution of the ablator is linked to the evolution of the electron collision frequency which is of crucial importance since it drives electron heating by laser energy absorption, and lattice-ion heating due to collisions between electrons and the lattice-ion system. Thus, an accurate description of electron collisions over the whole temperature range occurring in ICF, starting from a few tens of kelvins (solid state) up to tens of millions of kelvins (plasma state), is necessary. In this work, polystyrene ablator is considered and a model of chemical fragmentation is presented to describe the heated polystyrene evolution. Electron collisions are described by electron-phonon collisions in the solid state, and by electron-ion and electron-neutral collisions in plasma state. An effective electron collision frequency valid over the whole range of temperatures reached in ICF experiments is established and discussed. Thermal conductivity is also deduced from collisions in the plasma state and shows a good agreement with the one evaluated by ab initio calculations.

I. INTRODUCTION

To achieve controlled thermonuclear fusion, direct-drive inertial confinement fusion (ICF) has been intensively pursued over the past decades [1][2][3][4] . This process consists of focusing several intense laser beams onto a spherical target in order to compress it to reach temperature and density conditions that are required to generate fusion reactions. The target is usually composed of a cryogenic deuterium-tritium (DT) shell filled with gaseous DT. This cryo-DT shell is covered with an ablator that is usually made of polystyrene 5 . Upon the laser flux, the solid polystyrene layer absorbs the incident energy laser, is converted to a plasma, and is removed from the target leading to the DT compression through the so-called rocket effect. However, the overlapped laser beams can induce spatial variations of intensities on the ICF target surface bringing about a non-uniform ablation of the ablator. This is the so-called laserimprinting process 6 . If the target ablation is non-uniform, perturbations in density and pressure can be induced, becoming the seed of the Rayleigh-Taylor instability. As a consequence, ICF target performances are drastically reduced [7][8][9] because of these spatial perturbations. Despite this, the ablation process is partially described in the current ICF hydro-codes where materials are initially assumed in plasma state but the solid-toplasma transition is generally neglected. As a result, a factor of two was observed between experimental measurements and hydrodynamic simulations of shock velocity modulations 10 .

A first approach 5 modeling the evolution of the ablator properties during the solid-to-plasma transition has recently been developed. The developed model includes electron dynamic upon the laser radiation, and a two-temperature model a) Electronic mail: adrien.pineau@u-bordeaux.fr to describe the electron energy transfer to the lattice-ion system. One of the main conclusions is that the solid-to-plasma transition takes place over a timescale of a hundred of picoseconds, which is not small enough to be neglected in comparison to typical ICF laser pulses, of the order of a few nanoseconds. Among the various parameters involved in this model, the electron collision frequency appears to be one of the most important since it drives the laser absorption efficiency and energy transfer from electrons to the lattice-ion system. Within this first approach, a generic description of the electron collision frequency considering both solid and plasma state has been proposed. Electron-phonon collisions are considered in solid state and assumed to be proportional to the temperature of the lattice-ion system. Electron-ion collisions in plasma are modeled with the usual Spitzer's law. These two asymptotic regimes are then linked by the help of the mean free path. Although it led to some first conclusions, it has to be improved for more accurate ICF simulations, especially accounting for the ablator properties under consideration.

The goal of this work is to present an electron collision frequency model in polystyrene including its transformation toward a dense plasma. First of all, it is assumed that electronelectron collisions, which main influence is to make isotropic the electron distribution, can be neglected. This assumption is supported by the fact that the isotropisation takes place over ∼1 ps [START_REF] Spitzer | Physics of Fully Ionized Gases[END_REF] which is much shorter than the timescale of solidto-plasma transition (∼100 ps) 5 . A four-step scenario of the chemical fragmentation of the ablator with the lattice-ion temperature increase is developed (Sec. II). For each state, an accurate expression of the electron collision frequency is derived (Sec. III). Electron collisions with both acoustic and optical phonons are used for solid state. For the plasma state, standard electron-ion collision frequency is adapted for dense plasmas to account for screening effects. Collisions between electrons and neutrals are considered in states between solid and fully ionized plasma. An effective collision frequency valid over the whole temperature range reached in ICF experiments is finally presented in Section IV and used to calculate the thermal conductivity evolution with respect to the temperature of the ablator. A good agreement is obtained with thermal conductivity evolution deduced by ab initio calculations. Conclusions and outlooks are drawn in Section V.

II. SCENARIO OF POLYSTYRENE FRAGMENTATION

This section is devoted to the establishment of a scenario of polystyrene fragmentation in terms of the lattice-ion temperature. Polystyrene is a polymer molecule composed of carbon and hydrogen which the chemical formula is (C 8 H 8 ) n with C 8 H 8 the repeat unit and n their number. Each molecule of polystyrene consists of the repetition of several hundreds of this repeat unit. The structure of polystyrene is given in Fig 1. The proposed scenario of polystyrene fragmentation is given in Fig. 2. It presents the successive dissociation reactions of the different species occurring during the polystyrene heating. It has been obtained following an estimation of the chemical bond energies. This is a four-step scenario of fragmentation corresponding to three dissociations that are assumed instantaneous (∼ps) compared to the timescale of the solid-to-plasma transition. In addition, the solid-to-plasma transition is assumed to take place over a timescale shorter than the hydrodynamical timescale 5 . This allows to consider the fragmentation as an isochoric process.

C C C H C H C H C H C H H C H H                                          
At low temperature, the polystyrene is in a solid state. When the lattice temperature is about the melting point, it transforms into a liquid state, but the internal structure of the polystyrene molecules is not modified, polystyrene molecules do not break. As they are made of several hundreds of repeat units, it is assumed that it is still possible to define a local periodicity in the liquid state. Phonons still exist and it is assumed collisions are between electrons and phonons when polystyrene is in this state. Therefore, from a collisional point of view, there is no difference between polystyrene in solid and liquid state.

When the so-called ceiling temperature T c is reached, the depolymerization process (inverse process of polymerization) takes place and the polymer dissociates into monomers.

For polystyrene, T c =58 meV (670 K) [START_REF] Mark | Polymer Data Handbook[END_REF] and becomes styrene [START_REF] Feldmann | [END_REF][14][15][16] . The periodicity disappears and phonons are removed. Below 58 meV, matter is in the solid state and collisions are between electrons and phonons. Above 58 meV, matter is composed of atoms and molecules that will be gradually ionized as the temperature increases. In this regime, electrons collide with neutrals and ions.

When the temperature keeps increasing, styrene dissociates into smaller molecules that dissociate themselves into even smaller molecules up to a fully dissociated state composed only of partially ionized atoms. Each dissociation step is defined by a dissociation temperature. Regarding styrene, the main channel of dissociation gives one molecule of benzene and one molecule of acetylene 17 . This dissociation is expected in regards to the structure of the styrene (see Fig. 2). It is basically composed of a phenyl group and an alkyl group. The phenyl group is difficult to break due to its structure stability. Electrons involved in the carbon-carbon bonds are delocalized all over the cycle and the six bonds are equivalent with an energy between C -C and C = C bond energies. The alkyl group is composed of one C -C bond and one C = C bond. Due to a more important atomic orbital overlap in C = C than in C -C, it results the C = C bond is stronger than the C -C bond. Thus, the C -C alkyl bond appears as the weakest and then first breaks. Measurements of styrene fragmentation induced by 308nm and 355nm lasers showed that the dissociation takes place at 18 2.32±0.1 eV and 19 2.43±0.05 eV respectively. Therefore, it is assumed that the styrene dissociation happens between 2.32 eV and 2.43 eV. Between these two temperatures, both styrene and the benzene-acetylene mixture are considered.

Benzene and acetylene fragmentations should be then considered. Differently from the acetylene fragmentation which products can be easily identified due to its simple structure, the benzene fragmentation is more complicated to determine because of its stability. However, it is possible to assume the only bonds present in the fragments of benzene are C -C, C = C, C ≡ C and C -H. In addition, it is assumed that the aromatic carbon-carbon bond energy are between C -C and C = C bond energy. Table I proposes an overview of the different bond energies that can be found in the literature to our knowledge. It is observed that most of bond energies are between 4 and 8 eV. Therefore, a lot of molecules and radicals can appear and disappear on this small temperature range (in regards to the solid-plasma transition temperature range). In order not to deal with this very complex mixture evolving very fast, it is set that benzene and acetylene molecules dissociate directly into a carbon-hydrogen mixture between 4 and 8 eV. This scenario is supported by the fact that the target remains only a few picoseconds in this temperature range 5 . As for styrene dissociation, between 4 and 8 eV both the benzeneacetylene mixture and the carbon-hydrogen mixture are considered.

A summary of the various species that are considered in terms of the lattice-ion temperature T il is given in Fig. 3.

Considering polystyrene at standard density [START_REF] Bandrup | Polymer Handbook[END_REF]

ρ PS = C C CH CH H C C H C H H C H H                                       

III. ELECTRON COLLISIONS

It was established that two regimes can be exhibited in terms of the lattice-ion temperature. If it is smaller than 58 meV, collisions are between electrons and phonons. If it is higher, electrons collisions are with both neutrals and ions.

A. Electron-phonon collisions

Polystyrene is a dielectric material which is assumed to be ideal (no impurities nor defects) implying that only electron-phonon collisions take place. As the unit cell of the polystyrene is composed of several atoms, both acoustic and optical phonons have to be taken into account. The collision between electrons and phonons results in the absorption or emission of a phonon by an electron. The electron acoustic phonon (EAP) collision frequency and the electron optical phonon (EOP) collision frequency are given by Eqs. (1) and (2) respectively. The derivation of these collision frequencies is presented in Appendix A.

The (EAP) collision frequency is given by 28 :

ν A ( #» k ) = U 2 k B T il m * 2πρc 2 s h3 k k 2 + m * c s h 2 (1)
with k B the Boltzmann constant, h = h/2π the reduced Planck constant, the lattice-ion temperature T il , the material density ρ, the sound speed c s and U the deformation potential.

Electron optical phonon (EOP) collision frequency is given by [START_REF] Fischetti | [END_REF] :

ν LO ( #» k ) = e 2 m * ω LO h2 k 1 ε ∞ r - 1 ε 0 r × g LO ln 1 + 2m * ω LO /hk 2 + 1 1 + 2m * ω LO /hk 2 -1 + (g LO + 1) ln 1 + 1 -2m * ω LO /hk 2 1 -1 -2m * ω LO /hk 2 (2)
The first and second terms in the brackets correspond to the absorption and emission of an optical phonon, respectively. It is expressed in terms of the elementary charge e, the optical phonon energy hω LO , the dielectric constant ε 0 r and ε ∞ r at low and high frequency, respectively, and the phonon distribution function

g LO = [exp (hω LO /k B T il )-1] -1 .
For polystyrene, the dielectric constants are 12 ε 0 r = 2.5 and 29 ε ∞ r = 1.01. In order to get a macroscopic description for ICF purpose, these frequencies are averaged over the Brillouin zone weighted by the electron distribution function f ( #» k ) and the density of states g( #» k ):

ν = k BZ k min ν( #» k )g( #» k ) f ( #» k ) d #» k k BZ k min g( #» k ) f ( #» k ) d #» k (3)
For acoustic phonons, the different assumptions allowing to obtain Eq. ( 1) have led not to make differences between electron collisions by phonon emission or by phonon absorption. Therefore, one has k min = m * c s /h for EAP collision frequency because an electron cannot emit a phonon if its momentum is smaller than the phonon momentum. For EOP collision frequency, one has k min = 0 for absorption because an electron can absorb a phonon regardless of its momentum, and k min = k LO for emission in the same way as for acoustic phonons. Note that in practice Eq. ( 3) is split into two parts in order to calculate the emission part and the absorption part separately. The density of states g( #» k ) in Eq. ( 3) is taken in the free electron gas approximation g( #» k ) = 2V /8π 3 and the electron distribution function is given by the Fermi-Dirac distribution function f (

#» k ) = [exp ((E( #» k ) -µ)/k B T e ) + 1] -1
where the energy of the bottom of the conduction band has been set to zero. An approximation of the chemical potential µ is derived in Appendix C.

Two optical phonons are considered for polystyrene with hω LO,1 = 33.53 meV and hω LO,2 = 88.73 meV. The derivation of the relation dispersion leading to the phonon number and the optical phonon energies is given in Appendix B. Finally, the total electron collision is given by the Mathieusen rule [START_REF] Ashcroft | Physique des Solides[END_REF] :

ν e-ph = νA + νLO,1 + νLO,2 (4)
Figure 4 presents the total electron-phonon collision frequency given by Eq. ( 4) as a function of T e and T il , for T il ≤ 58 meV and for T e > T il because T e < T il is usually not observed in ICF experiments. It is set n e = 10 10 cm -3 corresponding to the number of electrons that are present in the conduction band of dielectrics [START_REF] Smith | [END_REF] . Note that the electronphonon collision frequency depends on the electron density due to the chemical potential. Within ICF conditions, values of both electron temperature and density are such that 5 the chemical potential is defined by Eq. (C6). Thus, the chemical potential exhibits a logarithmic dependence and electron density variations will slightly modify the collision frequency.

The increase in the collision frequency with the latticeion temperature is due to phonon population increasing. The higher the temperature T il , the larger the phonon density, and the larger the number of collisions. The evolution of the collision frequency as a function of T e is similar to the behavior of the collision frequency as a function of the electron energy E( #» k ), especially the one related to the emission of an optical phonon by an electron,that is defined by Eq. ( 2). Note that Eq. ( 2) can be written with respect to the electron energy by using E( #» k ) = h2 k 2 /2m * . This collision frequency first increases with the electron energy up to a maximum where the electron energy becomes of the order of magnitude of the phonon energy, and decreases for higher electron energies. Such a behavior is standard for electron optical phonon collision frequency 32 . In addition, values of the total collision frequency lie in the range of other dielectrics, as silicon dioxide 32 .

B. Electron collisions in the dense plasma state

In the plasma state, electrons can collide with ions and neutrals. These collisions depending on the particles (electrons, ions and neutrals) densities and the Coulomb logarithm, the latter are first evaluated before presenting the collision frequencies.

Ionization

This subsection is devoted to the presentation of the models used to describe ionization of styrene, benzene, acetylene, carbon and hydrogen. Thermal ionization models are considered.

Ionization of styrene, benzene and acetylene is described by the Saha equation. Because styrene exists for temperatures smaller than 2.43 eV (see Fig. 3), only the first ionization is considered. In the same way, benzene and acetylene exist for temperatures smaller than 8 eV and only the first ionization state are considered. Ionization energies for styrene, benzene and acetylene are 17,33 8.46 eV, 9.24 eV and 11.42 eV, respectively. Ionization of the CH mixture is obtained by ab initio calculations. As it exists starting from 4 eV, all of the ionization states are considered which ionization energies range from 11.26 eV to 490 eV. Saha equations and the model associated to ab initio calculations are detailed in Appendix D.

Figure 5 (a) shows the averaged ionization state < Z > as a function of temperature. Assuming that the plasma behavior, i.e. when density of charged particles in the medium becomes large enough so that its properties are modified, appears from an ionization degree 34 α = n i /n T ≈ 10 -3 , it can be deduced from Fig. 5(a) that the plasma behavior starts from around 650 meV. Note that for low temperatures, n e ≈ n i and α ≈ < Z >. Therefore, between 58 meV and 650 meV, styrene can be considered as a dense neutral gas.

Figure 5 (b) displays the density of electrons, neutrals and ions between 58 meV and 1 keV for T e = T il . It appears that neutrals are mainly present for temperatures smaller than around 4 eV. A drop in the neutral density is observed beyond this temperature because it has been defined as n n = n Tn i where n i = min (n e , n T ) for the CH mixture. Because the ionization model of CH gives only the electron density, it is impossible to access to the density of each ion population. Therefore, it is assumed that when electron density is higher than the total density, it means that all of the atoms are ionized so the density of neutral particles drops to zero. This may result in an overestimation of the ion density and an underestimation of the neutral density for temperature around 10 eV.

It is possible to deduce that electron-neutral collisions are the dominant process for temperatures up to 4 eV, there is a competition between electron-neutral and electron-ion for temperatures between 4 and 8 eV, and beyond, electron-ion collisions are the dominant process.

Usually, the Coulomb logarithm ln Λ appearing in the electron-ion collision frequency is obtained assuming particles interacting in a Coulomb potential, i.e. without screening effects. In such a potential, the Coulomb logarithm reads in the well-known form : with b min and b max the lower and upper cut-off lengths respectively. In dense plasmas, screening effects have to be taken into account and the previous definition is no longer valid.

(ln Λ) C = 1 2 ln 1 + b 2 max b 2 min (5)
Based on an evaluation of the stopping power, a definition of the Coulomb logarithm available for any potential U(r) has been proposed [START_REF] Ordonnez | [END_REF] :

ln Λ = 1 b 2 min b max 0 b sin 2 θ 2 db (6)
with b the impact parameter. The scattering angle θ of the electron is defined by [START_REF] Landau | Physique Théorique : Mécanique[END_REF] :

θ = π -2b +∞ r 0 1 r 2 1 - 2U(r) k B T e - b 2 r 2 -1/2 dr ( 7
)
where r 0 is the distance of closest approach and satisfies the relation 1 -2U(r 0 )/k B T eb2 /r 2 0 = 0. Interactions between ions of charge Ze and electrons are modeled by the Debye-Hückel potential U DH (r) which is the generalization of the Coulomb potential including screening effects [START_REF] Ramazanov | [END_REF] :

U DH (r) = Ze 2 r exp - r λ DH (8)
The Debye-Hückel length λ DH characterizing the screening effects is written in a two-temperature model 38 :

1

λ 2 DH = 4πn e e 2 k B T 2 e + T 2 F + 4πn i Z 2 e 2 k B T il (9)
where the Fermi temperature T F allows to take into account the degeneracy correction appearing in a dense plasma within a good approximation 39 . By using Eq. ( 8), the distance of closest approach for b = 0 is given by :

r ⊥ = λ DH Wp 2Ze 2 λ DH k B T e ( 10
)
where Wp is the principal branch of the Lambert Wfunction [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF] . Finally, the parameters b min and b max in Eq. ( 6) are defined as :

b min = 1 2 max [r ⊥ , λ B ] ; b max = λ 2 DH + a 2 i ( 11
)
with λ B = h/ √ m e k B T e the De Broglie length [START_REF] Huba | NRL Plasma Formulary[END_REF] ensuring to respect the Heisenberg uncertainty principle [START_REF] Landau | Quantum Mechanics[END_REF] and a i = (3/4πn i ) 1/3 the ion sphere radius allowing to take into account the effects of a high density plasma [START_REF] Gericke | [END_REF] . Note that the integration in Eq. ( 6) is done up to b = b max and not b → +∞. It has been assumed that U(r) = 0 for r > λ D leading to θ = 0. This ensures the convergence of the Coulomb logarithm while having a screened potential inside the Debye sphere and without double counting screening effects.

By modeling a CH plasma over a wide range of temperatures and densities occurring in ICF with ab initio calculations based on quantum molecular dynamics, a fitted formula for a generalized Coulomb logarithm (ln Λ) QMD has been obtained 44 . These ab initio calculations were performed using quantum molecular dynamics methods where the electronic potential is evaluated by the density functional theory. A CH plasma was modeled over a wide range of temperatures and densities relevant of ICF applications. The Coulomb logarithm was deduced from the thermal conductivity calculated in terms of the Onsager coefficients using the Kubo-Greenwood formalism. This allows to examine the validity of the previous two Coulomb logarithm evaluations designed to be introduced in models of laser-matter interaction.

Figure 6 compares the usual Coulomb logarithm (ln Λ) C defined by Eq. ( 5), the Coulomb logarithm (ln Λ) DH defined by Eq. ( 6)-( 7) using the Debye-Hückel potential, and the QMD Coulomb logarithm (ln Λ) QMD for T e = T il . The latter is plotted only for T > 8 eV since it has only been developed for a CH mixture. It is reminded a pure CH mixture is considered in this work only for T > 8 eV (see Fig. 3).

It appears that (ln Λ) C cannot be used for the range of temperature and density of ICF. Despite a reasonable order of magnitude for temperatures higher than 100 eV, it underestimates by a factor from 2 to 4 for temperatures between 8 and 100 eV. Moreover, it shows a non-physical behavior for temperatures smaller than 8 eV as it increases when the temperature decreases. This is due to the cut-off lengths. For (ln Λ) C , they are defined by Eq. ( 11), as for (ln Λ) DH , except that r ⊥ = Relatively good agreements are observed between (ln Λ) DH and (ln Λ) QMD . Despite differences of 40% can be observed for temperatures around 10 eV, they are quite similar for temperatures higher than 40 eV. For temperatures smaller than about 2 eV, the definition of (ln Λ) DH proposed in this work goes wrong as it leads to (ln Λ) DH 1. Thus, in that case, we impose (ln Λ) DH = 1. However, the consequences of not correctly evaluating ln Λ in this domain are relatively not important because the dominant processes are collisions between electrons and neutrals.

It can be noticed that the slope-changing points at 2 eV for ln Λ C and at 3 and 20 eV for ln Λ DH correspond to the temperatures where the switch between λ B and r 0 takes place in Eq. ( 11).

Collision frequencies

For electron-neutral collisions, a simple model based on the classical diffusion cross section σ is used [START_REF] Huba | NRL Plasma Formulary[END_REF] :

ν en = n n σ k B T e m e ( 12
)
with n n the density of neutral. The diffusion cross section is expressed by assuming neutrals as hard spheres and neglecting the deformation of the electronic cloud during the collision with an electron. Then 45 , σ = πa 2 with a the characteristic size of the neutral ranging from several tens to several hundreds of ångströms depending on if atoms or molecules are considered. Polarization effects due to the interaction between an electron and the electric dipole moment of neutral particles are indeed neglected in our expression of electronneutral collisions. This assumption is supported by the following considerations. Styrene, benzene and acetylene are nonpolar molecules, i.e. the total permanent dipole moment is equal to zero. Because permanent dipole moments of each C-H bond cancel each other out due to the molecule symmetry. In addition, models based on induced polarization are relevant for temperatures higher than several eV for which the plasma is almost fully ionized in our model and electronneutral collisions are negligible. Thus the polarization influence for electron scattering onto neutrals is expected to be negligible within the present conditions. For laser intensities in ICF applications, the electron quiver velocity is generally lower than the electron thermal velocity. Electron-ion collisions can thus be described in the classical regime [START_REF] Silin | [END_REF][47][48][49][50] . Since the solid-to-plasma transition takes place on a timescale smaller than hydrodynamic timescale 5 , the plasma density remains close to the solid density. The Spitzer's model is no longer valid and the Lee and More model is used 38 :

ν ei = 2 √ 2π 3
Z 2 e 4 n i √ m e (k B T e ) 3/2 (ln Λ) DH (1 + exp (-µ/k B T e ))F 1/2 (13)
where F 1/2 is the complete Fermi-Dirac integral:

F 1/2 µ k B T e = +∞ 0 √ t 1 + exp (t -µ/k B T e ) dt ( 14
)
and µ is defined by Eq. (C6)-(C8). The Coulomb logarithm (ln Λ) DH defined by Eq. ( 6)-( 8) is used to take into account the screening effects.

Figure 7 shows the total electron collision frequency ν ei + ν en as a function of electron and lattice-ion temperatures between 58 meV and 1 keV. As for the electron-phonon collision frequency, it is plotted for T e > T il since it is usually not observed T e < T il in ICF experiments.

The general behavior of the electron collision frequency is relatively standard. It increases with the lattice-ion temperature because when T il increases, the successive chemical dissociations of the ablator lead to a growing number of atoms and molecules. Concerning its evolution as a function of the electron temperature, it first increases between 58 meV and around a few eV. The electron temperature is not large enough to ionize molecules, so the density of neutrals remains quite constant and the electron collision frequency increases with the electron temperature (see Eq. ( 12)). For higher electron temperatures, atoms and molecules are ionized and the decrease in the electron collision frequency corresponds to the characteristic behavior of electron-ion collisions, i.e. the electron collision frequency evolves as T -3/2 e (see Eq. ( 13)). It appears electron-neutral collisions are the dominant process for electron temperature smaller than a few eV and beyond, electron-ion collisions are the dominant process.

IV. APPLICATIONS OF DERIVED COLLISION FREQUENCIES A. Eective electron collision frequency

A collision frequency ν e f f valid over the whole range of temperature occurring in inertial confinement fusion can be deduced from what has been presented previously. If T il ≤ T c , with T c = 58 meV, the ablator is in the solid state and electrons collide with phonons. Otherwise, it is a dense partially ionized gas and electrons collide with neutrals and ions. Thus ν e f f reads :

ν e f f = ν eph x PS + (ν ei + ν en )x AM ( 15
)
where ν eph , ν en and ν ei are defined by Eq. ( 4), Eq. ( 12) and Eq. ( 13), respectively. Fractions x PS and x AM , which are both function of the temperature, represent fractions of polymer and atoms/molecules densities respectively. They are the dominant species for temperatures smaller and higher than T c respectively. These fractions satisfy x PS → 1 and x AM → 0 when T T c , x PS → 0 and x AM → 1 when T T c , and x PS = x AM = 1/2 when T = T c . Variations of the species fractions around the transition are described by a second-order polynome. Such a definition allows to have an effective collision frequency that is continuous during the transition from the solid state to the plasma state taking place over a range δ T around T c . This temperature range being much smaller than the one considered to plot the effective collision frequency 51 in Fig. 8, this transition appears as a jump.

It can be seen that when the ablator is in the solid state, the collision frequency increases with the temperature. This is due to the excitation of phonon modes. At T = 58 meV, a drop by a factor 10 is induced by both phonon removal and small electron temperature in the plasma state. Then, due to the increase of electron density and temperature, the collision frequency increases up to around 10 16 s -1 for temperatures between 15 eV and 20 eV. Finally, beyond 20 eV, the collision frequency decreases with a typical plasma behavior.

The generic electron collision frequency ν c developed in 5 is also presented in Fig. 8. In that previous work, the collision frequency is defined as ν c = ν 0 T il /T 0 with ν 0 = 10 14 s -1 and

T 0 = 300 K if T il <0.1 eV. If T il >0.1 eV, ν -1 c = ν -1 mp f + ν -1 e-i which is defined in terms of ν m f p = v e n
1/3 a the mean free path with v e and n a = 3 × 10 22 cm -3 the electron thermal velocity and the atomic density respectively, and ν e-i = 3.10 -6 Zn e T -3/2 e ln Λ the usual Spitzer's electron-ion collision frequency with T e in eV and assuming Z ln Λ = 10, and n e = 3 × 10 22 cm -3 which corresponds to the order of magnitude of the critical density used in hydro-codes as initial condition. Note that this is a generic approach in the sense that it does not account for material specificities as its chemical fragmentation.

Large differences can be observed between the generic collision frequency and the one of this work in terms of behavior and amplitude depending on the temperature. At low temperatures, this is due to the fact that evaluating the electronphonon collision frequency as ν 0 T il /T 0 includes only collisions between electrons and acoustic phonons. Collisions between electrons and optical phonons are neglected. However, it has been shown 32 that in dielectrics, electron-optical phonon collisions are the dominant processes and electron-acoustic phonon collisions only gives a small correction to the total electron-phonon collision frequency. At high temperatures, the difference by a factor of 10 is explained by the difference by a factor of 10 between electron densities used in the various collision frequencies. As mentioned previously, a large jump around T c seems to appears. In addition to the fact that the energy range of Fig. 8 does not allow to see the continuous transition at T c , this jump also appears because the total collision frequency is plotted for T e = T il which does not represent relevant conditions of ICF. Within ICF conditions, when polystyrene depolymerizes at 0.058 eV, electron temperature is about a few eV 5 . Thus, having an electron temperature of a few eV instead of 0.058 eV leads to a decrease of the electronphonon collision (see Fig. 4) and an increase of the electronneutral collision frequency (see Fig. 7), so a reduction of the jump.

Despite such discrepancies, a relatively small difference in the solid-to-plasma transition timescale is expected because it mainly depends on the time required to reach the critical plasma density, the latter being mainly imposed by the photoionization rate 5 that is not considered in this work. To include photoionization processes, it is necessary to model the electron dynamics coupled to a description of the temporal evolution of both electron and ion-lattice temperatures. This is beyond the scope of the present work as done in 5 .

However a significant influence on both optical and thermal matter responses is expected. It appears from Fig. 8 that the effective collision frequency is mostly smaller than the laser pulsation considering a 351nm laser. This leads to a conductivity and a laser absorption proportional to the collision frequency. Thus, higher conductivity and laser absorption for smallest (solid state) and highest temperatures (full plasma state) are expected. In the transition zone for temperatures around a few eV, a smaller conductivity is expected. In addition, note that their value is also strongly driven by the free electron density that, in case where the laser interaction is included in the modeling, is driven by photoionization processes. Such a study is out of the scope of the present work and will be adressed in a forthcoming study.

B. Plasma thermal conductivity

In order to evaluate the thermal response of the ablator, this section is devoted to the evaluation of the thermal conductivity. Here, only the plasma state for temperatures higher than 1 eV is considered so that the electron density is large enough to ensure that the thermal conductivity is mainly due to electron collisions (which are well described in the present modeling). Below 1 eV, the ionized degree is very low (see Fig 5), the contribution of neutrals to the thermal conductivity is thus significant, and the present modeling no longer holds in that case.

Starting from the electron collision frequency in the plasma state, it is possible to deduce the thermal conductivity. This is a quantity of interest because it is involved (i) in the heat transport in plasmas, (ii) it is related to the evolution of hydrodynamic instabilities and (iii) it is one of the key parameter for target ablation when the electron density becomes higher than the critical density.

The thermal conductivity κ th in a dense plasma can be described by the Lee & More model 38 . However, electronelectron collisions are neglected in this model leading to an overestimation of the thermal conductivity. To include these collisions, the Lee & More model is weighted by a function R(Z) that can be defined as 52 R(Z) = (Z + 0.25)/(Z + 4.26). Thus, the thermal conductivity κ th is given by :

κ th = n e k B (k B T e ) m e ν e f f A µ k B T e R(Z) (16)
with the function A defined by : A good agreement is observed in Fig. 9 for temperatures higher than 8 eV. Differences for temperatures smaller than 8 eV come from the fact that the ab initio model considers a CH mixture for any temperatures and neglects the chemical fragmentation of the ablator below 8 eV whereas this fragmentation is taken into consideration in the present work (see Fig. 3). Below 8 eV, the ablator is mainly composed of benzene and acetylene, and below 2.3 eV, it is mainly composed of styrene. Thus, the total density is smaller leading to an electron density smaller than the one in ab initio calculations. It results κ th becomes smaller than κ QMD for temperatures below 8 eV. This origin of the differences is confirmed by comparing κ th,w/o F to the other two as it corresponds to the thermal conductivity defined by Eq. ( 16) but where the fragmentation processes are not taken into account and a plasma composed of carbon and hydrogen is considered.

A µ k B T e = 20 9
F 4 1 -16F 2 3 /15F 2 F 4 [1 + exp (-µ/k B T e )]
This smaller thermal conductivity will lead to a smaller heat flux and thus a slower heating of the ablator. Because at the same time, the ablator is ionized by the laser pulse, the conduction zone will be smaller when the critical density will be reached. In the end, thermal smoothing occurring in the conduction zone will be less efficient and the laser imprint will be more important as more perturbation modes will be able to be transmitted to the ablation front.

V. CONCLUSION

Two collision frequencies relevant from inertial confinement fusion conditions have been presented in this work. An electron-phonon collision frequency taking both acoustic and optical phonons into consideration has been derived in order to describe electron collisions when the ablator is in the solid state. Electron collisions when the ablator is in the plasma state are described considering both electron-neutral and electron-ion collisions. For the latter, a modified Coulomb logarithm has been used in order to account for screening effects, showing a good agreement when compared to ab initio calculations. By coupling these two collision frequencies with a four-step fragmentation model for polystyrene ablator, which describes the successive chemical dissociations of polystyrene leading to a carbon-hydrogen mixture, an effective collision frequency valid over the whole range of temperatures reached in inertial confinement fusion experiments has been obtained. This effective collision frequency has been compared to a generic collision frequency used to describe the solid-to-plasma transition of the ablator. Relatively large differences in behavior and amplitude depending on temperature are observed but the timescale of solid-to-plasma transition should not be strongly modified as it mainly depends on the photo-ionization rate. However, larger modifications of thermal and hydrodynamic responses are expected since they are driven by the electron collision frequency due to laser heating. In addition, this effective collision frequency has been used to evaluate the thermal conductivity for temperatures above 1 eV. When compared to a thermal conductivity obtained by ab initio calculations, good agreements have been observed for temperatures larger than 10 eV. But below 10 eV, a smaller conductivity was found, coming from the consideration of the chemical fragmentation of the ablator, which could lead to a more important laser imprint.

Such an accurate description of the collision frequency will be used in future studies dedicated to model the laser-induced solid-to-plasma transition. E(k) = h2 k 2 /2m * , with m * the effective electron mass assumed to be equal to the free electron mass. Electron-phonon collisions are assumed to be limited to the first Brillouin zone as electron energies do not exceed the energy E(k BZ ) associated to the Brillouin zone wave vector k BZ . This implies that umklapp processes, corresponding to collisions where the reciprocal lattice vector has to be considered in crystal momentum conservation [START_REF] Ziman | Electrons and Phonons[END_REF] , are not considered. Because the number of phonons and their respective energies are calculated in one dimension (see Appendix B), the Brillouin zone wave vector reads k BZ = 2π/a with a = 2d C-C the size of the unit cell and 54 d C-C = 1.53 Å the distance between two atoms of the carbon chain. This leads to E(k BZ ) ≈ 16 eV. Note that it has been assumed in Section II polystyrene is considered for T il 58 meV. In a one-temperature model, T e = T il and it clearly appears electrons have an energy much smaller than E(k BZ ). But, in the laser-matter interaction framework, electrons are heated first, well before the ion-lattice system. However, even in that case, it has been shown 5 that when T il ≈ 58 meV, electron energies are the order of magnitude of E(k BZ ). Therefore, both in one-temperature and two-temperature model, electron energies are not higher than E(k BZ ), thus confirming the description with one band in the first Brillouin zone.

Assuming the electron distribution function f e can be only modified by collisions, the Boltzmann equation reads in the relaxation time approximation

∂ t f e = -( f e -f 0 e )/τ( #» k ) 30 . f 0 e
is the unperturbed distribution function and τ( #» k ) is the relaxation time defined as:

1 τ #» k = ν( #» k ) = V 8π 3 S -#» k , #» k + S + #» k , #» k d #» k (A1) The term S ± ( #» k , # »
k ) corresponds to the electron transition rate from a state #» k to a state #» k by absorbing (-) or emitting (+) a phonon of momentum #» q . It is given by the Fermi golden rule [START_REF] Ziman | Electrons and Phonons[END_REF] :

S ± ( #» k , #» k ) = 2π h | #» k |H ± e-ph | #» k | 2 δ E #» k -E #» k ∓ hω β (q) (
A2) with E #» k and E #» k the energies of the states #» k and #» k respectively and ω β (q) the frequency of the β -mode phonon. If the umklapp processes are neglected, matrix elements read [START_REF] Ziman | Electrons and Phonons[END_REF] :

#» k |H ± e-ph | #» k = h 2ρV ω β (q) 1/2 g β + 1 2 ± 1 2 1/2 ×I #» q ,p ( #» k , #» k ) δ #» k - #» k ± #» q (A3)
which depends on the density ρ, the volume V and the lattice temperature T il through the Bose-Einstein distribution

g β = [exp hω β (q)/k B T il -1] -1 . δ #» k - #» k ± #» q repre-
sents the momentum conservation and I #» q ,p ( #» k , #» k ) the coupling between the lattice potential and the electronic states. It is calculated in the deformation potential approximation as a compromise between a simple model and good results leading to [START_REF] Ziman | Electrons and Phonons[END_REF] :

I #» q ,p ( #» k , #» k ) = i #» q • #» e #» q ,p U (q, ω β (q)) (A4)
The potential U (q, ω β (q)) depends on whether an acoustic or an optical phonon is considered. For an acoustic phonon, it is constant equal to a few eV in the limit of small momentum transfer and one has U = δ E /∆. It corresponds to a measure of the deformation potential δ E due to the lattice dilatation ∆ induced by the acoustic mode [START_REF] Ziman | Electrons and Phonons[END_REF] . Neither U nor δ E nor ∆ can be evaluated simply for dielectrics. Based on cross section evaluation, Fishcetti et al. [START_REF] Fischetti | [END_REF] set U ≈ 3.5 eV for SiO 2 , that is around 40% of the band gap energy E g . By analogy with metals where the deformation potential is proportional to the Fermi energy, it is assumed that the deformation potential of dielectric materials is proportional to the band gap. With E g = 4.05 eV 55 for polystyrene, this leads to U = 1.62 eV. This ensures to obtain the standard order of magnitude of a few eV for the deformation potential [START_REF] Ziman | Electrons and Phonons[END_REF] . For an optical phonon, it is given by the Fröhlich formalism [START_REF] Ziman | Electrons and Phonons[END_REF]56,57 :

U (q, ω β (q)) = -4πe 2 ρ ω 2 β ( #» q ) 1 ε ∞ r - 1 ε 0 r 1/2 1 #» q 2
(A5) where ε 0 r and ε ∞ r are the dielectric constants at low and high frequencies respectively. The momentum conservation equation giving

#» k 2 = #» k 2 + #» q 2 ± 2kq cos (θ )
, it is possible to write the term representing the energy conservation as:

δ E #» k -E #» k ∓ hω β (q) = m * h2 kq δ ± cos (θ ) + q 2k ∓ m * ω β (q) hkq (A6)
where the parabolic band approximation has been used E #» k = h2 #» k 2 /2m * with m * the effective electron mass assumed to be equal to the free electron mass m e . Finally:

S ± ( #» k , #» k ) = πqm * U 2 (q, ω β (q)) h2 ρV ω β (q) g β + 1 2 ± 1 2 × δ ± cos (θ ) + q 2k ∓ m * ω β (q) hkq (A7) It appears S ± ( #» k , # »
k ) depends only on q and cos θ and Eq. (A1) becomes

ν( #» k ) = V 4π 3 ∞ 0 1 -1 S -( #» k , #» k ) + S + ( #» k , #» k ) q 2 dqdcos (θ ) (
A8) Assuming q → 0, this leads to a linear relation dispersion for acoustic phonon: ω LA ( #» q ) ≈ c s q, and constant dispersion relation for optical phonon: ω LO ( #» q ) ≈ ω LO . The demonstration of these approximations is given in Appendix B. Finally, plugging Eq. (A7) and (A8), this leads to the following collision frequency for acoustic phonons [START_REF] Fischetti | [END_REF] :

ν AC ( #» k ) = U 2 k B T il m * 2πρc 2 s h3 k k 2 + m * c s h 2 (A9)
where it has been supposed the temperature is much higher than the phonon energy leading to a simplified phonon distribution: n AC ≈ n AC + 1 ≈ k B T il /hc s q. For optical phonons one it is possible to write Eq. (C1) in such a way that the Fermi-Dirac integral F 1/2 appears:

y = +∞ 0 √ t e t-x + 1 = F 1/2 (x) (C3)
As a consequence, looking for an expression of the chemical potential becomes equivalent to find the inverse function of F 1/2 . This is done by taking approximations of F 1/2 in the two asymptotic cases 61,62 :

x < 1.13 ⇔ y < 1.52 ⇒ F 1/2 (x) = √ π 0.5 + 2e -x (C4)

x > 1.13 ⇔ y > 1.52 ⇒ F 1/2 (x) = 2 3

x 2 + π 2 6

3/4 (C5)

Because these approximations lead to a discontinuity for x = 1.13, y = 1.52, it is possible to get a continuous approximation by weighting the previous two expressions with the logistic function V (y, y 0 , k) = [1 + e -a(y-y 0 ) ] -1 around these values. Thus, the following approximation of the chemical potential

Appendix D: Ionization

Saha equations are used in order to describe ionization of styrene, benzene and acetylene because it has already been shown valid for ICF conditions 44 . They were originally established assuming thermal equilibrium T e = T il but it is possible to establish modified Saha equations for T e = T il as follows [START_REF] Chen | [END_REF] :

n j+1 n e n j = 2 2πm e k B T e h 2 3/2 Q j+1 (T e )
Q j (T e ) exp -U j+1 k B T e (D1) where U j+1 is the ( j + 1)th ionization potential, n j and n j+1 are the densities of atoms ionized j and j + 1 times respectively, and Q j (T e ) and Q j+1 (T e ) are their partition functions reduced to the translational, rotational and vibrational contributions. As a first approximation, the rotational and vibrational parts can be neglected leading to the ratio Q j+1 (T e )/Q j (T e ) equal to unity. Because the mass of the atoms ionized j and j +1 times can be considered as the same, the translational partition functions cancel. Note that the prefactor 2 corresponds to the electron spin degeneracy. This equation involves three densities and the problem is closed with the two following relations: N is the number of ionized states taken into account. In an exact description, it is equal to the number of electrons of the particle. Z j = j and n j are the charge state and the density of particles ionized j times respectively. n n is the density of neutrals and n T is the total density. The averaged ionization state < Z > is given by:

< Z >= n e n T (D3)
The ionization model of the carbon hydrogen mixture is obtained by ab initio calculations and the averaged ionization state < Z > is given by 44 :

< Z >= ξ 2 1 + 4 ξ -1 Z max (D4)
with Z max = (Z H + Z C )/2 = 3.5 and:

ξ = α 0 n i 2πm e k B T h 2 3/2 exp - f Z k B T (D5)
while the average-ionization potential is given by:

f Z = α 1 + α 2 k B T 1 + 3Γ 0 1/4
-1 + α 3 (k B T ) 0.9

+ α 4 r 0 + α 5 r 2 0 + α 6 r 3 0 k B T (D6)
with r 0 = a i /a B the ion-ion interdistance defined as the ratio of the ion sphere radius a i = (3/4πn i ) 1/3 over the Bohr radius a B , Γ 0 = 1/r 0 k B T the ion-ion coupling parameter in atomic units, and α i fitting parameters: α 0 = 87.222, α 1 = 10.866, α 2 = -28.412, α 3 = 17.915, α 4 = -2.422, α 5 = 0.595 and α 6 = -2.369 × 10 -2 . Finally, the electron, ion and neutral densities n e , n i and n 0 are given by: n e =< Z > n T ; n i = min (n e , n T ) ; n 0 = n Tn i (D7) Note that defining n i by this way does not allow to access to densities of each ion existing n H + , n C + , n C 2+ , etc. As the rigorous relation linking the electron density and the density of each charge state is given by Eq. (D2), if n e is known, it is impossible to deduce each ion density n j .

nFIG. 1 .

 1 FIG. 1. Chemical structure of polystyrene. The circle indicates that the six C -C bonds in the phenyl group are equivalent.

FIG. 4 .

 4 FIG. 4. Electron-phonon collision frequency as a function of the electron and lattice-ion temperatures assuming a free electron density n e = 10 10 cm -3 . It is plotted for T e > T il which account for realistic physical conditions.

FIG. 5 .

 5 FIG. 5. (a) Averaged ionization state < Z > as a function of the temperature. (b) Density of species as a function of the temperature. The dotted curve represents the electron density n e , the dashed curve the total ion density n i , the dash-dot curve the total neutral density n n and the solid curve the total particles density n T = n n + n i . Vertical dashed lines indicate dissociation temperatures at 2.32, 2.43, 4 and 8 eV. Both (a) and (b) are plotted for T e = T il .

FIG. 6 .

 6 FIG. 6. Coulomb logarithm for various models as a function of the temperature. Calculations have been done with T e = T il . It is reminded vertical dashed lines indicate dissociation temperatures at 2.32, 2.43, 4 and 8 eV.

  2Ze 2 /k B T e . When T e becomes very small, one has b min ≈ λ B and b max ≈ λ DH . Thus, b 2 max /b 2 min ≈ m e k 2 B T 2 e /π h2 e 2 n e using the fact that λ DH ≈ k B T e /4πn e e 2 for small electron temperatures. Because a small variation of temperature leads to large variation of n e (see Fig. 5), the ratio b 2 max /b 2 min is driven by the behavior of n e . So b 2 max /b 2 min → +∞ and (ln Λ) C → +∞ when the temperature becomes very small.

FIG. 7 .

 7 FIG. 7. Total electron collision frequency ν ei + ν en in the plasma state as a function of electron and lattice-ion temperatures for T e > T il which account for realistic physical conditions.

FIG. 8 .

 8 FIG. 8. Effective collision frequency ν e f f developed in this work and the generic electron collision frequency ν c used in 5 as a function of the temperature assuming T e = T il .

FIG. 9 .

 9 FIG. 9. Analytical κ th , analytical κ th, w/o F without the fragmentation processes and ab initio κ QMD thermal conductivities as a function of the temperature with T e = T il .

FIG. 11 .

 11 FIG. 11. Polystyrene eigenmodes

1 <

 1 y < 2.5 : µ ≈ k B T e ln 4y/(2 √ πy) 1 + e a(y-y 0 ) + (3y/2) 4/3π 2 /61 + e -a(y-y 0 ) (C7) y > 2.5 : µ ≈ k B T e chosen as y 0 = 1.52 and a = 15 in order to have a good approximation in the transition zone.

Coulomb logarithm

VI. ACKNOWLEDGMENTS

The authors acknowledge Vladimir Tikhonchuk for fruitful discussions about Coulomb logarithm. Patrick Mora is also gratefully acknowledged for his input on the thermal conductivity. This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of "the Investments for the future" Programme IdEx Bordeaux -LAPHIA (ANR-10-IDEX-03-02). CEA is also acknowledged for supporting these studies.

VII. DATA AVAILABILITY STATEMENTS

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A: Electron-phonon collision frequency

Electrons in the conduction band are described under the parabolic band approximation and their energies are given by has:

Appendix B: Polystyrene dispersion relation

The polystyrene consists of a chain of the unit cell CH(Ph) -CH 2 where Ph refers to the phenyl group C 6 H 5 (see Fig. 1). Even if it is not always possible to define a periodicity from a macroscopic point view, it can be defined a local periodicity of the unit cell CH(PH) -CH 2 in order to get the phonon relation dispersion. However, the unit cell is here composed of 16 atoms leading to 48 eigenmodes in three dimensions. In order to get a problem easier to handle, only phonons propagating along the carbon chain, which is assumed linear, are considered in a first approximation. This allows us to consider a one dimension problem. In addition, the unit cell is also simplified by assuming it is composed of the phenyl group with a mass m Ph = 6m C + 5m H , the group CH with a mass m CH = m C + m H and the group CH 2 with a mass m CH 2 = m C + 2m H . This simplified form of the unit cell is given by Fig. 10. Within this simplification, one acoustic mode and two optical modes are expected. Finally, it is assumed only the closest neighbors interact among themselves and the interaction potential is taken in the harmonic approximation, i.e. U(x) = kx 2 /2 with k the interaction energy determined experimentally. This leads to :

The constants G and K characterize the CH -Ph and CH -CH 2 bonds respectively. They are defined in terms of the sound speed c s with G/K = 1.21375 58 . The constant a is the CH -CH 2 bond length, and quantities u X (x, y) represent the displacement of the group X having the initial position x and y.

The motion equations are given by the Euler-Lagrange equations :

By seeking for solutions in the form u X (x, y) = U X e i(k x x+k y y-ωt) with U X the amplitude, this leads to a linear system that has a solution if and only if :

= 0 (B8) It appears the term e ik y d 1 vanishes, because only longitudinal modes are considered. The nullity of the determinant leads to a 6th-order polynomial equation:

where the coefficients α, β , γ and δ are given by:

Finally, the eigenmodes of the polystyrene are given by :

Figure 11 plots these three modes. It is observed one acoustic mode and two optical modes as expected. The acoustic and optical modes are characterized in the limit of q → 0 by a linear and constant behavior respectively. It is found that phonon energies lie in the range of other dielectrics, as SiO 2 59,60 .


Pascal Aubertin

Julien Réthoré

René De Borst
email: r.d.borst@tue.nl

A coupled molecular dynamics and extended finite element method for dynamic crack propagation

Keywords: multiscale methods, molecular dynamics, extended finite element method, fracture, crack propagation

A multiscale method is presented which couples a molecular dynamics approach for describing fracture at the crack tip with an extended finite element method for discretizing the remainder of the domain. After recalling the basic equations of molecular dynamics and continuum mechanics, the discretization is discussed for the continuum subdomain where the partition-of-unity property of finite element shape functions is used, since in this fashion the crack in the wake of its tip is naturally modelled as a traction-free discontinuity. Next, the zonal coupling method between the atomistic and continuum models is recapitulated. Finally, it is discussed how the stress has been computed in the atomic subdomain, and a two-dimensional computation is presented of dynamic fracture using the coupled model. The result shows multiple branching, which is reminiscent of recent results from simulations on dynamic fracture using cohesive-zone models.

INTRODUCTION

Modern research into fracture commences with the seminal work of Griffith [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF]. Later, Irwin [START_REF] Irwin | Analysis of stresses and strains near the end of a crack traversing a plate[END_REF] and Rice [START_REF] Rice | A path independent integral and the approximate analysis of strain concentration by notches and cracks[END_REF] established the relation between the stress intensity factors and the energy release rate, and gave linear elastic fracture mechanics a firm basis. However, linear elastic fracture mechanics only applies to crack-like flaws in an otherwise linear elastic solid and when the singularity associated with that flaw is characterized by a non-vanishing energy release rate. The fracture and any dissipative processes must also be confined to a small region in the vicinity of the crack tip.

Linear elastic fracture mechanics provides a challenge to standard finite element approaches, since the polynomials that are conventionally applied in finite element methods cannot easily capture the stress singularity at the crack tip which is predicted in linear elastic fracture mechanics. However, methods have been developed to overcome this difficulty, e.g. the so-called quarterpoint elements [START_REF] Henshell | Crack tip finite elements are unnecessary[END_REF][START_REF] Barsoum | On the use of isoparametric finite elements in linear fracture mechanics[END_REF], and more recently, the advent of meshless methods [START_REF] Fleming | Enriched element-free Galerkin methods for crack tip fields[END_REF][START_REF] Krysl | The element-free Galerkin Method for dynamic propagation of arbitrary 3-D cracks[END_REF] and partition-ofunity based finite element methods [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Réthoré | A two-scale approach for fluid flow in fractured porous media[END_REF] have provided elegant solutions to incorporate stress singularities in domain-based discretization methods. On the other hand, boundary integral methods can naturally incorporate such singularities [START_REF] Albuquerque | Dual boundary element method for anisotropic dynamic fracture mechanics[END_REF].

When the region in which the separation and dissipative process take place is not small compared to a structural dimension, but any non-linearity is confined to a surface emanating from a classical crack tip, i.e. one with a non-vanishing energy release rate, cohesive zone models as introduced by Barenblatt [12] and Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] apply. The cohesive zone approach was extended by Hillerborg et al. [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF] and Needleman [START_REF] Needleman | A continuum model for void nucleation by inclusion debonding[END_REF] to circumstances where: (i) an initial crack-like flaw need not be present or, if one is present, it need not be associated with a non-vanishing energy release rate; and (ii) non-linear deformation behaviour may occur over an extended volume. Initially, cohesive-zone models were incorporated in finite element methods via special-purpose interface elements [START_REF] Rots | Smeared and discrete representations of localized fracture[END_REF][START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF], but more recently, partition-of-unity finite element methods have shown to be very amenable to the incorporation of cohesive-zone models, e.g. [START_REF] Wells | Discontinuous analysis of softening solids under impact loading[END_REF][START_REF] Wells | A consistent geometrically non-linear approach for delamination[END_REF][START_REF] De Borst | Numerical aspects of cohesive zone models[END_REF]. In particular, it naturally enables crack propagation, also in dynamics [START_REF] Duarte | A generalized finite element method for the simulation of three-dimensional crack propagation[END_REF][START_REF] Réthoré | An energy conserving scheme for dynamic crack growth with the extended finite element method[END_REF][START_REF] Menouillard | Efficient explicit time stepping for the extended finite element method[END_REF][START_REF] Remmers | A cohesive segments method for the simulation of crack growth[END_REF][START_REF] Remmers | The simulation of dynamic crack propagation using the cohesive segments method[END_REF] and in multi-phase continua [START_REF] Réthoré | A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks[END_REF].

In spite of the power of the cohesive-zone approach, and its wide applicability on a range of scales, it remains a phenomenological approach. Probably, quantum mechanics is physically the most appropriate theory to describe fracture, but the difficulties to relate quantum mechanics to continuum mechanics, e.g. via Density Functional Theory [START_REF] Parr | Local density functional theory of atoms and molecules[END_REF][START_REF] Car | Unified approach for molecular dynamics and density-functional theory[END_REF] presently seem insurmountable. One scale of observation higher is to use Molecular Dynamics to describe fracture processes from a more fundamental physics point of view. Indeed, researchers have recently used this approach to describe fracture, e.g. [START_REF] Abraham | Simulating materials failure by using up to one billion atoms and the world's fastest computer: work-hardening[END_REF][START_REF] Zhou | Three-dimensional fracture via large-scale molecular dynamics[END_REF][START_REF] Miller | Quasicontinuum models of fracture and plasticity[END_REF]. A disadvantage of the approach is that it is computationally demanding. For this reason multi-scale approaches have been introduced, in fracture [START_REF] Kolhoff | Crack propagation in bcc crystals studied with a combined finite element and atomistic model[END_REF], as well as in plasticity [START_REF] Shilkrot | Coupled atomistic and discrete dislocation plasticity[END_REF], in which only a part of the body is analysed using molecular dynamics, while the remaining part of the body is modelled using continuum mechanics and discretized using a finite element method. This manuscript furthers along this line and combines molecular dynamics for modelling the fracture process at the crack tip with an extended finite element method (XFEM), where the partition-of-unity property of the polynomial shape functions is exploited to model the crack in the wake of the tip as a traction-free discontinuity. It is noted that recently another approach has been published that couples atomistics and extended finite elements [START_REF] Gracie | Concurrently coupled atomistic and XFEM models for dislocations and cracks[END_REF], but the current paper makes a further advancement in that it includes dynamic loadings.

A major issue in multi-scale approaches as discussed above is the accurate coupling of both domains, especially when different descriptions are assumed on either domain. While the coupling can, in principle, either be achieved at a discrete interface, or on a zone of a finite size (overlap or zonal coupling), it is believed that zonal approaches, which include the Arlequin method [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Prudhomme | Computational analysis of modeling error for the coupling of particle and continuum models by the Arlequin method[END_REF], the bridging domain method [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF][START_REF] Zhang | A bridging domain and strain computation method for coupled atomistic-continuum modelling of solids[END_REF][START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Xu | Conservation properties of the bridging domain method fod coupled molecular/continuum dynamics[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF], discrete-to-continuum bridging [START_REF] Fish | Discrete-to-continuum bridging based on multigrid principles[END_REF], the discontinuous enrichment method [START_REF] Farhat | The discontinuous enrichment method for multiscale analysis[END_REF], and bridging scale decomposition [START_REF] Wagner | Coupling of atomistic and continuum simulations using a bridging scale decomposition[END_REF][START_REF] Farrell | Implementation aspects of the bridging scale method and application to intersonic crack propagation[END_REF] enable a more gradual transition from one domain to the other. The ability of a gradual transition is especially important for highly dissimilar domains and when wave propagation phenomena are considered, where preservation of the energy and avoiding spurious reflections when a wave exits one domain and enters the other can become an issue. Inspired by earlier work by Ben Dhia and Rateau [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF] and Xiao and Belytschko [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF] we have chosen a weak coupling between the models in the two adjacent domains.

This paper is organized as follows. First, we briefly list the equations of molecular dynamics and ways to ensure equilibrium of the atomistic domain before starting the computation that involves dynamic propagation of an existing, starter crack. This is followed by a succinct recapitulating of the governing equations of continuum mechanics, both in the strong and the weak forms. The discretization of the continuum subdomain is carried out using the extended finite element method, where the partition-of-unity property of finite element shape functions is used to model the traction-free discontinuity in the wake of the crack tip. Subsequently, it is recapitulated how both domains can be coupled, see also [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF] which also presents a detailed analysis of the energy conservation properties of the coupling scheme. The paper concludes with a full two-dimensional coupled analysis of dynamic crack propagation which shows multiple branching and suggests the formation of dislocations, which shows similarities with recent simulations on dynamic fracture using cohesive-zone models [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Remmers | The simulation of dynamic crack propagation using the cohesive segments method[END_REF].

MOLECULAR DYNAMICS

Governing equations in the atomic domain

For the discrete domain, i.e. m , we build a grid of N a atoms, and, accordingly, the initial value problem in this domain can be written as:

For 1 i N a (t) and t ∈[0; T ], given the initial conditions (d(0), ḋ(0)) find (d, f) ∈ D ad ×F ad such that:

m i di = f i (1)
with m i the mass of atom i and:

D ad ={d = (d i (t)) 1 i N a ∀t ∈[0, T ]} F ad ={f = (f i (t) =-∇ i U(d(t))) 1 i N a ∀t ∈[0, T ]} (2)
from where it transpires that the interatomic forces are derived from a potential energy U. d and f assemble the discrete displacements d i and forces f i of the individual atoms, respectively. The internal energy of the discrete domain can be viewed as the sum of each atomic contribution U j :

U = j U j (d) (3)
and the force f i acting on atom i can be written as the sum of all elementary forces:

f i =- *U *d i = j =i f ij (4)
In order to limit the cost of computing such a force, we reduce the summation by only including so-called 'nearest' neighbours, within a cut-off radius r c :

f i ≃ r ij <r c f ij (5)
where r ij is the interatomic distance. For details the reader is referred to the Appendix.

In the subdomain m the weak formulation becomes:

∀w * ∈ Ḋad,0 , given the initial conditions (d(0), ḋ(0)) find d ∈ D ad such that: a m (d, w * ) = 0 (6)
with w * the test function, and

a m (d, w * ) = N a i=1 m i di •w * i + N a i=1 ∇ i U(d)•w * i 2.

Equilibrium of the atomic domain

The atomistic domain has to be embedded in a continuum and we cannot consider periodic conditions. Thus, the zone of atoms has an internal boundary with the continuum zone, which has no physical character. Atoms that are on the boundary between the coupling domain and the continuum domain do not have neighbouring atoms in the continuum. This can produce nonphysical forces that break the equilibrium. In order to obtain a global stable system that is in equilibrium before loading the specimen, we therefore have to introduce so-called 'ghost' atoms in the continuum, which act as neighbours for the atoms on the boundaries. Their positions are directly derived from the continuum domain, and their velocities and accelerations do not have to be computed. Their creation merely aims to satisfy the equilibrium of the atomic domain.

For crack propagation we introduce a pre-crack or starter crack in the atomic domain. For this purpose we remove the atoms along a line that defines the position of starter crack. This initial crack crosses the continuum as well as the atomistic zone. The removal of atoms causes that the initial lattice is no longer at equilibrium and there is an additional boundary to consider. However, this boundary has a physical character because it stands for a crack surface.

As alluded to above, equilibrium in the atomistic zone has to be found prior to applying any load to the specimen. A widely used technique is the use of dynamic relaxation that stabilizes a given configuration with adaptative damping. As shown in [START_REF] Pan | The Meshless Dynamic Relaxation Technique for Simulating Atomic Structures of Materials[END_REF] the use of dynamic relaxation in molecular simulations can be effective in obtaining an equilibrium state. The solution of the problem is viewed as the steady-state solution of a damped wave equation. The classical equations of motion are augmented by a damping term to give:

Mẍ+Cẋ-∇ x U(x) = 0 (7)
When the steady-state solution is found, the acceleration and the velocity are zero. Therefore, the matrices M and C have no physical meaning and can be chosen such that the computation is accelerated in an optimal sense. At each time step of the resolution, the local atomistic stiffness is computed, using the potential energy function, and a critical time step is calculated. Then, either the mass or the damping can be adapted to obtain an optimal convergence. This method has been used for finding equilibrium in the atomistic zone and gives converged results within an acceptable computational time.

CONTINUUM MODEL

Assuming small deformation gradients for simplicity, the governing equations in the continuum subdomain M can be written in a standard manner as:

For x ∈ M (t) and t ∈[0; T ], given the initial conditions (u(x, 0), u(x, 0)) find (u, r) ∈ U ad ×S ad such that: ü = div r+g d [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF] with the mass density and u the continuum displacement vector, r the stress tensor, and g d the body force vector applied in M , subject to the boundary conditions

U ad ={u = u(x, t) ∈[H 1 ( M )] 3 ; u = u d on * 1 ∀t ∈[0, T ]} S ad ={r = K :∇u(x, t) ∈[L 2 ( M )] 6 ; r•n = F d on * 2 ∀t ∈[0, T ]} ( 9
)
where K is the fourth-order stiffness tensor, n the outward normal vector to * 2 ,a n du d and F d the prescribed displacements and tractions at * 1 and * 2 , respectively.

To allow for a discretization of the continuum subdomain we next specify the weak formulation:

∀v * ∈ Uad,0 , given the initial conditions (u(x, 0), u(x, 0)) find u ∈ U ad such that:

a M (u, v * ) =l M (v * ) (10)
with v * the test function, and

a M (u, v * ) = M ü•v * d + M e(u) : K : e(v * ) d (11) l M (v * ) = * 2 F d •v * dS+ M g d •v * d (12)
4. COUPLING METHOD

Coupling functions

In order to enable an efficient coupling between the two domains, a coupling zone is defined and a coupling function is computed which is used to obtain the global energy. First a coupling length L c is chosen such that will be the characteristic length of the coupling region. Subsequently, the patch of atoms, m , is included in the continuum at a given position. Atoms at a distance ℓ L c from the Molecular Dynamics Box (MD-Box) boundary are considered to be in the coupling where only the atomistic model applies, elements are removed. The resulting, discretized domain is shown in Figure 1. In principle, each atom must be coupled to the finite element discretization of the underlying continuum. It suffices to couple only a limited number of atoms to the underlying continuum without loss of accuracy. Since this directly affects the size of the coupling matrices that will be derived next, such a limited coupling markedly decreases the computational effort, and therefore, a bigger scale transition becomes possible at reasonable computational costs.

In the coupling zone c , we enforce a velocity coupling in a weak sense, and we distribute the energy between both models via a partition of unity [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF]. For this purpose we define the following functions, see Figure 2:

: M →[0, 1] (13)
: m →[0, 1] (14)
such that:

(x) = 1 for x ∈ M \ c (x) = 1 for x ∈ m \ c (x)+ (x) = 1 for x ∈ c ( 15
)
The displacement and velocity fields in the domains M and m have a different nature. In M we have a continuum field, while in m we have a discrete field, which is only defined at the geometrical points corresponding to the atoms. Therefore we construct a so-called 'mediator space', denoted by M, on which we project the fields u and ḋ in order to be able to compare them. The nature of M is constrained by the discrete character of the atomistic field. Indeed, its displacements cannot be extrapolated outside the atomic positions if we want to maintain a physical interpretation at the fine scale. Accordingly, M has to be a subspace of the physical atomistic space. More precisely, we project the velocities using an operator on a discrete subset c of the atomic positions included in c . Considering that M is built as a Hilbert space, we introduce a scalar product c from M×M onto R. With these definitions we formulate the velocity coupling operator as:

∀l * ∈ M, c(l * , u-ḋ) = * , u-ḋ M (16)
with c the classical scalar product on M. The global equations are coupled via Lagrange multipliers and can subsequently be written as:

∀(v * , w * , l * ) ∈ Uad,0 × Ḋad,0 ×M, given the initial conditions (u(x, 0), u(x, 0), d(0), ḋ(0))

find (u, d, k) ∈ U ad ×D ad ×M such that:

a ,M (u, v * )+a ,m (d, w * )+c(k, v * -w * )+c(l * , u-ḋ) =l ,M (v * ) (17)
The modified forms a ,M , a ,m and l ,M take into account the weighting functions (x) and (x),see [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF] for details. In the atomic subdomain the modified form a ,m that takes into account the distribution of the energy reads:

a ,m (d, w * ) = w * •m d-w * •f (18)
with i = (d i ), m =[ i i, j m i ], f =[f ,i ] (19)

Discretized problem

In a manner which is by now standard the interpolation of each component of the displacement field is enriched with discontinuous functions in order to properly capture the traction-free discontinuity in the wake of the crack tip:

∀x ∈ M , u h (x) = i∈N M N i (x) ūi + i∈N cut N i (x)H d ûi (20)
where N i are standard finite element shape functions supported by the set of nodes N M included in the discretized domain M . Nodes in N cut have their support cut by the discontinuity. They hold additional degrees of freedom ûi corresponding to the discontinuous function H d defined by:

H d (x) = x•n d x•n d ( 21
)
with n d the normal to the discontinuity . Symbolically, Equation ( 20) can be written as

∀x ∈ M , u h = N T U ( 22
)
where the matrix N contains the standard interpolation polynomials N i (x) as well as the discontinuous function H d , and the array U contains the displacement degrees-of-freedom ūi and ûi .The transition within the domain M between the subdomain where the nodes are 'enriched' and the part which has just the standard formulation does not affect the Molecular Dynamics computation other than through the coupling matrices.

With the latter symbolic notation the bilinear form a ,M and the linear form l ,M become:

a ,M (u h , v * h ) = V * T M Ü+V * T K U (23) l ,M (v * h ) = V * T F (24)
where the term that represents the body forces has been omitted for simplicity, and

M = M N T N d (25)
K = M ∇N T K∇N d ( 26
)
the mass and stiffness matrices, respectively. With the standard definition of the scalar product, the coupling term in the continuum can be discretized as follows:

c(k, v * h ) = V * T C M K = V * T F L M ( 27
)
with C M the continuum coupling matrix. The vector K contains the Lagrange multipliers and its size equals the c subset cardinal times the dimension of the space considered. F L M can be regarded as a fictitious force due to the coupling via the Lagrange multipliers. This force has a non-zero value only in the coupling zone c .

Using the Lagrange multipliers in the atomic domain similar to that in the continuum domain:

f L m = C m K ( 28
)
the weighted and coupled system ( 17) can be cast in a matrix-vector format:

V * T (M Ü+K U+C M K)+W * T (m d-f -C m K)+l * T (C T M U-C T m ḋ) = V * T F (29)
Since this set must hold for any admissible (V * , W * , l * ) we finally obtain:

M Ü+K U = F -C M K m d = f +C m K C T M U = C T m ḋ (30)
with (U, d, K) the set of unknowns. Details on the time integration scheme associated with this set of coupled ordinary differential equations are given in Reference [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF].

DYNAMIC FRACTURE

Mechanical quantities in the atomic domain

In order to extract mechanical quantities from the atomistic domain, we adopt a continuum mechanics point of view to derive classical stress quantities. The atomic stress tensor at an atom i is a measurement of the interatomic interactions of the atom with its neighbours. A widely used stress quantity defined on the atomistic domain is the virial stress, which takes into account the interactions and a kinetic energy contribution. Many formulations have been derived from this virial stress [START_REF] Basinski | Influence of shear stress on screw dislocations in a model sodium lattice[END_REF][START_REF] Lutsko | Stress and elastic constants in anisotropic solids: molecular dynamics techniques[END_REF][START_REF] Cheung | Atomic-level stress in an inhomogeneous system[END_REF], but, as pointed out by Zhou [START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF], these definitions, even perfectly correct in a statistical and thermodynamical sense, do not correspond to the Cauchy stress or any other mechanical stress. However, it can be shown that the interatomic interactions part of the virial stress reduces to the Cauchy stress with a physical meaning. We therefore adopt this definition for the stress tensor:

r i = 1 2V i r ij <r c f ji ⊗r ij ( 31
)
where V i is the volume of the atom i. Subsequently, the average of this atomistic stress tensor is computed over the volume around i within the cut-off radius. The average atomistic stress at the atom i thus reads:

r avg i = 1 N n,i N n,i j=1 r j (32)
Finally, the Von Mises atomic stress c at an atom i is defined as follows:

2( c,i ) 2 = ( avg xx - avg yy ) 2 +( avg yy - avg zz ) 2 +( avg xx - avg zz ) 2 +6( avg 2 xy + avg 2 yz + avg 2 xz ) (33)

Two-dimensional simulation

We now proceed with a two-dimensional simulation. A copper single crystal is considered in its (111) plane, so that the two-dimensional lattice is hexagonal. The Lennard-Jones potential is used in the molecular dynamics simulation with parameters from [START_REF] Agrawal | Predicting trends in rate parameters for self-diffusion on FCC metal surfaces[END_REF]: a = 0.415 eV and b = 0.2277 nm. The Young's modulus E and Poisson's ratio for the continuum are obtained following the procedure given in [START_REF] Aubertin | Coupling of atomistic and continuum models: dynamic crack propagation[END_REF]: E = 79.334 GPa and = 0.25. The copper atomic mass is taken as m = 0.105520602596×10 -24 kg (m Cu = (63.546/N A )g with N A the Avogadro number: N A = 6.02214179×10 23 mol -1 ), which corresponds to a mass density = 1865.250812586×10 3 kg m -3 .

In the present study, the temperature has not been taken into account, since the focus is on the coupling of Molecular Dynamics to an (extended) finite element method for crack propagation. When extending the methodology to explicitly include the temperature a 'thermal equilibrium' has to be achieved in addition to the mechanical equilibrium. This can for instance be done using the Nose-Hoover thermostat method [START_REF] Nose | Constant-temperature molecular dynamics[END_REF]. The computational domain has been plotted in Figure 3. It is 100 nm long and 77.5nm wide with an initial crack of 10 nm. A large MD-Box is considered in order to properly trace the crack propagation during dynamic loading. The finite element mesh consists of 1221 quadrilateral elements and 4868 nodes. The element size is about 10 times the interatomic distance. 36635 atoms are put in the vicinity of the crack tip. The width of the coupling domain is approximatively 3 nm and 33% of the atoms in this region hold Lagrange multipliers. For computational reasons, the results that are presented from now on have been obtained by only including the first neighbours in the atomistic interactions. Indeed, the equilibration techniques and the updates are expensive when we take into account many neighbours. Simulations on a smaller scale have indicated that the results are close to those obtained with more neighbours.

The test consists of applying a velocity on the top and bottom edges of the specimen. The prescribed velocity in this example is V p = 47.4ms -1 . The time step is 15.811388 fs (10 -15 s). Figure 4 gives the displacement field after the crack has propagated, at t = 33.209 ps. Close inspection reveals that, going from the atomistic zone to the continuum domain, no spurious wave reflections occur, which indicates that the coupling algorithm works properly, and is in agreement with results for linear elasticity presented in [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF]. We also see that, even though the prescribed loading corresponds to Mode-I and the atomistic lattice is a perfect single crystal, the crack path is not that which we would expect to obtain with classical continuum methods. We observe crack branching, the occurrence of dislocations, and locally, mixed-mode behaviour.

Figures 5 and6 show the evolution of the displacements and the Von Mises stress in the MD-Box during propagation. We recall that the stress representation is used to represent the behaviour of the interatomic forces qualitatively, and should not be considered as a quantitative result. The concept of stress at the atomistic level is discutable, and, moreover, as we simulate a two-dimensional lattice, the computed values are not representative of a real three-dimensional copper crystal. Yet, when considering Figure 6, the global behaviour of the stress field is consistent with classical results of dynamic crack propagation: A stress concentration forms, which is initially confined at the crack tip, and subsequently propagation occurs and stress waves are emitted from the tips. Moreover, the crack branching that is observed has been demonstrated also in experiments on fast crack propagation [START_REF] Coker | Dynamic crack growth along a polymer composite-homalite interface[END_REF], and in recent cohesive-zone simulations [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Remmers | The simulation of dynamic crack propagation using the cohesive segments method[END_REF].

CONCLUDING REMARKS

A numerical approach has been proposed for combining a molecular dynamics method and a finite element method that exploits the partition-of-unity property of finite element shape functions (extended finite element method). The aim is to simulate dynamic fracture in an efficient manner on basis of elementary physical principles. To this end the zone around the crack tip is modelled using molecular dynamics. Around this so-called Molecular Dynamics Box a continuum mechanics approach is adopted, with the finite element method used for discretization. The partition-of-unity property of the finite element shape functions is exploited to model the crack in the wake of its tip as a traction-free discontinuity. The coupling between the continuum and molecular dynamics zones has a zonal character where the energy is partitioned over both models and a weak velocity coupling is enforced. In this manner, spurious reflections are avoided and energy is conserved when a wave travels from one zone into another [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF].

Two-dimensional simulations show realistic crack patterns, including branching and local mixedmode behaviour. Also the stress concentrations at the crack tips, and the stress wave propagation show qualitatively realistic patterns. Quantitatively, less can be said as the concept of stress as it has been defined can be questioned at an atomistic level and as three-dimensional effects cannot be represented well in this two-dimensional computation. Furthermore, the computation is limited by the size of the MD-Box. When a crack tip reaches the boundary of the atomistic domain and intersects the coupling zone, the cracks are arrested and cannot propagate further, cf. The neighbours of each atom i are identified and stored in a table neighb i . This table contains N n,i pointers to the neighbouring atoms. A basic method would be to loop over all the atoms to identify a possible neighbour. This would result in a computational time that is O(N a 2 ) for N a atoms. To decrease the computational cost we therefore have implemented a cell method. A grid of regular cells is built in the wake of the atomistic domain, such that each atom i belongs to a cell C I . The cells dimensions are set such that all the neighbours of an atom are included in the neighbours cells of C I , Figure A1.

Such a method improves the computational efficiency regarding the neighbours list, but we have to determine the optimal number of cells. Indeed, if the grid is either too coarse, or too fine, the cell method becomes inefficient. In Figure A2, we have 7500 atoms in the lattice and we have plotted the CPU time needed for the list. In this case the optimal choice is a grid of ∼ 180 cells. Subsequently, the optimal number of cells has been determined for different lattices. This results in Figure A3 (left), in which the optimal number of cells has been plotted as a function of the number of atoms in the lattice. A logarithmic dependence can be identified that will be used in the sequel to build the grid in the computations. The CPU time spent for the computation of the neighbours list now shows a clear improvement, see Figure A3 (right). While the computation time of the basic method is O(N a 2 ), that of the cell method is O(N a √ N a ). For each atom we have built a table neighb i that contains its neighbours. During a computation of dynamic fracture, the neighbourhood of an atom, especially in the vicinity of the crack, will change, and a criterion has to be implemented to update the neighbours list. For this purpose, a Verlet list has been implemented. For each atom, more pointers are stored than just the neighbours. In addition to storing the atoms within the cut-off radius r c , all atoms within a larger radius, r v are stored. However, the interatomic force is computed only for the neighbours within an inner circle, that is for r <r c . Accordingly, if an atom which is originally located such that r c <r <r v , but which for which r <r c at a generic point in the computation, it is considered as an effective neighbour from that moment onwards, and is included in the computation of the interatomic forces. The criterion used to update the neighbours list is a maximum displacement criterion. At each time step the two maximum atomic displacements d 1,max and d 2,max are considered. If d 1,max +d 2,max r v -r c ,t h e list is updated, as, in this case, an atom can have displaced more than the distance allowed in the neighbours tables.

Figure 1 .

 1 Figure 1. Discretized domain with the coupling region. The dark, 'pear-shaped' area is the domain where an MD calculation is carried out ( m ). The coupling region c consists of elements that are surrounded by a bold line.

Figure 2 .

 2 Figure 2. Partition of unity for the energy distribution.

Figure 3 .

 3 Figure 3. Initial configuration.

Figure 4 .

 4 Figure 4. Deformed configuration at t = 33.209 ps.

Figure 5 .

 5 Figure 5. MD-Box: Displacement field around the crack tip during propagation: (a) displacement field at t = 11.070 ps; (b) displacement field at t = 17.395 ps; (c) displacement field at t = 23.721 ps; (d) displacement field at t = 30.046 ps; (e) displacement field at t = 36.372 ps; and (f) displacement field at t = 39.535 ps.

Figure 6 .

 6 Figure 6. MD-Box: Von Mises stress field around the crack tip during propagation: (a) Von Mises stress at t = 11.070 ps; (b) Von Mises stress at t = 17.395 ps; (c) Von Mises stress at t = 23.721 ps; (d) Von Mises stress at t = 30.046 ps; (e) Von Mises stress at t = 36.372 ps; and (f) Von Mises stress at t = 39.535 ps.

Figure A1 .

 A1 Figure A1. Atomic lattice and cells: the neighbourhood of an atom is included in the neighbours cells of C I .

Figure A2 .

 A2 Figure A2. CPU time needed for neighbours list calculation as a function of the number of cells in the grid. 7500 atoms.

Figure

  Figure A3. Left: optimal number of cells versus number of atoms in the lattice; Right: CPU time for neighbours list computation versus N a , number of atoms. The dashed line is the basic method and bold line represents the cell method.