Bruno Lombard
email: lombard@lma.cnrs-mrs.fr.

Denis Matignon
email: denis.matignon@isae.fr.

Yann Le

Yann Le Gorrec
email: yann.le.gorrec@ens2m.fr.

A fractional Burgers equation arising in nonlinear acoustics: theory and numerics ⋆

Keywords: fractional derivatives, Burgers equation, diffusive representation, numerical methods

teaching and research institutions in France or abroad, or from public or private research centers.

A fractional Burgers equation arising in nonlinear

acoustics: theory and numerics.

INTRODUCTION

The motivation of this work comes from musical acoustics, the underlying idea being to model wind instruments from the brass family, well-known for their brassy effects, which is a typically nonlinear phenomenon. Moreover, shock waves can be developped after some critical time, depending on the amplitude of the input or of the initial condition.

Such purely nonlinear effect can be obtained by means of solution to nonlinear PDEs, such as the inviscid Burgers equation: this equation has been widely studied, both from theoretical aspect (hyperbolic equation, conservation law) and numerical aspect (conservative schemes, such as Godunov). Now, a specific feature of wind instruments is the fact that viscous and thermal effects do play a role at the lateral walls of the duct, and contrarily to the open space case, the bounded space imposes a very specific behaviour, see [START_REF] Bruneau | General formulation of the dispersion equation in bounded visco-thermal fluid, and application to some simple geometries[END_REF] and Polak [1991]: evolution as the square root of the frequency, with much bigger physical coefficients, which translate into non negligible effects; to put it shortly, the presence of the boundary layer can almost be listened to! From a modelling point of view, the classical wave PDE is modified with a fractional derivative in time, this behaves as a damping term with long-memory decay, see e.g. [START_REF] Lokshin | Wave equation with singular retarded time[END_REF], [START_REF] Lokshin | Fundamental solutions of the wave equation with retarded time[END_REF].

⋆ Part of this work has been done within the context of the French National Research Agency sponsored project HAMECMOPSYS.

Further information is available at http://www.hamecmopsys.ens2m.fr/.

The one-way Burgers-Lokshin (BL) equation is the simplest model, that combines both these features (nonlinearity and fractional derivative), it has the following form:

∂ t u + c ∂ x u + ε ∂ α t u + b ∂ x (u 2 /2) = 0 for t > 0 , (1
) with compactly supported initial datum u(t = 0, x) = u 0 (x) at t = 0. The coefficients are c > 0 the sound speed, and ε > 0 which takes into account the specific length of both viscous and thermal effects and the radius of the duct; moreover the fractional order α ∈ (0, 1) is α = 1 2 ; b ≥ 0, or Burgers coefficient, quantifies the nonlinear effects.

Accounting for the existing literature on this sharp subject is not so easy, but we should mention first [START_REF] Makarov | Nonlinear and thermoviscous phenomena in acoustics[END_REF] for a general overview, then [START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF] for a case very close to ours, and finally [START_REF] Hélie | Volterra series for solving weakly nonlinear partial differential equations: Application to the Burgers equation with visco-thermal losses[END_REF], [START_REF] Hélie | Simulation of the weakly nonlinear propagation in a straight pipe: application to a realtime brassy audio effect[END_REF] for explicit solutions of this latter case, using Volterra series. The paper is organized as follows: in § 2, strong solutions and weak solutions are defined, and a useful energy balance is proved. Then in § 3, numerical methods are presented, that allow for efficient simulation of the oneway Burgers-Lokshin equation, taking advantage of the two different underlying structures: nonlinear hyperbolic equation, together with a linear diffusive system; in § 4, the results of numerical experiments are presented and carefully explained. Finally in § 5 a conclusion is drawn on this work, and many perspectives and open questions are proposed: one concerns the proof of existence of weak solutions to the original model, e.g. in the case of discontinuous inital data, or after a shock has appeared, i.e. when no more strong solutions can be expected. In order to be self-contained, Appendix A recalls definitions of fractional integrals and derivatives, together with their diffusive representations.

THEORETICAL RESULTS

Strong solutions

Let us first recall the following theorem, proved in [START_REF] Haddar | Well-posedness of nonlinear conservative systems when coupled with diffusive systems[END_REF]: Theorem 1. Let H be a Hilbert space, and A : D(A) ⊂ H → H be a maximal monotone operator, and F some non-linear function F : H → H, such that the semi-linear evolution problem

∂ t X + A X = F (X), t > 0, (2) with X(0) = X 0 ∈ D(A) is well-posed, for t ∈ [0, T max ),
meaning the existence and uniqueness of a strong solution

X ∈ C 1 ([0, T max ); H) ∩ C 0 ([0, T max ); D(A)).
Then, for any two positive pseudo-differential operators (PDOs) of diffusive type, such as the fractional integral I β in § A.1, or as the fractional derivative D α in § A.2, the pseudo-differential non-linear problem:

