Hideo Aochi

Ariane Ducellier

Fabrice Dupros

Mickael Delatre

Thomas Ulrich

Florent De Martin

Masayuki Yoshimi

Finite difference simulations

Introduction

The 2007 Mw 6.6 Niigata-ken Chuetsu-oki earthquake occurred on 16 July slightly off of the northwest coastline of the Japanese mainland (Figure 1). The hypocentral depth is about 10 km and the fault mechanism is reverse, causing no significant surface rupture.

Although modern structures incurred little significant damage, the event did shut down the Kashiwazaki-Kariwa nuclear power plant located above the inferred fault plane, which induced changes in the seismic hazard evaluation criteria applied to nuclear power plants. This earthquake's mechanism is quite complex and remains uncertain despite the dense observation network in operation (cf. [START_REF] Aoi | Source process of the 2007 Niigata-ken Chuetsu-oki earthquake derived from near-fault strong motion data[END_REF]. The main reason for this is that the subsurface structure is so complex that commonly used algorithms based on a 1D structure model are unable to accurately determine aftershock locations [START_REF] Kato | Kanazawa and Groupe for the aftershock observations of the 2007 Niigataken Chuetsu-oki earthquake, Imaging heterogeneous velocity structures and complex aftershock distributions in the source region of the 2007 Niigataken Chuetsu-oki Earthquake by a dense seismic observation[END_REF][START_REF] Kato | Reactivation of ancient rift systems triggers devastating intraplate earthquakes[END_REF][START_REF] Shinohara | Precise aftershock distribution of the 2007 Chuetsu-oki earthquake obtained by using an ocean bottom seismometer network[END_REF]. Immediately after the mainshock, high acceleration levels measured on seismic records along the coast tended to indicate that the fault orientation was northwest dipping, despite the fact that the geological structure would indicate a southeast dipping fault.

However on the basis of aftershock observations obtained from arrays deployed for one month after the mainshock, the interpretation that predominantly emerges is that the major part of faulting is southeast dipping [START_REF] Kato | Reactivation of ancient rift systems triggers devastating intraplate earthquakes[END_REF]. Multisegment models are inferred from the InSAR observations [START_REF] Aoki | Coseismic deformation due to the 2007 Chuetsu-oki earthquake (Mw=6.8)[END_REF][START_REF] Nishimura | Crustal deformation and a preliminary fault model of the 2007 Chuetsu-oki earthquake observed by GPS, InSAR, and levelling[END_REF], and the possibility of rupture transfer between multisegment models is also dynamically simulated [START_REF] Aochi | Dynamic rupture of crosscutting faults: A possible rupture process for the 2007 Mw 6.6 Niigata-ken Chuetsu-Oki earthquake[END_REF]. Thus understanding this earthquake is a seismologically challenging problem, and it is worth ascertaining how well we can reproduce the ground motion using the known seismological information for this earthquake.

