F Speisser

Francois Djeffal

Le Normand

N Boubiche

J El Hamouchi

J Hulik

M Abdesslam

C Speisser

F Le Normand

Dr F Le Normand

Kinetics of graphitization of thin diamond-like carbon (DLC) films catalyzed by transition metal

Keywords: diamond-like carbon, thin graphite films, metal catalyst, pulsed laser deposition, thermal treatment, Raman spectroscopy, optical transmission and absorption, electric transport measurements

come

Introduction

Graphene is a 2D planar material of hexagonal structure comparable to a honeycomb network [START_REF] Avouris | Graphene: Synthesis and Applications[END_REF]. This unique structure infers outstanding properties like very large electrical [START_REF] Neto | The electronic properties of graphene[END_REF] and thermal [START_REF] Jang | Thickness-Dependent Thermal Conductivity of Encased Graphene and Ultrathin Graphite[END_REF] conductivities, a high planar Young modulus combined with a high flexibility [START_REF] Lee | Measurement of the Elastic Properties and Intrinsic Strength of Monolayer Graphene[END_REF] and optical transparency [START_REF] Blake | Making graphene visible[END_REF]. In addition, ultra-thin graphitic films (UTGF) exhibit properties close to those of the graphene monolayer, at least with graphitic layers less than 5 monolayers [START_REF] Ferrari | Science and technology roadmap for graphene, related to-dimensional crystals, and hybrid systems[END_REF]. However, applications of UTGF as transparent electrodes require new elaboration technique with large size, low-temperature and thickness control. Moreover, it is also recommended to grow UTGF films directly on an appropriate substrate to avoid any transfer step. To reach these objectives at a moderate cost, the growth of UTGF by chemical vapor deposition (CVD) and related processes using metallic catalyst appears promising technique. However, the elaborated graphene films using this technique contain undesirable impurities like hydrogen and impose other transfer steps, which can degrade the grapheme performance [START_REF] Jia | Large Area, Few-Layer Graphene Films on Arbitrary Substrates by Chemical Vapor Deposition[END_REF]. Another solution to avoid the transfer step is the formation of graphitic layers from diamond-like carbon (DLC) films grown by pulsed laser deposition (PLD) at room temperature as well as by other PVD processes, followed by post-treatments. Under appropriate conditions, a tetrahedral amorphous carbon (ta-C) film is formed by PLD. These ta-C films exhibit many properties close to those of diamond, due to a high but tunable concentration of sp 3 -hybridized carbon [START_REF] Robertson | Diamond-like amorphous carbon[END_REF]. Moreover, they initially exhibit a sp 2 -enriched surface [START_REF] Fallon | Properties of filtered-ion-beam-deposited diamondlike carbon as a function of ion energy[END_REF][START_REF] Rey | Thermal stability of amorphous carbon films deposited by pulsed laser ablation[END_REF], that can be totally graphitized by subsequent thermal [START_REF] Orwa | Thermally-induced sp 2 clustering in tetrahedral amorphous carbon (ta-C) films[END_REF][START_REF] Chhowalla | Evolution of sp2 bonding with deposition temperature in tetrahedral amorphous carbon studied by Raman spectroscopy[END_REF][START_REF] Ferrari | Stress reduction and bond stability during thermal annealing of tetrahedral amorphous carbon[END_REF][START_REF] Sullivan | Stress relaxation and thermal evolution of film properties in amorphous carbon[END_REF] or laser [START_REF] Jung | The investigation of thermal effect on the properties of pulsed laser deposited diamond-like carbon films[END_REF] treatments. A double process has been described in the literature including: 1) the clusterization into graphitic cycles of sp 2 single or olefinic chain sites initially embedded in a large network of sp 3 carbon at rather low temperatures (low activation energy) and 2) the direct transformation of sp 3 carbon to sp 2 carbon at higher temperatures (high activation energy). These post-treatments not only infer a drastic stress reduction, but also the films exhibit low surface conductivity. A derived elaboration process includes the incorporation of a catalytic transition metal like Ni, Pt… during or after the DLC deposition. It is presumed these metallic nanoparticles would act as catalytic centers accelerating the graphitization at lower temperatures. This effect has been early described in the literature with a conductivity variation of several orders of magnitude, but this result has thereafter no more received much attention [START_REF] Dimigen | Tribological and electrical properties of metal-containing hydrogenated carbon films[END_REF]. The catalytic effect would indeed affect the graphitization process both on kinetic and thermodynamic viewpoints to elaborate high quality films, without structural or chemical modifications [START_REF] Sullivan | Stress relaxation and thermal evolution of film properties in amorphous carbon[END_REF]. To the best of our knowledge, few studies have investigated the kinetics of graphitization of thin diamond-like carbon (DLC) films at a given temperature [START_REF] Sullivan | The electronic transport mechanism in amorphous tetrahedrally-coordinated carbon films[END_REF]. In this perspective, in this paper we report both thermal and kinetic studies of the graphitization of Ni/ta-C/quartz films to elaborate few graphitic layers on top of the ta-C film. Raman investigations, optical transmission and conductivity measurements have been carried out to investigate the optical and electrical properties of the elaborated films for eventual applications as transparent electrodes, biosensors and electrochemical cells [START_REF] Tamulevičius | Diamond like carbon nanocomposites with embedded metallic nanoparticles[END_REF][START_REF] Bewilogua | History of diamond-like films: from first experiments to worldwide applications[END_REF].

Materials and Methods

Samples were prepared by a triple step process. First carbon was deposited on quartz (double-face polished) by PLD at room temperature with a vacuum base higher than 10 -3 Pa.

A KrF laser source (248 nm wavelength, 30 ns pulse frequency) was used with a nuclear-grade graphite target (purity 5N) at a substrate-target distance of 5 cm. The laser fluence was within the range 4-6 J/cm 2 , which corresponds to approximately 210 8 W/cm 2 , a power density for which we found a maximum of sp 3 carbon formation [START_REF] Stock | Diamond-Like Carbon layers obtained by Pulsed Laser Deposition in different conditions for conductive electrodes application[END_REF]. The deposited ta-C film thickness was around 15 nm. Further experimental details on the preparation, the determination of the sp 3 /sp 2 ratio and the carbon density of the films are reported in references [START_REF] Stock | Diamond-Like Carbon layers obtained by Pulsed Laser Deposition in different conditions for conductive electrodes application[END_REF][START_REF] Boubiche | Ultrathin DLC films coated with transition metals and its transformation into graphitic films by thermal or thermocatalytic treatments[END_REF]. Samples studied here have a density of about 3 g/cm 3 and 70% of sp 3 content.

Ni metal is subsequently deposited at room temperature by molecular beam epitaxy (MBE) on DLC/quartz samples at a minute amount of metal (around 0.15 nm of equivalent metal layer). The thickness of evaporated Ni was measured by Rutherford back scattering (RBS) as described in [START_REF] Bouanis | Diameter-controlled growth of SWCNTs using Ru as catalyst precursors coupled with hydrogen treatment[END_REF]. This is typically a deposition of less than one equivalent metallic monolayer, carried out 1) to avoid single layer formation and thus direct metallic surface conductivity and 2) to optimize the catalytic properties of small nanoparticles.

The samples Ni/ta-C/quartz were finally subjected to thermo-catalytic treatments in ultra-high vacuum (UHV) furnace (base pressure 10 -4 Pa). The samples were taken out of the UHV chamber after each treatment. Then new measurements (Raman spectroscopy, UV-vis transmission, conductivity) were carried out before a new thermo-catalytic treatment was implemented. The samples Ni-DLC-3 samples were treated from 573K up to 873K by steps of 100K (Table 1). After setting the samples into the UHV furnace, a heating rate of 10K/min was applied up to the desired temperatures. After treatment was achieved, samples were cooled down in UHV without instructions. Samples DLC-1 and DLC-2 were similarly treated at 773K from a first treatment of 30 min up to a cumulative time of 3765 minutes with variable increasing time sequence. Besides, a reference sample (DLC-4), similar to DLC-3 sample but without catalyst was investigated at 773K for 60 min. Table 1 summarizes the properties of the samples used in this work. Sample Ni-DLC-1 and Ni-DLC-2 exhibit similar characteristics, except the thickness values, which are slightly different. As it has been reported that the microstructure of ta-C films changes along the axis normal to the surface [START_REF] Siegal | Bonding topologies in diamondlike amorphous carbon films[END_REF], it is expected that the thickness could be a significant parameter in affecting the film performance. In addition, we check by ERDA that the samples are free from hydrogen impurities.

Sample