∂ t X + I β X + D α X + A X = F (X), t > 0, (3) with X(0) = X 0 ∈ D(A), has also a unique strong solution X ∈ C 1 ([0, T ′ max ); H) ∩ C 0 ([0, T ′ max ); D(A)) with T ′ max ≥ T max .
Remark 1. For instance, the well-posedness of problem ( 2) is garanteed in the case when F is locally Lipschitz on H, see [Pazy, 1983, Theorem 1.4

in ch. 6]. Remark 2. Even though u → F (u) = u 2 /2 is locally
Lipschitz, care must be taken not to make use the above theorem too rapidly: unfortunately, equation ( 1) does not fit in this framework, since the problem is nonlinear both in u and ∂ x u. Some adequate formulation is still to be found in order to give existence and uniqueness results on strong solutions to equation (1).

Weak solutions

Moreover, since it is known that, even with C ∞ initial datum u 0 , a shock (i.e. a discontinuity) may occur for t > 0 (typically at t = T max or t = T ′ max ), it is necessary to define generalized or weak solutions to (1).

Definition 1. A weak solution to (1) is a function u ∈ L 1 loc (R + × R), such that u 2 ∈ L 1 loc (R + × R), and ∀φ ∈ C 1 c ((-1, +∞) × R), the following identity holds: - R u 0 (x)φ(0, x) dx = ∞ 0 R (u ∂ t φ + cu ∂ x φ + b (u 2 /2)∂ x φ -ε u. ȟ1-α ⋆ ∂ t φ) dx dt .
where ȟ(t) = h(-t), the convolution on ∂ t φ becoming anticausal in the weak formulation. Proposition 1. Any strong solution is necessarily a weak solution (the proof goes classically by integration by parts), whereas the converse is true only under regularity assumptions on the solution. Remark 3. Existence of weak solution is still an open question. Also finding compatibility conditions, such as Rankine-Hugoniot conditions, is an open question.

Energy balance

Making use of appendix A.2, an energy density can be defined as follows: e ϕ (t, x) := 1 2 ∞ 0 ξ µ 1α (ξ) ϕ(ξ, t, x) 2 dξ. Then, with candidate energy function for equation (1)

E AB (t) := B A 1 2 u(t, x) 2 + ε e ϕ (t, x) dx ,
the following energy balance can be proved:

d dt E AB (t) = - B A ε ∞ 0 µ 1α (ξ) (∂ t ϕ(ξ, t, x)) 2 dξ dx + c u 2 (t, x) 2 + b u 3 (t, x) 6 B A .
Proposition 2. When u 0 is compactly supported and smooth enough, thanks to finite propagation speed, it can be inferred that equation ( 1) is dissipative on (A, B) = R.

NUMERICAL METHOD

Even though the theory for this equation seems not to be sufficiently mature so far, a numerical strategy has been adopted to solve (1), as follows:

(1) use a conservative numerical scheme for hyperbolic conservation laws, like Godunov, Total Variation Diminishing (TVD) or Weighted Essentially Non-Oscillatory (WENO), (2) use Diffusive Representation (DR), that help transform fractional derivatives and integral into a system of linear ODEs, (3) use splitting, before coupling with du dt ; this system is solved exactly in order to avoid numerical instabilities.

A semi-discretized version for simulation