The basic concept behind the seismic hazard evaluation is to be able to predict ground motion under a given situation (cf. [START_REF] Douglas | A survey of techniques for predicting earthquake ground motions for engineering purposes[END_REF]. As the time series analyses are increasingly called on in studying nonlinear soil-structure interaction in engineering seismology, the need is ever greater to provide the (input) ground motion quantitatively from seismological observations. Recent progress in numerical simulation and computational resources enables us to theoretically simulate ground motions at high frequencies up to several Hz by taking into account a complex source description and heterogeneous material parameters. However the frequency limit still stands at around 0.5 -1 Hz when comparing with the observed data. Well-known, successful examples can be cited for the 1995 Hyogoken Nanbu (Kobe), Japan, earthquake or the 2004 Parkfield, California, earthquake (cf. [START_REF] Pitarka | Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogoken-nanbu, Japan, earthquake[END_REF][START_REF] Sasetyan | Caractérisation du spectre de réponse en champ proche d'une faille par simulation des mouvements forts du sol[END_REF]. Both examples feature very shallow strike-slip faulting, relatively better studied than the fault mechanism for the 2007 Chuetsu-Oki earthquake we are dealing with in this paper. In any case, we believe it important to demonstrate how seismological insights serve to reproduce the ground motions in such difficult and complex cases.

In this paper, to study the validity of 3D subsurface structure models, we begin with simulations of the aftershocks, considered as point sources. We thence move on to the mainshock simulation using different finite source models provided by some seismological inversion studies. We compare the simulated ground motions with the records obtained on the permanent networks in terms of the waveforms and response spectra. In this paper, our aim is not to tune up the parameters to obtain the best model, but rather to analyze each model in detail and discuss how we can procure reliable input ground motion in the near field.

Numerical Simulation Methodology

All the simulations in this paper are carried out using the finite difference method based on the staggered grid with 4 th order in space and 2 nd order in time [START_REF] Aochi | The 1999 Izmit, Turkey, earthquake: Nonplanar fault structure, dynamic rupture process, and strong ground motion[END_REF], Dupros et al., 2008, references therein). The finite fault is approximated by a series of point sources in space (Olsen, 1994;Graves, 1996) and an arbitrary source time function can be considered at each point. How material heterogeneity dealt with is considered a crucial issue in the finite difference scheme (cf. [START_REF] Moczo | 3D heterogeneous staggered-grid finite-difference modelling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities[END_REF][START_REF] Sasetyan | Caractérisation du spectre de réponse en champ proche d'une faille par simulation des mouvements forts du sol[END_REF]. In this study, we follow Graves (1996) by averaging the material parameters, since the material interface is not always precisely defined at the FD grids we used.

We define our physical model volume as 110 km (EW) x 120 km (NS) x 30 km (depth). We use a grid spacing of ∆s = 200 m and a time step of ∆t = 0.01 sec for most of the simulations; in some cases, we calculate with ∆s = 100 m and ∆t = 0.005 sec, and ∆s = 75 m and ∆t = 0.0035 sec are used to check the convergence of simulations. For all cases, ten finite difference grids are added at each edge as a Perfectly Matched Layer absorbing boundary [START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Komatitsch | An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[END_REF]). Thus, the numerical dimension is 570

x 620 x 160 = 56 million grids for ∆s = 200 m, 1120 x 1220 x 310 = 423 million for ∆s = 100 m and 1356 x 1716 x 429 = 1 billion for ∆s = 75 m. The calculation duration is 60 seconds.

The upper frequency limit in the simulations is usually estimated as max min / (5 ) f V s = ∆ , where min V denotes the minimum wave velocity in the structural model being used [START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF].

Based on several tests using different grid spacings and filter frequency ranges, we find it empirically to be a quite useful parameter.

Our simulation procedure is illustrated schematically in Figure 2. The input files both for the structure and the source need to be handled with care. For the structure, we read the original files provided by different studies and assign the material property at each grid in our simulation program. For the source models, however, we format the source files prior to the simulations after carefully ascertaining how the inversions have been achieved. Since the synthetic near-field ground motion is quite sensitive to the source description in time and space, in order to reproduce the ground motions using the source models inverted from the observations, we must understand exactly how they are solved in terms of the Green's function being used and detailed source time description (cf. Aochi et al., 2011).

The numerical dimension is not overly large, compared to some advanced simulations of wave propagation using finite difference methods (cf. [START_REF] Olsen | Ground Motion Estimates Using an Ensemble of Large Earthquakes on the Southern San Andreas Fault With Spontaneous Rupture Propagation[END_REF][START_REF] Furumura | An integrated simulation of ground motion and tsunami for the 1944 Tonankai earthquake using high-performance super computers[END_REF]. In order to be able to repeat several simulations for verification and calibration purposes, we improved time performance through a hybrid implementation using MPI

3D Structural Models Geological Features

In this section, we will be explaining three 3D structure models we use in this study.

Their key features are summarized in Table 1, and the cross-sections of each model are shown on Figure 3. The flat part of this region and the folded hills in between, as well, are composed of thick sedimentary layers. Conversely, the surrounding mountains to the south and east are characterized by rocks. The complex lateral variations in basement structure along the NE-SE fault striking are ascribed to Miocene rifting during the opening stage of the Sea of Japan and the subsequent shortening of the crust [START_REF] Okamura | Rifting and basin inversion in the eastern margin of the Japan Sea[END_REF][START_REF] Kato | Reactivation of ancient rift systems triggers devastating intraplate earthquakes[END_REF]. Because the Mw 6.6 Chuetsu earthquake in 2004 occurred 30 km southeast of the area ruptured in 2007, this region is recognized as an active fracture zone or strain concentration zone, which had been identified through structural geology [START_REF] Okamura | Rifting and basin inversion in the eastern margin of the Japan Sea[END_REF] and GPS observations [START_REF] Sagiya | Continuous GPS array and present-day crustal deformation of Japan[END_REF]. Because this area exhibits complex 3D structures, the velocity structure must be correctly reproduced if we are to generate sound synthetic ground motions.

The ERI Model

The

The NIED Model

The second model (hereinafter referred to as the "NIED" model) is taken from tune up the structure along cross-sections between source and receiver, particularly for the aftershocks of the 2004 Chuetsu earthquake (2D tuning). The file (J-SHIS) provided is described in terms of depth to the layer boundaries at mesh nodes of 45″ in longitude and 30″ in latitude (the 3 rd mesh for land planning in Japan), namely about a 1 km x 1 km mesh. It contains 32 layers between the surface and the basement characterized by an S-wave velocity (Vs) of 3300 m/s, but many of these are not present in the region of concern. The minimum

Vs is set at 350 m/s.

The GSJ Model

The third model, referred to in this paper as the "GSJ" model, is presented in [START_REF] Sekiguchi | 3D subsurface structure model of the Niigata sedimentary basin[END_REF]. They improved the earlier NIED model [START_REF] Fujiwara | A subsurface structure modelling of whole of Japan for strong-motion evaluation[END_REF]

Comparison

The original features of each model are provided in the corresponding references.

Figure 3 shows the cross-sections for these three models derived from the numerical simulations. We note that while there are many layers in the NIED and GSJ models, all of them are not always present at the finite difference grid point. We then interpolate the material parameters by averaging (Graves, 1996) instead of estimating precisely the interface plane between the finite difference grids [START_REF] Moczo | 3D heterogeneous staggered-grid finite-difference modelling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities[END_REF]. While uncertainties linked to the numerical processing remain on the scale of a calculation grid, the models used do adequately represent characteristic features in Figure 3. The three models are visually similar in terms of the shape of the basin structure. Although the tomography study (the ERI model) is carried out independently, it does detect the basin structure (low velocity zone) well. As the ERI model is inverted from dense observation data mainly above the aftershock area, the resolution toward the north (Y=20 km) is not good enough to be comparable with the other two models. The NIED and GSJ models are similar, but the bedrock depth seems different:

the NIED model is generally deeper.

Aftershock simulations

Model settings

In order to examine the structural models in detail, we begin by simulating some small earthquakes. We chose two aftershocks of the 2007 Niigata-ken Chuetsu-oki earthquake that were well located by [START_REF] Kato | Kanazawa and Groupe for the aftershock observations of the 2007 Niigataken Chuetsu-oki earthquake, Imaging heterogeneous velocity structures and complex aftershock distributions in the source region of the 2007 Niigataken Chuetsu-oki Earthquake by a dense seismic observation[END_REF] and [START_REF] Shinohara | Precise aftershock distribution of the 2007 Chuetsu-oki earthquake obtained by using an ocean bottom seismometer network[END_REF] and also have mechanisms obtained by the F-net broadband seismograph network (http://www.fnet.bosai.go.jp). Table 2 summarizes the source parameters. The hypocental depths given by the two catalogues differ by a few kilometres, although the epicenters are relatively well determined. The seismicity distribution in this area is often shallower than that routinely obtained (cf., the progress report of the Japanese national project on the observation along the strained zone, MEXT, 2010).

We impose a smooth bell-shaped source time function (cubic B-spline function) of a 0.5second duration at the hypocenter (Figure 4), which is reasonable for events of Mw 4.4. We simulate these sources for each of the three geological models.

Before discussing the simulation results, a question may be raised concerning the quality of finite difference simulations. In Appendix, we have provided our synthetic comparisons between the finite difference and spectral element methods using the ERI model for the first aftershock. It is thus found that the numerical simulations are reliable enough to provide the ground motion under the given model and conditions.

Near-Field Ground Motion

Figure 4 shows the comparison at station KZK (F-net) using the three structural models, all of which are simulated with s ∆ = 100 m. This station is located at a depth of 60 m on terrain characterized as rocky, not only locally but also regionally to a certain extent. In the simulations, the event origin times are set at 0. We align the observations for the origin time reported by the relocated earthquake catalogue (Table 2). The seismograms are not filtered at this station. First we observe that the P-wave arrives simultaneously for all the cases (UD component), but the S-wave arrival does differ somewhat. The ERI model, for which the earthquakes are relocated, yields the best arrival time, whereas the other two models, NIED and GSJ, display a delay of a few seconds (EW and NS components). Then, the two models predict much longer and stronger ground shaking than observations in the later phases. This would tend to indicate that their S-wave structure might need to be better tuned up, supposing that the relocated earthquake parameters are correct here. For this KZK station, the ERI model appears to be better constrained, presumably thanks to their 3D tomography coverage from the source area to the station. However, we find that the EW component of the first aftershock (07/16) displays a clear discrepancy between the synthetics and the observation, even with ERI model, as the later phases are those that become dominant in terms of amplitude. This issue is not perfectly resolved and represents a common tendency for this aftershock along the coastline (NIG018 and KK as well). The synthetics do not properly account for the later phases, possibly due to the more complex layer that is shallower than the given resolution (3 km in depth, see Table 1).

The Regional Wave Field

Figures 5 and6 show the comparison at different stations for different structural models for Aftershocks 1 and 2 respectively. Each seismogram is aligned as was described above. The seismograms are filtered with a 0.1 -0.5 Hz band-pass filter and only the EW component is shown; the response spectra for 5 % damping are also compared. There are still difficulties in reproducing both arrivals and later phases in a satisfactory manner, since they must exhibit a correct velocity profile on the direct path (first arrivals) as well as in a larger surrounding volume for the later phases. In this sense, it is not possible to determine whether one model is better than another at this stage insofar as none of them are perfect for a particular study.

However, some remarks can be made, also briefly summarized in Table 3. First, some models are seen to generate later phases that are too large at some stations, for example, the GSJ model at NIGH12. The first phase, on the other hand, is fairly accurately rendered by some models at most stations because the source and structure parameters are relatively well established. We note some time shift in phase among the simulations and/or with respect to the observation, at NIG004 and NIG025 among others, for example. This indicates that further improvement may be achieved both in source and structural parameters, but this is not this study's objective. Generally in the near field (cf., NIG018, NIGH11, NIG019, NIG017, NIG016, counter-clockwise), where the propagation path is relatively simple, the reproduction is quite good (especially Figure 6

Discussion

As already stated, the simulations are sometimes not as satisfactory for Aftershock 1.

Reasons for this could be that the source parameters are not well constrained or the local structure for this earthquake is not precise enough. As the arrival times of the main phases are fairly good in different directions, the earthquake location appears reliable. Our next step, then, was to invert Aftershock 1's focal mechanism using the ERI model and apply the neighbourhood algorithm so as to generate a wide range of parameters without prior information. After several trials using different stations and different components of the seismograms, we always obtain a focal mechanism of approximately (strike, dip, rake) = (181°, 65°, 84°), which differs from the reference solution only by 10 degrees. We may accordingly conclude that the focal mechanism obtained is stable. The remaining open question is: how is it possible to have only one component that does not fit well?

Mainshock simulations

Finite source models and simulations in 1D layered models

In this section, we chose two finite source models, both involving a southeast dipping reverse fault. Figure 7 shows the projection of the fault plane and the final slip obtained and Note that the data at NIG018 seems to have been distorted by the strong liquefaction at the site, and the record at KZK is saturated for the mainshock. In Figure 8, we show the source time function for each model, which is not retrieved directly from their provided file but rather is based on our input file after changing the file format to make it compatible with our code (e.g. Figure 2). This is important so as to verify our simulations and also for our subsequent discussion of how the strong ground motion is radiated.

We first entered these input files in the simple 1D structure and verified that the source model was properly taken into account. [START_REF] Aoi | Source process of the 2007 Niigata-ken Chuetsu-oki earthquake derived from near-fault strong motion data[END_REF] introduced the moving source effect when calculating their Green's functions semi-theoretically. We approximated this effect numerically by means of 5 x 5 sub-sub-point sources distributed on each sub-fault in our finite difference simulation, as we already demonstrated in other earthquakes (e.g. [START_REF] Aochi | Dynamic rupture of crosscutting faults: A possible rupture process for the 2007 Mw 6.6 Niigata-ken Chuetsu-Oki earthquake[END_REF]. For this earthquake, we observe that the moving source effect is also visible at NIG004.

Table 5 gives the 1D velocity structure models derived from the inversions. These are representative of their reference velocity models to a certain extent, but it should be borne in mind that the two inversions adopt different 1D velocity models independently calibrated at every station, that is, the wave field we simulate here does not simultaneously represent all the stations. We chose NIG004 and NIG016 for, as seen on Figure 9, the 1D structure from These verifications not only confirm our numerical implementations but also reflect some of the inversion procedures.

We first carried out simulations in 1D layered models. The reason for this is that most seismological and seismic hazard applications still call on this simple structure except for a few well-known areas in the world, such as southern California, southeastern Japan and Taiwan. In practice, the two source models that were adopted are deduced using this approximation. It accordingly is worthwhile to observe the effects of a 3D structure with respect to the 1D layered models.

Simulations in 3D structure models

It is not at all evident that the combination of any 3D structure model with any finite source model obtained within a 1D layer model can coherently reproduce the observation, although this procedure is recommended to predict the ground motion for scenario earthquakes intended for use in quantitative seismic hazard studies. It is thus useful now to seek to obtain the characteristic ground motions at different stations and discuss the capability and limitations of the given models. ). The ERI structure, for example, produces much larger later phases at NIG004, NIG025 and NIGH07, not noticeable in earlier comparisons of aftershocks (Figures 5 and6), implying that the validity of structure should be evaluated taking into account the frequency content of the source time function. Furthermore, it is noted that the observed motion is much larger (spectra over all the frequency range) than any synthetics at NIG010 and NIG011, and this discrepancy is more obvious for the mainshock than for the aftershock. This indicates that the mainshock generates much larger surface waves than the aftershock.

In this section, we combine different source and structure models to simulate the ground motions. For the same reason already cited in the previous section for the aftershock simulations, it is impossible to determine whether one model is better than the others, for this depends on which aspect one is looking at, since no combination is currently perfect.

Discussion

The numerical simulations reflect the current level of seismological knowledge concerning this earthquake. For practical reasons, it is thus important to study the variations in the calculated result from the standpoint of seismic hazard studies. It was predictable that the simulations for the mainshock would not fit the data well in terms of waveforms because the fault models are derived supposing some 1D structures: this may be the limit of our seismological knowledge. However, for engineering purposes, the variation of the simulated ground motions is meaningful. Figures 12 and13 show the PGV (Peak Ground Velocity) map (different frequency ranges up to 0.5 Hz) for the three structural and two source models. Each source model is also calculated with the reference 1D velocity models (Table 5). As observed in the seismograms (Figures 10 and11), the source model from Hikima & Koketsu yields stronger ground motion near the coastline rather than offshore; furthermore, the PGV with three asperities and obtained the comparable synthetic ground motions (frequency range 0.05 -1.6 Hz), strongly affected by these asperities as well as by their local 3D structure. The JNES model may be more precise around this point of interest but remain essentially local in scope. Thus for the purposes of a regional discussion as provided in this paper, this model should be further compiled with other models.

We are also interested in other engineering parameters which characterize the ground motion. Figure 14 shows the comparison of PGV values (up to 0.5 Hz) at K-net and KiK-net stations. In terms of PGV (one of the most simple engineering parameters), we can confirm that the simulations are globally consistent with observations within this frequency range except for NIG018, and probably for other stations where later phases are not well modelled (e.g. NIG010 in the ERI structure). This is because the PGV is less sensitive to the details of the rupture process, being mainly affected by macroscopic parameters like fault location, fault geometry, rupture directivity, rupture velocity, and final magnitude. In this sense, both source models are suitable for simulating regional ground motion around the fault.

As already seen, this earthquake is an example in which it is difficult to quantitatively reproduce the waveforms. One possible reason is that the complex 3D structure masks the real features of this earthquake. Many kinematic inversions adopt the southeast dipping fault geometry, but this orientation is not the only one sustainable. [START_REF] Aochi | Dynamic rupture of crosscutting faults: A possible rupture process for the 2007 Mw 6.6 Niigata-ken Chuetsu-Oki earthquake[END_REF] demonstrate the possibility of a dynamic rupture transfer from a northwest dipping sub-fault to another southeast dipping sub-fault. The two models we used in this study differ in fault orientation, in the position of asperities and in rupture timing, and actually do not fit the same stations. Our study reported in this paper does not aim to calibrate the parameters, but the computing performance we have achieved will allow us to investigate further the source parameters in a 3D structure as well as to refine the structural models. We see that even the subsurface structure obtained through fine tomography might not be sufficient to reproduce observed strong motions, nor those created through the use of both geologic and geophysical data. For some stations suffered from the strong ground motion (c.f. NIG018), it will be necessary to study by coupling with the nonlinear site effects at local level. We may conclude that a combination of tomography and data compiling in velocity structure modelling are called for.

Conclusion

Finite difference simulations of seismic wave propagation are carried out in the Niigata area, Japan, for the 2007 Mw 6.6 Niigata-ken Chuetsu-Oki earthquake at low frequencies. Some of the calculations are extended up to 1 Hz considering the minimum velocity in the medium of 350 m/s. However we limit our discussion to 0.5 Hz throughout this study because of lack of precision in the model. We test three 3D structure models, all of which are built differently in various studies. From the aftershock simulations, it is seen that none of the models are uniform in their resolution and precision in this region, as the later phases are poorly reproduced in some more distant stations. The model based on 3D tomography (ERI model) is good enough for the near field in terms of body waves (arrival time), but its precision for shallow sediment layers is insufficient to reproduce the later phases properly, while the other models (NIED and GSJ models), based on various 1D/2D profiles and geological interpretation, work well for the site response, but do sometimes cause a time shift in phases.

For the mainshock simulations, we adopt two different finite source models [START_REF] Aoi | Source process of the 2007 Niigata-ken Chuetsu-oki earthquake derived from near-fault strong motion data[END_REF]Hikima and Koketsu, 2008), which differ in the near-field ground motion, especially above the fault plane (but under the sea) and also along the coastline. It is found that each model is calibrated differently for the given stations. For engineering purposes, the variation observed in simulated ground motion is significant, but for seismological purposes, further parameter calibration is desirable and possible for such a complex 3D case using current high performance computing.

Figure A1 shows the comparison of the seismograms at NIG004 filtered in the frequency range between 0.1 and 0.5 Hz using the ERI model for Aftershock 1. Despite the completely different numerical procedures, the overall waveforms are well enough reproduced to allow the characteristics of the ground motion to be discussed. Figure A2 shows the goodness-of-fit (GOF; [START_REF] Kristekova | Time-frequency misfit and goodness-of-fit criteria for quantitative comparison of time signals[END_REF] for the surrounding K-net stations as a function of epicentral distance. This criterion provides a succinct evaluation of how similar two signals are by a score assigned on a scale of 0 (poor) to 10 (excellent), and ranges higher than 6.5 and 8.5 are considered as good and excellent, respectively. It is observed that the GOD is better at closer distance because of the reduced influence of numerical dispersion and of little impact from absorbing conditions. In fact, the worst station, NIG002, which is qualified as fair (range 4.5 to 6.5), is very close to the model edge. Otherwise, it is confirmed that the simulations are fair enough both in amplitude and phase for the purpose of our discussions.

Tables Table 1: 2) at station KZK using three different structural models. All the seismograms are not filtered. The earthquake location and the given moment release function are shown at the top.

Response spectra are calculated at the bottom for the waveforms up to 0.5 Hz. 2) at different stations for each structure model. The ground motion is calculated over 60 seconds, then filtered between 0.1 and 0.5 seconds. The response spectra are also shown. 2) at different stations for each structure model.

