Fabien Brieuc

Yael Bronstein

Hichem Dammak
email: hichem.dammak@centralesupelec.fr

Philippe Depondt

Fabio Finocchi
email: fabio.finocchi@insp.upmc.fr

Marc Hayoun

Zero-point energy leakage in Quantum Thermal Bath molecular dynamics simulations

The quantum thermal bath (QTB) has been presented as an alternative to pathintegral based methods to introduce nuclear quantum effects in molecular dynamics simulations. The method has proved to be efficient, yielding accurate results for various systems. However, the QTB method is prone to zero-point energy leakage (ZPEL) in highly anharmonic systems. This is a well known problem in methods based on classical trajectories where part of the energy of the high frequency modes is transferred to the low frequency modes leading to a wrong energy distribution. In some cases, the ZPEL can have dramatic consequences on the properties of the system. Thus, we investigate the ZPEL by testing the QTB method on selected systems with increasing complexity 1 in order to study the conditions and the parameters that influence the leakage. We also analyze the consequences of the ZPEL on the structural and vibrational properties of the system. We find that the leakage is particularly dependent on the damping coefficient and that increasing its value can reduce, and in some cases, completely remove the ZPEL. When using sufficiently high values for the damping coefficient, the expected energy distribution among the vibrational modes is ensured. In this case, the QTB method gives very encouraging results. In particular, the structural properties are well reproduced. The dynamical properties should be regarded with caution although valuable information can still be extracted from the vibrational spectrum, even for large values of the damping term.

Introduction

The molecular dynamics (MD) simulation technique is a powerful tool to investigate the properties of complex atomic systems. At low temperature and/or in systems containing light elements such as hydrogen, nuclear quantum effects can play a major role on the behavior of the system. However, it is currently a computational challenge to account for the quantum nature of nuclei in MD simulations.

Over the past years, several techniques have been proposed to deal with this issue. Among them, the most common methods are based on the path integral formalism such as path integral molecular dynamics (PIMD). In this formalism, each quantum nucleus is described by a ring of classical monomers (or beads) connected through harmonic springs [1][2][START_REF] Tuckerman | Statistical mechanics: Theory and Molecular Simulation[END_REF] . When the number of beads is large enough, the statistical averages converge towards the exact quantum result. However, in order to compute time correlation functions, approximate methods such as centroid MD [START_REF] Cao | [END_REF] or ring-polymer MD 5,6 are needed. These path-integral methods are computationally demanding when the number of beads increases, in particular at low temperature.

Recently, alternative methods based on a modified Langevin equation have been pro-posed 7,8 . Among them, the quantum thermal bath (QTB) 7 is an approximate yet efficient method to include nuclear quantum effects in MD simulations. Although exact only in the case of a system of harmonic oscillators, the QTB provides satisfactory results in many anharmonic systems. 9-14 A first advantage is its implementation without any additional computational cost compared to standard MD. Hence, large and complex systems can be in principle treated by QTB-MD. Moreover, the method can give information about dynamical properties of the system. Finally, its formulation is not system-dependent, in particular, no knowledge of the system's vibrational density of states is needed beforehand. However, the QTB method has several drawbacks. First, the method can fail when dealing with highly anharmonic systems. 15 Second, the QTB technique is prone to zero-point energy leakage (ZPEL), like any other method based on classical trajectories. 16 The ZPEL is a known problem where a part of the energy of the high-frequency modes is transferred to the low-frequency ones which is due to the classical nature of MD trajectories.

The ZPEL was observed in many different systems (water clusters and liquid water, Lennard-Jones systems, ...) [17][18][19][20][21] ; in particular, in the case of the QTB, the ZPEL has been recently pointed out by Bedoya-Martinez and coworkers. 20 However, no systematic or general study of ZPEL within the QTB framework has been done up to now.

Several solutions to the ZPEL problem within QTB-MD simulations have recently been suggested. Bedoya-Martinez and coworkers tried to modify the noise power spectrum in order to obtain the expected energy distribution. However, this solution is system-dependent and only worked for weakly anharmonic systems. 20 Ganeshan and coworkers proposed a deterministic approach to suppress ZPEL, which unfortunately requires the knowledge of the vibration normal coordinates prior to the simulation. 19 Here, we investigate the conditions leading to the ZPEL within QTB-MD simulations in various systems in order to get a better understanding of the validity of the QTB method.

More precisely, we focus on the conditions and the parameters that influence the ZPEL and on the consequences for the system's properties. After a brief presentation of the QTB method, we study selected anharmonic systems with increasing complexity. First, we investigate two simple models: two coupled harmonic oscillators and a one-dimensional chain of atoms.

Then, we focus on more realistic systems, a Lennard-Jones aluminium crystal and the phase transitions in BaTiO 3 . In the last section, we discuss our results and their implications.

The Quantum Thermal Bath method

The QTB method is based on a modification of the Langevin thermostat in order to include nuclear quantum effects in MD simulations. Both in the standard (i.e. classical) Langevin thermostat and in the QTB method, the equation of motion for one degree of freedom x of mass m and submitted to the internal force f (x) reads:

7 mẍ = f (x) -mγ ẋ + R(t) (1)
The last two terms correspond to the friction and stochastic forces of the thermostat, respectively.

The random force is described by a stationary stochastic process R(t) whose distribution is Gaussian with zero mean:

< R(t) > = 0 (2) < R(t)R(t + τ ) > = +∞ -∞ I R (ω, T ) e -iωτ dω 2π . ( 3
)
Equation 3 is the Wiener-Khinchin theorem, which relates the autocorrelation function < R(t)R(t + τ ) > of the stochastic process to its power spectral density (PSD) I R (ω, T ) at temperature T . The dynamical properties obtained using eq. 1 are directly related to this PSD. The closely related PSD of the position, I x , is obtained from the fluctuation-dissipation theorem 22 , which reads in the classical case:

χ′′ (ω) = ω 2k B T I x (ω, T ) (4)
with χ′′ (ω) the imaginary part of the susceptibility χ(ω) that connects the Fourier transform of the position x(ω) to the Fourier transform of the random force R(ω) within the linear response theory:

x(ω) = χ(ω) R(ω) (5)
From this expression, we obtain a linear relation between the PSD of the position, I x , and the PSD of the stochastic force, I R :

I x (ω, T ) = | χ(ω)| 2 I R (ω, T ) (6)
and the fluctuation-dissipation theorem can be rewritten as follows:

I R (ω, T ) = 2k B T ω χ′′ (ω) | χ(ω)| 2 . ( 7
)
In the case of an harmonic oscillator with an angular frequency ω 0 , using eq. 1 and 5 the susceptibility writes

χ(ω) = 1 m [ω 2 0 -ω 2 + iγω] . (8)
By introducing this expression in eq. 7, the PSD of the random force is obtained in the classical case as a white noise:

I R (ω, T ) = 2mγ k B T ∀ω. (9)
By using this expression and the Wiener-Khinchin theorem (eq. 3) in eq. 1, the standard Langevin dynamics is obtained. In this case, the equipartition of the energy is ensured, and all harmonic vibrational modes have the same average energy (k B T ), which is independent of the angular frequency ω.

In the quantum case, the average energy of a vibrational mode is given by

θ(ω, T ) = ω   1 2 + 1 exp ω k B T -1   ( 10
)
in the harmonic approximation. The main idea of the QTB method is to replace the PSD of the classical random force ( 9) by the one corresponding to the energy distribution of eq.

10. This is done, in practice, by using the quantum version of the fluctuation-dissipation theorem as developed by Callen and Welton 23 and reviewed by Kubo 22 , which gives:

χ′′ (ω) = ω 2θ(ω, T ) I x (ω, T ) (11)
and, through the Wiener-Khinchin theorem, leads to the PSD of the colored noise R(t) as used in the QTB method:

I R (ω, T ) = 2mγ θ(ω, T ) (12)
In contrast to the Langevin thermostat, I R is ω dependent and the random force R(t) is obtained using the procedure 24,25 described in Appendix A.

The use of an angular frequency cut-off ω cut is necessary during the generation of these random forces 26 because the average energy of a harmonic oscillator diverges at high frequencies. Thus, the QTB method contains two free parameters : the friction coefficient γ and the angular frequency cut-off ω cut . The values of these parameters must be carefully chosen.

When using the Langevin thermostat, it is generally assumed that the friction coefficient γ has to be small enough so that the forces associated with the thermostat do not significantly perturb the dynamics of the system. 27 Moreover, as already stated by Barrat and Rodney 26 , ω cut must be chosen of the order of a few times the highest angular frequency observed in the system to prevent the inclusion of non-physical high frequency modes. Too high values for ω cut and γ could lead to the divergence of the total energy. 26 In the simulations that are presented in this work, we found that a reasonable value for ω cut is approximately 2ω max with ω max being the highest angular frequency in the system. In the following, we focus on the results of QTB-MD simulations when increasing the friction coefficient γ in eq.1 and, consistently, in the power spectrum of the stochastic force (eqs. 9 and 12).

3 Model systems

Coupled harmonic oscillators

In this section, we study the behaviour of QTB-MD on a simple model consisting of two coupled one-dimensional harmonic oscillators. Thanks to the small number of degrees of freedom, we can directly compare the QTB-MD results with the numerical solution of the time-independent Schrödinger equation here. The system is described by the Hamiltonian