Laser fluence (J/cm the green laser at 532 nm. The light was filtered to get a laser power of 1.5 mW on the sample to avoid damages in the carbon layer, using a 50 µm hole (confocal configuration), a monochromator of 2400 lines/mm and a magnitude lens of 100. Raman spectra were recorded in the 300-3000 cm -1 range. For Raman mapping, a set of 2020 spectra distant each other by 0.5 m was recorded. The spectra were then treated by removing the background and normalizing the signal to the intensity of the G signal. For some films, the light was focused either on the surface or at the interface with quartz. In the last case, the sample was turned with the light crossing the transparent quartz down to the interface. This was possible as the quartz samples were also polished on bottom face (RMS in the range 0.5-1.0 nm). Moreover, we tuned the intensity ratio of the G line with the Si line at 520 cm -1 to check the accuracy of the focalization. The conductivity measurements were performed using Hall effect technique (ECOPIA HMS-5000 set up) in the Van der Pauw configuration. Electric conductivity, in the planar orientation, as well as the nature, the mobility and the density of the carriers were determined. Knowing the thickness t (cm) of the conductive layer, the conductivity  (S/cm) of the film can be expressed as:

 = w / (t  L Rs) (1)
where w and L are the width and the length of the sample, respectively, and Rs (/) denotes the sheet resistance.

The optical transmission of the films has been measured using a NIR Perkin ELMER spectrophotometer Lambda 19. The spectra were recorded from 250 to 850 nm, using a sweep rate of 120 nm/min. The slit width was fixed to 2 nm. Moreover, the recorded spectra were analyzed according to the approach described in [START_REF] Changshi | Optical gap determination[END_REF], derived from Tauc's analysis, which can separately determine the nature of the optical transition and deduce the optical band gap Eg.

Results and discussion

A) Structural characterization A B C D

Fig.1 A) and B) SEM images of the Ni-DLC-1 sample after 3765 min of cumulative thermocatalytic treatment at two magnitudes; C) Histogram of the size distribution of set of nanoparticles displayed in Fig.1A on more than 2 500 nanoparticles using ImageJ software; D) AFM image of sample DLC-1 (without nickel).

A SEM image of the Ni-DLC-1 sample after 3765 min of cumulative thermocatalytic treatment at 773K is displayed at two magnitudes in Fig. 1A andB, respectively. The low magnitude spectrum shows a large and homogeneous spread of particles on the substrate.

However, images recorded with a higher magnitude reveal that these spots are split into many single nanoparticles with a mean size of about 20-30 nm. This splitting can be the result of a possible carburization of the Ni particles. Ni in the presence of carbon can form a carbide phase Ni3C which is thermodynamically stable within the temperature range 473-673K [START_REF] Sinharoy | The formation and decomposition of nickel carbide in evaporated nickel films on graphite[END_REF].

The decomposition of unstable Ni3C particles above 673K might lead to such configuration for the nanoparticles. The Ni nature of the nanoparticles is validated by energy-dispersive Xray spectroscopy. The size distribution of the set of nanoparticles displayed in Fig. 1C and obtained over 2500 particles over a surface area of 208 μm 2 is randomly distributed with a maximum around 0.3 μm (Fig. 1C). The surface roughness was measured from AFM to 0.6 nm on sample DLC-1 (without nickel) (Fig. 1D).

B) Temperature treatment by Raman spectroscopy

We first investigated by Raman spectroscopy the effect of UHV thermal posttreatments. The sequence of Raman spectra (Ni-DLC-3 sample) are plotted on the surface and at the interface, respectively, as initially prepared and then after subsequent thermocatalytic treatments from 573 to 873K (Fig. 2A andB). From these figures, a broad complex band within the 1100-1650 cm -1 range characteristic of carbon contributions is recorded. These Raman spectra are roughly similar with previous reports on ta-C films prepared by PLD using the same nanosecond KrF source [START_REF] Fallon | Properties of filtered-ion-beam-deposited diamondlike carbon as a function of ion energy[END_REF] or under similar conditions [START_REF] Kalish | Thermal stability and relaxation in diamond-like-carbon. A Raman study of films with different sp 2 /sp 3 fractions (ta-C to a-C)[END_REF]. This broad band points out a D band centered around 1360 cm -1 and a G band centered at 1560-1590 cm -1 . The intense G band is due to the stretching 6E2g mode of sp 2 -hybridized carbon whatever the nature (single, chain, cycle, …) of sp 2 carbon sites. By contrast, the D band is due to a multiphotonic process of sp 2 carbons and involves defects like edge sites of graphitic cycles. Thus, the ratio ID/IG of the intensities of these bands is a measurement of the cyclisation of the sp 2 carbon of nanometer size when starting from a highly defective carbon material like ta-C prepared at room temperature [START_REF] Ferrari | Interpretation of Raman spectra of disordered and amorphous carbon[END_REF]. A careful analysis at the initial spectra focused on the surface shows that a D contribution is initially present at a low concentration. Thus, very small graphitic-like fragments are formed in the initial DLC structure. Nevertheless, it should be highlighted that the DLC structure consists in a sp 3 -rich carbon network surrounding small sp 2 -rich carbon fragments [START_REF] Robertson | Diamond-like amorphous carbon[END_REF]. Therefore, other very weak contributions arise at around 1060 cm -1 (T mode), 1190 cm -1 (1) and 1480 cm -1 . (3) carbon vibration modes at the interface between the graphitic clusters in growth or in formation and the carbon sp3 network. The two latter lines are characteristic of disordered sp 3 hybridization state of carbon [START_REF] Robertson | Diamond-like amorphous carbon[END_REF][START_REF] Chu | Characterization of amorphous and nanocrystalline carbon films[END_REF], like sp3-hybridized carbons in close vicinity of the sp2 clusters. It is generally reported that the Raman ratio between sp3 and sp2 carbon is very weak [START_REF] Wada | Raman efficiency measurements of graphite[END_REF], but in our case with a laser source at 2.33 eV and band gap around 1 eV, we may have some sp 3 carbon vibration modes.

Finally, the attribution of the T mode at 1060 cm -1 is more discussed, as being due either to carbyne chains or to diamond nanocrystalline form [START_REF] Chu | Characterization of amorphous and nanocrystalline carbon films[END_REF]. Besides, these lines are much more intense on the spectra recorded at the interface, which means that a higher level of sp 3 carbon is formed at the interface (Fig 2B). This suggests a gradient of sp 2 /sp 3 carbon across the film, in agreement with the subplantation model of carbon DLC film growth [START_REF] Kenzie | Compressive-stress-induced formation of thin-film tetrahedral amorphous carbon[END_REF][START_REF] Siegal | Ultrahard carbon nanocomposite films[END_REF][START_REF] Siegal | Amorphous-tetrahedral diamondlike carbon layered structures resulting from film growth energetics[END_REF]. By contrast the assignment of these lines to some Si-C bonds can be excluded [START_REF] Inoue | Raman spectra of amorphous SiC[END_REF].

A B C Figure 2: Raman spectra of Ni-DLC-3 sample after consecutive thermocatalytic treatments from 573 to 873K (100K stepwise, 60 min). A) surface and B) Interface analysis, respectively; C) Ni-DLC-3 (black) and DLC-4 samples (red) at 873K, 60 min.

Notable changes occur after successive thermocatalytic treatments from 573 to 873K. The intensity of the D band increases, the position of the G band shifts towards higher wavenumbers, from 1570 cm -1 to 1575-1580 cm -1 , and a narrowing of both D and G (from 170cm -1 to 130 cm - 1 ) bands occurs. By comparison, the ratio ID/IG is around 0.6 on the sample DLC-4 (without Ni) instead of 0.75 on the sample Ni-DLC-3 after 60 min at 873K (Figure 2-C). This is an indication of the catalytic effect of the Ni particle for sp 2 carbon cyclization. These two effects are explained by the process of cyclization of carbon sp 2 fragments aggregating carbon sp 2 sites.

Also, the arising at high temperature (T = 873K) of a weak 2D band at around 2700 cm -1 can be noticed (Fig. 2 A). This 2D band is a multiphotonic process, and its occurrence can be considered as a signature of some weak organization of the graphitic fragments in the direction normal to the planar graphene plane. In contrast, sp 3 carbon bands do not change in intensity and remains weak. The lower wavenumber lines at around 500 and 800 cm -1 are due to Si-O and Si-Si vibrations and do not evaluate in intensity.

However, it is important to note that the DL films are not in an equilibrium state, and thus evolution of the carbon structure is not only temperature-dependent, but also may be timedependent. This is why we performed a comprehensive kinetics study at a given temperature equals to 773K.

C) Kinetics of graphitization 1) Raman spectroscopy