Figure 7: The two finite source models used in this study. Also see Table 4 for their characteristics.

  for the Niigata area by carefully calibrating layer depths and rock parameters based on seismological observations and taking into account the regional dependency of material parameters. The file (appended CD-ROM of[START_REF] Sekiguchi | 3D subsurface structure model of the Niigata sedimentary basin[END_REF] provides the depth of layer boundaries at mesh nodes of 0.00625 degrees longitude and 0.004167 degrees latitude, i.e., about a 0.5 x 0.5 km mesh. Model 2 designated in their CD-ROM contains 50 layers above the Moho. Because the parameters beneath the basement (Vs = 3300 m/s) are not supplied, we have used those from their previous study in other regions, the central part of Japan[START_REF] Sekiguchi | Broadband ground motion reconstruction for the Kanto Basin during the 1923 Kanto earthquake[END_REF].

  for Aftershock 2). The response spectra indicate the validity of the average response around the given stations. The simulation of Aftershock 2 is found to be more consistent with the observation than that of Aftershock 1, for which the simulations underestimate the long-period potion of the spectra at the nearest station along the coast. This once again indicates that Aftershock 1 could afford further calibration, probably of the focal depth and/or local structure around the source. In the Kashiwazaki area in the simulation of Aftershock 2, station KZK (a rock site) is well reproduced by the ERI model, while NIG018 (a soil site) is rendered best by the GSJ model, reflecting local geological conditions.

Aoi

  et al. corresponds better to NIG004 and that from Hikima & Koketsu to NIG016. The compared ground motions were filtered between 0.1 and 0.5 Hz and calculated with s ∆ = 200 m. The origin time is set at 10:13:22.16 (JST) for the observation, and the simulations start at t=0. Our synthetic ground motion for the Aoi et al. model is consistent with the observation at NIG004. Note that their inversion uses only 14 seconds, consisting mainly of the S-wave. This confirms our implementation of the finite source models. At NIG016, the synthetic motion from Hikima & Koketsu reproduces the characteristic waveforms, especially for the vertical component, with time shift of a few seconds. Such a time shift is common in inversions, as the location of the hypocenter and the origin time are not always the same.

Figures 10 and 11

 11 Figures 10 and 11 show the comparison of a number of synthetic ground motions and

  generated from the model by Aoi et al. is centred along the blind fault trace. This latter feature is common for thrust fault earthquakes due to the geometry. In other words, it is confirmed that the model from Hikima & Koketsu integrates some complexity in the rupture process that is acquired particularly from the very near-source records such as the station at the Kashiwazaki-Kariwa nuclear power plant. On the other hand, different 3D structure models briefly give the similar characteristic. The expansion of the green area is consistent with the existence of low-velocity sediments beneath the Sea of Japan between Sado Island and the mainland as well as in the Niigata plain (see topographical map in Figure3), which cannot be seen in 1D structure modelling. The amplification under the sea is stronger for the ERI structural model, which is the only model constrained from the OBS stations, although the shallow sediment layer under the sea is difficult to evaluate. Cirella et al. (2008) also obtained the finite source model by the joint inversion of strong motion and GPS data and show their forward modelling result in the 1D structure. Concerning the ground motion pattern above the fault plane, the model by Cirella et al. (2008) is closer to the one by Hikima & Koketsu used here. Kawabe and Kamae (2010) simulate the wave propagation in a 3D structure model provided by JNES (Japan Nuclear Energy Safety Organization, internal report, 2005) focusing on the Kashiwazaki-Kariwa nuclear plant. They used a characterized, simple source model

Figure 1 :

 1 Figure 1: Map of the area surrounding the 2007 Niigata Chuetsu Oki earthquake. The star

Figure 2 :

 2 Figure 2: Schematic representation of the simulation procedure.

Figure 3 :

 3 Figure 3: Cross-sections of different 3D structure models in the numerical simulations. See

Figure 4 :

 4 Figure 4: Comparison of synthetic and observed ground motion for the selected aftershocks