H: H = 1 2 m ẋ2 1 + 1 2 mω 2 1 x 2 1 + 1 2 m ẋ2 2 + 1 2 mω 2 2 x 2 2 + C 3 (x 1 -x 2 ) 3 + C 4 (x 1 -x 2 ) 4 (13)
where x 1 and x 2 are the positions of the two oscillators, ω 1 and ω 2 are their angular frequencies, m is their mass and C 3 and C 4 are coupling constants. The Hamiltonian H can be written in a dimensionless form, H = H/ ω 1 , so that:

H = q2 1 2 + q 2 1 2 + q2 2 2 + Ω 2 q 2 2 2 + c 3 (q 1 -q 2 ) 3 + c 4 (q 1 -q 2 ) 4 (14)
where the following variables are used:

Ω = ω 2 ω 1 , ξ = mω 1 , q i = x i ξ , c 3 = C 3 ξ 3 ω 1 , c 4 = C 4 ξ 4 ω 1 , t * = ω 1 t, qi = dq i dt * . ( 15
)
q 1 and q 2 are the reduced positions of the two oscillators and Ω is the ratio of the frequencies of the two oscillators (we set ω 1 > ω 2 ). The non-linear coupling terms introduce a controllable degree of anharmonicity in the system which in turn leads to a clear illustration of the ZPEL within the QTB method and allows for the analysis of the conditions leading to this phenomenon. In particular, we study here the influence of cubic and quartic coupling terms on the energies of the oscillators ǫ 1 and ǫ 2 , that are:

ǫ 1 = q2 1 2 + q 2 1 2 , ǫ 2 = q2 2 2 + Ω 2 q 2 2 2 . ( 16
)
The QTB-MD simulations were performed with a friction coefficient γ = 4 × 10 -4 ω 1 , a cut-off frequency ω cut = 2ω 1 and a time step δt = 0.05ω -1 1 . Average values are computed using at least 30 independent trajectories that are 10 7 time steps long each. The ratio Ω is varied in the 0.05-0.8 range and the parameters c 3 and c 4 are varied in the 0-25 × 10 -4 and 0-40 × 10 -4 ranges respectively, so that we cover a large range of coupling energies (figure 1). The temperature is set to k B T = 0.03 ω 1 (e.g. T ∼ 60 K if ω 1 = 2π × 40 THz) so that the thermal energy contribution to the energies of the oscillators is negligible with respect to their zero-point energies.

The exact quantum calculation shows that the energies of the oscillators are almost independent of the anharmonic coupling intensities for the range of coupling values studied Figure 2: Zero-point energy leakage quantified through the deviation factor ζ (eq. 17) as a function of the ratio of frequencies Ω (eq. ( 15)). Top panel: cubic coupling (c 3 = 0, c 4 = 0). Lower panel: quartic coupling (c 4 = 0, c 3 = 0). here and are equal to their zero-point energies; hence, in reduced units, ǫ 1 = 0.5 and ǫ 2 = Ω/2.

Figure 1 shows the average energies obtained with the QTB method in two distinct cases: Ω = 0.5 with only a cubic coupling and Ω = 0.25 with only a quartic coupling. As expected, in the uncoupled case, i.e. c 3 = 0 and c 4 = 0, the QTB method gives the expected quantum energies for the two oscillators, corresponding to their zero-point energies. In contrast, when the coupling constants c 3 or c 4 are increased, the QTB-MD energies diverge from the exact results: part of the energy of oscillator 1 is transferred into oscillator 2, hence, ZPEL occurs.

In the following, we investigate how the ZPEL depends on the three parameters (Ω,c 3 ,c 4 ) that define the Hamiltonian.

First, we adopt the following deviation factor ζ in order to quantify the ZPEL:

ζ = ∆ǫ exact -∆ǫ QTB ∆ǫ exact = (ǫ exact 1 -ǫ exact 2 ) -(ǫ QTB 1 -ǫ QTB 2 ) ǫ exact 1 -ǫ exact 2 (17)
With this definition, the leakage is maximum when ζ = 1, i.e. when the system has reached an equipartition of the energy:

ǫ QTB 1 = ǫ QTB 2 .
In contrast, there is no leakage when ζ = 0, i.e. when ǫ QTB 1,2 = ǫ exact 1,2 . In figure 2, the results obtained for ζ as a function of Ω for different values of c 3 and c 4 are presented. One can note that the ZPEL strongly depends on the ratio of frequencies and is present only for certain values of Ω. In the cubic case, it occurs only near Ω = 0.5 (figure 2.a). Indeed, cubic terms in the potential are known to be responsible for frequency doubling, that is the second harmonic generation (2ω). This is confirmed by the vibrational spectrum of the two oscillators computed from QTB-MD in the cubic case (figure 3.a): harmonics at 2ω 2 , ω 1 -ω 2 , and ω 1 + ω 2 are visible. Therefore, at Ω = 0.5, there is a resonance between the couple of modes (ω 1 ; 2ω 2 ) and (ω 2 ; ω 1 -ω 2 ). Similarly, the quartic terms are responsible for the generation of modes with frequency 3ω; ZPEL is indeed observed near the resonance at Ω = 1/3 (figure 2.b). With increasing quartic coupling, significant ZPEL also occurs for smaller values of Ω.

0 0.5 1 1.5 2 arbitrary units γ=4x10 -4 ω 1 ω 1 -ω 2 2ω 2 ω 1 +ω 2 (a) 0 0.5 1 1.5 2 γ=4x10 -3 ω 1 (c) 0 0.5 1 1.5 2 arbitrary units ω/ω 1 γ=4x10 -4 ω 1 (b) 0 0.5 1 1.5 2 ω/ω 1 γ=4x10 -3 ω 1 (d)

Influence of the friction coefficient γ

The damping coefficient is now varied from 4 × 10 -4 ω 1 to 2 × 10 -2 ω 1 . We focus on the frequency range where the ZPEL is important: Ω = 0.5 for cubic coupling and Ω = 0.25 for quartic coupling (see figure 2). Figure 4 shows that the ZPEL strongly depends on γ. Increasing γ can limit the leakage and even practically remove it in the case of the cubic coupling. In particular, for c 3 = 2.4 × 10 -4 , a value of γ equal to 4 × 10 -3 ω 1 is sufficient to remove the ZPEL (ζ = 0.08). Figure 3.c

shows the vibrational spectra obtained in this case with the larger γ: while the ZPEL has been suppressed, the peaks corresponding to the resonances (2ω 2 , ω 1 -ω 2 and ω 1 + ω 2 )

have disappeared. This further illustrates the relation between the mode resonances and the ZPEL. Moreover, increasing γ also leads to a broadening of the peaks of the oscillators in the spectra, consistently with the fact that the full width at half maximum in the case of a harmonic oscillator is γ/2π in a Langevin dynamics and for a spectrum in frequency. The case of the quartic coupling is more complicated and even for large values of γ, the ZPEL is not completely suppressed (figure 4). Figure 3 also shows that increasing γ in the case of a quartic coupling with Ω = 0.2 and c 4 = 15.5 × 10 -4 is not sufficient to suppress all of the resonances between the different modes.

In order to estimate the characteristic time t tr of the energy transfer between the two oscillators, we performed NVE calculations where only oscillator 1 is initially excited. t tr can then be roughly estimated by calculating the typical time at which oscillator 2 starts to get excited. Figure 5 shows the evolution of t tr for Ω = 0.5 as a function of the cubic coupling constant c 3 . As expected, the characteristic time for transfer is directly related to the strength of the coupling. To remove the ZPEL, we need to choose a value for γ that is greater than the typical transfer frequency ω leakage = 1/t tr . For example, in the case of

c 3 = 4 × 10 -4
, we find that t tr ∼ 400ω -1

1
and thus ω leakage ∼ 2.5 × 10 -3 ω 1 . Accordingly, figure 4 shows that a value of γ = 10 -2 ω 1 or higher is necessary to remove the leakage i.e.

the ZPEL is removed if γ ≫ ω leakage . In conclusion, this simple model raises several important issues: the role of resonances and the possibility to remove or at least significantly reduce the effects of the ZPEL by increasing γ beyond the typical frequencies for energy transfer between the modes. We now address these issues on a more complex model.

One-dimensional chain of atoms

We consider a one-dimensional chain of atoms, consisting of 3 oxygen atoms interspaced with 3 hydrogen atoms, with periodic boundary conditions. The interactions between the atoms are described by two interatomic potentials. On the one hand, the O-H interaction is a Morse-type potential derived by Johannsen for hydrogen-bonded systems: 28

V OH (r) = u 0 a + be a(r-r 0 ) a e -b(r-r 0 ) -1 + b e a(r-r 0 ) -1 -u 0 ( 18
)
where r is the O-H distance, u 0 is the height of the potential barrier, r 0 the equilibrium O-H distance, a and b are two parameters. The values of the parameters are set so that: r 0 = 0.96 Å (which corresponds to the length of the covalent bond in the OH -ion), a ≃ 7.11 Å -1 , b ≃ 2.00 Å -1 and u 0 = 2.73 eV so that the O-H stretching frequency (ν OH ) in the harmonic approximation of the potential V OH approximately equals 100 THz. On the other hand, the O-O interaction is described by a standard Morse potential:

V OO (R) = C 0 1 -e -α 0 (R-R 0 ) 2 -C 0 ( 19
)
where C 0 and α 0 are the depth and the width of the potential respectively and R 0 the O-O equilibrium distance. The parameters are the following: C 0 = 3.81 eV, R 0 = 2.88 Å and α 0 varies so that the value of the O-O frequency (ν OO ) lies between 10 and 60 THz. The QTB-MD simulations are performed with a 0.1 fs time step and equilibrium averaged values were obtained using 12 independent trajectories of 3 ns each.