The sequences of Raman spectra, recorded at both the surface and at the interface, are displayed as a function of the cumulative time of thermocatalytic treatment (from 30 to 3765 min) in Fig. 3A andB, respectively. From these figures, we observe a progressive increase of the D line compared to the G line. At the interface, we note also some weak but significant carbon sp 3 and T bands (with some dispersion of the intensities due to the error in the focalization at the interface), whereas the D and G bands characteristics of sp 2 carbon remains unchanged. We have calculated the ratio of intensities (surface areas) ID/IG according to a process analysis that has been described in more details in literature [START_REF] Boubiche | Ultrathin DLC films coated with transition metals and its transformation into graphitic films by thermal or thermocatalytic treatments[END_REF] and which is summarized here. The whole Raman lines in the range 1000-1650 cm -1 is fitted with five gaussian -shape contributions which are the G, D, 3, 1 and T lines around 1560-1590 cm -1 , 1360 cm -1 , 1470-1495 cm -1 , 1180 cm -1 and 1060 cm -1 , respectively. The baseline was first subtracted from Raman spectra before starting the analysis. Fig. 3 is an example of deconvolution of the Raman spectrum of the sample DLC-04 (773K-60 min). This analysis provides important information about each peak (FWHM, wavenumber, area and peak intensities). The intensity ratio of the two peaks D and G is used to evaluate the crystallite size of the sp2-hybridized carbons [START_REF] Ferrari | Interpretation of Raman spectra of disordered and amorphous carbon[END_REF]. The surface ID/IG curves are plotted in Figure 3A for the Ni-DLC-1 and Ni-DLC-2 samples. For both samples, Ni-DLC-1 and Ni-DLC-2, this ratio starts with a slow raise, and then the increase becomes more pronounced above 200-300 minutes. The saturation apparently does not occur, even after near 4 000 minutes of thermal treatment. As these samples have rather similar characteristics, both curves give a good indication of the experimental dispersion of the results.

A B C D

Fig 3 A) Sequence of Raman spectra as a function of cumulative thermocatalytic treatments at 773K in the range time 0-4000 minutes. Ni-DLC-1 sample. A) on the surface; B) at the interface; C) Fitting of DLC-4 at surface, 773K, 60 min and D) Variation of the Raman intensity (surface areas) ratio ID/IG as a function of cumulative thermocatalytic treatment at 773K. Ni-DLC-1 and Ni-DLC-2 samples.

If we consider that the variations of the D band are mainly due to the nucleation and growth of graphitic cycles, the increase of the ID/IG ratio is an indication of the formation and growth of graphitic cycles in the course of the thermocatalytic treatments. This formation of graphitic cycles clearly takes a long time before completing stabilization and is not even achieved at 773K after 4 000 minutes of cumulative thermal treatments. Meanwhile, the G band remains stable around 1575 cm -1 and the FWHM of the G band decreases from 164 cm - 1 to 137 cm -1 . It is expected however that the process involves not only the formation and growth of graphitic rings, but also ordering of the graphitic cycles and a reorientation of them in a direction parallel to the surface [START_REF] Chu | Characterization of amorphous and nanocrystalline carbon films[END_REF]. The occurrence of a weak 2D band at 2700 cm -1 above around 700 min of cumulative treatment, even if it remains invariably low, is an indication of a starting ordering of the graphitic layers.

2) Raman mapping

The Figs.

A B C D

Fig. 4: A) and B): Raman images with a surface area of 100 μm² of the Ni-DLC-1 sample showing the ratio of normalized intensities (ID / IG) and (I2D / IG), respectively, after a thermocatalytic treatment of 3765 min at 773K. C) optical image of the same area. D) (ID/ IG) image of the same sample after thermocatalytic treatment of 180 min at 773K with the same amplitude (different area).

3) UV-vis transmission

The Then, some saturation is attended after approximately 1000 min (Fig. 5C in which the

A B C D

Fig.5: A) sequence of optical transmission on samples Ni-DLC-1 with cumulative thermocatalytic treatment. The transmission curve of the quartz is also reported. Scan rate of 120 nm/min; B) UV-Vis spectra of Ni-DLC3 (red) and DLC4 (black) at 873K, 60 min; C) Transmission (%) at 550 nm as a function of time (corrected from the quartz absorption); D) Band gap as a function of time (extracted from the transmission curves according to [24].

transmissions at 550 nm are corrected for the quartz absorption). These results are in good agreement with the progressive formation of graphitic aromatic clusters. However, a close inspection in the near-UV range reveals some slight changes in the transmittance behavior. We finally compare the samples Ni-DLC-3 and DLC-4 (without Ni) of same thickness and both treated 60 min at 873K (Fig. 5B). The sample DLC-4 exhibits a transmission of around 90% at 550 nm, much higher than the transmission on Ni-DLC-3 sample around 65%. This is a clear evidence of the catalytic behavior of Ni which accelerates the formation of aromatic domains, highly absorbing the light.

4) Transport measurements

We have first checked that the untreated Ni-DLC-01 and Ni-DLC-02 samples exhibit no conductivity within the detection limits (less than 10 -3 S/cm). The same behavior is observed for the DLC-4 sample without catalyst annealed at 773K for 60 min. increase with a linear slope of 0.32, the second one occurs in the range 200-1500 min with a steeper increase (linear slope of 0.7) and finally a saturation or even a slight decrease is observed above 1000 min. The overall increase is larger than one order of magnitude (Fig. 6A), from initially around 8.10 2 S/cm to 7.10 3 S/cm after 1245 min of cumulative thermocatalytic treatment. The two samples yield very similar conductivities, except for the third regime. The main conclusion is that the transport conductivity of Ni-DLC samples is never stabilized at 773K of annealing temperature over a large range time.

Fig. 6: Logarithmic plots of A) Conductivity (S/cm) measured in the Van der Pauw configuration (the linear correlations display the two first regimes as discussed in the text); B) Carrier mobility as a function of the cumulative annealing time (Ni-DLC-2 sample).

We compared the conductivities of Ni-DLC-3 and DLC-4 samples treated under the same conditions (873K, 30 min). We find 1.9 10 3 S/cm for the metal-loaded sample and we do not observe any conductivity on the unloaded sample. The catalytic effect of the metal on conductivity is clearly observed.

The mobilities of the carriers are reported in Fig. 6B for Ni-DLC-2 sample. The Hall effect measurements indicate that conductivity of thin layers graphite was both due to positive (holes) and negative charge carriers (electrons). This affects the determination of the Hall constant of the carrier mobility and of the carrier concentration.

The Hall constant RH is then expressed as

RH = (1/e)(p -nb 2 ) / (p + (nb)) 2 (2)
where p and n are the positive and negative carrier concentrations, respectively, and b = e/h is the ratio of the mobility values. The total density of charge carriers n+p can then be expressed as:

n+p = (1+R)  (1-Rb 2 ) / eRH  (1+Rb) 2 (3)
The electron and hole mobilities e and h, respectively, defined as  = (een + hep), can be expressed as:

e = RH b / (1-Rb 2 ) (4)
and

h = e / b = RH / (1-Rb 2 ) (5)
where R = n/p is the ratio of charge densities. Neglecting again the contribution of impurities, we take this ratio as R = 1.33 obtained on pyrolytic graphite [34]. To support this assumption, it has been checked by a statistical analysis (more than 100 independent Hall measurements) that the ratio of major n type over p type conductivity is 1.24. Fig. 6B reveals that the electron mobility increases with time to reach a maximum mobility value of 20 cm 2 /(Vs) for a thermocatalytic time equals to 2000 min. The carrier density given by equation ( 3) is quoted to around 10 20 -10 21 cm -3 . It must be noted that the transport properties obtained, with mobility of around 20 cm 2 /V.s at the best, are quite comparable with a recent publication on ta-C:H samples annealed at higher temperatures [START_REF] Zhai | Evolution of structural and electrical properties of carbon films from amorphous carbon to nanocrystalline graphene on quartz glass by HFCVD[END_REF].

Discussion

It is worthy to mention that the kinetics of surface conductivity is rather complex and its interpretation needs further investigation. However, it is possible to rely such behavior with the Raman analysis.

The relation between the conductivity and the nanostructure has been often studied in ta-C films, but only the thermal effect has been investigated on ta-C [START_REF] Sullivan | The electronic transport mechanism in amorphous tetrahedrally-coordinated carbon films[END_REF][START_REF] Ilie | Electronic transport, photoconductivity and photoluminescence in amorphous carbon[END_REF][START_REF] Ilie | Role of sp 2 phase in field emission from nanostructured carbons[END_REF] and ta-C-H films [START_REF] Dasgupta | Electrical conductivity of amorphous carbon and hydrogenated amorphous carbon[END_REF]. Two processes occur that contributes in improving the conductivity with time or temperature: 1) clustering of the sp 2 phase into sp 2 cyclic domains and 2) carbon sp 3 to carbon sp 2 conversion. However, owing to high activation energy we can rule out this last process that occurs only at high temperatures. Moreover, the sp 2 cluster orientation is also important in the surface conductivity, with an optimum configuration for surface conductivity when the cluster domains are oriented parallel to the surface. At room temperature, transport in ta-C films can be described by hopping between near neighbor sp 2 clusters and the conductivity . The latter is expressed as [START_REF] Ilie | Electronic transport, photoconductivity and photoluminescence in amorphous carbon[END_REF][START_REF] Ilie | Role of sp 2 phase in field emission from nanostructured carbons[END_REF]:

 = K(d 2 /T)  exp{-[E / (kT)] -(2d)/a]} ( 6
)
where a is the localization radius of the conducting states, d is the distance between the near neighbor sp 2 clusters, E expresses the activation energy of the transport process and K is a constant. Assuming a random and homogeneous distribution of graphitic spherical clusters, the localization radius a can be assigned to the cluster radius. The hopping distance d can be expressed as d = D -2a, where D is the mean distance between two near-neighbor sp 2 clusters.