Figure 5 :

 5 Figure 5: Comparison of synthetic and observed ground motion for the first aftershock (Table

Figure 6 : 5 :

 65 Figure 6: Figure 5: Comparison of synthetic and observed ground motion for the second

Figure 8 :

 8 Figure 8: Source time function prepared as the input file for the numerical simulations. Each

Figure 9 :

 9 Figure 9: Synthetic ground motions at NIG004 and NIG016 for the source models from Aoi et

Figure 10 :

 10 Figure 10: Comparison of the synthetic ground motions at different stations for the Aoi source

Figure 11 :

 11 Figure 11: Comparison of the synthetic ground motions at different stations for the Hikima &

Figure 12 :

 12 Figure 12: Map of peak ground velocity (PGV) for each component (x: east-west, y: north-

Figure 13 :

 13 Figure 13: PGV map for the source model from Hikima & Koketsu.

Figure 14 :

 14 Figure 14: Comparison of horizontal (geometric mean) PGV for six simulations (see also

Figure A1 :

 A1 Figure A1: Comparison of synthetic seismograms at NIG004 calculated using the finite

Figure A2 :

 A2 FigureA2: The goodness-of-fit between the two simulations at the surrounding K-net stations.

Figure 3 : 39 Figure 4 :

 3394 Figure 3:

Figure 5 :

 5 Figure 5:

Table 4

 4

summarizes their characteristics. We used model B (southeast dipping) drawn from

[START_REF] Aoi | Source process of the 2007 Niigata-ken Chuetsu-oki earthquake derived from near-fault strong motion data[END_REF]

. The model initially obtained by

Hikima and Koketsu (2007)

was recently integrated into the work of

[START_REF] Miyake | Source fault of the 2007 Chuetsu-oki, Japan, earthquake[END_REF]

. We found hypocental locations only about 1 km apart; however, due to differences in strike and dip, the disrepancy between the two fault planes is greater. Both models are analysed using available near-field strong ground motion.

The only significant difference is that Hikima and Koketsu (2007) used the nearest station data (KK) at the Kashiwazaki-Kariwa nuclear power plant, while

[START_REF] Aoi | Source process of the 2007 Niigata-ken Chuetsu-oki earthquake derived from near-fault strong motion data[END_REF]

do not.

Table 2 :

 2 Characteristics of the three different structural models. See the references in the text. Source parameters of the aftershocks used in the simulations. The relocated hypocentral parameters are based on Kato et al. (2008) and Shinohara et al. (2008). The focal mechanisms are obtained routinely by F-net of NIED.

		ERI model	NIED model	GSJ model
	Area covered	Around the area of the 2004	All of Japan	Niigata area
		and 2007 earthquakes		measuring about 200
				km x 200 km
	Principal data	P-and S-wave tomography	Seismic reflection,	NIED model,
	included	of the aftershock	geological	geological
			interpretation, 1D	interpretation,
			tuning, 2D tuning	geophysical borehole
				data, 1D tuning
	Given	3 km (N35°E) x 5 km	45″ in longitude and	0.00625 degrees in
	resolution	(N125°E) x 3 km at least	30″ in latitude	longitude and
				0.004167 degrees in
				latitude
	Vs minimum	866 m/s	350 m/s	400 m/s
	Given	Vp, Vs by point	Vp, Vs, ρ by layer	Vp, Vs, ρ by layer
	property

Table 3 :

 3 Brief summary of characteristics of the three 3D structural models from the aftershock simulations.

	Model	Positive characteristics	Aspects to be improved
	ERI	Synthetics at near-field rock site,	Finer shallow structure, larger area
		especially in the south
	NIED	Synthetics at near-field soil site,	Parameters in shallower soft layers,
		good site response generally	deep structure
	GSJ	Synthetics at near-field soil site,	Same as above
		good site response generally,
		similar to NIED

Table 4 :

 4 Characteristics of the finite fault models used in this study. Some parameters were not specified by Hikima and Koketsu (2007), but we have supplied an estimated condition indicated in the blank.

	Model	Aoi et al.	Hikima & Koketsu
	Reference	Model B from Aoi et al. (2008) Hikima & Koketsu (2007)
	Hypocentral location	37.54°N, 138.61°E	37.54014°, 138.62653°E
	Focal depth	8.9 km	10 km
	Fault geometry	Strike N49°E, dip 42°	Strike N38°E, dip 34°
	Total seismic moment	1.62 x 10 19 N.m	Not specified (1.02 x 10 19 N.m)
	Sub-fault number	15 (strike) x 12 (dip)	15 (strike) x 9 (dip)
	Sub-fault size	2 km x 2 km	2 km x 2 km
	Source time discretization 6 multi-windows	7 multi-windows
	Rake	Free in each window	Free in each window
	Source time function	A smoothed ramp function	Unknown (Ramp function)
	Moving rupture effect	Yes	Unknown (No)

Acknowledgments

This work is carried out mainly in the framework of the French national ANR DEBATE project (DEvelopment of Broadband Accleration Timeseries for Engineers; 2009-2012). We thank all the colleagues for making available their numerical results: Drs. A. Kato, S. Aoi, H. Sekiguchi, K. Hikima and K. Koketsu. Many fruitful discussions were held with Profs. Raul Madariaga, Dr. Fabian Bonilla and Prof. Hiroshi Kawase, and Dr. Ken'ichi Tsuda from the DEBATE project. We thank Prof. Peter Moczo and Dr. Eiichi Fukuyama for improving the manuscript. We also express our appreciation for the data provided by K-net, KiK-net, and Fnet of the National Research Institute for Earth Science and Disaster Prevention. HA thanks NIED for his invitation to Japan in October 2010. MY thanks AIST-GSJ for funding his visit to BRGM in 2009-2010. Part of the simulations were carried out at the super computing center (CINES) in France.

Appendix: Comparison between finite difference and spectral element methods

Finite difference simulations are widely used for practical applications, but their quality is sometimes debated especially in the framework of synthetic benchmark tests (e.g. [START_REF] Day | Tests of 3D elastodynamic codes: Final report for Lifelines project 1A01[END_REF][START_REF] Chaljub | Quantitative comparison of four numerical predictions of 3D ground motion in the Grenoble Valley, France[END_REF]. As pointed out by many researchers, the finite difference scheme does not represent any interface correctly (free surface, faulting and material interface) but always gives an approximate solution. In this study, we do not deal explicitly with the interface plane in the numerical simulations even though the models initially supplied (the NIED and GSJ models) are defined by layers. As none of the three models are defined at the finite difference grids of calculation, estimating the material interfaces at given points requires further assumptions.

In this appendix, we are describing a comparative test using the finite difference as well as the spectral element methods (De Martin, 2011). We adopt the first example, namely the ERI structure model for the first aftershock. In the spectral element scheme, the material properties are interpolated at each of the GLL (Gauss-Lobatto-Legendre) points. The entire domain is uniformly meshed with hexahedra, constrained by the lowest S-wave velocity. The mesh is composed of 208,936 geometrical nodes, 196,425 hexahedra 1200 m in size in all directions and 16,757 quadrangles for absorbing boundaries by paraxial approximation different from our finite difference simulation. The polynomial order of the basic functions is N = 4, and the total number of degrees of freedom is 46 million (including the redundant GLL at the interfaces between the CPUs). Also, the finite difference simulation here is carried out with a grid spacing of 200 m.


Erwan Faou

F R Éd

Éric Rousset

LANDAU DAMPING IN SOBOLEV SPACES FOR THE VLASOV-HMF MODEL

Keywords: 1991 Mathematics Subject Classification. 35Q83, 35P25 Vlasov equations, Scattering, Landau damping, HMF model. 1

We consider the Vlasov-HMF (Hamiltonian Mean-Field) model. We consider solutions starting in a small Sobolev neighborhood of a spatially homogeneous state satisfying a linearized stability criterion (Penrose criterion). We prove that these solutions exhibit a scattering behavior to a modified state, which implies a nonlinear Landau damping effect with polynomial rate of damping.

Introduction

In this paper we consider the Vlasov-HMF model. This model has received much interest in the physics literature for many reasons: It is a simple ideal toy model that keeps several features of the long range interactions, it is a simplification of physical systems like charged or gravitational sheet models and it is rather easy to make numerical simulations on it. We refer for example to [START_REF] Barré | The Vlasov equation and the Hamiltonian Mean-Field model[END_REF], [START_REF] Yamaguchi | Stability criteria of the Vlasov equation and quasi stationary states of the HMF model[END_REF], [START_REF] Barré | Algebraic damping in the one-dimensional Vlasov equation[END_REF], [START_REF] Caglioti | Long time estimates in the mean field limit[END_REF], [START_REF] Caglioti | Quasi-stationary states for particle systems in the mean-field limit[END_REF] for more details.

We shall study the long time behavior of solutions to this model for initial data that are small perturbations in a weighted Sobolev space to a spatially homogeneous stationary state satisfying a Penrose type stability condition. We shall prove that the solution scatters when times goes to infinity towards a modified state close to the initial data in a Sobolev space of lower order. This result implies a nonlinear Landau damping effect with polynomial rate for the solution which converges weakly towards a modified spatially homogeneous state. In the case of analytic or Gevrey regularity, this result has been shown to hold for a large class of Vlasov equations that contains the Vlasov-Poisson system by Mouhot and Villani [START_REF] Mouhot | On Landau damping[END_REF] (see also the recent simplified proof [START_REF] Bedrossian | Landau damping: paraproducts and Gevrey regularity[END_REF]). Some earlier partial results were obtained in [START_REF] Caglioti | Time asymptotics for solutions of Vlasov-Poisson equation in a circle[END_REF], [START_REF] Hwang | On the existence of exponentially decreasing solutions of the nonlinear Landau damping problem[END_REF]. The related problem of the stability of the Couette flow in the two-dimensional Euler equation has been also studied recently [START_REF] Bedrossian | Inviscid damping and the asymptotic stability of planar shear flows in the 2D Euler equations[END_REF]. The question left open in these papers is the possibility of nonlinear Landau damping for Sobolev perturbations. In this case, one cannot hope for an exponential damping, but we can wonder if it could occur by allowing polynomial rates. For the Vlasov-Poisson system, it was proven in [START_REF] Lin | Small BGK waves and nonlinear Landau damping[END_REF] that this is false with rough Sobolev regularity due to the presence of arbitrarily close travelling BGK states. Nevertheless, this obstruction disappears for sufficiently high Sobolev regularity as also proven in [START_REF] Lin | Small BGK waves and nonlinear Landau damping[END_REF]. These arguments can also be extended to a large class of Vlasov equations and in particular the HMF model. Note that the regularity of the interaction kernel does not play an essential part in the argument of [START_REF] Mouhot | On Landau damping[END_REF] (though the decay in Fourier space provided by the Coulomb interaction seems critical), it is more the nonlinear "plasma-echo" effect which is crucial to handle. Besides the physical interest of the HMF model, it is thus mathematically interesting to study the possibility of nonlinear Landau-damping in Sobolev spaces for this simple model.

∂ t f (t, x, v) + v∂ x f (t, x, v) = ∂ x R×T P (x -y)f (t, y, u)dudy ∂ v f (t, x, v),
where (x, v) ∈ T × R and the kernel P (x) is given by P (x) = cos(x). Note that the main difference with the Vlasov-Poisson equation is the regularity of the kernel: in this latter case, P (x) = k≥0 k -2 cos(kx) is the kernel associated with the inverse of the Laplace operator.

The HMF model is thus the simplest nonlinear model with the structure (1.1). We consider initial data under the form f 0 (x, v) = η(v) + εr 0 (x, v) where ε is a small parameter and r 0 is of size one (in a suitable functional space). This means that we study small perturbations of a stationary solution η(v). We shall thus write the solution at time t under the form

f (t, x, v) = η(v) + εr(t, x, v).
We are interested in the study of the behavior of f when time goes to infinity. To filter the effect of the free transport, it is convenient to introduce (as in [START_REF] Mouhot | On Landau damping[END_REF], [START_REF] Bedrossian | Landau damping: paraproducts and Gevrey regularity[END_REF]) the unknown g(t, x, v) = r(t, x + tv, v) that is solution of the equation

(1.2) ∂ t g = {φ(t, g), η} + ε{φ(t, g), g}.
where

(1.3) φ(t, g)(x, v) = R×T (cos(x -y + t(v -u)))g(t, y, u)dudy and {f, g} = ∂ x f ∂ v g -∂ v f ∂ x g is the usual microcanonical Poisson bracket.
We shall usually write φ(t) when the dependence in g is clear.

We shall work in the following weighted Sobolev spaces, for m 0 > 1/2 be given, we set

(1.4) f 2 H n = |p|+|q| n T×R (1 + |v| 2 ) m 0 |∂ p x ∂ q v f | 2 dxdv,
and we shall denote by H n the corresponding function space. Note that compared to the usual Sobolev space H n , there is also a fixed weight (1 + |v| 2 ) m 0 in physical space. The interest of the weight is that functions in H 0 are in L 1 and thus it allows to get a pointwise control in Fourier (see Lemma 2.1 below). We do not include the dependence in m 0 in the notation since m 0 will be fixed. We shall denote by • or F the Fourier transform on T × R given by fk

(ξ) = 1 2π T×R f (x, v)e -ikx-iξv dxdv.
Note that due to the regularity of the interaction kernel and the conservation of the L p norms, it is very easy to prove the global well-posedness of the Vlasov-HMF model in H s for every s ≥ 0. Nevertheless, in order to study the asymptotic behaviour of g, the regularity of the kernel is not of obvious help. Indeed, when performing energy estimate on (1.2), it costs one positive power of t each time one puts a v derivative on the kernel. 1.2. The Penrose criterion. We shall need a stability property of the reference state η in order to control the linear part of the Vlasov equation (1.2). Let us denote by η, the spatially homogeneous stationary state and let us define the function

K(n, t) = -np n nt η0 (nt)1 t≥0 , t ∈ R, n ∈ Z,
where (p k ) k∈Z are the Fourier coefficients of the kernel P (x). We shall denote by K(n, τ ) = R e -iτ t K(n, t) dt the Fourier transform of K(n, •). We shall assume that η satisfies the following condition

(H) (1 + v 2 )η(v) ∈ H 5 and ∃ κ > 0, inf Im τ ≤0 |1 -K(n, τ )| ≥ κ, n = ±1.
Note that thanks to the localization property of η in the first part of the assumption, the Fourier transform of K can be indeed continued in the half plane Im τ ≤ 0. Here, the assumption is particularly simple due to the fact that for our kernel, there are only two non-zero Fourier modes. This assumption is very similar to the one used in [START_REF] Mouhot | On Landau damping[END_REF], [START_REF] Bedrossian | Landau damping: paraproducts and Gevrey regularity[END_REF] and can be related to the standard statement of the Penrose criterion. In particular it is verified for the states η(v) = ρ(|v|) with ρ non-increasing which are also known to be Lyapounov stable for the nonlinear equation (see [START_REF] Marchioro | A note on the nonlinear stability of a spatially symmetric Vlasov-Poisson flow[END_REF]).

1.3. Main result. In the evolution of the solution g(t, x, v) of (1.2), an important role is played by the quantity

(1.5) ζ k (t) = ĝk (t, kt), k ∈ {±1}, such that φ(t, g) = 1 2 k∈{±1} e ikx e iktv ζ k (t). Note that for k = 0, ζ k (t) is the Fourier coefficient in x of the density ρ(t, x) = R f (t, x, v) dv.
This quantity also plays a key part in the analysis of [START_REF] Mouhot | On Landau damping[END_REF], [START_REF] Bedrossian | Landau damping: paraproducts and Gevrey regularity[END_REF]. Note that here we need only to control two Fourier modes due to our simple interaction kernel. Let us define for every s ≥ 4 and T ≥ 0 the weighted norm

Q T,s (g) = sup t∈[0,T ] g(t) H s t 3 + sup t∈[0,T ] sup k∈{±1} t s-1 |ζ k (t)| + sup t∈[0,T ] g(t) H s-4 .
Our main result is:

Theorem 1.1. Let us fix s ≥ 7 and R 0 > 0 such that Q 0,s (g) ≤ R 0 and assume that η ∈ H s+4 satisfies the assumption (H). Then there exists R > 0 and ε 0 > 0 such that for every ε ∈ (0, ε 0 ] and for every T ≥ 0, we have the estimate

Q T,s (g) ≤ R.
As a consequence, we obtain the following scattering result:

Corollary 1.2. Under the assumption of Theorem 1.1, there exists a constant C and g ∞ (x, v) ∈ H s-4 such that for all r ≤ s -4 and r ≥ 1,

(1.6) ∀ t ≥ 0, g(t, x, v) -g ∞ (x, v) H r ≤ C t s-r-3 .
The consequence of such results is the following nonlinear Landau damping effect: as

g(t, x, v) is bounded in H s-4 , the solution f (t, x, v) = η(v)+εr(t, x, v) = η(v)+εg(t, x-tv, v) satisfies ∀ n ∈ Z * , ∀ξ ∈ R, ∀ α + β = s -4, | fn (t, ξ)| = ε|ĝ n (t, ξ + nt)| ≤ Cε ξ + nt α n β .
The last estimate being a consequence of the elementary embedding Lemma 2.1. This yields that for every n = 0, fn (t, ξ) tends to zero with a polynomial rate.

Moreover, by setting

η ∞ (v) := η(v) + ε 2π T g ∞ (x, v)dx,
we have by the previous corollary (and again Lemma 2.1) that for r ≤ s -4,

∀ ξ ∈ R, | f0 (ξ) -η∞ 0 (ξ)| ≤ C ξ r t s-r-3
In other words, f (t, x, v) converges weakly towards η ∞ (v).

The remaining of the paper is devoted to the proof of Theorem 1.1. We shall obtain Corollary 1.2 in section (4) as an easy consequence. As pointed out in [START_REF] Mouhot | On Landau damping[END_REF] the control of the "plasma echoes" that can be seen as kind of resonances is crucial to prove nonlinear Landau damping. These resonances occur when nt = kσ in the last integral term of (2.6). The main structural property of the Vlasov-HMF model that makes possible the following short proof of nonlinear Landau damping in Sobolev spaces is that the resonances are easy to analyze, the only possibility is when n = k = ±1 and t = ±σ. Moreover, the structure of the nonlinearity then allows to control it without loss of decay. Making an analogy with dispersive equations (see [START_REF] Klainerman | The null condition and global existence to nonlinear wave equations[END_REF] for example), the nonlinearity of the Vlasov-HMF model could be thought as a nonlinearity with null structure. Our approach also allows to handle the case of a kernel with a finite number of modes which allows more resonances, we briefly sketch the modification in section 5. Nevertheless, it is still unclear if this can be done for the Vlasov-Poisson equation.

A priori estimates

In this section, we shall study a priori estimates for the solution of (1.2). Let us fix s ≥ 7 and introduce the weighted norms:

(2.1) N T,s (g) = sup t∈[0,T ] g(t) H s t 3 , M T,γ (ζ) = sup t∈[0,T ] sup k∈{±1} t γ |ζ k (t)| so that (2.2) Q T,s (g) = N T,s (g) + M T,s-1 (ζ) + sup [0,T ] g(t) H s-4 .
Let us take R 0 > 0 such that Q 0,s (g) ≤ R 0 . Our aim is to prove that when ε is sufficiently small we can choose R so that we have

Q T,s (g) ≤ R
for every T ≥ 0.

In the following a priori estimates, C stands for a number which may change from line to line and which is independent of R 0 , R, ε and T .

We shall make constant use of the following elementary lemma.

Lemma 2.1. For every α, β, n ∈ N with α + β = n we have the following inequality:

(2.3) ∀ k ∈ Z, ∀ ξ ∈ R, | fk (ξ)| 2 n/2 C(m 0 ) k -α ξ -β f H n ,
where C(m 0 ) depends only on m 0 and where x = (1 + |x| 2 ) 1/2 for x ∈ R.

Proof. We have by using the Cauchy-Schwarz inequality that

k α ξ β fk (ξ) = 1 2π T×R ∂ α x ∂ β v f (x, v)e -ikx e -ivξ dxdv C f H n R (1 + |v| 2 ) -m 0 dv 1/2
.

The previous inequality with α = β = 0 yields the result when k = 0 or |ξ| 1 and we conclude by using x 2 α/2 |x| 2 for |x| > 1.

2

(2.4) M T,s-1 (ζ) ≤ C R 0 + εR 2 .
Proof. We first note that in Fourier space, the equation (1.2) can be written after integration in time,

(2.5) ĝn (t, ξ) = ĝn (0, ξ)

+ t 0 p n ζ n (σ)η 0 (ξ -nσ)(n 2 σ -nξ)dσ + ε k∈Z p k t 0 ζ k (σ)ĝ n-k (σ, ξ -kσ)(nkσ -kξ)dσ,
for all (n, ξ) ∈ Z × R, with p k = 1 2 for k ∈ {±1} and p k = 0 for k = ±1, and where the ζ k (t) are defined by (1.5). Setting ξ = nt in (2.5), the equation satisfied by (ζ n (t)) n=±1 can be written under the almost closed form

(2.6) ζ n (t) = ĝn (0, nt) - t 0 p n ζ n (σ)η 0 (n(t -σ))n 2 (t -σ)dσ -ε k∈{±1} p k t 0 ζ k (σ)ĝ n-k (σ, nt -kσ)kn(t -σ)dσ.
To study the equation (2.6), we shall first consider the corresponding linear equation, that is to say that we shall first see

(2.7) F n (t) := ĝn (0, nt) -ε k∈{±1} p k t 0 ζ k (σ)ĝ n-k (σ, nt -kσ)kn(t -σ)dσ
as a given source term and we shall study the linear integral equation

(2.8) ζ n (t) = t 0 K(n, t -σ)ζ n (σ) dσ + F n (t) n = ±1
where the kernel K(n, t) has been introduced in section 1.2. For this linear equation, we have the estimate:

Lemma 2.3. Let γ ≥ 0, and assume that η ∈ H γ+3 satisfies (H). Then, there exists C > 0 such for every T ≥ 0, we have

M T,γ (ζ) ≤ CM T,γ (F ).
Let us postpone the proof of the Lemma and finish the proof of Proposition 2.2.

From the previous Lemma and (2.3), we first get that (2.9)

M T,s-1 (ζ) ≤ C g(0) s + εM T,s-1 (F 1 ) + εM T,s-1 (F 2 )
with where F 1 corresponds to the term with k = -n in (2.7) and F 2 corresponds to the term with k = n, hence

F 1 n (t) = -n 2 p -n t 0 ζ -n (σ)ĝ 2n (σ, n(t + σ))(t -σ) dσ, n = ±1, F 2 n (t) = n 2 p n t 0 ζ n (σ)ĝ 0 (σ, n(t -σ))(t -σ) dσ, n = ±1.
Let us estimate F 1 n , by using again (2.3) and the definition (2.1) of N σ,s , we get that

|F 1 n (t)| ≤ C t 0 (t -σ) σ 3 M σ,s-1 (ζ)N σ,s (g) σ s-1 t + σ s dσ ≤ C R 2 t s-1 +∞ 0 1 σ s-4 dσ ≤ C R 2 t s-1 provided s ≥ 6. This yields that for all T ≥ 0 M T,s-1 (F 1 ) ≤ CR 2 .
To estimate F 2 n , we split the integral into two parts: we write

F 2 n (t) = I 1 n (t) + I 2 n (t) with I 1 n (t) = n 2 p n t 2 0 ζ n (σ)ĝ 0 (σ, n(t -σ))(t -σ) dσ, n = ±1, I 2 n (t) = n 2 p n t t 2 ζ n (σ)ĝ 0 (σ, n(t -σ))(t -σ) dσ, n = ±1.
For I 1 n , we proceed as previously, 4 dσ and hence since s ≥ 6, we have M T,s-1 (I 1 ) ≤ CR 2 . To estimate I 2 n , we shall rather use the last factor in the definition of Q s,T in (2.2). By using again (2.3), we write

|I 1 n (t)| ≤ CR 2 t 2 0 σ 3 (t -σ) σ s-1 t -σ s dσ ≤ CR 2 t s-1 +∞ 0 1 σ s-
|I 2 n (t)| ≤ t t 2 M σ,s-1 (ζ) σ s-1 g(σ) H s-4 t -σ s-5 dσ ≤ CR 2 t s-1 +∞ 0 1 τ s-5 dσ ≤ CR 2 t s-1
and hence since s ≥ 7, we find again

M T,s-1 (I 2 ) ≤ CR 2 .
By combining the last estimates and (2.9), we thus obtain (2.4). This ends the proof of Proposition 2.2.

It remains to prove Lemma 2.3.

Proof of Lemma 2.3. Let us take T > 0, and let us set for the purpose of the proof K(t) = K(n, t), F (t) = F n (t)1 0≤t≤T . Since we only consider the cases n = ±1, we do not write down anymore explicitly the dependence in n. We consider the equation (2.10)

y(t) = K * y(t) + F (t), t ∈ R
setting y(t) = 0 for t ≤ 0. Note that the solution of this equation coincides with ζ n (t) on [0, T ] since the modification of the source term for t ≥ T does not affect the past. By taking the Fourier transform in t (that we still denote by • ), we obtain

(2.11) ŷ(τ ) = K(τ )ŷ(τ ) + F (τ ), τ ∈ R,
with K(τ ) = K(n, τ ). Under the assumption (H), the solution of (2.11) is given explicitely by the formula

(2.12) ŷ(τ ) = F (τ ) 1 -K(τ ) .
Let us observe that since (1 + v 2 )η 0 ∈ H 5 , we have by (2.3) that for α ≤ 2 and for t > 0 (2.13)

|∂ α t K(t)| ≤ C t 4 ∈ L 1 (R + ).
Note that by definition of K(t), the function K(t) is continuous in t = 0, but not C 1 . Using an integration by parts on the definition of the Fourier transform, we then get that (2.14)

|∂ α τ K(τ )| ≤ C τ 2 , α ≤ 2.
To get this, we have used that the function t η0 (t) vanishes at zero. By using this estimate on K, (H) and that F (τ ) ∈ H 1 τ (the Sobolev space in τ ) since F is compactly supported in time, we easily get that y defined via its Fourier transform by (2.12) belongs to H 1 τ . This implies that t y ∈ L 2 and thus that y ∈ L 1 t . These remarks justify the use of the Fourier transform and that the function y defined through its Fourier transform via (2.12) is a solution of (2.10). Moreover, thanks to (2.12) and (H), we get that ŷ can be continued as an holomorphic function in Im τ ≤ 0 and thanks to a Paley Wiener type argument, that y vanishes for t ≤ 0. We have thus obtained an L 1 solution of (2.10) that vanishes for t ≤ 0. By a Gronwall type argument, we easily get that there is a unique solution in this class of (2.10) and thus we have obtained the expression of the unique solution.

We can thus now focus on the proof of the estimate stated in Lemma 2.3. Note that a L 2 -based version of this estimate would be very easily obtained. The difficulty here is to get the uniform L ∞ in time estimate we want to prove.

We shall first prove the estimate for γ = 0. Let us take χ(τ ) ∈ [0, 1] a smooth compactly supported function that vanishes for |τ | ≥ 1 and which is equal to one for |τ | ≤ 1/2. We define χ R (τ ) = χ(τ /R) and χ R (∂ t ) the corresponding operator in t variable corresponding to the convolution with the inverse Fourier transform of χ R (τ ). Thanks to (2.14), we have that for R large

t 2 |(1 -χ R (∂ t ))K(t)| ≤ C α≤2 ∂ α τ ((1 -χ R (τ )) K(τ )) L 1 (R) ≤ C |τ |≥R/2 1 τ 2 ≤ C R and hence (2.15) (1 -χ R (∂ t ))K(t) L 1 (R) ≤ C R ≤ 1 2
for R sufficiently large. This choice fixes R.

To estimate the solution y of (2.10), we shall write that

y = χ 2R (∂ t )y + (1 -χ 2R (∂ t ))y =: y l + y h .
By applying (1 -χ 2R (∂ t )) to (2.10), we get that

y h = K * y h + (1 -χ 2R (∂ t ))F = (1 -χ R (∂ t )K * y h + (1 -χ 2R (∂ t ))F since (1 -χ R ) = 1 on the support of 1 -χ 2R .
Therefore, we obtain thanks to (2.15) and the fact that χ 2R (∂ t ) is a convolution operator with a L 1 function, that

y h L ∞ ≤ 1 2 y h L ∞ + C F L ∞
and hence

y h L ∞ ≤ 2C F L ∞ .
For the low frequencies, we can use directly the form (2.10) of the equation: We can write that

ŷl (τ ) = χ 2R (τ ) 1 -K(τ ) χ R (τ ) F (τ ).
Since the denominator does not vanish thanks to (H), we obtain again that y l can be written as the convolution of an L 1 function -which is the inverse Fourier transform of χ 2R (τ )/(1 -K(τ )) -by the function χ R (∂ t )F which is a convolution of F by a smooth function. Thus we obtain by using again the Young inequality that

y l L ∞ ≤ C F L ∞ . Since y L ∞ ≤ y l L ∞ + y h L ∞ ,
we get the desired estimate for γ = 0. To get the estimate for arbitrary γ, we can proceed by induction. We observe that ty(t) = K * (ty) + F 1 with F 1 = (tK) * y + tF . Using the result γ = 0, we obtain that ty L ∞ ≤ C F 1 L ∞ . Now since η 0 ∈ H γ+3 , for γ = 1, we obtain that tK ∈ L 1 and thus

F 1 L ∞ ≤ C tF L ∞ + y L ∞ ) ≤ C (1 + t)F L ∞ .
The higher order estimates follow easily in the same way.

Estimate of N T,s (g).

Proposition 2.4. Assuming that η ∈ H s+2 verifies the assumption (H), then there exists C > 0 such that for every T > 0, every solution of (1.2)

such that Q T,s (g) ≤ R enjoys the estimate N T,s (g) ≤ C(R 0 + εR 2 )(1 + εR)e CεR .
Proof. To prove Proposition 2.4, we shall use energy estimates. We set L t [g] the operator

L t [g]f = {φ(t, g), f } such that g solves the equation ∂ t g = L t [g](η + εg).
For any linear operator D, we thus have by standard manipulations that

d dt Dg(t) 2 L 2 = 2ε Dg(t), D(L t [g]g(t)) L 2 + 2 Dg(t), D(L t [g](η)) L 2 = 2ε Dg(t), L t [g]Dg(t) L 2 + 2ε Dg(t), [D, L t [g]]g(t) L 2 +2 Dg(t), D(L t [g](η)) L 2 ,
where [D, L t ] denotes the commutator between the two operators D and L t . The first term in the previous equality vanishes since L t [g] is the transport operator associated with a divergence free Hamiltonian vector field. Consequently, we get that

(2.16) d dt Dg(t) 2 L 2 2ε Dg(t) L 2 [D, L t [g]]g(t) L 2 + 2 t 0 Dg(t) L 2 D(L t [g](η)) L 2 .
To get the estimates of Proposition 2.4, we shall use the previous estimates with the operator D = D m,p,q defined as the Fourier multiplier by k p ξ q ∂ m ξ for (m, p, q) ∈ N 3d such that p+q s, m m 0 and the definition (1.4) of the H s norm. To evaluate the right hand-side of (2.16), we shall use Lemma 2.5. For p + q γ and m m 0 , and functions h(t) and g(t), we have the estimates

D m,p,q , L σ [g] h(σ) L 2 ≤ C m σ,γ+1 (ζ) h(σ) H 1 + m σ,2 (ζ) h(σ) H γ , (2.17) D m,p,q L σ [g] h(σ) L 2 ≤ C m σ,γ+1 (ζ) h(σ) H 1 + m σ,2 (ζ) h(σ) H γ+1 , (2.18)
for all σ, where ζ is still defined by ζ k (t) = ĝk (t, kt), k ∈ { ± 1}, and where

m σ,γ (ζ) = σ γ sup k∈{±1} |ζ k (σ)| ,
with a constant C depending only on γ, and in particular, does not depend on σ.

Let us finish first the proof of Proposition 2.4. By using the previous lemma with γ = s and (2.17) with h = g and (2.18) with h = η, we obtain from (2.16) that

d dt g(t) 2 H s ≤ t 2 m t,s-1 (ζ) η H 1 + ε g(t) H 1 g(t) H s + 1 t s-3 m t,s-1 (ζ) η H s+1 g(t) H s + ε t s-3 m t,s-1 (ζ) g(t) 2 H s .
This yields using the fact that

M t,γ (ζ) = sup σ∈[0,t] m σ,γ (ζ), g(t) H s ≤ g(0) H s + t 3 M t,s-1 (ζ) η H s+1 + εR + εR t 0 1 σ s-3 g(σ) s dσ for t ∈ [0, T ].
From the Gronwall inequality, we thus obtain

g(t) H s ≤ g(0) H s + t 3 M t,s-1 (ζ) η H s+1 + εR e εR +∞ 0 dσ σ s-3 .
By using Proposition 2.2, this yields

N T,s (g) ≤ R 0 + (R 0 + εR 2 )(C + εR) e CεR .
This ends the proof of Proposition 2.4.

Let us give the proof of Lemma 2.5.

Proof of Lemma 2.5. We give the proof of (2.17), the proof of the second estimate being slightly easier. In the Fourier side, we have for L σ [g](h) the expression

(F L σ [g]h) n (ξ) = k∈{±1} kp k ζ k (σ) ĥn-k (σ, ξ -kσ)(nσ -ξ).
Consequently, we obtain that

F ([D m,p,q , L σ [g]h) ) n (ξ) = k∈{±1} kp k ζ k (σ) n p ξ q ∂ m ξ ĥn-k (σ, ξ -kσ)(nσ -ξ) - (n -k) p (ξ -σ) q ∂ m ξ ĥn-k (σ, ξ -kσ)(nσ -ξ) .
For k = ±1, we can thus expand the above expression into a finite sum of terms under the form

I k n (σ, ξ) = kp k ζ k (σ)k p 1 (n -k) p-p 1 +α kσ q 1 +α (ξ -kσ) q-q 1 +β ∂ m 1 ξ ĥn-k (σ, ξ -kσ) where 0 ≤ p 1 ≤ p, 0 ≤ q 1 ≤ q, m -1 ≤ m 1 ≤ m, α + β = m 1 -m + 1, α, β ≥ 0.
Moreover, if m 1 = m, then we have p 1 + q 1 > 0.

We have to estimate n ξ | k∈±1 I k n (σ, ξ)| 2 dξ by isometry of the Fourier transform. We note that for a fixed k ∈ {±1} then for |n -k| + |ξ -kσ| ≤ |k|σ, we have

|I k n (σ, ξ)| ≤ Cσ p+q+1 |ζ k (σ)||n -k||∂ m 1 ξ ĥn-k (σ, ξ -kσ)| whereas for |n -k| + |ξ -kσ| ≥ |k|σ, we have |I k n (σ, ξ)| ≤ C σ 2 |ζ k (σ)|(|n -k| + |ξ -kσ|) γ |∂ m 1 ξ ĥn-k (σ, ξ -kσ)|.
Consequently by taking the L 2 norm, we find that k∈±1

I k n (σ, ξ) L 2 ≤ C m σ,γ+1 (ζ) h(σ) H 1 + m σ,2 (ζ) h(σ) H m .
This ends the proof of the Lemma.

2.3.

Estimate of g H s-4 . To close the argument, it only remains to estimate g H s-4 .

Proposition 2.6. Assuming that η ∈ H s+2 verifies the assumption (H), then there exists C > 0 such that for every T > 0, every solution of (1.2) such that Q T,s (g) ≤ R enjoys the estimate g(t)

H s-4 ≤ C R 0 + εR 2 )e CεR , ∀t ∈ [0, T ].
Proof. We use again (2.16) with D = D m,p,q but now with p + q ≤ s -4. By using Lemma 2.5, we find

(2.19) d dt g(t) 2 H s-4 ≤ m t,s-3 (ζ) η H s-3 g(t) H s-4 + ε g(t) 2 H s-4 . This yields g(t) H s-4 ≤ g(0) H s-4 + η H s-3 M t,s-1 (ζ) t 0 1 σ 2 dσ + εM t,s-1 (ζ) t 0 1 σ 2 g(σ) H s-4 dσ.
By using Proposition 2.2, we thus get

g(t) H s-4 ≤ C R 0 + εR 2 ) + εR t 0 1 σ 2 g(σ) H s-4 dσ.
From the Gronwall inequality, we finally find

g(t) H s-4 ≤ C R 0 + εR 2 )e CεR .
This ends the proof of Proposition 2.6.

Proof of Theorem 1.1

The proof of Theorem 1.1 follows from the a priori estimates in Propositions 2.2, 2.4 and 2.6 and a continuation argument. Indeed, by combining the estimates of these three propositions, we get that

Q T,s (g) ≤ C(R 0 + εR 2 )(1 + εR)e CεR
assuming that Q T,s (g) ≤ R. Consequently, let us choose R such that R > CR 0 , then for ε sufficiently small we have R > C(R 0 + εR 2 )(1 + εR) e CεR and hence by usual continuation argument, we obtain that the estimate Q T,s (g) ≤ R is valid for all times.

Proof of Corollary 1.2

In view of (2.5), let us define g ∞ (x, v) by

g ∞ (x, v) = g(0, x, v) + +∞ 0 {φ(σ, g), η + εg(σ)} dσ.
Note that the integral is convergent in H s-4 since thanks to (2.18), we have

{φ(σ, g), η + εg(σ)} H s-4 ≤ C(R) 1 σ 2 + σ 3 4
σ s-3 . Note that for the last estimate, we have used that by interpolation

g H s-3 ≤ C g 3 4 H s-4 g 1 4 H s ≤ C(R) σ 3 4 .
From the same arguments, we also find that

g(t) -g ∞ H s-4 ≤ C(R) +∞ t 1 σ 2 + 1 σ s-3-3 4 dσ ≤ C(R) t .
In a similar way, by using again (2.18), we have for r ≤ s -4 and r ≥ 1,

g(t) -g∞ H r ≤ C(R) +∞ t 1 σ s-r-2 + 1 σ s-3 dσ ≤ C(R) 1 t s-r-3 + 1 t s-4 ≤ C(R) t s-r-3 .

The case of a kernel with a finite number of modes

In this section, we briefly indicate the modifications in the case that in (1.1), the kernel P is defined by

P (x) = M k=1 p k cos(kx),
for some p k ∈ R and for a fixed M. For the Penrose criterion (H) , it suffices to consider that it holds for any n, |n| ≤ M, n = 0.

We can use the weighted norms

Q T,s (g) = sup t∈[0,T ] g(t) H s t 2M +1 + sup t∈[0,T ] sup |k|≤M, k =0 t s+1-2k |ζ k (t)| + sup t∈[0,T ] g(t) H s-2M -2 .
One can then obtain that Theorem 5.1. Let us fix s ≥ max(4M + 2, 2M + 5) and R 0 > 0 such that Q 0,s (g) ≤ R 0 and assume that η ∈ H s+4 satisfies the assumption (H). Then there exists R > 0 and ε 0 > 0 such that for every ε ∈ (0, ε 0 ] and for every T ≥ 0, we have the estimate

Q T,s (g) ≤ R.
It is then easy to get from this result a nonlinear damping effect as previously.

The proof of this result follows exactly the same lines as the proof of Theorem 1.1. The estimates for sup t∈[0,T ] g(t) H s t 2M +1 and for sup t∈[0,T ] g(t) H s-2M -2 can be obtained exactly in the same way as in Proposition 2.4 and Proposition 2. If k and n have opposite sign, then, we can proceed as in the estimate of F

1. 1 .

 1 The Vlasov-HMF model. The Vlasov-HMF model reads(1.1)

6 .

 6 The only technical difference is that in Lemma 2.5, we definem σ,γ (ζ) = σ γ sup |k|≤M, k =0 |ζ k (σ)|.The only part were we need to be careful is to estimatesup t∈[0,T ] sup |k|≤M, k =0 t s+1-2k |ζ k (t)| as in Proposition 2.2 since more resonances are possible in the integral equation (2.6). By using the Volterra equation (2.8) for |n| ≤ M, n = 0, we still get that sup t∈[0,T ] t s+1-2n |ζ n (t)| ≤ sup t∈[0,T ] t s+1-2n |F n (t)| and we only need to estimate the right hand side.The only difficulty is to estimate the contribution of the integral termsJ n = sup t∈[0,T ] t s+1-2n |k|≤M, k =0 |p k |t 0 |ζ k (σ)| |ĝ n-k (σ, nt -kσ)||kn|(t -σ)dσ for |n| ≤ M.

  .1. Estimate of M T,s-1 (ζ). Towards the proof of Theorem 1.1, we shall first estimate ζ k (t), k = ±1.Proposition 2.2. Assuming that η ∈ H s+2 verifies the assumption (H), then there exists C > 0 such that for every T > 0, every solution of (1.2) such that Q T,s (g) ≤ R enjoys the estimate

  1 n in the proof of Proposition 2.2, we find (σ)| |ĝ n-k (σ, nt -kσ)||kn|(t -σ)dσ bounded since s ≥ 4M + 2 and |n| ≥ 1. Now let us assume that k and n have the same sign. We can assume that k ≥ 1 and n ≥ 1, the other situation being similar. If n > k, then we have that nt -kσ ≥ (k + 1)t -kσ ≥ t and hence we can use the same bound as above. If n = k, we can proceed exactly as for the term F 2 n in the proof of Proposition 2.2. It remains to handle the case n < k which is new. For this one, we split the time integral in the region σ ≤ n 2k t and the region σ ≥ n 2k t. For the first region we have nt -kσ ≥ nt/2 and hence this part of the integral can be handled as previously. For the region σ ≥ n 2k t, we estimate it by This term is uniformly bounded since 1 ≤ n ≤ k -1.

		t
	t s+1-2n	0	|ζ k ≤ CQ t,s (g) 2 t s+1-2n	0	t	σ 1+2M (t -σ) σ s+1-2k t s dσ
						≤ CQ t,s (g) 2 t 2-2n	0	+∞	1 σ s-4M dσ
	t which is uniformly t s+1-2n 0 |ζ n 2k t	t σ s+1-2k dσ ≤ C t 2n-2k+2 .

k (σ)| |ĝ n-k (σ, nt -kσ)||kn|(t -σ)dσ ≤ CQ t,s (g) 2 t s+1-2n t

INRIA & ENS Cachan Bretagne, Avenue Robert Schumann F-35170 Bruz, France. E-mail address: Erwan.Faou@inria.fr Laboratoire de Mathématiques d'Orsay (UMR 8628) Université Paris-Sud et Institut Universitaire de France E-mail address: frederic.rousset@math.u-psud.fr


Marcin Golabczak
email: marcin.golabczak@p.lodz.pl

Philippe Jacquet
email: bphilippe.jacquet@ecam.fr

Corinne Nouveau
email: ccorinne.nouveau@ensam.eu

Romain Fliti
email: romain.fliti-4@etudiants.ensam.eu

Tribological investigations of TiC+a-C:H Coatings Manufactured on X38CrMoV5-1 Steel Using PVD Technology

Keywords: PVD method, Cathodic Arc Evaporation, TiC+a-C:H coatings, Tribological investigations, Pin-on-disc test, Wear resistance

come   L'archive ouverte pluridisciplinaire

Introduction

X38CrMoV5-1 steel (ASTM designation) is a very commercial hot working steel characterized by very good mechanical properties (high tensile strength, hardness, high abrasion resistance, high yield stress -2200 MPa). For this reason X38CrMoV5-1 steel is widely used in various branches of industry e.g. in die casting and plastic forming, cf. [1]. Moreover, its abrasive wear resistance and hardness (especially in the surface layer) can be significantly improved through thermochemical treatments or by application of surface coatings (nitriding, CVD and PVD coatings) as well as by laser alloying of the surface with tungsten carbide [2]. Steel producers try to outdo one another in developing specialized grades of steel and treatments for them. In this paper, the technology of manufacturing of TiC+a-C:H protective coatings on X38CrMoV5-1 steel using PAPVD (Plasma Activated Physical Vapour Deposition) method [3][4][START_REF] Wendler | Functional Coatings by PVD or CVD Methods[END_REF] is proposed. Moreover, the tribological investigation concerning measurement of friction coefficient and wear resistance of manufactured coatings is presented. For the investigations the X38CrMoV5-1 steel samples from LaBoMaP-Centre Arts et Métiers de Cluny-France were used. The chemical composition (in wt%) of investigated X38CrMoV5-1 steel is: 0.42% C, 0.5% Mn, 1.2% Si, 5.5% Cr, 1.4% Mo, 0.3% V, Fe balance. The series of samples of 30 mm diameter and 2 mm thickness were made. The main requirement of the preparation process of samples for deposition TiC+a-C:H protective coatings is to obtain the solid surface layer with low roughness and devoid of all impurities. To this end, the technological process of preparation of the samples is elaborated by using a precision grinderpolisher machine.

Preparation of X38CrMoV5-1 steel samples for investigations

The polishing process of samples made of X38CrMoV5-1 steel was carried out using Phoenix Beta 2 (Buehler-Germany) dual platen grinder-polisher machine equipped with a Vector power head and applying Buehler abrasive materials and lubricants. Process stages of realized operations contained:

 two-stage grinding of samples using BuehlerMet silicon carbide (SiC) abrasive paper P600 (grits size ø 26 µm) and P1200 (grits size ø 15 µm);  lapping of samples using medium hard woven silk cloth VerduTex and monocrystalline diamond suspension MetaDi (grains size ø 3 µm);  polishing of samples using soft synthetic pad ChemoMet and aluminium oxide (Al 2 O 3 ) final polishing suspension MasterPrep (grains size ø 0,05 µm). The elaborated technological process ensures suitable preparation of the steel samples in the range of required roughness parameters. It also ensures proper purity of the samples. The surface of polished samples had a silver glossy colour and invisible tool marks.

Experimental apparatus and procedure

For deposition of TiC+a-C:H coatings on X38CrMoV5-1 circular steel samples the modified hybrid PVD Cathodic Arc Evaporation (PVD CAE) system was used. The process of manufacturing of these coatings was realized in the Institute of Materials Science at Technical University of Lodz by using multisource, hybrid factory-scale equipment of type URM 079. This equipment allows for deposition of coatings by a physical method, cf. [3][4][START_REF] Wendler | Functional Coatings by PVD or CVD Methods[END_REF]. It has two independent continuous arc power supplies with magnetic filtration of microdrop for production of metallic plasma (maximum discharge current 180 V), two pulsating carbon plasma power supplies (average pulse power 1 kW) and magnetron source (power 10 kW) with a flat circular target of 100 mm of dimension. This equipment allows for plasma cleaning of the surface by four independent ion guns (4 kW/100 mA) in argon atmosphere. In the range of high voltages the polarization voltage of substrate was controlled from 800-2500 V. In the range of low voltages it was controlled from 0-200 V. The nominal vacuum was below 1,0x10 -4 Pa. The scheme of the vacuum chamber of PVD stand and its equipment is presented in Figure 1.

Ar, N2, C2H2, CH4, O2

To the pump unit The results obtained here show that the surface layer of X38CrMoV5-1 steel samples with TiC+a-C:H coatings have a regular, smooth structure and lower average roughness parameters (Ra=0.08 μm) than uncoated samples (Ra=0.087 μm).

Gas flow controllers

Continuous arc power supply

Ionic gun power supply

Ionic gun power supply

Pulsating carbon plasma power supply

Continuous arc power supply

Ionic gun power supply

Ionic gun power supply

Pulsating carbon plasma power supply

Magnetron power supply

Sample bias control unit

Sample

Tribological investigations of TiC+a-C:H coatings manufactured on X38CrMoV5-1 steel samples

Tribological investigations of TiC+a-C:H coatings manufactured on X38CrMoV5-1 steel samples were done by using computer-controlled pin-on-disk tribometer (CSM Instruments-Switzerland) with a heated sample holder and with applying steel (100Cr6) and corundum (Al 2 O 3 ) balls as a counter-samples. Tribological tests were conducted not only at ambient conditions (25°C) but also at a temperature 200°C because in some cases, moulds are heated in this range of temperature. Parameters of the friction coefficient test were as follows: load -5N, frequency -5 Hz, counter-sample ball radius -6 mm, average velocity -10 cm/s, distance 200-800 m. The results of friction coefficient measurements of uncoated X38CrMoV5-1 steel samples at 25ºC and 200ºC with using counter-samples made of different materials are presented in Figs. 34. Comparison of friction measurements of X38CrMoV5-1 steel with TiC-a-C:H coating at 25ºC and 200ºC with using counter-samples made of different materials are shown in Fig. 5 and Fig. 6. From the results presented above it follows that the measured friction coefficient of X38CrMoV5-1 steel samples with manufactured TiC-a-C:H coatings is significantly lower as compared with the steel samples without these coatings, especially at ambient conditions (25ºC). Moreover, it can be observed that the lowest values of friction coefficient both at 25ºC and 200ºC take place for corundum (Al 2 O 3 ) ball as a counter-sample and even after a 800 meter test the friction coefficient does not increase (Fig. 6). It means that TiC-a-C:H coatings manufactured on X38CrMoV5-1 are characterized by very high wear resistance.

X38CrMoV5-1 steel covered with TiC+a-C:H coating Results of the friction measurements performed with using

X38CrMoV5-1 steel covered with TiC+a-C:H coating Results of the friction measurements performed with using

Wear rate of TiC+a-C:H coatings manufactured on X38CrMoV5-1 steel

For investigations of wear tracks of TiC+a-C:H coatings manufactured on X38CrMoV5-1 steel samples the optical profilometry was done by using the WYKO NT 1100 optical profiler (Veeco Instruments-USA). It allows for obtaining the profile of the track. The coatings wear tracks look like circular grooves. With the cross-sectional area of the track taken on several points it is possible to calculate the worn volume of the coatings and then to calculate the wear rate. In order to obtain a representative average, for measurements we have chosen 8 sections, one at every 45°. One multiplies the average section by the perimeter to have the volume: (1) where V track is the worn volume (m 3 ), r is the track radius (mm) and S is the average section (mm²). The wear rates of the coating and ball are calculated according to the equation [START_REF] Holmberg | [END_REF]: (2) where: V is the worn volume (m 3 ), w is the load (N) and s is the sliding distance (m).

Obtained results are depicted in Figures 78. However pictures of wear tracks of TiC+a-C:H coatings using optical microscopy are shown in Figures 910. The wear rate of TiC+a-C:H coating with using 100Cr6 ball as a counter-sample at the temperature 25°C and after 800 meter test is very low (0.085x10 -6 mm 3 /Nm). However, at the temperature 200°C the wear rate rapidly increases to 2.899x10 -5 mm/Nm after the 450 meter test (Fig. 7 and Fig. 9). The wear rate of TiC+a-C:H coating with using Al 2 O 3 ball as a counter-sample at the temperature 25°C and 200°C as well as after the 800 meter test is very low and its averaged values are respectively 0.167x10 -7 mm 3 /Nm and 0.252x10 -7 mm 3 /Nm (Fig. 8 and Fig. 10).

Wear rate of the counter-samples made of 100Cr6 and Al 2 O 3 balls

The investigations of the wear rate of the counter-samples were done basing on the morphology of the balls scars with using optical microscope. The morphology of the balls scars gives information about the friction mode and the wear mode. The size of the scar allows for calculating the ball wear rate. The worn material of the balls is similar to a spherical cap.

The worn volume of the balls is calculated according to the equations [START_REF] Whitenton | Diffusion in Solids and Liquids VIII[END_REF]:

(3)

( ) or
( ) (4)
where: V cap is the worn volume (m 3 ), h cap is the height of the cap (mm), R ball is the radius of the ball (mm) and r scar is the radius of the scar (mm). The wear rate of 100Cr6 ball used as a counter-sample in friction test of TiC+a:-C:H coating at temperature 25°C is very small after the 200 meter test (6.206x10 -8 mm 3 /Nm) and after the 800 meter test decreases even to 1.777x10 -8 mm 3 /Nm. It means that the 100Cr6 ball probably wears during the running-in and acting at the same time like a polishing leveled the roughness of the coating. Since the surface of the ball and the coating are very smooth, the friction coefficient decreases and the wear slows down. However at the temperature 200°C the wear rate of 100Cr6 ball rapidly increases to 4.080x10 -6 mm/Nm after 450 meter test (Fig. 11 and Fig. 13). The wear rate of Al 2 O 3 ball used as a counter-sample at the temperature 25°C is also very small after the 200 meter test (2.628x10 -8 mm 3 /Nm) and after the 800 meter test it decreases to 6.790x10 -9 mm 3 /Nm. At the temperature 200°C the wear rate of Al 2 O 3 ball also decreases. After the 200 meter test it is equal to 1.870x10 -8 mm 3 /Nm, however after the 800 meter test it is equal to 7.870x10 -9 mm 3 /Nm (Fig. 12 and Fig. 14).

Final remarks and conclusions

Investigations carried out in this paper have confirmed that PVD CAE method enables deposition of protective TiC+a-C:H coatings on X38CrMoV5-1 steel samples. Elaborated technological polishing process of X38CrMoV 5-1 steel samples ensures suitable preparation of their surfaces in the range of required roughness parameters. It also ensures proper purity of the samples. From the results it evidently follows that TiC+a-C:H coatings deposited on X38CrMoV5-1 steel samples considerably improve their tribological properties. It was shown that the measured friction coefficient of steel samples with PVD coatings is significantly lower than the friction coefficient of uncoated steel especially with using Al 2 O 3 corundum ball as a counter-sample. Wear tracks of TiC+a-C:H coating sliding against corundum ball are very small even after a 800 meter test and the wear of the coating is not observed. It means that TiC+a-C:H coatings are characterized by a very high wear resistance.

temperature control unit Fig. 1 .

 unit1 Fig. 1. Scheme of arrangement of separate parts in the vacuum chamber of PVD Cathodic Arc Evaporation stand used for deposition of TiC+a-C:H coatings.

Fig. 2 .

 2 Fig. 2. Microstructure images of X38CrMoV5-1 steel samples: a) uncoated after polishing, b) with TiC+a-C:H coatings.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Friction measurements of uncoated X38CrMoV5-1 steel conducted at the temperature 25ºC and 200ºC using 100Cr6 ball as a counter-sample.