The potential energy of an hydrogen atom is given by V OH (r) + V OH (R -r) which is a double-well potential. Within this model, we can define short "covalent" O-H bonds (∼ 1 Å) and longer "hydrogen bonds" (∼ 1.9 Å). Although this model cannot represent a real physical system such as an ice cluster, it is characterized by realistic O-H frequencies and mode couplings and is useful to assess the nature and effects of ZPEL in realistic hydrogen-bonded systems. A normal mode analysis of the system yields one low-frequency, ν 2 , corresponding to the O-O lattice mode, and two very similar optical high-frequencies, ν 1 , corresponding to the O-H stretching modes. In analogy with the previous model, the O-H stretching modes roughly play the role of the high-frequency oscillator while the O-O lattice mode corresponds to the low frequency oscillator. In the following, we show the influence of the parameter Ω = ν 2 /ν 1 and the friction coefficient γ on the ZPEL at T = 600 K. The frequency ν 2 is varied through the parameter α 0 (eq. ( 19)) while ν OH is fixed at 100 THz (thus, the frequency ν 1 is almost constant). The QTB-MD results are compared with those from PIMD simulations, using a Trotter number P = 20 which ensures a good convergence of all the physical quantities in all cases studied here. For each QTB-MD simulation, we checked that the total energy of the system, as well as the kinetic and potential energies, are in good agreement with the reference values given by PIMD.

In order to evaluate the leakage, we compare the kinetic energy of the light atoms, significantly involved in the high-frequency modes, to that of the heavier atoms, mainly involved in the low-frequency modes. Thus, the effective temperatures T H and T O of H and O atoms are defined from the kinetic energies:

k B T H 2 = 1 N H N H i=1 E (i) k , k B T O 2 = 1 N O N O i=1 E (i) k (20)
where N H = 3 and N O = 3 are the numbers of H and O atoms respectively, and

E (i) k
the average kinetic energy of atom i. In a classical system, equipartition ensures that the kinetic energy is equally distributed among all degrees of freedom: they all have the same effective temperature. This is not true in the quantum case: high-frequency modes have more kinetic energy and their effective temperature is therefore greater. This is the case for the QTB method and for PIMD, which serves as a reference here. From figure 6, one sees that, as expected, the leakage tends to increase the effective temperature of light atoms and decrease that of heavier atoms. In this case, the ZPEL can be quantified through the deviation factor:

ζ = (T H -T O ) (PIMD) -(T H -T O ) (QTB) (T H -T O ) (PIMD) . (21)
ζ = 0 if there is no leakage and 0 < ζ < 1 if leakage occurs and its dependence on Ω is shown in figure 7. Similarly to the coupled harmonic oscillators' model (section 3), ZPEL occurs mostly for Ω ∼ 1/2. Figure 7 also shows that the ZPEL can be substantially decreased by increasing γ, as in the previous model. On the other hand, important ZPEL is observed for Ω < 0.2: this corresponds to a highly anharmonic regime where a structural transition occurs and therefore corresponds to a different physical situation than the other values of Ω. increasing γ, the QTB-MD d OO distribution coincides with the PIMD one.

ZPEL effects on structural properties

ZPEL effects on vibrational properties

We have seen in the case of the two coupled harmonic oscillators that increasing γ has consequences on the vibrational spectrum of the system; in particular, the peaks are broadened and the peaks corresponding to the mode resonances disappear when γ is large enough (see figure 3). Figure 9 shows the vibrational spectrum of the one-dimensional chain of atoms for Ω = 0.5 and for two different values of γ (0.2 and 10 THz). For γ = 0.2 THz, ZPEL occurs while for γ = 10 THz, the ZPEL is almost fully removed (see figure 7). We can see that increasing the friction coefficient leads to broader peaks as expected. However, the positions of these peaks hence the mode frequencies, are not modified by the large value of γ. Even with a large damping term, the vibrational spectrum still yields useful information about the mode frequencies in this case.

In conclusion, in a system containing different chemical elements, the kinetic energy ratio between them can be used as an indicator of the ZPEL. As in the case of two coupled harmonic oscillators, the ZPEL is intrinsically related to resonances between vibrational modes. Correlatively, increasing the friction coefficient allows to remove the ZPEL. In this case, the quantum structural properties are well reproduced; the dynamical properties should be regarded with caution but the vibrational spectrum still contains useful information.

Applications to realistic systems

We now investigate the effect of the friction coefficient γ on the ZPEL for systems with many degrees of freedom.

Lennard-Jones Aluminium

Using QTB-MD simulations, Bedoya-Martínez et al. 20 have evidenced the ZPEL at T = 10 K in a crystal of aluminium modeled by a Lennard-Jones potential (ǫ/k B = 1450.6 K, σ = 2.54 Å, cutoff = 1.37σ = 3.49 Å). In their paper, they showed that the energy is transferred from the high-energy modes to the low-energy modes because the QTB method is unable to fully counterbalance this leakage. We carried out QTB-MD simulations using a 1 fs time step and different values of γ. We confirm that, with γ = 0.9 THz, QTB-MD fails to give the correct quantum energy distribution, as illustrated by the full circles in figure 10. Indeed, the resulting distribution is intermediate between the quantum and the classical homogeneous distributions. However, with a higher value of γ (10 THz), the energy distribution from QTB-MD is very close to the expected quantum distribution θ(ν, T ) (eq. 10), as given by the open circles in figure 10. Therefore, for large enough damping, the ZPEL is neutralized by the QTB. The inset of figure 10 provides the evolution of the slope of the energy distribution, normalized by that of the quantum distribution, as a function of γ. The larger the friction coefficient, the lower the ZPEL, up to γ = 9 THz for which a plateau value is reached. From this value of γ upwards, the leakage is satisfactorily reduced and the energy distribution obtained by QTB-MD is the one initially introduced in the colored noise.

The disadvantage of any thermostat involving a damping term, as in the QTB or the Langevin thermostat, is the possible broadening of the vibrational peaks and the possible occurrence of a spurious high-frequency tail in the phonon density of states (DOS). For small values of γ, i.e. when γ is lower than the full width at half maximum (∆ω) of the peaks of the DOS, increasing γ does not significantly perturb the spectrum. Conversely, for large values of γ, the broadening induced by the damping term is of the order of ∆ν = ∆ω/2π = γ/2π. This slope / QD slope (THz)

Figure 10: Scaled kinetic energy distribution of aluminum at T = 10 K as a function of the frequency of the modes, obtained from QTB-MD simulation with γ values of 0.9 and 10 THz, while the frequency cut-off (20 THz) is chosen equal to twice the highest frequency of the system. The solid line corresponds to the quantum distribution (QD) θ(ν, T ) (eq. 10 with ν = ω/2π). Inset: evolution as a function of γ of the slope of the energy distribution, normalized according to the quantum distribution.

issue is shown in figure 11. In the case of standard MD with a Langevin thermostat, the DOS is obtained by normalizing the Fourier transform of the velocity autocorrelation function by k B T . In the case of QTB-MD, k B T must be replaced by θ(ω, T ) (eq. 10). Figure 11.a) shows that, when the ZPEL is removed, the DOS obtained from QTB-MD trajectories is close to that derived from standard-MD. In contrast, figure 11.b) shows that, when ZPEL occurs, the DOS cannot be obtained from the QTB-MD trajectories, since the number of high-frequency modes or low-frequency modes are underestimated and overestimated, respectively.

Barium titanate

BaTiO 3 (BTO) is a strongly anharmonic ferroelectric crystal characterized by a complex energy landscape. Moreover, quantum effects have been shown to influence its structural properties. 29,30 It undergoes a complex sequence of structural phase-transitions 31 as tem-

Conclusion and practical consequences

We have performed a systematic and quantitative study of zero-point energy leakage (ZPEL) in QTB-MD simulations. The aim here is to assess the reliability of the QTB method on various systems with different degrees of complexity. We have found that the ZPEL is intrinsically related to resonances between vibrational modes and, as in realistic systems many modes can resonate, this is hardly avoidable in practice. However, increasing the damping term γ significantly reduces the leakage and can even, in some cases, remove it entirely. A striking example is provided by our results on BTO, as with small damping term, the phase diagram obtained by the QTB method is wrong, while for larger damping, one recovers the complete sequence of phase transitions at the correct temperatures.

This effect can be explained as follows. The QTB method connects a classical system to a thermal bath which follows the quantum fluctuation-dissipation theorem. Therefore, there is no equipartition of the energy since the QTB pumps more energy into high-frequency modes than low-frequency ones. The ZPEL results from the transfer of energy from highfrequency to low-frequency modes: the obtained energy distribution is therefore the result of the balance between QTB pumping and damping on the one hand and energy transfer within the system on the other hand. Increasing the damping term will increase the pumping rate with respect to the internal equilibration and the QTB energy distribution becomes closer to the quantum one. Moreover, when γ is larger than the characteristic frequency of the energy transfer between vibrational modes, the effects of the resonances between these modes are hindered. Hence, a simple and effective way to prevent ZPEL to occur in QTB-MD simulations is to increase the damping term γ.