Then D can be expressed as a function of the cluster density N as D = k / N 1/2 , where N represents the cluster density. Thus, the hopping distance d can be expressed as a function of the cluster density N and the cluster mean size as:

d = (k / N 1/2 ) -2a (7)
at a fixed temperature and from expressions ( 6) and ( 7), the conductivity can be given by:

 = K'(T)  [(k / N 1/2 ) -2a] 2  exp{-[2k / (aN 1/2 )] (8)
where K'(T) is a constant that involves all the temperature-dependent parameters. In the absence of any more data about the independent time variations of the localization radius a(t)

and the cluster density N(t), it is difficult to go further in the analysis without supplementary assumptions. However, the localization radius is correlated to the cluster size La in a hopping transport model as:

a = C  La / [Eπ(La) -E] 1/2 (9)
where Eπ(La) and E are the electronic gaps between π-π* and -* states, respectively. Eπ(La) is strictly depending on the cluster size as reported by Robertson [START_REF] Robertson | Electronic and atomic structure of amorphous carbon[END_REF]. For large clusters however, like here as ascertained by large Raman ID/IG ratio, Eπ(La) is constant, so a  C1  La with C1 = C / [Eπ -E] 1/2 . In the phase of graphitic clusters growth with increasing ID/IG, ID/IG  La 2 [START_REF] Ferrari | Interpretation of Raman spectra of disordered and amorphous carbon[END_REF]. Moreover, we can assume that the density N(t) is rapidly constant as the metal particle acts as a nucleation center for graphitic clusters, then N(t)  N0. Then equation {8} reduces to:

Ln  = C(T) -(D / (ID/IG) 1/2 ) (10)
A plot of Ln  = versus (ID/IG) -1/2 is displayed in Figure 7 for both samples Ni-DLC-1

and Ni-DLC 2 samples. It shows a reasonable agreement with a negative linear slope of -1/2, considering the error on ID/IG.

Another explanation can be found in close insight on the variations of the Raman ratio ID/IG with time. This ratio is expected to be dependent on both the density and on the size of the cluster radius, and clearly supports the increase of the conductivity with time. Moreover, the power laws with time distinguish two regimes. We propose that the first regime is a nucleation kinetics domain, where the conductivity is mainly governed by the variations of N 1/2 , whereas the second regime is a growth kinetics regime which is mainly governed by the variations of the mean size a. This could be related to the slope ratio between the first and the second regime around 0.32 / 0.7  0.46, not far from ½ as expected from equation [START_REF] Robertson | Diamond-like amorphous carbon[END_REF]. Moreover, the improvement of the conductivity at a given temperature may be also the result of a reorientation of the sp 2 clusters from the normal direction to the surface to an orientation parallel to the substrate. In the PLD process under our conditions it is believed that carbon ions impinge the surface with kinetic energy around 100 eV or more. Within these kinetic energies sp 2 carbon may be formed with a preferential orientation along the pathway of carbon ions inside the ta-C film, thus rather perpendicular to the surface. This has been shown for example for the formation of graphitic clusters normal to the surface during the diamond nucleation assisted by carbon ion bombardment [START_REF] Lifshitz | Visualization of Diamond Nucleation and Growth from Energetic Species[END_REF]. It is expected that the graphitic clusters will tend to reorient themselves with time or temperature in the direction parallel to the surface or the interface with substrate to minimize the bi-axial stress of the graphitic cluster sheets.

Fig. 7: Correlation between the conductivity variation and the Raman ID/IG ratio according to the equation (10).

Therefore, the improvement of the conductivity at a given temperature may be also the result of a reorientation of the sp 2 clusters from a direction normal to the surface to an orientation parallel to the substrate. This is supported by the literature report where conductivities were recorded in both direction, perpendicular and parallel directions to the surface as a function of thermal treatment. A cross is clearly observed between two measurements [START_REF] Ilie | Role of sp 2 phase in field emission from nanostructured carbons[END_REF]. This process is expected to be a slow process as it involves the movement of the carbon atoms of the graphitic cluster. Finally, the last regime is a saturation/degradation of the conductivity. It is proposed that in this regime the density and/or the size of the graphitic clusters become so high that their coalescence starts.

Conclusions

In

  4A and B display the Raman surface mapping of the normalized intensities ratio ID/IG and I2D/IG obtained on the sample Ni-DLC-1 after thermocatalytic treatment of 3765 min at 773K, respectively, associated with the optical image on the same area (Fig.4C). The two ID/IG images (Fig.4A and D) are displayed with the same color amplitude of 0.3 (0.5-0.8 and 0.9-1.2 range, respectively, with mean values 0.65 and 0.95, in agreement with the ID/IG intensities reported in Fig.3Dof 0.75 and 0.95 at 180 min and 3765 minutes, respectively).Clearly, the last image displays bright areas contrasted with darker zones. The limited resolution of the camera however precludes to assign these areas to catalytic particles which are in the nanometer scale. The size of the catalytic domains (Fig1B) is however of the same order of magnitude than the bright spots delivered by the optical image in Fig.4C. Therefore, it is tempting to assign the small dark red domains with high ID/IG ratio values (1.15-1.20 in the intensity scale of Fig.4A) to the Ni large domains observed in Fig.1. These areas can be assigned to preferential graphitic areas. In contrast, the areas with the lower ID/IG ratio values are considered as carbon sp 3 -enriched domains. The mapping at 180 min exhibits much less contrast indicating that between 180 and 3765 min both nucleation and growth of domains with high ID/IG occur. Reversely, we fail to find a matching between the (I2D/IG) and the (ID/IG) mappings (Fig.4A and 4B, respectively), or a matching between the (I2D/IG) mapping and the optical image (Fig 4Band 4C, respectively). The reason may be due to the very weak intensity of the 2D line, even after a long treatment at 773K. Moreover, it is likely that local conditions like roughness, orientation of the sp 2 cycles are more important at this level.

  Fig 5A plots some optical transmission curves of the Ni-DLC-1 sample with successive thermocatalytic post-treatments at 773K in the UV-visible range 250-800 nm. The transmission of the quartz alone remains constant, at around 93%. In this figure, it is clearly shown the transmittance is decreased to 65% at 550 nm for Ni-DLC-1 sample. The influence of the metal deposition on the transmission spectrum of the untreated sample (not shown) is negligible. After a first treatment, at 773K for 30 min, the transmittance is slightly decreased to 55% at 550 nm.

Fig

  5D plots the variation of optical band gap Eg as function of time extracted according to[START_REF] Changshi | Optical gap determination[END_REF]. Eg is initially measured on DLC-01 and DLC-02 samples at 1.1 eV. It drops to 0.95 eV with deposition of nickel. Starting the thermal treatment at 773K, it drops again to slightly less than 0.8 eV. It is shown in Fig.5Dthat the optical band gap increases first with time to reach a maximum value of 0.85-088 eV for an annealing time equals to 60 min. After 60 min, the optical band gap is decreased to near 06-065 eV with some trend to saturation. This signifies that the cumulative annealing time has a significant impact on the DLC film electrical behavior.

Fig

  Fig 6A illustrates the variation of the electrical conductivity, , versus the cumulative annealing time. From this figure, it is clearly shown that the samples Ni-DLC-1 and Ni-DLC-2 exhibit high conductivity values in the 10 3 -10 4 S/cm range. The conductivity behavior shows three power regimes in a logarithmic scale. The first one, up to 200 min, displays a small

  this work, we have studied the kinetics of graphitization at 773K of thin diamondlike carbon (DLC) films coated with Ni metallic nanoparticles. DLC films are deposited at room temperature by pulse laser deposition (PLD) on a transparent quartz substrate. Moreover, a Ni metal (less than one equivalent monolayer) is deposited on the surface of DLC using molecular beam epitaxy technique at room temperature in order to create Ni nanoparticles. Ultra-high vacuum annealing treatments (range 573-873K with 60 min) show a transformation of carbon sp 2 single or chain sites into cyclic aromatic domains which is probed by the Raman ID/IG ratio.Studied for the first time to our knowledge, a comprehensive kinetic study of thermocatalytic treatments on Ni-DLC samples in the range 30-3760 min at 773K shows a continuous process of sp 2 clusters formation and growth using Raman spectroscopy, optical transmission and electrical conductivity measurements. The kinetic behavior is analyzed in terms of nucleation and growth of graphitic sp 2 clusters, probably encompassed by a long-range kinetics of reorientation of the graphene layers of the graphitic clusters moving from normal to parallel direction to the surface and the interface. A correlation is found between the Raman ID/IG ratio and the conductivity in agreement with the nucleation and growth of graphitic clusters. Further studies are however required to monitor independently the kinetic of graphitic cluster formation, growth, orientation and coalescence. High surface electrical conductivities are reached accordingly, at more than 610 3 S/cm, with carrier mobility equals to 20 V/cm 2 s and a high density of carriers which are almost equally positive and negative. These characteristics make the elaborated DLC film potentially suitable for developing new high performance conductive electrodes and optoelectronic sensors.