For 0 < α < 1, using appendix A for classical diffusive representations (see e.g. [START_REF] Hélie | Representations with poles and cuts for the time-domain simulation of fractional systems and irrational transfer functions[END_REF]), or following [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputotype fractional derivatives[END_REF], upon setting ξ = θ 2 , another exact diffusive representation of ∂ α t u = I 1-α (∂ t u) can be given by:

∂ α t u(t) = +∞ 0 φ(θ, t) dθ d dt φ(θ, t) = -θ 2 φ + γ α θ 2α-1 ∂u ∂t . φ(θ, 0) = 0, (4)
where we have set

γ α := 2 sin πα π . Then, the continuum (4) is discretized into                d α t u ≈ N ℓ=1 µ ℓ φ(θ ℓ , t) ≡ N ℓ=1 µ ℓ φ ℓ , dφ j dt = -θ 2 j φ j + γ α θ 2α-1 j ∂u ∂t , j = 1, • • • N, φ j (0) = 0.
(5)

The number N of relaxation mechanisms plays a role both on the accuracy of the approximation and also on the computing cost. The determination of weights µ ℓ and quadrature nodes θ ℓ is discussed further. Combining (1) and ( 5) yields the first-order system without fractional derivatives

                   ∂u ∂t + c ∂u ∂x + b ∂ ∂x u 2 2 = -ε L ℓ=1 µ ℓ φ ℓ , ∂φ j ∂t + c α θ 2α-1 j c ∂u ∂x + b ∂ ∂x u 2 2 = -θ 2 j φ j -c α θ 2α-1 j ε L ℓ=1 µ ℓ φ ℓ . (6)
Taking the vector of (N + 1) unknowns

U = ( u φ 1 . . . φ N )
T , (7) the system (6) can be put in the form

∂ ∂t U + ∂ ∂x F(U) = S U, (8)
where

F = (F 1 F 2 . . . F N +1
) T is the nonlinear flux function

F 1 = c u + b u 2 2 , F j = γ α θ 2α-1 j F 1 , j = 2 • • • N + 1, ( 9
) and S is the (N + 1) × (N + 1) diffusive matrix

S = -          0 ε µ 1 • • • ε µ N 0 θ 2 1 + Ω µ 1 • • • Ω µ N . . . . . . . . . . . . 0 Ω µ 1 • • • θ 2 N + Ω µ N          , (10)
with Ω = γ α ε. The eigenvalues of the Jacobian matrix J = F ′ in (9) are real: c + b u, and 0 with multiplicity 2 N . Remark 4. These eigenvalues do not depend on the coefficients of the diffusive representation.

Numerical Modelling

Numerical Scheme In order to integrate the system (8), one introduces a uniform mesh size ∆x and a variable time step ∆t. The approximation of the exact solution U(x j = j ∆ x, t n = t n-1 + ∆t) is denoted by U n j . Unsplit integration of (8) is not optimal, because the time step stability condition typically implies

∆t ≤ min ∆x c max , 2 ̺(S) , (11)
where c max := c + b max j (u n j ) is the maximum numerical velocity at time t n , and S is the spectral radius of S which grows in an unbounded way with N , thus penalizing the standard CFL condition. Moreover, it requires to build an adequate scheme for the coupled system.

A more efficient strategy is adopted here, which consists in splitting the original system (8) into a propagative part (12a) and a diffusive part (12b)

       ∂ ∂t U + ∂ ∂x F(U) = 0, (12a)
∂ ∂t U = S U. (12b)
The discrete operators associated with steps (12a) and (12b) are denoted by H a and H b , respectively. The secondorder Strang splitting (see e.g. [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF], [START_REF] Holden | Operator splitting for the KDV equation[END_REF]) is then used between t n ant t n+1 , leading to the time-marching scheme:

• U (1) j = H b ( ∆t 2 ) U n j , • U (2) j = H a (∆t) U (1) j , • U n+1 j = H b ( ∆t 2 ) U
(2) j .

(13)

The propagative part (12a) is solved by any standard scheme for nonlinear hyperbolic PDE:

u n+1 j = u n j - ∆t ∆x F 1 j+1/2 -F 1 j-1/2 , φ n+1 j,ℓ = φ n j,ℓ -γ α θ 2α-1 ℓ ∆t ∆x F 1 j+1/2 -F 1 j-1/2 , ( 14
)
where F 1 j±1/2 is the numerical flux function of the advection-Burgers part in (9). In practice, a second-order TVD scheme with MC-limiter is used in our numerical experiments, see LeVeque [2002]. Stability analysis of ( 14) yields the optimal CFL condition

Υ = c max ∆t ∆x ≤ 1. ( 15
)
Since the physical parameters do not vary with time, the diffusive part (12b) can be solved exactly. This gives

H b ∆ t 2 U j = e S ∆ t 2 U j . (16)
It is computed numerically using a (6, 6) Padé approximation in the "scaling and squaring method", see [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[END_REF]. If the physical parameters are constant in space, the computation is done only once at each time step, leading to a negligible computational cost. This part of the splitting is unconditionally stable, so that the global stability requirement is (15) and is not penalized by the diffusive part. In other words, the time step only depends on the advection and Burgers coefficients in (1). In particular, ∆t does not depend on the coefficients of the diffusive representation.

Coefficients of the Diffusive Representation It remains to determine the 2 N coefficients of the diffusive representation µ ℓ and θ ℓ in (10). This issue is crucial both for the accuracy of the modeling and for the computational efficiency of the method, see [START_REF] Hélie | Representations with poles and cuts for the time-domain simulation of fractional systems and irrational transfer functions[END_REF], and [START_REF] Deü | Simulation of fractionally damped mechanical systems by means of a Newmarkdiffusive scheme[END_REF]. Two strategies exist for this purpose.

The first one relies on Gaussian quadratures of improper integrals (4). A Laguerre quadrature formula can be used for this purpose [START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF]. This technique is frequently used in the literature, see e.g. [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputotype fractional derivatives[END_REF], [START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF], ensuring µ ℓ > 0 and θ ℓ > 0. We have tried this approach, but the convergence rate is very poor. As can be shown in the limit-case of the linear regime (b = 0), a large number of memory variables is required to get acceptable accuracy, which greatly increases the computational cost, see Appendix B.

When a characteristic length of the waves exists, a more efficient strategy can be proposed based on the Fourier relation (B.2). The original problem (1) and the first-order system (8) differ only in their symbol χ(ω): meaningful, and has proven its efficiency in a previous work about poroelastic waves, see [START_REF] Blanc | Biot-JKD model: simulation of 1D transient poroelastic waves with fractional derivatives[END_REF].

(B.

NUMERICAL EXPERIMENTS

A 20 m long domain is discretized on 1000 grid nodes. The celerity of advection is c = 300 m/s, and the nonlinear coefficient is b = 1 in (1). The CFL number is Υ = 0.95. One takes N = 6 memory variables only in (7). The computations are initialized by a Gaussian pulse

u 0 (x) = u m exp - x -x 0 σ 2 , ( 17
)
with amplitude of velocity u m = 100 m/s, spatial shift x 0 = 1 m, wavelength λ = 1, and standard deviation σ = λ 2 √ ln 100 . In the inviscid case, the smooth Gaussian pulse develops a shock in finite time, yielding a decay of energy. The initial data give the time break

t * = e 2 σ b u m = 0.021 s. ( 18
)
Simulations are performed up to t = 0.04 s > t * , which corresponds to 769 steps in the inviscid case.

Figure 1 displays the simulations with the Gaussian pulse.

In the inviscid case ε = 0 (a), one shows the initial and the final solutions. Classically, the right part of the pulse has developped a shock. Dissipative models are examined in subfigures (b-e), for various values of ε and α in (1); the numerical value of the inviscid case is shown in red line for comparison. As ε and α increase, the amplitude of waves decreases, and the wave slows down. The shock is also smeared. For ε = 10 and α = 0.5, the shock totally disappears. In subfigure (f), one shows the time evolution of the energy E n := j (u n j ) 2 , and the time break ( 18) is denoted by a vertical dotted line. In the inviscid case, the energy E n is constant up to t * , and then it decreases. On the contrary, dissipation occurs at all times when ε = 0.

Another numerical experiment, made with a door pulse as initial condition, would show the following: in the inviscid case, the door pulse splits into two parts: a rarefaction fan and a right-going shock, which decreases the energy. These two waves interact at a specific time. Observations similar to the Gaussian case could be done: decrease of amplitude and slowdown of waves when ε and α increase, smearing of shock.

CONCLUSION AND PERSPECTIVES

In this work, the following questions have been addressed: a fractional nonlinear PDE arising in musical acoustics has been studied, weak solutions have been properly defined, an energy balance has been proved. A numerical method to solve advection-Burgers equation with time fractional derivative has been proposed. Three parts are involved: a classical TVD scheme for nonlinear hyperbolic equations, a diffusive representation of the fractional derivatives to avoid storing past values of the solution, and a splitting to couple the previous tools in an efficient manner, see e.g. [START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube with an array of Helmholtz resonators[END_REF].

So far, many perspectives are left open:

• find an appropriate framework for existence and uniqueness of strong solutions. • prove existence of weak solutions to the original model (in the sense of definition 1); it should be useful in at least two situations: in the case of discontinuous inital data, or after a shock has appeared, i.e. when no more strong solutions can be expected.

• recast this coupled model in the setting of infinitedimensional port-Hamiltonian systems (pHs) including dissipation. Based on the energy balance derived in § 2.3, and following e.g. [START_REF] Maschke | On alternative Poisson brackets for fluid mechanical systems and their extension to Stokes-Dirac structures[END_REF], one can choose as Hamiltonian functional H 0 (u) := 1 0 (c u 2 2 + b u 3 6 ) dz for the Burgers part; and using e.g. [START_REF] Gorrec | Diffusive systems coupled to an oscillator: a Hamiltonian formulation[END_REF], one can choose as Hamiltonian functional H Φ := 1 2 ∞ 0 µ 1-α (ξ) ξ ϕ(ξ) 2 dξ for the diffusive part. So far, the port-Hamiltonian formulation of the infinitedimensional coupled system does not prove straightforward.

• ensure positivity of weights µ ℓ in the linear optimization procedure, since it is crucial for well-possedness of the diffusive model (6). To ensure µ ℓ > 0, we will implement optimization with constraint based on Shor's algorithm, see [START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF]. • perform a rigorous analysis of the shock zones, based e.g. on a matched asymptotic analysis as in [START_REF] Sugimoto | Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[END_REF], to predict the width of the boundary layer. Indeed, depending on the parameters of the fractional terms, shocks are smeared or may even disappear, as shown by the numerical simulations

  Fig. 1. Gaussian pulse. Snapshots of u in the inviscid case and for various dissipative parameters ε and α (a-e). The subfigure (f) shows the time evolution of the energy; the vertical dotted line denotes the time of break t * (18).

ACKNOWLEDGEMENTS

The authors would like to thank the three anonymous reviewers, who carefully read the original submission, made useful comments, and greatly helped enhance the final version of the paper; they are gratefully acknowledged.

Appendix A. FRACTIONAL DERIVATIVES AND DIFFUSIVE REPRESENTATIONS

A.1 Fractional integral

Let u ∈ L 2 (0, T ), for any β ∈ (0, 1), define the causal convolution kernel h β (t) := 1 Γ(β) t β-1 for t > 0, then the fractional integral of order β ∈ (0, 1) of u is defined by

, the fractional integral can be reformulated by the following input-output representation :

with e ξ (t) := e -ξ t , and [e ξ ⋆ u](t) = t 0 e -ξ (t-τ ) u(τ ) dτ . The following infinite-dimensional dynamical system can be seen as a state-space realization of the fractional intergral of order β: ∂ t ϕ(ξ, t) = -ξ ϕ(ξ, t) + u(t), ϕ(ξ, 0) = 0 , (A.1)

(A.2)

A.2 Fractional derivative

For short, fractional derivative of order α ∈ (0, 1) of u in a generalized sense is defined by

and a careful computation shows that the following inputoutput representation holds:

The following infinite-dimensional dynamical system can be seen as a state-space realization of the fractional derivative of order α:

FOURIER ANALYSIS OF THE LINEAR MODEL

B.1 Full model

The Fourier transforms in time and space are denoted by

where ω is the angular frequency and k is the wavenumber. Applying these transforms to the fractional equation (1) provides the nonlinear equation

with the fractional derivative coefficient

3) When c = 0 and b = 0, one obtains the dispersion relation

It follows the phase velocity υ := ω / ℜe(k) and the attenuation β a := -ℑm(k):

(B.5) One deduces the elementary properties:


Michael J Walock

Yujiao Zou

Denis Lagadrillere

Corinne Nouveau

Andrei V Stanishevsky

Cr-W-C-N nanocomposite coatings via reactive magnetron sputtering

While binary tungsten carbide can form smooth, hard films, these suffer from low fracture toughness. Tungsten nitride films are frequently harder, but are more brittle. Chromium nitride has excellent wear and oxidation resistance, but films often form with low hardness. Composites of these binary compounds offer a possibility to tailor the material for a desired combination of properties. To this end, we have used reactive RF-magnetron sputtering with Cr and WC targets to form quaternary composites, with nitrogen as the reactive gas. The coatings were deposited on Si, Ti, and steel substrates. The nitrogen partial pressure was varied to investigate the relationship between the film properties and the deposition conditions. Energy dispersive spectroscopy showed changes in the chemical composition as a result of the change in nitrogen partial pressure. X-ray diffraction illuminated the structure as either a solid solution with a B1 NaCl structure, or a nanocomposite with the average crystallite size under 11 nm. Optical interferometer revealed low compressive stresses. And nanoindentation established that the films are hard and adherent.

I. INTRODUCTION

Tungsten carbide tools have played a prominent role since the beginning of the last century. In part, this is due to the material's high melting point, high hardness, and high oxidation resistance [START_REF] Warren | Oxidation of tungsten and tungsten carbide in dry and humid atmospheres[END_REF]. However, it's most stable form, h-WC, tends to be brittle and is easily damaged [START_REF] Li | Theoretical study on the stability, elasticity, hardness and electronic structures of WC binary compounds[END_REF]. To counter this, a composite of metal (usually, cobalt) and tungsten carbide was formulated in the 1920s. This cemented-WC has excellent hardness, strength, and fracture toughness. Unfortunately, this improvement comes at the price of thermal stability. At temperatures as low as 600ºC, cemented-WC begins to oxidize and its strength begins to degrade [START_REF] Acchar | Strength Degradation of a Tungsten Carbide-Cobalt Composite at Elevated Temperatures[END_REF]. In addition, the prominent wear mechanism seems to be preferential removal of the metal binder [START_REF] Blombery | Abrasive wear of tungsten carbide-cobalt composites. I. Wear mechanisms[END_REF][START_REF] Bailey | Wear of some cemented tungsten carbide tools in machining oak[END_REF][START_REF] Gant | Wear of tungsten carbide-cobalt hardmetals and hot isostatically pressed high speed steels under dry abrasive conditions[END_REF]. Today's tools are typically coated with titanium-based compounds to improve the wear and oxidation resistance.

With its excellent mechanical, chemical, physical, and electronic properties, tungsten nitride has found use in a wide variety of applications. It has proven to be an effective diffusion barrier in microelectronics , preventing the formation of copper silicides up to 600ºC [START_REF] Becker | Diffusion barrier properties of tungsten nitride films grown by atomic layer deposition from[END_REF]. Tungsten nitride has shown to be useful in the catalytic conversion of nitrous oxide [START_REF] Shi | Catalytic activities of tungsten nitride for NO dissociation and reduction with hydrogen[END_REF]. As a protective coating, tungsten nitride has excellent wear properties up to 300ºC [START_REF] Polcar | Structural and tribological characterization of tungsten nitride coatings at elevated temperature[END_REF]. Unfortunately, this limits the use of tungsten nitride in mechanical applications.

The greatest advantages of CrN is its chemical and oxidation resistance. It is a very good barrier to corrosion [START_REF] Lin | Corrosion of CrN and CrN/TiN coated heat-resistant steels in molten A356 aluminum alloy[END_REF][START_REF] Paulauskas | Corrosion behavior of CrN, Cr2N and [pi] phase surfaces on nitrided Ni-50Cr for proton exchange membrane fuel cell bipolar plates[END_REF]. In addition, CrN prevents oxidation up to 700ºC. However, these films have the drawback of forming porous, columnar structures [START_REF] Hurkmans | Chromium nitride coatings grown by unbalanced magnetron (UBM) and combined arc/unbalanced magnetron (ABS (TM)) deposition techniques[END_REF][START_REF] Bertrand | Properties of reactively RF magnetron-sputtered chromium nitride coatings[END_REF][START_REF] Han | The cavitation-erosion phenomenon of chromium nitride coatings deposited using cathodic arc plasma deposition on steel[END_REF][START_REF] Hones | Structural and mechanical properties of chromium nitride, molybdenum nitride, and tungsten nitride thin films[END_REF]. This tends to lower Paper Number: T-2 their apparent hardness, and reduces their corrosion and oxidation resistance. One solution is the application of a metal buffer layer between the substrate and the CrN film. Both Cr and Ni interlayers have shown improved oxidation and corrosion resistance over single CrN films [START_REF] Brandl | Corrosion behaviour of hybrid coatings[END_REF][START_REF] Creus | Improvement of the corrosion resistance of CrN coated steel by an interlayer[END_REF][START_REF] Han | Corrosion and tribological studies of chromium nitride coated on steel with an interlayer of electroplated chromium[END_REF].

There are clear limits to the use of binary compounds in modern applications. One approach to improving these compounds is the addition of a third material. For example, adding aluminum to chromium nitride, forming CrAlN, improves the hardness, oxidation resistance, and corrosion resistance [START_REF] Barshilia | A comparative study of reactive direct current magnetron sputtered CrAlN and CrN coatings[END_REF]. Superlattices of binary compounds may also improve performance. For example, CrN/WN superlattices show improved hardness over single layer CrN or WN films [START_REF] Wu | Microstructure evaluation and mechanical properties of nanolayered chromium nitride/tungsten nitride coating[END_REF].

However, the most successful improvements seem to be from the formation of nanocomposites [START_REF] Musil | Hard and superhard nanocomposite coatings[END_REF]. A true nanocomposite consists of at least two distinct phases, with an average crystallite size less than 10 nm. This can lead to some extraordinary improvements without the disadvantages of superlattice structures. For example, nanocrystalline (nc-) TiN combined with an amorphous (a-) Si3N4 matrix has reached hardness greater than diamond [START_REF] Veprek | Composition, nanostructure and origin of the ultrahardness in nc-TiN/a-Si3N4/a-and nc-TiSi2 nanocomposites with HV = 80 to ≥ 105 GPa[END_REF][START_REF] Veprek | A concept for the design of novel superhard coatings[END_REF][START_REF] Brazhkin | Harder than diamond: dreams and reality[END_REF][START_REF] Veprek | Different approaches to superhard coatings and nanocomposites[END_REF][START_REF] Veprek | Design of ultrahard materials: Go nano![END_REF]. In addition to improving hardness, toughness can also be significantly improved over traditional materials. For example, nc-TiC/a-C [START_REF] Voevodin | Nanocrystalline carbide/amorphous carbon composites[END_REF] and nc-WC/a-C [START_REF] Voevodin | Tribological performance and tribochemistry of nanocrystalline WC/amorphous diamond-like carbon composites[END_REF] form hard, but supertough nanocomposites.

While the research of nc-TmN/a-nitride (such as the nc-TiN/a-Si3N4) and nc-TmC/a-C (such as nc-TiC/a-C) nanocomposites has become well established, less information is available on other quaternary systems. For example, the Cr-W-C-N system may provide a separate path to improved hardness, toughness, and oxidation-resistance. This study will report one aspect of the investigation into the Cr-W-C-N system: the effect of changing the N2 partial pressure.

II. EXPERIMENTAL A. Sample Preparation

Samples were prepared by RF magnetron sputtering in a modified commercial sputtering system (Nordiko 3500). This system consists of two 4 inch magnetron sputtering guns powered by individual 1.25 kW RF generators (Nordiko SG-1250), operating at 13.56 MHz. The guns are arranged confocally, at ± 45º with respect to the sample normal. Substrate to target distance is 8 cm. The two targets were Cr (Neyco, 99.95 %) and WC (Neyco, 99.5 %). The base pressure is 6×10 -5 Pa. For improved adhesion, the substrates were cleaned and a CrWC buffer layer was deposited. Cleaning was provided by a dc pulsed power supply operating at 12 kV (25 µs on, 50 ms off) for 5 minutes, in an Ar atmosphere at 2 Pa. The CrWC layer was deposited with the Cr and WC targets at -400 V and -600 V, respectively, and an Ar pressure of 0.6 Pa. The target thickness for this layer was 50 nm. For the CrWCN layer, the total pressure was maintained at 1 Pa; the Ar and N2 pressures were adjusted to obtain Ar:N2 ratios of 80:20, 60:40, and 40:60. The Cr and WC targets were held at -600 V and -600 V, respectively. The deposition time was varied to obtain a thickness of 2 µm for the CrWCN layer. The substrates were grounded and the substrate temperature was less than 200ºC during deposition.

B. Characterization Paper Number: T-2

Samples deposited on Si (100), Ti, and carbon steel (AFNOR XC100) were subjected to various characterization techniques. The surface morphology was investigated with scanning electron microscopy (SEM, JEOL JSM-5900LV) and atomic force microscopy (AFM, Veeco Topometrix Explorer).

Chemical composition was determined via energy dispersive spectroscopy (EDS, Oxford INCA x-act). Secondary electron images were captured at 20 kV accelerating voltage, in a chamber at 1×10 -4 Pa. For EDS, the accelerating voltage was 5 kV. Surface roughness was explored (across a 5 x 5 µm 2 area) with the AFM in contact mode, using v-shaped, high resonance frequency SiN cantilevers with a pyramidal tip of 50 nm radius and a force constant of 0.032 N/m. The thickness was measured with an optical profilometer, across an artificial step created on the sample during deposition. The crystal structure was investigated by x-ray diffraction (XRD, Philips X'pert MPD) with a Cu K-α tube (λ = 0.15406 nm), operating at 45 kV, 40 mA. The detector was scanned between 20º and 70º, with a constant take off angle of 5º. Structure was determined by comparing the diffractograms with the ICDD powder diffraction file [START_REF]Powder Diffraction File (Database)[END_REF]. The crystallite size was estimated with the Debye-Scherrer equation.

Stress was determined by the wafer curvature method. A 380 µm thick Si (100) substrate is placed in an optical interferometry, with a sodium lamp (λ = 589.3 nm). Newton's rings are created in each sample. The diameters of the resulting dark interference fringes were measured with image analysis software, and the radius of curvature was determined by linear regression. A modified version of Stoney's equation [START_REF] Janssen | Stress and strain in polycrystalline thin films[END_REF], accounting for the elastic anisotropy of Si (100) substrates, is used to quantify the stress.

Mechanical properties were also investigated. Using the Oliver and Pharr method [START_REF] Oliver | Improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF], hardness and Young's modulus are calculated from the loading-unloading curves of a nanoindenter (MTS Nanoindenter XP). The tool uses a Berkovich tip, and is calibrated with silica samples. The maximum indentation depth was 500 nm. The maximum indentation depth was 500 nm. The values for hardness and Young's modulus were taken at 100 nm indentation depth where the surface effects should be minimized. Although the maximum average hardness is observed at the 100 nm indentation depth, we assume that there is already the substrate effect (soft commercially pure titanium) even at this depth [START_REF] Chen | Numerical study on the measurement of thin film mechanical properties by means of nanoindentation[END_REF]. Hold times of 10 s (at maximum load) and 50 s (at 10 % maximum load, during unloading) were used to minimize thermal drift and creep effects.

III. RESULTS AND DISCUSSION

A. Chemical Composition

The CrWC underlayer was targeted with a 2:1 ratio of Cr:W. From EDS, the actual composition is 31 at % C, 45 at % Cr, 24 at % W; the resulting Cr:W ratio is 1.875:1. Any O2 present is below the detection limits (~ 0.1 at %) of the analysis system. For the CrWCN samples, the atomic concentrations are summarized in Table 1. The target metal ratio for all three samples was 2:1. However, the actual ratios for Cr:W are 1.71, 1.36, and 2.13 for the samples deposited in Ar/N2 ratios of 80/20, 60/40, and 40/60, respectively. There may be several factors at work. For example, CrN formation is more energetically favorable over WCN. In addition, the addition of N2 will more strongly affect the WC sputtering rate. While there is no clear trend for the Cr:W ratio, there is a clear trend with regards to the incorporation of N2 into the coatings. As N2 concentration increased, the films incorporated more N2 into the structure. While assuming most of this is bound as a nitride, the EDS results indicate a significant amount of free N2 may be present. This gas may be trapped in grain boundaries and voids created during deposition.

B. Surface morphology

The surfaces of samples are shown in the secondary electron and AFM images of Fig. 1 (a)-(c). For the most part, the surfaces are relatively clean. However, there are few large spherical drops on each surface. These are more than likely areas of overgrowth due to defects/ contaminates on the Si substrate surface, such as dust particles. Neither the SEM nor SFM images show any significant differences between the films at this scale. In Fig. 2, the roughness is plotted versus the N2 concentration. There is a clear trend to a reduced surface roughness with increased N2 concentration.

C. Structural Analysis

The structures of the three samples was studied with XRD. The peak parameters were determined with the CMPR program [START_REF] Toby | CMPR-a powder diffraction toolkit[END_REF]. Sample A, Fig. 3(a), has peaks at 36.9º, 42.7º and 62.3º representing the (111), (200), and (220) orientations. The preferred orientation is (200), with a lattice parameter of (0.42 ± 0.01) nm. In sample B, Fig. 3(b), the same orientations are present at 36.9º, 42.8º, and 62.5º, respectively, with no observed change in the lattice parameter. However, there is a change in the preferred orientation from (200) to (111). When the nitrogen content is further increased (as in sample C, Fig. 3(c)), this preference remains. In addition, the (111), (200), and (220) peaks shift to higher diffraction angles; this reduces the lattice parameter to (0.417 ± 0.001) nm. Additionally, these patterns are typical of the B1 (NaCl) structure [START_REF] Hones | Mechanical properties of hard chromium tungsten nitride coatings[END_REF][START_REF] Lousa | Multilayered chromium/chromium nitride coatings for use in pressure die-casting[END_REF][START_REF] Yau | Tungsten doped chromium nitride coatings[END_REF].

where λ is the x-ray wavelength, K is a shape factor (assumed to be 0.94, due to predominately cubic crystal shape), βsample is the FWHM (in radians) of the sample peak (preferred orientation), βSi is the instrument broadening contribution, and θ is the preferred peak location (in radians). From fitting, the (200) peak of sample A has a FWHM of 1.15º. The (111) peak of sample B has a FWHM of 0.91º, and the (111) peak of sample C yields a FWHM of 0.97º. The instrumental broadening is 0.31º. Inserting these numbers into the Equation 1, the average crystallite size of sample A is (8.56 ± 0.02) nm, sample B is (10.9 ± 0.2) nm, and sample C is (10.1 ± 0.2) nm. There are differences in the average crystallite size, but the only correlation is oxygen contamination. According to the EDS results, sample B is the most contaminated, followed by sample C and finally, by sample A. There are differences in the average crystallite size, but these variations are small to be statistically important. While the crystallite size follows somewhat the level of oxygen contamination in the film, the role of the impurities in the film growth in not obvious in this case.

D. Stress Measurements

From the linear regression of the measured dark fringe interference rings, the radius of curvature were calculated for each Si sample. This value was placed into a modified Stoney's equation, Equation 2where ts is the substrate thickness, ts is the film thickness, R is the radius of curvature, and the (s11+s12) -1 factor accounts for the elastic anisotropy of Si. From optical profilometry, the coatings were approximately 2.5 µm thick. The estimated stress in the samples is compressive, with values of -(0.14 ± 0.06) GPa, -(0.32 ± 0.13) GPa, and -(0.13 ± 0.05) GPa for samples A, B, and C, respectively. There is no regular pattern to the stress with regards to N2 partial pressure, but there is a correlation with WC content. Sample B has more WC present than either sample A or sample C. Reduced WC and increased Cr content seems to play a role in reducing the average stress of the coating.

E. Nanoindentation

Nanoindentation is a very powerful tool to elicit the mechanical properties of thin film coatings. Specifically, the hardness and Young's modulus of each coating can be extracted from the loading-unloading curves, Fig. 4(a). As shown in Fig. 4(b), the apparent hardness values are 15.4 GPa, 14.2 GPa, and 8.9 GPa for sample A, B, and C, respectively. Similarly, in Fig. 4(c), the calculated Young's moduli are 232 GPa, 221 GPa, and 168 GPa, for samples A, B, and C, respectively. There may be a correlation between the composition and the mechanical properties. While there is no correlation with the N2 partial pressure, a relationship seems to exist between the mechanical properties and WC concentration. Increased amounts of WC appear to increase both the hardness and modulus of the sample.

However, these numbers are lower than a previously reported study of this system [START_REF] Walock | Sputtered tungsten-based ternary and quaternary layers for nanocrystalline diamond deposition[END_REF]. This is probably due to a combination of low stress, microporosity, substrate effects, and oxygen contamination in the current samples. In many materials, high stress states will magnify the elastic response of the material [START_REF] Lu | Recent advances on understanding the origin of superhardness in nanocomposite coatings: A critical review[END_REF]. Another factor is the addition of Cr. Cr tends towards columnar grain growth. This may induce voids and trap N2. For example, ceramic coatings with only 10 % porosity show 2-3x lower hardness than their nonporous counterparts [START_REF] Luo | Porosity-dependence of elastic moduli and hardness of 3Y-TZP ceramics[END_REF][START_REF] Milman | Temperature dependence of hardness in silicon-carbide ceramics with different porosity[END_REF]. And while the hardness and Young's modulus were evaluated at an indentation depth of 100 nm, far less than the standard 10 % of film thickness, there may be substrate effects [START_REF] He | Finite element modeling of indentation into superhard coatings[END_REF]. In addition, there are O2 impurities. While there is no correlation in the spread of the values with the O2 contamination, other systems are known to suffer with this type of impurity [START_REF] Veprek | Degradation of superhard nanocomposites by built-in impurities[END_REF][START_REF] Hao | Superhard nitride-based nanocomposites: role of interfaces and effect of impurities[END_REF].

IV. CONCLUSION

CrWCN samples were successfully prepared with RF magnetron sputtering. From optical profilometry, the coating thicknesses were estimated at 5 µm. From the diffraction data, the samples may be a solid solution with a B1 NaCl fcc structure, or a possible nanocomposite. In either case, the average unit cell is 0.42 nm for samples A and B, and 0.417 nm for sample C. While there is no significant deviation in the lattice parameters, there is variation in the average crystallite size. This appears to correlate with the level of O2 contamination. In addition, the combination of this contamination with low stress, porosity, and possible substrate effects may be factors in the reported hardness values.

However, these films are adherent, with low surface roughness. These properties may lead to a variety of industrial uses. Future work will further evaluate the wear and tribological properties of these films. Future depositions of the CrWCN system will seek to minimize this contamination, and produce optimized coatings for industrial use.

Figure 1 .

 1 Figure 1. Secondary electron and AFM (insets) images of the CrWCN samples produced in an Ar/N2 environment of (a) 80/20, (b) 60/40, and (c) 40/60.

Figure 2 .

 2 Figure 2. Roughness parameter as a function of N2 concentration.

Figure 3 .

 3 Figure 3. X-ray diffractogram for CrWCN samples. Depositions were carried out in an Ar/N2 environment with ratios of (a) 80/20, (b) 60/40, and (c) 40/60. Note the change in preferred orientation from (200) to (111) with decreasing N2 content.

Figure 4 .

 4 Figure 4. (a) The loading-unloading curve of sample C. The (b) hardness and (c) Young's modulus of the CrWCN samples. To avoid substrate effects, the values are taken at an indentation depth of 100 nm (dashed line).

Table 1 .

 1 Summary of the EDS results for CrWCN sample prepared in different gas environments.

		Sample Sample			Composition Composition Composition Composition Composition
	Name	Ar/N2 ratio	Cr	W	C	N	O
	A	80/20	27	15.8	16.9	36.9	3.6

By fitting the various peaks, the FWHM of each peak can be accurately determined. With this information, the the Debye-Scherrer equation is used to estimate the average crystallite size, Equation1Paper Number: T-2

ACKNOWLEDGMENTS

The authors would like to thank Romaric Masset for fabricating the steel substrates. This work was supported by the U.S. National Science Foundation (DMR-0806521) and the Regional Council of Burgundy, France.