This raises the issue that within the frame of a Langevin simulation, one should decrease, not increase, the damping term in order not to alter the dynamics of the system too dramatically. Careful analysis of the effect of damping on both structural and dynamical properties tends to show that this question should be addressed with care in each specific case, but that QTB-MD simulations turn out to be relatively robust and yield excellent results as long as one keeps in mind the physics of the problem. Indeed, we have seen that the mode frequencies obtained from QTB-MD vibrational spectra are not dramatically altered by the increase of γ, even though a large γ implies a broadening of the peaks. This allows us, for example, to study the O-H stretching or bending modes in hydrogen-bonded materials since the corresponding frequencies are usually much larger than γ. On the contrary, we expect the low frequency part of the spectrum to be substantially affected by a large damping term.

Therefore, the QTB method is an efficient tool to study a large variety of anharmonic systems provided that the value of the friction coefficient is large enough to ensure that the ZPEL remains negligible. In this case, the QTB method presents several advantages compared to path integral methods : its computational cost is similar to that of standard MD simulations, enabling the study of large and complex systems, and dynamical properties are directly accessible making possible the confrontation of QTB-MD results to spectroscopic measurements for example.

A Generation of the random force

This part presents the technique used to generate the fluctuating force R(t) which is a random gaussian variable. This technique has been proposed by Dammak et al. 7 and is based on a procedure proposed by Maradudin et al. to generate random surfaces of specific roughness 24 .

Here, we want to generate the stationary Gaussian process R(t) with the following properties:

R(t) = 0 (22) R(t)R(t + τ ) = +∞ -∞ I R (ω, T )e -iωτ dω 2π (23)
The second equation is the Wiener-Khinchin theorem that relates the autocorrelation of R(t)

to its power spectral density, I R (ω, T ). In QTB-MD simulations, I R is given by eq. 12. The value of the random noise at a time t n = nδt is a Gaussian random variable that can be written as a sum of independent Gaussian random variables X j with weights W j as

R n ≡ R(t n ) = +∞ j=-∞ W j X j+n (24)
The variables X j have zero mean and a standard deviation of unity. From eqs. 23 and 24, the weights W j obey the following relation:

+∞ j=-∞ W j W j-l = +∞ -∞ I R (ω, T )e -iωt l dω 2π (25)
with t l = lδt. W (ω) is defined as the Fourier transform of W (t):

W j ≡ W (t j ) = +∞ -∞ W (ω)e -iωt j dω 2π (26)
Using expression (26), in the continuous limit (δt → 0):

+∞ j=-∞ W j W j-l = 1 δt ∞ -∞ W (ω) W (-ω)e -iωt l dω 2π (27)
Considering W (t) as even and real, the function W (ω), also even and real, is obtained using eq. 25 and 27:

W (ω) = δtI R (ω, T ) (28)
and:

W j = √ δt +∞ -∞ I R (ω, T )e -iωt j dω 2π (29)
In practice, the time and the pulsations are discretized. The Fourier transform are expressed so that

W j = 1 √ N N/2 l=-N/2+1
Wl e -i2πjl/N (30)

X j = 1 √ N N/2 l=-N/2+1
Xl e -i2πjl/N (31) with N the total number of MD steps. By comparing eq. 29 with the integral form of eq. 30:

W j = 1 √ Nδω +∞ -∞
W (ω)e -iωt j dω (32) and using Nδtδω = 2π, we obtain that

Wl = 1 √ Nδt I R (ω l ) (33)
with ω l = lδω. From eq. 30 and 31, R n writes

R n = N/2 l=-N/2+1
W-l Xl e -i2πnl/N (34) and using eq. 33, R n finally is: Rl e -i2πnl/N (36)

and then obtain that

Rl = I R (ω l ) δt Xl (37)
the Gaussian random variables Xl can be rewritten as

Xl = Ml + i Ñl √ 2 (38)
with Ml and Ñl independent Gaussian random variables with zero mean and a standard deviation of unity. Moreover, to ensure that the variables Xl are real, Ml = M-l and Ñl = -Ñ-l are required. Finally:

Rl = I R (ω l ) 2δt Ml + i Ñl (39)
In practice, the random forces R n are obtained using the following steps:

B Langevin barostat

This section gives the equations of the Langevin barostat used to fix the hydrostatic pressure in BTO. The extension of the Langevin method to the isothermal-isobaric ensemble has been achieved by Quigley and Probert 34,35 , giving rise to an algorithm in which random and friction forces are applied, not only on the atomic coordinates, but also on the supercell vectors. In the following expressions, second-rank tensors are written in bold. The equations of motion on the local mode i (with mass m) using the Langevin barostat are:

d p i dt = f i -γ p i + R i - p G W g p i - 1 N f . T r(p G ) W g p i (40)
with f i = -∇ u i Φ( u 1 , ..., u N ) the internal force. The terms -γ p i and R i correspond to the friction and the random forces of the thermostat (Langevin or QTB). The momentum p i is related to the position u i by

d u i dt = p i m + p G W g u i (41)
while the matrix of the supercell vectors h and its conjugate momentum p G evolve according

to dh dt = p G h W g (42)
and

dp G dt = V (t)(X -P ext Id) + 1 N f i p 2 i m Id -γ G p G + L G ( 43
)
in which V (t) is the supercell volume, W g is the "mass" associated to the barostat, N f is the number of degrees of freedom, P ext is the external pressure, Id is the identity tensor and X is the internal pressure tensor 34 . In the right member of Eq. 43, one recognizes a friction force on the supercell -γ G p G (γ G is a friction coefficient for the barostat) and a random force L G , a 3 × 3 matrix whose components are randomly drawn at each time step in a gaussian with variance 2γ G Wgk B T δt . This random force on the barostat is symmetrized at each time step to avoid global rotation of the supercell during the simulation.

Figure 1 :

 1 Figure 1: Average energies, ǫ 1 and ǫ 2 , of the two oscillators, and average coupling energy ǫ c computed by QTB-MD as a function of the intensity of the coupling constants c 3 and c 4 . Top panel: cubic coupling (c 3 = 0, c 4 = 0) with Ω = 0.5. Lower panel: quartic coupling (c 4 = 0, c 3 = 0) with Ω = 0.25. By symmetry, ǫ c = 0 in the cubic case.

oscillator 1 oscillator 2 Figure 3 :

 23 Figure 3: Vibrational spectra (in logarithmic scale) of oscillators 1 and 2 obtained by QTB-MD simulation in the case of a cubic coupling (top panel) with Ω = 0.5 and c 3 = 2.4 × 10 -4 , and in the case of a quartic coupling (lower panel) with Ω = 0.2 and c 4 = 15.4 × 10 -4 . The spectra are computed for two selected values of the friction coefficient: γ = 4 × 10 -4 ω 1 and γ = 4 × 10 -3 ω 1 .

Figure 3 .

 3 b shows, in the case of Ω = 0.2 and c 4 = 15.4 × 10 -4 , that many other modes than ω 1 and ω 2 also appear in the spectrum. Hence, multiple resonances are likely to occur leading to significant ZPEL for

Figure 4 :

 4 Figure 4: Effect of the damping coefficient γ (given in units of ω 1 ) on the energies of the two oscillators (eq. 16) as a function of the coupling constants. Top panel: cubic coupling (c 3 = 0, c 4 = 0) and Ω = 0.5. Lower panel: quartic coupling (c 4 = 0, c 3 = 0) and Ω = 0.25. The solid line and symbols represent the results obtained from QTB-MD and the grey dashed lines represent the exact results.

Figure 5 :

 5 Figure 5: Energy transfer time between the two oscillators t tr (in units of ω -1 1 ) as a function of the cubic coupling constant c 3 for Ω = 0.5 estimated from NVE simulations where only one oscillator is initially excited. Here, c 4 = 0.

Figure 6 :

 6 Figure 6: Effective temperatures T O and T H of oxygen and hydrogen atoms calculated with QTB-MD and PIMD at T = 600K. Here Ω = 0.5 and γ = 0.2 THz. The arrows indicate the temperature shifts due to the ZPEL.

Figure 8 Figure 7 :

 87 Figure 7: ZPEL quantified by the deviation factor ζ using the definition in eq. (21) as a function of the frequency ratio Ω.

Figure 8 :

 8 Figure 8: Probability distribution of the O-H (a) and the O-O (b) distances computed from standard-MD, PIMD and QTB-MD, for selected values of the friction coefficient γ (0.2 and 20 THz) and for Ω = 0.5.

Figure 9 :

 9 Figure 9: Vibrational spectra obtained by QTB-MD simulation for Ω = 0.5 and two values of the friction coefficient: γ = 0.2 THz (blue dashed line) and γ = 10 THz (red full line).

Figure 11 :Figure 12 :

 1112 Figure 11: Vibrational density of states (DOS) as a function of the frequency for different values of γ. They have been computed through the Fourier transform of the velocity autocorrelation function which is normalized by k B T in the case of Langevin MD or by θ(ν, T ) (eq. 10) in the case of QTB-MD. Two values of the friction coefficient γ are used: (a) 10 THz and (b) 2 THz.

I

  R (ω l ) Xl e -i2πnl/N(35) Let us define the discrete Fourier transform Rl such that

1 . 1 2. 3 .

 113 Generation of independent Gaussian random numbers Ml and Ñl for l = 1, .., N/2 -Computation of Rl using expression (39) Symmetrization of Rl : Rl = RN-l for l = N/2 + 1, .., N -1 4. Cancellation ( Rl = 0) for l = 0 and N/2 5. Computation of R n from eq. 36

Acknowledgement