2 ) DLC thickness (nm) Ni thickness (nm) Annealing time (min) at 773K

						Annealing
						temperature T
						at 60 min
	Ni-DLC-1	5	16	0.16	30-3765	773
	Ni-DLC-2	5	13	0.13	30-3765	773
	Ni-DLC-3	4.1	9.5	0.13	60	573-873K
	DLC-4	5.6	9.5	-	60	873

Table 1 : Main characteristics of the Ni/DLC/quartz samples

 1

	SEM images were obtained on a Zeiss GeminiSEM 500 working at 3 keV with secondary
	electrons emission and total yield detection modes, to contrast the topographic from the
	chemical information. Atomic force microscopy (AFM) both on a NT-MDT Stand Alone
	SMENA set-up is used in the resonant semicontact mode. Raman spectroscopy and mapping
	measurements have been carried out using a LabRam/Aramis Jobin/Yvon spectrometer with

Acknowledgements:

The authors would like to thank: T. Fix (Hall measurements), F. Stock (optical transmission analysis); N. Roques, D. Muller (NRA and RBS) at ICube/MaCEPV, Strasbourg; J. Faerber (SEM) and J. Arabski (MBE) at IPCMS, Strasbourg, are acknowledged for their contributions. Part of this research was supported by Exceptional National Program (178 / P.N.E), funded by the Ministry of Higher Education and Scientific Research of Algeria, which is thanked for this fellowship.


Saeed Varasteh Yazdi

Ahlame Douzal-Chouakria

Patrick Gallinari
email: patrick.gallinari@lip6.fr

Manuel Moussallam
email: manuel.moussallam@deezer.com

Time warp invariant dictionary learning for time series clustering: application to music data stream analysis

Keywords: Time series clustering, dictionary learning, sparse coding

This work proposes a time warp invariant sparse coding and dictionary learning framework for time series clustering, where both input samples and atoms define time series of different lengths that involve variable delays. For that, first an l0 sparse coding problem is formalised and a time warp invariant orthogonal matching pursuit based on a new cosine maximisation time warp operator is proposed. A dictionary learning under time warp is then formalised and a gradient descent solution is developed. Lastly, a time series clustering based on the time warp sparse coding and dictionary learning is presented. The proposed approach is evaluated and compared to major alternative methods on several public datasets, with an application to deezer music data stream clustering.

Introduction

Sparse coding and dictionary learning become popular methods in machine learning and pattern recognition for a variety of tasks as feature extraction, reconstruction and classification. The aim of sparse coding methods is to represent input samples as a linear combination of few basis functions called atoms composing a given dictionary. Sparse coding problem is formalised basically as an optimisation problem that minimises the error of the reconstruction under l 0 or l 1 sparsity constraints. The l 0 constraint leads to a non convex and NP-hard problem, that can be solved efficiently by pursuit methods such as orthogonal matching pursuit omp [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation. Information Theory[END_REF]. Relaxing the sparsity constraint from l 0 to l 1 norm yields a convex sparse coding problem, also known as lasso problem [START_REF] Tibshirani | The lasso problem and uniqueness[END_REF]. The dictionary for the sparse representation can be selected among pre-specified family of basis functions (e.g. gabor) or be learned from the training data to sparse represent the input data. A two-step strategy is commonly used: 1) keep the dictionary fixed and find a sparse representation by using a sparse approximation, 2) keep the representation fixed and update the dictionary, either all the atoms at once [START_REF] Engan | Method of optimal directions for frame design[END_REF] or one atom at a time [START_REF] Aharon | k-svd: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. In the context of classification, the dictionary is generally learned from the labeled training dataset and the sparse codes of the testing samples are used for their classification [START_REF] Lee | Efficient sparse coding algorithms[END_REF][START_REF] Mairal | Discriminative learned dictionaries for local image analysis[END_REF][START_REF] Yazdi | Time warp invariant ksvd: Sparse coding and dictionary learning for time series under time warp[END_REF]. In clustering setting, we distinguish two principal sparse coding and dictionary learning approaches. The first category of approaches assumes the data structured into a union of subspaces [START_REF] Bradley | K-plane clustering[END_REF][START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF][START_REF] Tseng | Nearest q-flat to m points[END_REF][START_REF] You | Scalable sparse subspace clustering by orthogonal matching pursuit[END_REF] where each sample may be represented as a linear combination of some input samples that ideally belong to the same subspace. Several of these approaches are related to sparse subspace clustering where samples are first sparse coded based on the input set as a dictionary [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF][START_REF] You | Scalable sparse subspace clustering by orthogonal matching pursuit[END_REF], then a spectral clustering [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] of the sparse codes is used to cluster the data. The number of subspaces as well as their dimension may be fixed beforehand or induced from the affinity graph. In the second category of approaches, a sparse coding and dictionary learning framework is proposed to simultaneously learn a set of dictionaries, one for each cluster, to sparse represent and cluster the data [START_REF] Ramirez | Classification and clustering via dictionary learning with structured incoherence and shared features[END_REF].

For temporal data analysis, sparse coding and dictionary learning are especially effective to extract class specific latent temporal features, reveal salient primitives and sparsely represent complex temporal features. However, what makes temporal data particularly challenging is that salient events may arise with varying delays, be related to a part of observations that may appear at different time stamps. This work addresses the problem of time series clustering under sparse coding and dictionary learning framework, where both input samples and atoms define time series that may involve varying delays and be of different lengths. For that, in the first part, an l 0 sparse coding problem is formalised and a time warp invariant orthogonal matching pursuit based on a new cosine maximisation time warp operator is proposed. Subsequently, the dictionary learning under time warp is formalised and a gradient descent solution is developed. In the second part, a time series clustering approach based on the time warp invariant sparse coding and dictionary learning is proposed. The main contributions of the paper are:

1. We propose a time series clustering approach under sparse coding and dictionary learning setting. 2. We propose a tractable solution for time warp invariant orthogonal matching pursuit based on a new cosine maximisation time warp operator. 3. We provide a sparse representation of the clustered time series and learn, for each cluster, a sub-dictionary composed of the most discriminative primitives. 4. We conduct experiments on several public and real datasets to compare the proposed approach to the major alternative approaches, with an application to deezer music data stream clustering.

The reminder of the paper is organised as follows. Section 2 formalises the time series clustering problem under sparse coding and dictionary learning setting. Section 3 proposes a solution for sparse coding and dictionary learning under time warp, then presents the time series clustering method. Finally, Section 4

presents the experiments and discusses the results obtained.

Problem statement

This Section formalises the time warp invariant sparse coding and dictionary learning for time series clustering. In the following, bold lower case letters are used for vectors and upper case letters for matrices. Let X = {x i } N i=1 be a set of N input time series x i = (x i1 , ..., x iqi ) t ∈ R qi of length q i . We formalise the problem of time series clustering under the sparse coding and dictionary learning setting as the estimation of: a) the partition

C = {C l } K l=1 of X into K clusters and b) the K sub-dictionaries {D l } K l=1
, to minimise the inertia goodness criterion (i.e., the error of reconstruction) as:

min C,D K l=1 xi∈C l E(x i , D l ) (1)
where D l = {d l j } K l j=1 the sub-dictionary of C l is composed of K l time series atoms d l j ∈ R pj . Note that, both input samples x i and atoms d l j define time series of different lengths that may involve varying delays. E(x i , D l ) the error of reconstruction, under time warp, of x i based on the sub dictionary D l = {d l j } K l j=1 is formalised as:

E(x i , D l ) = min αi x i -F i (D l )α l i 2 2 s.t. α l i 0 ≤ τ. (2)
where

F i (D l ) = [f i (d l 1 ), ..., f i (d l K l )] ∈ R qi×K l
is the transformation of D l to a new dictionary composed of warped atoms f i (d l j ) ∈ R qi aligned to x i to resorb the involved delays w.r.t x i . α l i = (α l 1i , ..., α l K l i ) t is the sparse codes of x i under D l and τ the sparsity factor under the l 0 norm.

Proposed solution