Fig. 6 .

 6 Fig. 6. Comparison of friction measurements of X38CrMoV5-1 steel with TiC-a-C:H coating conducted at the temperature 25ºC and 200ºC using Al 2 O 3 corundum ball as a counter-sample

Fig. 7 .

 7 Fig. 7. The wear rate of the TiC+a-C:H coating manufactured on X38CrMoV5-1 steel with using 100Cr6 ball as a counter-sample: a) coating wear rate as a function of temperature, b) coating wear rate as a function of sliding distance.

Fig. 9 .

 9 Fig. 9. Wear tracks of TiC+a-C:H coating manufactured on X38CrMoV5-1 steel sliding against 100Cr6 ball: a) after 800 meter test at 25°C, b) after 450 meter test at 200ºC.

Fig. 10 .

 10 Fig. 10. Wear tracks of TiC+a-C:H coating manufactured on X38CrMoV5-1 steel sliding against Al 2 O 3 ball: a) after 800 meter test at 25°C, b) after 800 meter test at 200ºC.

  Fig. 11. The wear rate of 100Cr6 ball used as a counter-sample: a) ball wear rate as a function of temperature, b) ball wear rate as a function of sliding distance.

Fig. 12 .Fig. 14 .

 1214 Fig. 12. The wear rate of Al 2 O 3 ball used as a counter-sample: a) ball wear rate as a function of temperature, b) ball wear rate as a function of sliding distance.