We gratefully acknowledge important contributions to this work by Grégory Geneste, who provided us the code to perform the calculations on BTO and enlightened many issues about ferroelectric transitions. We also thank Jean-Louis Barrat for providing the values of the LJ parameters they used in their paper to describe the Aluminium crystal. Y.B. acknowledges financial support from the Conseil Régional d'Île-de-France through the DIM Oxymore. This work was performed using HPC resources from the "mésocentre" computing center of CentraleSupélec, which is supported by CentraleSupélec and CNRS.


M Thomas
email: marc.thomas@onera.fr

T Malot

P Aubry

C Colin

T Vilaro

Pierre Bertrand

The prospects for additive manufacturing of bulk TiAl alloy

Keywords: gamma titanium aluminide, selective laser melting, laser metal deposition, electron beam melting, heat treatments

  L'archive ouverte pluridisciplinaire

Introduction

Regarding the development of advanced materials, the manufacturing techniques cover the compositional and microstructural aspects to satisfy the industrial requirements, but also the shape building aspect to achieve the geometrical, technical and functional properties of components. Thus, attempts are currently made in the aerospace sector to meet such property requirements by using near net or net shape processing.

Due to the rapid progress made in the field of Additive Manufacturing (AM) these last years, such technologies seem suitable to manufacture net shape and fully dense complex metallic components.

The attractiveness of these direct manufacturing technologies is that it reduces the time to market of products with no specific tooling required. Practically, most of the AM processes start from a precursor powder that is consolidated layer-by-layer using directly a three dimensional computer aided design (CAD) data file. In the case of aeroengine components, AM is highly appealing for several reasons. First of all, it allows a better buy-to-fly ratio owing to a low material waste factor.

Secondly, AM shortens the lead time by reducing machining and joining steps and by avoiding time consuming tooling reworking. Thirdly, it allows weight saving and additional functionalities by adding supplementary fonctions such as cooling channels. Hence, AM already appears as a costeffective processing route for a number of existing aeroengine components.

Available techniques for metals are laser metal deposition (LMD), electron beam melting (EBM), and selective laser melting (SLM). For laser metal deposition (LMD) process, the powder is injected through a nozzle onto the substrate to be molten layer by layer. The other two processes (SLM and EBM) consist of creating a powder bed where the powder is deposited through a roller or a blade into a compacted layer which is then selectively melted by the scanned laser or electron beam, with the part being built up layer by layer. Now, for an improved energy efficiency and power in gas turbine engines, an increase of their maximum work temperature is required. Therefore, the development of advanced materials with high temperature resistance such as gamma titanium aluminides attracted considerable attention in the last two decades. Due to their ordered structure and to their lightweight, this family of intermetallic alloys exhibits a good combination of properties, such as a high melting point, low density and strength at temperatures up to ~750°C [START_REF] Liu | Tensile properties and fracture toughness of TiAl alloys with controlled microstructures[END_REF] . Recently, the Ti-(47-48)Al-2Cr-2Nb alloy has been introduced in the low pressure (LP) turbine of the GEnX jet engine (General Electric) which powers Boeing787 and 747-8 [START_REF] Clark | Propulsion System[END_REF] .

It is recalled that the conventional processing route for LP turbine blades is either gravity or centrifugal casting. The latter is preferred to limit the formation of gas porosities and of mis-runs.

However, cast processing generally yields a relatively coarse microstructure with metallurgical defects such as porosity shrinkage. This leads to relatively poor room-temperature mechanical properties, in particular low ductility and toughness. The lack of ductility also renders machining difficult. Moreover, wide property variability is generally observed for different thickness of cast parts and location of test specimens. In a general manner, these casting related features might restrict the use of γ-TiAl alloys to simple shape and moderate dimension. Therefore, alternative processing routes, including AM technologies, have been investigated with the aim to produce complex shaped TiAl components with a more uniform microstructure. A number of experimental studies already deal with the process of γ-TiAl by means of Additive Manufacturing (AM).

1° LMD or cladding process