To resolve the clustering problem defined in Eq. 1, we use a two steps iterative refinement process, as in standard kmeans clustering. In the cluster assignment step, D l 's are assumed fixed and the problem remains to resolve the sparse coding based on the warped dictionary F i (D l ) defined in Eq. 2. The cluster assignment is then obtained by assigning each x i to the cluster C l whose sub dictionary D l minimises the reconstruction error. In the dictionary update step, the learned sparse codes and the clusters C l are that time fixed and the problem in Eq. 1 defines a dictionary learning problem to minimise the clustering inertia criterion and represent sparsely samples within clusters. For the cluster assignment, we propose in Section 3.1, a time warp invariant orthogonal matching pursuit based on a new cosine maximisation time warp operator. In Section 3.2, a gradient descent solution for dictionary learning under time warp is developed, then the clustering algorithm for time series under sparse coding and dictionary learning setting is given in Section 3.3.

Time warp invariant sparse coding

For the sparse coding under time warp problem given in Eq. 2, we define F i (D l ) as a linear transformation of D l based on the warping function f i (d l j ) = ∆ l ij d l j , where the projector ∆ l ij ∈ {0, 1} qi×pj specifies the temporal alignment that resorbs the delays between x i and d l j . The problem given in Eq. 2 is then formalised as:

min αi, ∆ l i x i - K l j=1 ∆ l ij d l j α l ji 2 2 (3)
s.t. α l i 0 ≤ τ, ∆ l ij ∈ {0, 1} qi×pj , ∆ l ij 1 pj = 1 qi .
with

∆ l i = {∆ l ij } K l j=1 .
The last constraint is a row normalisation of the estimated ∆ l ij to ensure for x i equally weighted time stamps. To resolve this problem, we propose an extended variant of omp that can be mainly summarised in the following steps:

1. For each d l j , estimate ∆ l ij by dynamic programming to maximise the cosine between x i and d l j . 2. Use the projector ∆ l ij to align d l j to x i . 3. Estimate the sparse codes α l i based on the aligned atoms. For that and to estimate ∆ l i , we propose a new operator costw to estimate the cosine between two time series under time warp. To the best of our knowledge, it is the first time that the cosine operator is generalised to time series under time warp. Then, we present a time warp invariant omp (twi-omp), that extends the standard omp approach to sparse code time series under non linear time warping transformations.

Cosine maximisation time warp (costw):

The problem of estimating the cosine between two time series comes to find an alignment between two time series that maximises their cosine. Let x = (x 1 , ..., x qx ), y = (y 1 , ..., y qy ) be two time series of length q x and q y . An alignment π of length |π| = m between x and y is defined as the set of m increasing couples:

π = ((π1(1), π2(1)), (π1(2), π2(2)), ..., (π1(m), π2(m)))
where the applications π 1 and π 2 defined from {1, ..., m} to {1, .., q x } and {1, .., q y } respectively obey the following boundary and monotonicity conditions:

1 = π 1 (1) ≤ π 1 (2) ≤ ... ≤ π 1 (m) = q x 1 = π 2 (1) ≤ π 2 (2) ≤ ... ≤ π 2 (m) = q y and ∀ l ∈ {1, ..., m}, π 1 (l + 1) ≤ π 1 (l) + 1 , π 2 (l + 1) ≤ π 2 (l) + 1, (π 1 (l + 1) -π 1 (l)) + (π 2 (l + 1) -π 2 (l)) ≥ 1 Algorithm 1 M axT riplet(u, v, z) Input: u, v and z. 1: if f (u) ≥ f (v) and f (u) ≥ f (z) then 2: return u; 3: else if f (v) ≥ f (u) and f (v) ≥ f (z) then 4:
return v; 5: else 6:

return z; 7: end if

Intuitively, an alignment π between x and y describes a way to associate each element of x to one or more elements of y and vice-versa. Such an alignment can be conveniently represented by a path in the q x × q y grid, where the above monotonicity conditions ensure that the path is neither going back nor jumping. We will denote A as the set of all alignments between two time series. The cosine maximisation time warp can be formalised as:

costw(x, y) = s(π * ) (4) π * = arg max π∈A s(π) s(π) = |π| i=1 x π1(i) y π2(i) |π| i=1 x 2 π1(i) |π| i=1 y 2 π2(i)
where s(π) is the cost function of the alignment π. The solution for costw is obtained by dynamic programming thanks to the recurrence relation detailed here after.

Let x qx-1 = (x 1 , ..., x qx-1 ), y qy-1 = (y 1 , ..., y qy-1 ) be two sub-time series composed of the q x -1 and q y -1 first elements of x and y, respectively. In the case of aligned time series, that do not include delays and with the same length (i.e., q x = q y ) the following incremental property of the standard cosine can be established:

cos(xq x -1, y qy -1 ) = f (< xq x-1 , y qy -1 >, xq x-1 2 , y qy -1 2 ) cos(x, y) = f ((< xq x-1 , y qy -1 >, xq x -1 2 , y qy -1 2 ) ⊕ (xq x , yq y )) (5)
where

f is a real function defined as f (a, b, c) = a √ b √ c with (b, c ∈ R * + )
and ⊕ is an operator that associates to a triplet (a, b, c) and a couple (u, v) a new triplet as:

(a, b, c) ⊕ (u, v) = (a + uv, b + u 2 , c + v 2 )
For time series including delays and based on the incremental property given in Eq. 5, let us introduce the computation and recurrence relation that allows to estimate the alignment π * that maximises costw(x, y) in Eq. 4.

Computation and recurrence relation: Let us define M ∈ R qx×qy the matrix mapping x and y of general term M i,j = (a i,j , b i,j , c i,j ). Based on the Algorithm 2 twi-omp(x, D, τ )

Input: x, D = {dj} K j=1 , τ Output: α, ∆ 1: r = x, Ω = {φ} 2: while |Ω| ≤ τ do 3:
Select the atom dj (j / ∈ Ω) that maximizes |costw(r, dj)| 4:

Update the set of selected atoms Ω = Ω ∪ {j} and SΩ = [∆jdj]j∈Ω 5:

Update the coefficients: αΩ = (S T Ω SΩ) -1 (S T Ω x) 6:

Estimate the residual: r = x -SΩαΩ 7: end while incremental property established in Eq. 5, computing recursively for (i, j) ∈ {1, ..., q x } × {1, ..., q y } the terms M i,j as:

∀ i ≥ 2, j = 1 M i,1 = (a i-1,1 , b i-1,1 , c i-1,1 ) ⊕ (x i , y 1 ) ∀ j ≥ 2, i = 1 M 1,j = (a 1,j-1 , b 1,j-1 , c 1,j-1 ) ⊕ (x 1 , y j ) and ∀ i ≥ 2, j ≥ 2 M i,j = M axT riplet      (a i,j-1 , b i,j-1 , c i,j-1 ) ⊕ (x i , y j ) (a i-1,j , b i-1,j , c i-1,j ) ⊕ (x i , y j ) (a i-1,j-1 , b i-1,j-1 , c i-1,j-1 ) ⊕ (x i , y j ) and M 1,1 = (x 1 y 1 , x 2 1 , y 2
1 ), we obtain costw(x, y) = f (M qx,qy ) with a quadratic complexity of O(q x q y ). The two first equations give the first row and column updates, the third equation gives the recurrence formula that ensures the cosine maximisation at each M i,j cell and M axT riplet function (Algorithm 1) retains the triplet that maximises the cosine at M i,j .

Time warp invariant OMP (twi-omp):

Based on the defined costw, let us present the time warp invariant omp (twi-omp) to sparse code a given time series x based on a dictionary D = {d j } K j=1 under time warp conditions (Algorithm 2). The proposed twi-omp follows the three steps given in the previous section. First, perform a costw between x and each d j to estimate ∆ = {∆ j } K j=1 and find the atom d j that maximises costw(x, d j ) (line 3 in Algorithm 2). Then, update the set Ω of the yet selected projected atoms and the dictionary S Ω = [∆ j d j ] j∈Ω of the yet selected warped atoms (line 4). The updated S Ω is then used to estimate the sparse coefficients of x (line 5-6). The process is reiterated on the residuals of x until the sparsity factor τ is reached.

Time warp invariant dictionary learning

For the dictionary learning step, the problem in Eq. 1 becomes to learn the dictionary D under time warp where, that time, the sparse codes α l i and ∆ l i are assumed fixed as:

min D K l=1 xi∈C l x i - K l j=1 ∆ l ij d l j α l ji 2 2 s.t. d l j 2 = 1. ( 6
)
This problem is then resolved as K single dictionary learning problems to learn each sub-dictionary D l that minimises the inertia of the cluster C l :

J l = min D l xi∈C l x i - K l j=1 ∆ l ij d l j α l ji 2 2 (7)
which is equivalent to

J l = min D l xi∈C l qi t=1 (x it - K l j=1 α l ji (t,t )∈π * ij d l jt ) 2 (8)
s.t. d l j 2 = 1
where x it is the t th time instant of x i and π * ij denotes the optimal alignment path between x i and d l j . To resolve the Eq. 8, we propose a gradient descend method based on the following update rule at iteration m for the atom d l j :

d l(m+1) jt = d l(m) jt -η m ∂J l ∂d l(m) jt (9)
d l(m+1) j = d l(m+1) j d l(m+1) j 2
with,

∂J l ∂d l jt = x i ∈C l q i t=1 -2α l ji (xit -α l ji d l jt -α l ji (t,t )∈π * ij (t =t ) d l jt ( 10
) - j =j α l j i (t,t )∈π * ij d l j t )
where η is the learning rate. In the following section, we show how the time warp invariant omp and dictionary learning are involved for time series clustering.

Time warp invariant dictionary learning for time series clustering

For time series clustering, the clustering criterion given in Eq. 1 is minimised by an iterative process involving, respectively, time warp invariant sparse coding (twi-omp) and dictionary learning for cluster assignments and dictionary Algorithm 3 twi-dlclust(X, K, τ )

Input: X = {xi} N i=1 , K, τ . Output: {C1, ..., CK }, {D1, ..., DK } 1: {Clustering Initialisation:} 2: Define the affinity matrix S ∈ R N ×N of general term: 3:

s ii =costw(xi, x i ) 4: Apply the affinity propagation (or spectral clustering) to cluster S into 5: K clusters: C1, ..., CK 6: {Sub-dictionary initialisation:} 7: for l = 1, ..., K do 8:

Initialise D l randomly from C l 9:

repeat 10: Sparse code each xi ∈ C l : [α l i , ∆ l i ] =twi-omp(xi, D l , τ ) 11:
Update each d l j ∈ D l by using Eq. 9 and 10. 12:

until Convergence (stopping rule) 13: end for 14: repeat 15:

{Cluster assignment:} 16:

Sparse code each xi ∈ X based on each D l (l = 1, ..., K):

17: [α l i , ∆ l i ] =twi-omp(xi, D l , τ ) 18:
Assign xi to the cluster C l whose D l minimises E(xi, D l ): 19:

C l = {xi / l = min l xi - K l l =1 ∆ l ij d l j α l ji 2 2 } 20:
{Dictionaries update:} 21:

for l = 1, ..., K do 22:

Update each d l j ∈ D l by using Eq. 9 and 10.

23:

end for 24: until Convergence (no changes in cluster assignments) update steps (Algorithm 3). In the initialisation step, a clustering (e.g., spectral clustering, affinity propagation) is performed on the costw matrix S to determine an initial partition {C l } K l=1 of X (line 1-5). A sparse coding and a dictionary learning are then performed on the samples of each cluster to initialise the sub-dictionaries {D l } K l=1 (line 6-13). Based on the initial partition {C l } K l=1 and sub-dictionaries {D l } K l=1 , the cluster assignment step consists to perform a sparse coding of each input sample based on each dictionary D l , then to assign it the cluster whose dictionary minimises its reconstruction error (line 15-19). Subsequently, in the dictionary update step, the atoms d l j of each dictionary are updated by using the formula given in Eq. 9 and 10 (line 20-23).

Experimental study

In this section, we evaluate the proposed time series clustering under dictionary learning setting (twi-dlclust) on several synthetic and real datasets, including multivariate and univariate time series, that may involve varying delays and be of different or equal lengths. The proposed twi-dlclust clustering method is compared to two major alternative approaches, the subspace sparse clustering (ssc) [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF] and the Dictionary Learning with Structured Incoherence (dlsi) [START_REF] Ramirez | Classification and clustering via dictionary learning with structured incoherence and shared features[END_REF]. For ssc, two variants ssc-bp [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF] and ssc-omp [START_REF] You | Scalable sparse subspace clustering by orthogonal matching pursuit[END_REF] are studied for a sparse coding under l 0 and l 1 norms, where an orthogonal matching pursuit and a basis pursuit methods are used respectively. For dlsi, both sample-based and atom-based affinity matrix initialisations proposed in [START_REF] Ramirez | Classification and clustering via dictionary learning with structured incoherence and shared features[END_REF] are studied. The Matlab codes of these methods are available online4 .

Data description

We have considered in Table 1 two groups of datasets. The first group is composed of the top 12 datasets for which the ground truth clustering is given. The four first datasets are composed of public multivariate time series that have different lengths and involve varying delays. In particular, digits, lower, and upper datasets give the description of 2-D air-handwritten motion gesture of digits, upper and lower case letters performed on a Nintendo (R) Wii device by several writers [START_REF] Chen | 6dmg: A new 6d motion gesture database[END_REF]. The char-traj dataset gives the 2-dimensional handwritten character trajectory performed on a Wacom tablet by the same user [START_REF] Frank | Uci machine learning repository[END_REF]. The ecg-mit dataset was obtained from the mit-bih Arrhythmia [START_REF] Goldberger | Physiobank, physiotoolkit, and physionet[END_REF] database where the heartbeats represented by qrs complexes. The 7 remaining datasets are composed of univariate time series of the same lengths that involve significant delays [START_REF] Keogh | The ucr time series data mining archive[END_REF]. The last two datasets are provided by deezer5 , the online music streaming service that offers access to the music content of nearly 40 million licensed tracks. deezer data, for which we have no ground truth, give the description of streaming data of music albums, randomly selected among 10 5 French user streams and recorded from October 2016 to September 2017. They are composed of univariate time series that give the daily total number of streams per album from its release date to September 2017; this study consider only the streams of a duration ≥ 30 seconds. In particular, deezer15 and deezer30 are provided for the streams analysis over the crucial early period after the album release date. They give the description of the prefix time series on the early period covering a cumulative number of 10 3 streams (in red in Figure 1). In addition, for the pertinence of the analysis, the prefix time series of length < 7 days are extended to 15 days in deezer15 and to 30 days in deezer30. Table 1 gives some characteristics of the studied datasets: the size of the clusters when the ground truth is known, the size of the validation and evaluation sets and the length of the time series that may be variable or fixed.

Validation protocol

For the top 12 datasets in Table 1, for which the ground truth partition is known, the proposed method twi-dlclust as well as the alternative clustering approaches are applied to cluster the data. For alternative methods, time series of different lengths are zero padded beforehand. The adjusted Rand index [START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF] is then used to evaluate the goodness of the obtained clusterings. The Rand index lies between 0 and 1, it measures the agreement between the obtained clusters and the ground truth ones. The higher the index, the better the agreement is.

In particular, the maximum value "1" of the Rand index is reached when the obtained partition and the ground truth one are identical. For deezer datasets, the ground truth being unknown, a dtw-based within-class W r ratio 6 is used.

The lower the within-class ratio W r , the better the clustering is. W r is as well used to select the optimal number of clusters. Finally, the parameters related to each studied method, indicated in Table 2, are learned by line/grid search on the validation set, the best parameters are then used to perform the clustering

6 Wr = K l=1
x,y∈C l dtw(x,y) x,y∈X dtw(x,y) λ ∈ {0.001, 0.01}, lag of 0.01 l1 sparsity regularisation dlsi λ ∈ {0.001, 0.01}, lag of 0.01 l1 sparsity regularisation η ∈ {0, 0.1, 0.01} dictionary incoherence regularisation K l = 5, ∀ l ∈ {1, ..., K} Sub-dictionary D l size twi-dlclust sc ∈ [0, 100], lag of 10 Sakoe-Chiba band width τ ∈ {1, 2, 3, 4, 5} l0 sparsity threshold

K l = 5, ∀ l ∈ {1, ..., K} Sub-dictionary D l size
on the evaluation set. The process is iterated over 10 runs and the averaged performances are reported in Tables 3 and4.

Results and discussion

Table 3 gives for the top 12 datasets the obtained adjusted Rand index values. The best values are indicated in bold, the non significantly different ones from the best (t-test at 5% risk) are in italic and the remaining results are significantly different from the bold values. For the two l 0 sparse coding methods ssc-omp and twi-dlclust, the learned sparsity coefficient τ is given between brackets. The two last rows give, over all the datasets, the number of times a method reaches the best value as well as its average ranking. From Table 3, we can see that the proposed twi-dlclust reaches the best clustering results with 9 times (9 out of 12) as the best values, 2 times as significantly non different from the best and obtained the lowest average ranking. The second best results are obtained by ssc-bp and dlsi-s, followed by ssc-omp. Although the l 1 sparse coding models (here ssc-bp and dlsi-s) are known to be more efficient than the l 0 models, twi-dlclust even involving an l 0 sparse coding leads to the best results. While twi-dlclust and dlsi-s involve smaller size sub-dictionaries (K l = 5), ssc-omp and ssc-bp are based on larger dictionary of the size of the evaluation set (Table 1). Finally, by comparing the two l 0 sparse coding methods ssc-omp and twi-dlclust, we can see that twi-dlclust leads for all datasets to sparser solutions with a lower or equal sparsity coefficient τ than ssc-omp.