In the mid-nineties, methods for joining and for surface modification were investigated by laser welding for a Ti-46Al-2Mo alloy [START_REF] Uenishi | Processing of Intermetallic Compounds for Structural Applications at High Temperature[END_REF] . Crack-free welds could be obtained only for low scan rates and with pre-heating temperatures above 400°C. First laser cladding trials that have been carried out on titanium aluminide alloys range back to the late nineties [START_REF] Griffith | Proceedings of the International Congress on Applications of Lasers and Electro-Optics (ICALEO'98)[END_REF] . In a general manner, laser cladding operations were aimed at overlaying a base material, typically TA6V, with TiAl to form a protective layer against oxidation. Dilution rate is an important factor since some dilution is required to guarantee a sound metallurgical bond, but in a limited extent to reduce the problems associated with distortion and a brittle heat-affected zone. For instance, in a more recent work of laser cladding coatings of Ti-48Al-2Cr-2Nb alloy on Ti6Al4V, some dilution with base material that may degrade coating properties, was observed since the coating composition (Ti-46Al-2Cr-2Nb at.%.) was slightly richer in Ti than the original powder (Ti-48Al-2Cr-2Nb at.%.) [START_REF] Cárcel | Laser cladding of TiAl intermetallic alloy on Ti6Al4V. Process optimization and properties[END_REF] . However, most of their work was based on the important cracking sensitivity of this intermetallic material. By using penetrating fluid tests for single tracks, uncracked conditions were identified by using lower powder feed rate and lower scanning speed. This was related to the decreased cooling rate for lower scanning speed which prevents cracking. In addition, the authors reported that heating previously and during the process tends to improve the results in terms of cracking.

The LMD process was also evaluated to process bulk Ti-47Al-2Cr-2Nb samples [START_REF] Moll | Laser Forming of Gamma Titanium Aluminide[END_REF] . Provided a continuous deposition was ensured, thermal stresses were reduced and the deposited material contained no cracks. However, not much detail is given about the process parameters. Subsequent heat treatments with super-transus conditions, i.e. in the single α phase field, were then performed to homogenize the microstructure. However, this did not seem to improve tensile properties with respect to the as-built condition. In the meantime, Srivastava et al have demonstrated that the heterogeneous microstructure of a laser fabricated Ti-48Al-2Mn-2Nb alloy was associated with the many thermal excursions as the successive layers are deposited [START_REF] Srivastava | The influence of thermal processing route on the microstructure of some TiAl-based alloys[END_REF] . For this particular alloy, the processing parameters that have significant effects on the quality of the component have been identified [START_REF] Srivastava | The optimisation of processing parameters and characterisation of microstructure of direct laser fabricated TiAl alloy components[END_REF] . The amount of laser energy and the volume of powder were found to determine the stability of the build-up rate and therefore the uniformity of build geometry and quality [START_REF] Srivastava | The effect of process parameters and heat treatment on the microstructure of direct laser fabricated TiAl alloy samples[END_REF] . However, crack-free samples and a good quality surface finish were not achieved in their work.

According to Zhang et al, the focus of the powders must be below the deposition plane for an optimum performance in terms of build-up rate and coupling of the deposit to the substrate [START_REF] Zhang | Characterization of Laser-Deposited TiAl Alloys[END_REF] . It was also found in this study for the Ti-48Al-2Cr-2Nb alloy that each deposited layer consists of a metastable microstructure with the massive γ phase on the top of the layer while the high temperature α-phase, which is then ordered to α 2 -Ti 3 Al, was retained on the bottom of the layer. Much higher cooling rate was obtained at the bottom of the layer because of a rapid rate of heat extraction.

The directionality of thermal gradients along the Z-built axis tends to favour texture effects. For the stoichiometric binary TiAl alloy, sequential laser remelting tends to create the directional growth of metastable γ phase with (111) plane perpendicular to the thermal gradient direction [START_REF] Liu | Directional growth of metastable phase γ in laser-remelted Ti-Al[END_REF] . Qu et al investigated the phase transformation for the Ti-47Al-2.5V-1Cr alloy manufactured by LMD [START_REF] Qu | Microstructure and mechanical properties of laser melting deposited γ-TiAl intermetallic alloys[END_REF] . A directionally solidified columnar grain structure was formed in the deposition direction due to the high temperature gradient. Subsequent solid state phase transformation took place in the heat affected zone of the deposited material and led to the formation of a fully lamellar microstructure with colony size in the range of 50-100 μm and lamellar spacing of 0.3-0.5 μm. At room temperature, a lack of tensile ductility (0.6%) was observed in both longitudinal and transverse sections of the as-deposited samples. The authors further investigated the effect of heat treatments using the same Ti-47Al-2.5V-1Cr alloy [START_REF] Qu | The effects of heat treatment on the microstructure and mechanical property of laser melting deposition γ-TiAl intermetallic alloys[END_REF] . However, the compositional inhomogeneity of the asdeposited samples were not completely removed after heat treating at 1125°C for 30 min. So, the tensile ductility of test samples was not really improved by this heat treatment at low temperature, with maximum values of 0.6%.

2° The EBM process

Regarding now the development of the EBM process, first attempts to consolidate Ti-47Al-2Cr-2Nb powder were made by Cormier [START_REF] Cormier | Electron beam melting of gamma titanium aluminide[END_REF][START_REF] Cormier | Freeform fabrication of titanium aluminide via electron beam melting using prealloyed and blended powders[END_REF] . This earlier work revealed a dramatic decrease of aluminium from 46 at.% to 39 at.% by comparing the original powder with the EBM-fabricated test sample. This is related to the fact that aluminium is a volatile element with high vapor pressure, which is evaporated from the surface due to the vacuum exposure. However, it is difficult to identify a clear trend in their work about the Al loss conditions because of the scarce information given on the process. melting (EBM) [START_REF] Murr | Characterization of titanium aluminide alloy components fabricated by additive manufacturing using electron beam melting[END_REF] . Residual porosity was attributed both to non-optimized build conditions and to Ar bubbles from the atomised powder. Surprisingly, the authors did not observe any changes in the chemistry of the coupons compared to the original powder, in particular no aluminium loss. This is presumably due to the fact that the EBM machine did operate without multiple melt passes required for optimized settings.

Rastkar et al evaluated the possibility of making hard surface layers by using EBM with a Ti-45Al-2Nb-2Mn-1B alloy [START_REF] Rastkar | Surface transformation of Ti-45Al-2Nb-2Mn-1B titanium aluminide by electron beam melting[END_REF] . By using high power and high scanning speed, the cooling rate is fast enough to largely suppress α to α 2 phase transformation and promote a microstructure with the hard ordered B2 phase. On the other hand, a decreased scanning speed induces higher local superheat, thereby increasing the evaporation of aluminium [START_REF] Semiatin | Diffusion Models for Evaporation Losses during Electron-Beam Melting of Alpha/Beta-Titanium Alloys[END_REF] . Ge et al also reported the strong influence of energy input on Al evaporation for Ti-47Al-2Cr-2Nb powders [START_REF] Ge | Effect of process parameters on microstructure of TiAl alloy produced by electron beam selective melting[END_REF] . An Al loss up to 15 at.% was observed at the highest energy input. As a consequence, the phase transformation in samples was really different depending on the energy input.

Biamano et al investigated the possibility to produce fully dense near-net-shape parts made of the Ti-48Al-2Cr-2Nb alloy composition [START_REF] Biamino | Electron beam melting of Ti-48Al-2Cr-2Nb alloy: Microstructure and mechanical properties investigation[END_REF] . A process optimization was needed to completely eliminate the residual porosity located between two subsequent layers. However, smaller spherical pores were found to remain in the powder particles due to argon gas entrapment during the atomization process. In parallel, the Al loss was reported to be approximately 1 at.%. The as-built samples were then hot isostatically pressed (HIPed) at 1260°C under a pressure of 170 MPa for 4 h, thus reducing the residual porosity down to 1%. A beneficial effect of HIP and subsequent heat treatment at 1320°C for 2 hours was to significantly reduce the data scattering in tensile properties, i.e. down to 6% for YS at room temperature.

Schwerdtfeger et al made a comprehensive study of the effect of EBM process parameters on the resulting microstructures and on the Al loss of the Ti-48Al-2Cr-2Nb alloy [START_REF] Schwerdtfeger | Selective electron beam melting of Ti-48Al-2Cr-2Nb alloy: Microstructure and aluminium loss[END_REF] . A too high energy input should not be recommended, first because this results in insufficient interlayer connection due to significant differences in powder layer thickness deposited in the successive layers, and second because of a higher Al evaporation. Attempts were then made to decrease overheating in the melt pool for an expected reduced Al evaporation. Under optimized processing conditions, an Al loss as low as 0.5 at.% was then achieved. In more details, a decrease in overheating was obtained by lowering layer thickness and by increasing scan speed. However, this limited Al loss is obtained at the expense of the microstructure because the preferred lamellar microstructure with small colony sizes tends to develop for low scan speeds. Furthermore, because the Al evaporation takes place in the top region of the melt pool, an inhomogeneous distribution of Al element cannot be prevented through the sample height. Local variations in the Al concentration can then affect the solidification pathway and subsequently the size and morphology of the microstructure.

A high niobium content alloy was also investigated for blade production and mechanical testing [START_REF] Sabbadini | Third generation gamma TiAl alloys: recent developement at Avio[END_REF] .

In the as-built condition, the Ti-45Al-8Nb-2Cr alloy develops a fully lamellar microstructure with lamellar grains of about 100 μm. A higher oxidation resistance was also reported for this alloy with respect to Ti-48Al-2Cr-2Nb. Test results on blades show that FE prediction of critical location is correct and fatigue limit is comparable with the values obtained on test samples.

3° The SLM process

Whereas the SLM process is already used for a number of metallic materials (steels, titanium alloys, nickel base superalloys, cobalt chromium and aluminum alloys), it has not been demonstrated so far as a suitable process for γ-TiAl alloys. In 2011, Löber et al made some attempts of SLM processing of a Ti-48Al-2Cr-2Nb alloy which resulted in multi-cracks samples [START_REF] Löber | Selective Laser Melting of Titaniumaluminides[END_REF] . More recently, Löber et al used a TNM alloy composition which is also susceptible to cracking because of the small volume fraction of β-phase as the starting powder [START_REF] Löber | 'Selective laser melting of a beta-solidifying TNM-B1 titaniumaluminide alloy[END_REF] . A very limited fraction size of powder (45-63 μm) with a layer thickness of 75 μm was used. Laser power and scanning speed were adapted to stabilize the process and to obtain the best morphology for single beads. Fragmentation of the beads, the so-called balling effect, was reported at low laser power (50 W) and low scanning speeds (50-100 mm/s), but also at intermediate laser power (100 W) and very high scanning speeds (350-2100 mm/s). Pronounced cracking which was related to high cooling rate during solidification, was also observed at low power and low scanning speeds (150-300 mm/s). A set of process parameters in an intermediate range (scanning speed: 150-300 mm/s; laser power: 100-250 W) was identified to generate stable beads without fragmentation or crack formation. Now, because of a non-uniform cooling rate during the SLM process, the microstructure was found to be inhomogeneous with coarser and finer regions next to each other. A two-step heat treatment caused a microstructural homogenization in the first step at 1230°C in the α/β-region and allowed the precipitation of γ lamellae in the second step at 950°C. Subsequent to this heat treatment, a duplex microstructure consisting of β and α 2 grains with some γ lamellae in α 2 /γ colonies is obtained.

In these prior investigations, one can note that no much work was devoted to AM of bulk TiAl materials. In a recent paper, it has been demonstrated that optimal process parameters can be determined for LMD process to prevent cracking due to built-up residual stresses during fast cooling [START_REF] Thomas | Laser Metal Deposition of the Intermetallic TiAl Alloy[END_REF] . Post heat-treatments were then used to restore homogeneous lamellar or duplex microstructures. Subsequent room temperature tensile tests provided a good validation of the TiAl material soundness and homogeneity. Based on this promising result with the LMD process, the present work is aimed at investigating the influence of the SLM process parameters for the commercial Ti-47Al-2Cr-2Nb (at.%) alloy in order to produce fully dense parts without cracking.

Conventional post-heat treatments were also given to the as-built material in order to homogenize the microstructure.

Experimental procedure

The atomised powder was sieved to collect powder particles smaller than 50 µm. Chemical analysis of the selected powder which was performed by plasma emission (ICP) is given in Table I. The SLM tests were carried out by means of a Phenix PM100 machine. Successive powder layers of 60 μm high were deposited without any additional densification of the layer. To prevent any melt pool oxidation, a neutral atmosphere was maintained in the process chamber with an argon flow of 50 l / min. A series of experiments have been performed to identify the effect of process parameters such as the scan speed, the diameter of the laser beam, the overlapping ratio and the preheating temperature of the TA6V substrate. Our purpose was to identify the range of values of individual processing parameters for which stable build-up rate can be obtained. For microstructural observations, samples were cut by electro discharge machining, then embedded in a conductive resin for subsequent mechanical polishing until a mirror finished-surface was obtained. SEM observations of un-etched specimens were performed on a Zeiss 962 DSM equipment, by using either back-scattered or secondary electron conditions. For optical observations, some samples were etched with the Kroll reagent to reveal the microstructure as well as the metallurgical defects.

Experiments were also conducted on the possibility of recycling the powder by using the unmelted powder during the built. The oxygen content and the morphology of recycled powder particles were analyzed. For this investigation, conductive pads were used to seize the powder on the substrate after each laser melt.

Experimental results

The first process parameter to be considered was the speed of the laser beam, which can affect the morphology of the molten bead. Different scanning speeds ranging from 10 mm/s to 160mm/s have been compared in Figure 1. Irregular beads are produced with low speeds as illustrated in Figure 1a.

An increasing scanning speed tends to reduce the apparent width of the bead. At intermediate speeds, the bead formation remains stable as revealed in Figure 1b. Scanning speeds between 30 and 60 mm/s yield the most stable bead. On the other hand, higher scanning speeds induce a fragmentation effect which is characterized by elongated molten discontinuous materials as shown in Figure 1d. This instability of melt pool is probably due to surface tension, which favors the formation of a fragmented bead. It is interesting to note that these different beads showed no visible crack using optical microscopy. In a second stage, XY 10 * 10 mm surface layers have been produced with different overlapping ratios according to the favorable scanning speeds. The overlapping ratio cannot be optimized independently of the scanning speed since a lower speed results in wider beads. The optimised overlapping ratio was determined by comparing different overlapping ratios from 0 to 30% between adjacent beads. A too high overlapping ratio contributes to significantly increase the lead time.

Conversely, if the overlapping is not sufficient, then this can lead to the formation of elongated cavities in the solid sample, as observed in Figure 2. The last stage which corresponds to the manufacturing of bulk cubic samples, was conducted for an overlapping ratio of 25%. To overcome the problem of porosities for low overlapping ratio, a "double melting" strategy was investigated. This consists of a second laser scanning over the entire surface before adding the next powder layer in order to smooth the surface. Indeed, melting the surface by laser is a method applied industrially for improving the roughness. However, this built strategy was found to enhance the risk of interrupted laser due to the high reflectivity of smooth surfaces. Consequently, a "cross strategy" was adopted, even though the surface roughness was not so good as illustrated in Figure 3.

Moreover, the process chamber was pre-heated at 200°C to allow the residual stresses to be relieved during the process due to the cracking sensitivity of the TiAl material. taken as a criteria in order to identify the best processing conditions. It appears that the distance between cracks is reduced, i.e. the number of cracks is increased for higher speeds. This is consistent with the fact that the residual stresses are higher at high speeds, which logically generate more cracks. On the other hand, the average distance between the cracks does not clearly change for different overlapping ratios.

Even if the speed of the laser beam is the predominant parameter for the total energy absorbed by the material, the microstructural evolutions for different scanning speeds are relatively low. Quite expectedly, a higher speed produces a slightly finer microstructure.

As emphasized previously in some published works, the intermetallic TiAl alloy tends to crack during multiple heating and cooling cycles due to fast cooling rate. Therefore, processing conditions that induces a lower temperature gradient of the solidification front should be favorable to prevent cracking. Thus, an extended time of interaction between the laser beam and the powder should decrease the cooling rate. However, the cooling rate is higher by SLM than by LMD. For instance, by SLM it seems that a delay time of half a second before each new laser melting does not eventually lead to a more complete cooling. With or without delay time, SLM exhibits fast cooling. This is also the reason why the resulting microstructure is much finer than that obtained by LMD.

Moreover, due to the cumulative heating and cooling cycles, the microstructure is quite heterogeneous. This is the reason why some heat treatments were applied to homogenize the microstructure.

A first heat treatment was carried out at 900°C for 415 hours under argon. In this case, the initial dendritic structure is decomposed in γ and α 2 grains of approximately 5 μm in size (Figure 4a). A few lamellar colonies can also be detected. The second heat treatment consists of an annealing at 1250°C for 4h under argon followed by air cooling and a further ageing at 900°C for 4h under argon also followed by air cooling, giving rise to a duplex microstructure with lamellar γ+α 2 colonies as well as γ grains (Figure 4b). A few Nb-rich precipitates are also present at some triple boundaries. Right in the middle of Figure 4b, a large α phase is characterized by some discontinuous coarsening.

Conclusion

A large number of SLM experiments were conducted under controlled atmosphere by changing the processing parameters to manufacture a series of TiAl beads, XY surface layers and cubes. Whereas LMD process parameters have been identified to eliminate the cracking effect for the TiAl intermetallic alloy, the SLM process induces higher cooling rates that makes difficult to alleviate the cracking effect. By changing the process parameters in the sense that more heat generates lower cooling rate, crack-free samples should be produced. In the present study, the range of process parameters did not allow to completely suppress the cracking effect. Further investigations by changing the process parameters in a larger range and by introducing a preheating source should be able to produce sound bulk TiAl samples.

Conventional post-heat treatment can be successfully used to obtain a more uniform microstructure, which is important to guarantee minimum values for mechanical properties. In that sense, industrial components produced by additive manufacturing appear to be feasible for this intermetallic TiAl material.

Table captions

Table I : Chemical analyses of the initial powder

Figure 1 :

 1 Figure 1: Optical microscopy of beads manufactured by SLM with different speeds: (a): 10 mm / s; (b): 50 mm / s; (c): 70 mm / s; (d): 160 mm / s

Figure 2 :

 2 Figure 2: Optical microscopy of XY surface layers using the following parameters: power 50 W, speed 60 mm/s, distance between two laser lines 200 μm

Figure 3 :

 3 Figure 3 : SEM micrographs using secondary electrons of the surface of SLM samples using the following parameters: power 50 W, speed 50 mm/s, distance between beads 160 μm: (a) : crossstrategy ;(b) : double melting strategy

Figure 4 :

 4 Figure 4 : SEM micrographs using back-scattered electrons of SLM samples after heat treatments: (a): 415h at 900°C; (b): 4h at 1250°C + 4h at 900°C

Figure 5 :

 5 Figure 5 : SEM micrographs of the TiAl powder; (a): fresh powder; (b) used powder

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 12345 Figure 1: Optical microscopy of beads manufactured by SLM with different speeds: (a): 10 mm / s; (b): 50 mm / s; (c): 70 mm / s; (d): 160 mm / s

Table I :

 I Chemical analyses of the initial powder

	Ti (at%)	Al (at%)	Cr (at%)	Nb (at%)	Fe (at%)	C / ppm	O / ppm
	48.5	47.21	2.12	2.11	0.03	105	745

Acknowledgements

The authors would like to acknowledge the FRAE (Fondation de Recherche en Aéronautique et Espace) institution for their financial support.


Tuan Anh Pham
email: tuan-anh.pham@inria.fr

Nhan Le Thanh
email: nhan.le-thanh@inria.fr

Checking the Compliance of Business Process in Business Process Life Cycle

Keywords: Business Process Compliance, OWL 2, SWRL, Reasoner

Business process compliance has become more and more important function for business process management (BPM). One of challenges in this area is to check the business process compliance in the business process life-cycle (design-time, run-time). In this paper, we propose a description logic-based approach for business process compliance checking during two phases of the business process life-cycle. In our approach, business process and the set of regulations are represented in a machine readable form. And we use that knowledge bases to check the compliance between them.

Introduction

Business process compliance (BPC) checking in companies is a crucial feature for a BPM system. In essence, many approaches are developed to formally and automatically prove that business processes comply with relevant constraints like rules, laws. The requirement for business process compliance checking based on a set of constraints might emerge in different phases of the process life-cycle. During design time, the compliance of a process model with a set of constraints is checked for detecting the structure error (infinite loop, i.e.) and the compliance violation in the early step. At runtime, the progress of a potentially large number of process instances is monitored to detect or even predict compliance violations. After studying most important publications on this topic, we focus on some following research questions:

 What phase of the business process life-cycle our research approach can apply?  What kind of rule our approach can cover?  How to adapt our solution to existing business process definition language and system?

After providing our analysis on some related works, each research question will be answered in section 4.

Our contributions in this paper are:

 Proposing a method for representing and integrating business process and regulatory, business rules or laws.  Building a method for business process compliance checking. This paper is organized as follows, section 2 reviews the related work for compliance checking of BPM. Section 3 provides our research methodology. Section 4, we present our solution. Finally, Section 6 concludes this paper and indicates next steps of our work.

2

Related Work

Checking the compliance of business processes is a challenging task: the number and complexity of business rules is increasing and the rules are subject to constant change. Becker et al. [START_REF] Becker | Generalizability and Applicability of Model-Based Business Process Compliance-Checking Approaches -A State-of-the-Art Analysis and Research Roadmap[END_REF] define business process compliance management (BPCM) as "the steady modeling, refinement, and analysis of business processes regarding the fulfillment of regulatory compliance". Ramezani et al. outline the interdependencies between BPM and Compliance Management (CM) and describes CM as a "methodology to elicit, specify and formalize, implement, check and analyze, and optimize compliance requirements in organizations" [START_REF] Ramezani | Separating Compliance Man gement and Business Process Management[END_REF]. Works on Business process compliance have focused on examining whether a given process model is compliant with a certain reference model/pattern. On the technical aspect, the business process pattern initiative has identified various patterns for the specification of control-flow [START_REF] Van Der Aalst | Workflow patterns[END_REF], data-flow [START_REF] Russell | Workflow data patterns: Identification, representation and tool support[END_REF], and resources [START_REF] Russell | Workflow data patterns: Identification, representation and tool support[END_REF] in business process management systems. The work in [START_REF] Tran | Broadening the use of process patterns for modeling processes[END_REF] deals also with the planning layer by formalizing process patterns using UML concepts. These compliance works have focused on the structural level of process models, while another line of works focuses on the combination of data and structure [START_REF] Awad | Towards Resolving Compliance Violations in Business Process Models[END_REF][START_REF] Weber | Semantic business process validation[END_REF][START_REF] Governatori | Compliance checking between business processes and business contracts[END_REF][START_REF] Thom | Inventing Less, Reusing More, and Adding Intelligence to Business Process Modeling[END_REF][START_REF] Kharbili | Towards a framework for semantic business process compliance management[END_REF]. The frameworks in [START_REF] Arbab | Towards using reo for compliance-aware business process modeling[END_REF][START_REF] Kharbili | Towards a framework for semantic business process compliance management[END_REF][START_REF] Gschwind | Applying patterns during business process modeling[END_REF], for example, provide general compliance criteria for assessing the compliance of processes with semantic constraints. In addition, some compliance works aimed at supporting specific purposes, for example: correcting process models at design time [START_REF] Arbab | Towards using reo for compliance-aware business process modeling[END_REF], verifying changes in existing models [START_REF] Arbab | Towards using reo for compliance-aware business process modeling[END_REF], identifying compliance in the context of process mining [START_REF] Gschwind | Applying patterns during business process modeling[END_REF], and identifying violations of execution order compliance rules [START_REF] Awad | Towards Resolving Compliance Violations in Business Process Models[END_REF].

In our approach, we take advantage of Color Petri Net (CPN)'s [START_REF] Jensen | Coloured Petri Nets[END_REF] color set for checking the compliance of processed data in data-flow with the constraints in the set of business rules. A color set can be defined in many types: int, string or object. When firing a transition, the value of each color can be changed, the new value must respect to the constraint in the set of predefined rules. This work will be explained in more detail in section 4.

On the other hand, for representing the business rules, some works use an ontologybased approach [START_REF] Reynares | An automatable approach for SBVR to OWL 2 mappings[END_REF] [9] [10] [START_REF] Bernotaityte | Developing SBVR Vocabularies and Business Rules from OWL2 Ontologies[END_REF]. They translate the Semantic of Business Vocabulary and Rule (SVBR) [16] vocabulary to OWL [29] and Semantic Web Rule Language (SWRL) [START_REF] Horrocks | Swrl: A semantic web rule language combining owl and ruleml[END_REF], they provide the mapping or the rule in order to translate each property of SVBR (the definition of OMG) to a set of axioms in an ontology.

In our approach, we classify the business rule into five main type of business rules, and we use also Attempto Controlled English (ACE) [START_REF] Auer | OntoWiki -A Tool for Social, Semantic Collaboration[END_REF] for defining the business rule, each ACE phrase will be represented by an axiom in business rule ontology. We consider not only the business rule representation aspect but also the compliance of business process with a set of business rules.

The advantage of our approach is to ensure that the business process is well-defined at design time and executed correctly during runtime. We must consider this aspect because a business process must always respect a set of predefined rules during the business process life-cycle. Therefore, the difference between our work and the previous related works is that our solution allows the system to check not only the consistency of business process and business rule but also the consistency of the integration between them at design-time and run-time. And the other difference is to use the reasoning for preventing the potential sematic error of a business process instance which occurs on the deduced knowledge or the generated data during both phrases of business process life-cycle (design-time, run-time), this kind of error cannot be detected by some other business process compliance checking approach (query-based approach, i.e.).

Research Methodology

Fig. 1. Research Methodology

Initially, we started with requirements analysis based on case studies. We also formulated a state-of-the-art of business process compliance. Based on the state of the art, we choose a theoretical basis and propose a solution that is also implemented in prototypes.

The solutions do not necessarily address all aspects of the requirements analysis at once but may also focus on certain aspects. Using the developed prototypes, we are able to analyze and evaluate developed solutions using data from practical applications. This may lead to further development and implementation iterations (e.g., in case the developed concepts do not yet cover all relevant aspects or do not yet yield adequate solutions). The evaluation of developed solutions may also result in a completely new iteration leading to modifications or refinements of the solution when studies reveal additional requirements.

Requirement analysis

State of the art

Theory and Solution proposing

Design and implementati on

Evaluation

In this section, we will answer the mentioned research questions above. Fig 2 illustrates the sketch of our approach.

Fig. 2. Sketch of our solution

In Fig 2, a business process can be defined by some graphic design language (BPMN, UML, i.e.). Depending on the selected design language, the user's defined model will be translated into a Color Petri Net [START_REF] Jensen | Coloured Petri Nets[END_REF] (CPN) graph for checking the correctness of the model at the design time. In order to do that automatically, the CPN graphs are represented in a machine readable form. Description logic and first order logic are chosen to describe the model and the constraint for ensuring the structure correctness of the business process model. An inference engine (reasoner) is used to reason and verify the consistency of the model's knowledge base. If the knowledge base is consistent, the business process is validated. On the other hand, a set of business rules is represented in a business rule knowledge base. The inference engine is used to check the consistency of this knowledge base for ensuring that there are no conflicts inside the set of business rules. After having two consistent knowledge bases, we combine them into one knowledge base. The combined knowledge base is used for checking the compliance of a business process instance with a set of business rules inside business rule knowledge base.

Representation of Business Process Knowledge Base

As mentioned above, we use Color Petri Nets for checking the model of a business processes; in this section, we introduce the method of building business process, this method helps the user to represent a business process by an ontology. The advantage of this method is to allow the user to check the consistency business process automatically by the reasoning. The TBox of business process ontology is defined as follow: CPN is the concept for representing all CPN graphs. A CPN graph is well-defined if and only if it has at least one place, one transition, one input arc and output arc. We define in BPO following classes: CPN, Place, Transition, OutputArc, In-putArc and some properties which define the relations between them, hasPlace, hasTransition, hasInputArc, hasOutputArc. Place represents the properties of place, we define a concept Place. A place may have a token or not, it has also at least one InputArc or one OutputArc. The concept Transition is defined for all transitions. A transition must have at least one InputArc and one OutputArc. It's one of the minimum conditions for having a well-defined CPN graph. A transition may have only one guard function or not. The concept InputArc defined for all input arcs. An input arc has only one source place and one target transition. It may be marked by only one expression or not. An OutputArc has only one source transition and only one target place. It also may have only one expression or not. A business process instance is a set of individuals inside the business process ontology (BPO).

Representation of Regulatory

The different structural categories of business rules are (Wagner 2005):

1. Integrity (or constraints); For example: Each company must have one and only one director. 2. Derivation (conditions resulting in conclusions); For example: Platinum customers receive a 5% discount. John Doe is a platinum customer. As a conclusion, John Doe receives a 5% discount. 3. Reaction (Event, Condition, Action, Alternative action, Post-condition); For example: An invoice is received. If the invoice amount is more than $2,000 then a supervisor must approve it. 4. Production (condition, action); For example: If there are no defects in the last batch of cars then the batch is approved. 5. Transformation (change of state); For example: A man's age can change from 28 to 29, but not from 29 to 28.

Business Rules Ontology.

In this section, we introduce the method of building a business rule ontology (BRO). As mentioned above, there are five type of rules. For each type of rules, we create a set of axioms in the BRO. We also introduce some transitive rules for allowing the reasoner to reason on BRO and BPO to detect the potential semantic error automatically.

Integrity Rule.

The integrity rule have the same meaning with a constraint in the relational database. In table 1, we define the cardinality rules. It will be translated into a set of cardinality axiom inside BRO.

Example 1:

The rule "Someone that owns at least 2 cars" is represented by an OWL 2 syntax as follow:

ObjectMinCardinality(2 hasCar SomeOne)

Derivation Rule. This kind of rule allows the system to deduce a new knowledge. If a set of facts satisfies the derivation rule, the reasoner will deduce a new fact from the existing facts. We use SWRL rule to represent this kind of rule.

Example 2:

A rule "Platinum customers receive a 5% discount" is represented by a SWRL rule as follow:

PlatinumCustomer(w)->hasDiscount(x,5)

Reaction Rule.

One of the important rule is the reaction rule which allows the user to define the relationship between a set of actions in a specific domain. We propose six kind of relationships between the tasks: dependency, parallel execution, choice execution, sequential exlusion, parallel exclusion and choice exclusion. These relationships are represented by an OWL object property.

Example 3:

A rule "Task A is depended on task B" is represented by an object property as follow:

hasDependencyTask(B,A)

Production Rule and Transformation Rule. This kind of business rules is represented in the form "IF something DO something". For representing this form with OWL language, we use SWRL [30].

Example 4: a rule "if a customer is a VIP member, they do not have to provide more information" is represented by a SWRL rule as follow:

isVIPMember(cus)-> not provideInfor(cus) Transitive Rule.
We define a set of transitive rules inside the business rule ontology. This kind of error occurs after a sequence of activities. At design-time, when a user designs a business process instance, the user's defined business process will be combined with the set of transitive rules for detecting the potential semantic error when they use these chosen activities in a business process. At run-time, if there are some changes on the business process instance, transitive rule will be used for detecting the error which can occur with this change. Each rule is represented by a set of axioms and SWRL rule inside BRO. In order to check the compliance of business process with business rules, we merge BRO and BPO into one ontology; two concepts Transition in BPO and Task in BRO are defined as two equivalence concepts. The business term individuals can be used as a color and a token in CPN graph (business process). During the execution of business process, the value of individual can be changed but the change must respect the constraints inside BRO (TBox and Properties).

At design time, when a user defines a business process, the business term will be used to name an item. Each transition individual in BPO is equivalent to an action individual in BRO. Depending on the user's given order, BPO editor will generate a set axioms inside BPO.

For example: there are two tasks inside BRO, which is defined that b depends on a as follow:

ObjectPropertyAssertion(:hasDependencyTask ∶b∶a)
It means that a must be executed before b, but at design time a user define that a executes after b, and the rule is generated as follow:

ObjectPropertyAssertion(:hasDependencyTask ∶a∶b) Two rules above are opposite, so the merged ontology of BRO and BPO will be inconsistent. It can be checked by a reasoner (Pellet, Hermit). Because the property hasDependencyTask is defined as a TransitiveObjecProperty in the reaction rule 1 in table 3.

At runtime, we use the same approach to verify the consistency of merged ontology. If a user modifies a business process, the modification will be generated and insert into BPO; for each modification, the reasoner will check the consistency of merged ontology and notify the result to the user automatically.

Conclusion

In this paper, an ontology-based approach for business process compliance checking is proposed. It takes important features of the ontology which are the reasoning capabilities, the possibility to express complex actions, and its declarative semantics to validate not only the consistency of business rules and business process but also the compliance of business process with a set of business rules. The advantage of this approach is to allow the system to detect the semantic flaws of business process automatically at design time and run-time. Nevertheless, by using this approach, if BRO has many concepts and properties, the reasoning may take long time for checking the consistency of BPO and BRO ontology. According to that, future theoretic works involve three main issues. The first one is to focus on the distributed reasoning. The second one will be achieved by selecting the related rule of an action for the validation. And the last goal is to consider the business process execution and work with a data source.

Fig. 3 . 3 ,

 33 Fig. 3. Integration of two ontologies

Table 1 . Transitive Rule 4.3 Checking the Business Process Compliance Using Reasoning

 1