For deezer data we have performed each clustering method for several number of clusters and the within-class ratio of the obtained partitions reported in Table 4. For simplicity, the dlsi approach is conducted only with dlsi-s variant, dlsi-a being highly equivalent in Table 3. We can see easily that, for both datasets and almost all the number of clusters, the best values are reached by twi-dlclust, followed by the l 1 sparse code approaches ssc-bp and dlsi-s, then by ssc-omp. Finally, note that from both Tables 3 and4, ssc-omp and ssc-bp lead to the lowest performances with a slightly better results for ssc-bp as using an l 1 norm sparse coding. These results may be partly explained by the fact that both ssc-omp and ssc-bp are purely sparse coding methods based on one global dictionary fixed beforehand, unlike dlsi and twi-dlclust that learn one sub-dictionary per cluster.

In the second study, we analyse more closely the obtained clusterings. For instance, based on Figure 2 that displays the progression of the within-class ratio w.r.t the number of clusters, a partitioning into four clusters is performed on deezer30. Accordingly, Figure 3, shows for each of the four clusters (each row), the profile of the medoids (in the first column), the closest albums to the medoid in the second column and at the third column, the atom that most contributes to sparse represent the cluster's samples. deezer data provide additional album descriptive features as "Full" (composed of several tracks) or "Single" (composed of one track), if it is a "Deluxe" edition, namely a re-edition of the album featuring extra contents related to the album, as well as the artist popularity before and after the release. The analysis of some album characteristics brings meaningful interpretation of the extracted clusters (Figure 3).

The first cluster is composed of 71% of "Full" albums and 15% of "Deluxe" editions. It corresponds to album releases with flat stream profiles. Such behaviour usually occurs when the content has already been published ("Deluxe" versions) or for lesser-known artists, as assessed by the cluster medoid "Empereur du Sale" album of the rapper "Lorenzo" that released several singles a few weeks before the album release date and although not highly popular has still a steady fan base.

In the second cluster, 75% of the albums are "Single". The fast decrease stream profile just after the release date is not surprising for short albums (composed of 1-4 tracks). Indeed, a "Full" album is generally released shortly after the "Single" release, inducing a decrease of streams for the "Single" few weeks after its release. The cluster medoid "Ethologie" is produced by the rapper "Dehmo" that has not released albums since a long period, that may explain the burst of streams for the new content just after its release.

The cluster 3 is composed of 69% of "Full" albums mainly produced by artists that became popular after their album release. This is reflected by the stream profiles that initially evolve at low level then increase significantly several days/weeks after the album release. This is confirmed by the medoid album "Be Mine" a single produced by "Ofenbach" that was in fact revealed to the public with that album.

Finally, the cluster 4 comprises a majority of "Single" albums (84%) produced by very popular artists with a huge fan base and immediate success. The medoid album "Divide" produced by "Ed Sheeran" was one of the biggest hits of 2017. Although the stream profiles of the clusters 4 and 2 seem similar, albums of cluster 4 concern more established artists in their second/third albums while cluster 2 is more related to emerging works and first successes.

The aim of the last study is to analyse the pertinence of the learned subdictionaries {D 1 , ..., D K } for both dlsi and twi-dlclust; the dictionary for the other methods ssc-omp and ssc-bp is not learned but fixed beforehand. For that, for each method dlsi and twi-dlclust, the atoms of the learned subdictionaries are gathered together to built one global dictionary ∪ K l=1 D l . Let us denote D G1 and D G2 the global dictionaries obtained for dlsi and twi-dlclust, respectively. Subsequently, the samples in X are sparse coded, by using first an l 1 norm regularisation based on D G1 , then a twi-omp based on D G2 . For instance, for deezer30, Figure 4 shows for the 278 albums the learned sparse codes based on D G1 (on left) and on D G2 (on right), organised for interpretation purpose per cluster {C 1 , ..., C 4 } and per sub-dictionary {D 1 , ..., D 4 }. It emerges from Figure 4, that sparse codes based on D G2 highlight clearly a block structure that reflects the discriminative performance of the sub-dictionaries composing D G2 (learned by twi-dlclust). Indeed, sparse codes show that each sub-dictionary D l is mainly involved to reconstruct samples of C l . On the other hand, the structure of the sparse codes based on D G1 (learned by dlsi) seems much less sparser and less discriminative. We can note, in particular, that the atoms d

Conclusion

This work proposes a time warp invariant sparse coding and dictionary learning for time series clustering where both input samples and atoms define time series that may have different lengths and involve varying delays. For that, first a time warp invariant orthogonal matching pursuit based on a new cosine maximisation time warp operator is proposed. Then, a dictionary learning approach under time warp is formalised and a gradient descent solution is developed. The proposed time series clustering allows to sparse represent the clustered time series and learn, for each cluster, a sub-dictionary composed of the most discriminative primitives. The conducted experiments show that although twi-dlclust involves an l 0 sparse coding approach based on a very small size sub-dictionaries, it leads to the sparser and the best clustering results, while revealing atoms with a good discriminative capacity to represent the time series of each cluster.

Fig. 1 .

 1 Fig. 1. An album streaming time series, in red the prefix time series covering a cumulative number of 10 3 streams.

Fig. 2 . 1 Fig. 3 .

 213 Fig. 2. Number of clusters K vs. Within-class ratio Wr.

Fig. 4 .

 4 Fig. 4. Sparse representations based on: DG1 learned by dlsi (left) and DG2 learned by twi-dlclust (right).

  primitives involved to reconstruct the samples of all the clusters.

Table 1 .

 1 Data description

	Dataset	Nb. class Valid. set Eval. set Length range
			size	size
	digits	10	100	100	29∼218
	lower	26	130	260	27∼163
	upper	26	130	260	27∼412
	char-traj	20	100	200	109∼205
	ecg-mit	4	40	160	541
	cbf	3	30	900	128
	facefour	4	24	88	350
	lightning2	2	60	61	637
	lightning7	7	70	73	319
	cc	6	300	300	60
	trace	4	100	100	275
	ecg200	2	100	100	96
	deezer15	-	-	281	15∼301
	deezer30	-	-	278	30∼301

Table 2 .

 2 Parameter Line/Grid values

	Method	Line / Grid values	Desc.
	ssc-omp	τ ∈ {1, 2, 3, 4, 5}	l0 sparsity threshold
	ssc-bp

Table 3 .

 3 Adjusted Rand index

	Dataset	ssc-omp (τ ) ssc-bp dlsi-s dlsi-a twi-dlclust (τ )
	digits	0.839 (2)	0.856	0.854	0.841	0.940 (1)
	lower	0.935 (3)	0.943	0.937	0.934	0.970 (1)
	upper	0.940 (2)	0.942 0.940	0.938	0.942 (1)
	char-traj	0.947 (5)	0.977	0.978 0.971	0.965 (3)
	ecg-mit	0.327 (2)	0.789	0.772	0.773	0.792 (2)
	cbf	0.558 (2)	0.668	0.599	0.601	0.770 (2)
	facefour	0.810 (5)	0.722	0.767	0.769	0.776 (3)
	lightning2	0.559 (2)	0.559 0.559 0.519	0.559 (2)
	lightning7	0.793 (2)	0.808	0.724	0.747	0.814 (3)
	cc	0.736 (5)	0.630	0.813	0.791	0.910 (1)
	trace	0.680 (5)	0.752	0.755	0.753	0.805 (1)
	ecg200	0.547 (4)	0.631 0.689 0.664	0.653 (3)
	Nb. Best	2	2	3	0	9
	Avg. Rank	4.00	2.83	2.92	3.58	1.67

Table 4 .

 4 Within-class ratio Wr per number of clusters K.

	Dataset	K ssc-omp ssc-bp dlsi-s twi-dlclust
			(τ = 5)			(τ = 2)
	deezer15 2	0.266	0.310 0.245	0.262
		3	0.201	0.253 0.177	0.145
		4	0.212	0.146 0.114	0.112
		5	0.188	0.118 0.112	0.096
		6	0.133	0.106 0.074	0.069
	deezer30 2	0.339	0.348 0.322	0.303
		3	0.226	0.292 0.273	0.173
		4	0.175	0.241 0.133	0.127
		5	0.154	0.119 0.110	0.096
		6	0.100	0.080 0.085	0.085
	Nb. Best		0	1	1	8
	Avg. Rank		3.40	3.30	2.05	1.25

ssc-omp: https://goo.gl/E6khsq, ssc-bp: https://goo.gl/719pvx and dlsi: https://goo.gl/X5nZgE.

https://www.deezer.com/fr/

Acknowledgment

This work is supported by the French National Research Agency (ANR-Locust project) and Bpifrance funds in the frame of the French National PIA Program.