Fatemehsadat Moosavi

Alex Lemarchand

Cyrille Bazin

Maria Konstantakopoulou

Mathieu Frégnaux

Chris- Tophe Colbeau-Justin

Touraj Tavakoli Gheinani

Andrei Kanaev
email: andrei.kanaev@lspm.cnrs.fr

Mehrdad Nikravech

Photocatalytic nanocomposite anatase-rutile TiO 2 coating

Keywords: nanocomposite TiO2 coating, anatase, rutile, adherence, photocatalysis

Nanocomposite TiO2 coatings of mixed anatase-rutile crystalline polymorphs with strong adherence to glass substrates were prepared via spray-plasma deposition method. The rutile content was varied between 0 and 30 mol%. The coating showed an enhanced photocatalytic activity towards cefixime antibiotic decomposition in aqueous solutions under UVA (=365 nm) light illumination. The activity significantly increased after the rutile phase formation, suggesting a synergetic effect of the two-component solids.

Introduction

Titanium dioxide (TiO2) has deserved large interest in photocatalysis, becoming nowadays a key material in the environmental protection [START_REF] Fujishima | TiO2 photocatalysis: Fundamentals and applications[END_REF][START_REF] Hashimoto | TiO2 photocatalysis: A historical overview and future prospects[END_REF][START_REF] Serpone | Semiconductor photocatalysispast, present, and future outlook[END_REF]. However, studies of its optimal composition under UV and sun illuminations still continue. Advances in preparation methods of photcatalytic titania coatings have been recently reviewed [START_REF] Dell'edera | Photocatalytic TiO2-based coatings for environmental application[END_REF]. Between different polymorphs, anatase TiO2 with the bandgap of 3.2 eV generally showed a higher activity and rutile TiO2 with the bandgap energy of 3.0 eV possesses a lower activity, which has been mainly related to the radicals generation at respective (101) and (110) thermodynamically stable crystalline facets [START_REF] Liu | Nanosized anatase TiO2 single crystals for enhanced photocatalytic activity[END_REF][START_REF] Sun | Evaluating intrinsic photocatalytic activities of anatase and rutile TiO2 for organic degradation in water[END_REF][START_REF] Ahmed | Photocatalytic activities of different well-defined single crystal TiO2 surfaces: anatase versus rutile[END_REF]. Attempts to increase of the adsorption capacity of rutile TiO2 have been undertaken, conducting to a significant enhancement of the photocatalytic activity [START_REF] Holm | Impact of rutile and anatase phase on the photocatalytic decomposition of lactic acid[END_REF]. In the same time, one of the best photocatalytic performances showed a commercial composite TiO2, Degussa P25, which up to now is considered as the reference photocatalytic material. This material appeared in a powder form, composed of agglomerated nanoparticles of 15-30 nm size of mixed anatase (~70 %) and rutile (~30 %) polymorphs with a small addition of an amorphous material [START_REF] Bickley | A structural investigation of titanium dioxide photocatalysts[END_REF][START_REF] Ohtani | What is Degussa (Evonik) P25? Crystalline composition analysis, reconstruction from isolated pure particles and photocatalytic activity test[END_REF]. Previous studies [START_REF] Hurum | Explaining the enhanced photocatalytic activity of Degussa P25 mixed-phase TiO2 using EPR[END_REF][START_REF] Carneiro | How phase composition influences optoelectronic and photocatalytic properties of TiO2[END_REF] have revealed a significant inhibition of the electron-hole recombination process in this material, which has a consequence of the photocata-lytic activity enhancement [START_REF] He | Facile formation of anatase/rutile TiO2 nanocomposites with enhanced photocatalytic activity[END_REF][START_REF] Liu | Engineering of anatase/rutile TiO2 heterophase junction via in-situ phase transformation for enhanced photocatalytic hydrogen evolution[END_REF][START_REF] Zhang | Enhanced photocatalytic degradation of gaseous toluene and liquidus tetracycline by anatase/rutile titanium dioxide with heterophase junction derived from materials of Institut Lavoisier-125(Ti): Degradation pathway and mechanism studies[END_REF]. It has been also suggested that small particles of size ~10 nm possess a higher activity compared to larger particles [START_REF] Jang | Effect of particle size and phase composition of titanium dioxide nanoparticles on the photocatalytic properties[END_REF][START_REF] Marugan | Optical density and photonic efficiency of silicasupported TiO2 photocatalysts[END_REF][START_REF] Tieng | New homogeneously doped Fe(III)-TiO2 photocatalyst for gaseous pollutant degradation[END_REF], which corresponds to Degussa P25 TiO2 nanopowders.

Despite of much progress in understanding of the functional properties, large scale applications of this composite-phase photocatalyst is rather limited because of the severe problems related to environmental pollution by photocatalytic nanoparticles. The filtering of effluents via filtration from suspended nanoparticles in a liquid phase is generally expensive and of insufficient effectiveness. Furthermore, a tendency of semiconducting nanoparticles to aggregation decrease both light coupling efficiency and charge transfer efficiency to active surface sites, which finally worsen their functional response. The depositions of Degussa P25 TiO2 nanopowders on appropriate supports, though being successful [START_REF] Fanou | Elaboration of novel nanoparticulate TiO2-P25@n-TiO2 composite for photocatalysis[END_REF], have not resulted in mechanically strong coatings. Therefore, elaboration of photocatalytic nanocoatings with the structural composition equivalent to that of Degussa P25 remains a challenging task.

In this work we developed a preparation method of photocatalytically active coatings with crystalline grain size ~10 nm and unique anatase-rutile phase composition, which can be adjusted. These coatings were mechanically stable, forming strong covalent bonds with supports, which

Experimental

The TiO2 films were deposited on glass beads by using the low-pressure spray plasma method [START_REF] Baba | Growth of ZnO thin films by spray plasma technique: Correlation between spectroscopic measurements and film properties[END_REF][START_REF] Nikravech | Spray plasma device, a new method to process nanostructured layers. Application to deposit ZnO thin layers[END_REF]. Before operation, the reaction chamber was maintained under vacuum of 20 Pa by a continuous pumping. Microscopic droplets of about 0.5 m size of a titanium tetra isopropoxide (Sigma-Aldrich) solution in 2-propanol (Sigma-Aldrich) with volume ratio 1:1 were produced by an ultrasonic nebulizer Mist Maker DK24 (ServoVendi) operating at 1 MHz frequency and injected into the reaction chamber via argon carrier gas at a flow rate of 300 ml/min. A special aerosol conditioner was used to introduce aerosol from atmospheric pressure into the low pressure plasma reactor. The oxygen gas as a reactive component was admixed to the argon gas flow at a rate varying between 0 and 25 ml/min. A plasma discharge was initiated in the reaction chamber by using an inductive radiofrequency generator (Hüttinger Trumpf Qinto 3013, 3 kW, 13.5 MHz) of 300 W power. The glass beads, preliminarily cleaned in sulfuric acid and distilled water, were positioned in a fixed bed inside the plasma reactor chamber and the deposition time was fixed to 60 min. In order to improve the crystallinity, the coatings were post-treated thermally in a furnace RHTC-80230115 (Nobertherm) at 450 °C during 4 hours.

The x-ray diffraction (XRD) patterns were collected at a grazing incidence using Panalytical Empyrean diffractometer equipped with Cu-K source (λ=1.541 Å). The xray incident beam was collimated with 1/16° slits coupled with a 5 mm mask to ensure that the beam footprint did not exceed the sample surface. The diffractometer was equipped with PIXcel1D 255-channels detector operating in scanning line mode coupled with a parallel plate collimator (acceptance angle of 0.18°). The patterns were collected with an incident angle of 0.5° in the range from 15° to 70° (2θ) with a step of 0.03° and counting time per step of 75 s. The phase quantification was estimated from the Grazing Incidence XRD (GIXRD) patterns with rutile/anatase ratio calculated from the respective (110) and (101) peak intensities via equation Fr=1-(1+1.265•Ir110/Ia101) -1 [START_REF] Spurr | Quantitative analysis of anataserutile mixtures with an x-ray diffractometer[END_REF]. The high-resolution transmission electron microscopy (HRTEM) was performed using JEOL 2011 equipment operating at 200 kV with LaB6 electron emission source. The images were recorded using Gatan MultiScan CCD camera positioned at the output of Gatan Imaging Filter system. Samples were prepared by tearing off some pieces of the thin film and depositing on holey carbon film supported 400 mesh copper grids. The interplanar lattice spacing of TEM patterns was analyzed with ImageJ 1.53k software [23]. The XPS analysis was realized on Thermofisher Scientific Escalab 250xi equip-ment with spot size of 650 µm and flood gun assistance (Al-K, =1486.6 eV, Pass Energy (survey/HR) 100 eV / 20 eV, dwell time 100 ms, step (survey/HR) 1 eV / 0.1 eV). Spectra were then charge corrected by shifting all peaks to the adventitious C1s spectral component (C-C, C-H) binding energy set to 284.8 eV. Quantification is performed based on the photopeak areas after a Shirley type background subtraction using the Thermofisher Scientific Avantage© software and its "ALTHERMO1" library as sensitivity factor collection.

The time-resolved microwave conductivity (TRMC) measurements providing lifetime of photoinduced charges were realised at the experimental installation described in Ref. [START_REF] Colbeau-Justin | Structural influence on charge-carrier lifetimes in TiO2 powders studied by microwave absorption[END_REF]. The method is based on a variation of the microwave power reflected by a semiconductor irradiated with a pulsed laser, inducing changes of the electrical conductivity. The samples were irradiated with 360-nm light pulses of 8 ns duration and 1.3 mJ energy delivered by cw OPO Nd:Yag laser (EKSPLA) operating at 10 Hz. The irradiated area on a sample was 5 mm in diameter. The measurements were conducted at room temperature under ambient atmosphere. The recorded signal was accumulated over 200 laser pulses.

The photocatalytic activity of the prepared materials was evaluated by the degradation of cefixime antibiotic in aqueous solutions under UV-lamp illumination in a continuous-flow fixed-bed reactor described in Ref. [START_REF] Mendez | Photocatalytic activity of nanocoatings based on mixed oxide V-TiO2 nanoparticles with controlled composition and size[END_REF]. In brief, 10 mg/l cefixime solution in distilled water was placed in a reservoir and pushed by a peristaltic pump into the 6-mm gap between two coaxial quartz reactor tubes, filled with 100 g of the coated glass beads of 3-mm diameter, forming a packing bed of about 6 cm height. The light source, UVA lamp Phillips TL BLB with 8 W power emitting at 365 nm with full-width at half maximum 1/2=16 nm, of a cylindrical geometry was inset in the inner reactor tube maintained at 25 °C by circulating the gas flow through a heat exchanger. The pollutant concentration was periodically measured by absorption at 286 nm, using a compact optical fiber coupled PC-plugged AvaSpec-Mini UV/VIS/NIR Spectrometer (spectral resolution of 0.5 nm). Before illumination, the pollutant flow was passed through the reaction medium during 4 hours in order to attain the adsorption-desorption equilibrium of the pollutant molecules on the photocatalyst. Then the UVA lamp was set "on" and the pollutant concentration measurements were continued under illumination. Blank tests were realized in the empty reactor in order to evaluate the photolytic decomposition of pollutant by lamp photons. The reproducibility of the experimental series was verified. We notice that the applied monitoring degradation kinetic was used to evaluate the photocatalyst activity and does not indicate that cefixime molecules were totally oxidized into small inorganic species such as H2O, CO2, etc., which can be inferred using Total Organci Carbon (TOC) analysis. The intermediate cefixime oxidation products can persist in our case after the irradiation period. Such intermediate products of cefixime photodegradation have been recently reported [START_REF] Zhao | Green and efficient degradation of cefixime by 3D flower-like BiOBr: Performance and degradation pathway[END_REF].

Results and discussion

Table 1 summarizes the preparation conditions and resulting crystalline composition of the prepared coatings. The respective structural data were obtained from the analysis of GIXRD patterns, some of which are shown in Fig. 1. The oxygen gas flow in the plasma appeared to be one of key factors in the elaboration of composite anatase-rutile TiO2 coatings. Indeed, in the absence of oxygen almost pure anatase phase was observed in sample D1 (Fig. 1c), while an increase of the oxygen flow resulted in a strong increase of the rutile content up to 16 mol% in sample D2 (Fig. 1b) and 27 mol% in samples D3 (Fig. 1c). At even higher doubled oxygen content, the rutile phase insignificantly increased to ~30 mol%. The characteristic size of the crystalline domains, of both anatase and rutile polymorphs, determined by Scherrer relation was very similar of 91 nm.

The HRTEM images confirm these results, showing well crystallised materials with the mean grain size about 10 nm: anatase and rutile domains coexist in the appropriate samples (example of D4 is shown in Fig. 2). The appearance of rutile phase after heat treatment at 450 °C was unexpected. In fact, previous studies have shown that rutile appeared in single TiO2 nanoparticles of size larger than 14 nm [START_REF] Zhang | Thermodynamic analysis of phase stability of nanocrystalline titania[END_REF][START_REF] Bouslama | New nanoparticulate TiO2-Al2O3 photocatalytic media: Effect of particle size and polymorphism on photocatalytic activity[END_REF] and in agglomerated nanopowders at temperatures above 550 °C (with the heat treatment duration of 4 h) [START_REF] Bouslama | New nanoparticulate TiO2-Al2O3 photocatalytic media: Effect of particle size and polymorphism on photocatalytic activity[END_REF][START_REF] Zhang | UV Raman spectroscopic study on TiO2. I. Phase transformation at the surface and in the bulk[END_REF]. The rutile phase appearance is suggested to occur through the interface nucleation at a contact area between two anatase nanoparticles in a particular orientation [START_REF] Zhang | New kinetic model for the nanocrystalline anatase-to-rutile transformation revealing rate dependence on number of particles[END_REF]. More favourably, it is initiated at the twin (112) anatase interface, which involves structural elements common to rutile [START_REF] Penn | Formation of rutile nuclei at anatase 112 twin interfaces and the phase transformation mechanism in nanocrystalline titania[END_REF]; this transformation requires lower activation energy compared to that at the surface or in the bulk. On the other hand, the mixed anatase and rutile phase composition has been reported in films deposited with pulsed micro-plasma cluster source [START_REF] Barborini | Engineering the nanocrystalline structure of TiO2 films by aerodynamically filtered cluster deposition[END_REF]. A variation of the stoichiometry has been disregarded in these experiments because of the strictly homogeneous elemental condition O:Ti=2:1. As a result, these authors have suggested that the structure of the very small titania clusters is predisposed for the rutile formation while that of larger particles (of size above 5 nm) converges to anatase crystallites. A control of the crystalline anatase/rutile composition could be therefore achieved via selection of cluster during the deposition process.

Despite of the oxygen gas insertion in the plasma chamber, its impact on the films stoichiometry can be neglected in our experiments. Indeed, the XPS elemental analysis of the film composition of D1 and D3 samples resulted in the same ratio of O/Ti=1.97, which can be considered as stoichiometric. Most probably, and in agreement with earlier studies [START_REF] Barborini | Engineering the nanocrystalline structure of TiO2 films by aerodynamically filtered cluster deposition[END_REF][START_REF] Nikravech | Role of reactive species in processing materials at laboratory temperature by spray plasma devices[END_REF], the oxygen content affects the size and surface state of the condensed species coming from plasma and impinging substrates. The influence of plasma oxygen on physicochemical properties of oxide coatings has been previously reported [START_REF] Nikravech | Spray plasma device, a new method to process nanostructured layers. Application to deposit ZnO thin layers[END_REF][START_REF] Zhang | UV Raman spectroscopic study on TiO2. II. Effect of nanoparticle size on the outer/inner phase transformations[END_REF]. The polydispersity and/or inhomogeneity of the condensed species controlled by the oxygen content apparently define the anatase/rutile ratio of prepared films.

An interesting observation is however the formation of similar sizes of anatase and rutile crystallites evidenced by GIXRD and TEM data. This makes doubtful the proposed cluster size effect on TiO2 polymorph formation [START_REF] Barborini | Engineering the nanocrystalline structure of TiO2 films by aerodynamically filtered cluster deposition[END_REF], since smaller clusters and bigger nanoparticles are not expected to result in the same size entities after crystallisation. We believe that the surface state is responsible for the rutile phase appearance. In fact, the standard free energy change is expressed by

G 0 =fG 0 (T)2-fG 0 (T)1+A22-A11 (1)
where fG 0 is the standard free energy of two phases formation and A and  are respectively surface area and surface free energy of the particles. The condition of G 0 =0 enables the phase transformation 12 (anataserutile) and it obviously depends on the surface state and size of a particle and temperature T [START_REF] Zhang | Thermodynamic analysis of phase stability of nanocrystalline titania[END_REF]. Low pressure plasmas, under which the layers were elaborated, are known as thermodynamically nonequilibrium media. It is noteworthy that in these conditions oxygen molecules undergo dissociative reactions leading to generation of reactive oxygen radicals [START_REF] Baba | Growth of ZnO thin films by spray plasma technique: Correlation between spectroscopic measurements and film properties[END_REF][START_REF] Nikravech | Role of reactive species in processing materials at laboratory temperature by spray plasma devices[END_REF]. The exothermic recombination reactions, particularly those of oxygen radicals on the surface of deposit, could also enhance local temperature leading to the rutile phase formation on some sites while the surroundings remain cold. More oxygen radicals are produced in plasma, more recombination reactions take place and, consequently, more rutile crystallites could be produced, in accordance with the correlation between rutile contents and oxygen gas flow observed in our experiment (Table 1).

The Ti2p, O1s and Si2p XPS spectra of coatings, one example of which is given in Fig. 3, were quite similar. The crystalline composition of the prepared coatings is summarized in Table 1. The Ti2p core level peaks (23.3 at.%) with binding energies of 458.9 eV (2p3/2) and 464.6 eV (2p1/2) and the satellite peak related to surface plasmon at 472.2 eV are characteristic of Ti 4+ cation in TiO2 [START_REF] Iatsunskyi | Structural and XPS studies of PSi/TiO2 nanocomposites prepared by ALD and Agassisted chemical etching[END_REF]. No significant signal of partially reduced Ti n+ (n<4) to the lower energies from the principal Ti2p peaks [START_REF] Benito | Mixed Ti-O-Si oxide films formation by oxidation oftitanium-silicon interfaces[END_REF] have been ob-served in our spectra, which confirms stoichiometry of the produced TiO2 films. However, a weak Ti2p3/2 contribution suggests a minor (~2 at.%) appearance of suboxide Ti 3+ at 457.4 eV, which nature will be discussed below. The O1s spectra were de-convoluted into three peaks at 530.2 eV (46.3 at.%) 531.2 eV (3.4 at%) and 532 eV (8.9 at.%). The lower energy one originates from Ti-O-Ti bonds, while the less intense one at higher binding energy was attributed to the physisorbed species (Ophs), typically oxidised carbonaceous species, supported by observation of four C1s contributions in the vicinity of 285 eV with the total content of 15.5 at% (not shown). The O1s intermediate energy chemical environment could belong to O in cross-linking bonding Ti-O-Si [START_REF] Iatsunskyi | Structural and XPS studies of PSi/TiO2 nanocomposites prepared by ALD and Agassisted chemical etching[END_REF][START_REF] Zhong | Mechanism of the performance improvement of TiO2-x-based field-effect transistor using SiO2 as gate insulator[END_REF][START_REF] Hu | Tunable optical properties of co-sputtered Ti-SiO2 nanocomposite thin films[END_REF]. We notice that a weak peak associated with Si-O-Si bonds at about 533.0 eV could not be identified in the intense spectrum of O1s because of the screening of pure silica substrate by TiO2 coating. This O1s analysis might indicate (O)3Ti-O-Si(O)3 bonds formed between the coating and substrate. However, Ti-O-H bonds [START_REF] Iatsunskyi | Structural and XPS studies of PSi/TiO2 nanocomposites prepared by ALD and Agassisted chemical etching[END_REF] could also contribute at this energy that may compromise the assignment. The Si2p spectra were therefore recorded to verify this point.

The Si2p spectrum in Fig. 3 with the maximum at 102.3 eV and full-width at half maximum E1/2=1.8 eV (Fig. 3) can be due to Si 3+ suboxide [START_REF] Mehonic | Structural changes and conductance thresholds in metalfree intrinsic SiOx resistive random access memory[END_REF][START_REF] Brassard | High-k titanium silicate thin films grown by reactive magnetron sputtering for complementary metal-oxide-semiconductor applications[END_REF], which chemical shift of 1.05 eV per oxidation state of Si n+ has been earlier established [START_REF] Hollinger | Probing the transition layer at the SiO2-Si interface using core level photoemission[END_REF]. Its intensity was significantly weaker compared with Ti2p and O1s spectra, owing to the screening effect of the TiO2 film. Surprisingly, the signal of Si 4+ at 103.8 eV characteristic of silica was missing, which is expected to come from silica substrate. We notice in this respect that the peak at ~102-103 eV is a common feature of silicon exposed to plasmas and according to Cardinaud et al. [START_REF] Ch | Analyse XPS des surfaces de Si et SiO2 exposées aux plasmas de CHF3 et CHF3-C2F6. Polymérisation et gravure[END_REF] could only be indicative of the Si n+ oxidation state and belong to an ultrathin surface SiO2 layer. Although Si-O-Si moieties make surface non-reactive, they are expected to break upon thermal activation into fragments Si-O + Si, forming thus Si 3+ point defects that promote the film adherence (see e.g. [START_REF] Yates | Water interactions with silica surfaces: A big role for surface structure[END_REF][START_REF] Wendt | The interaction of water with silica thin films grown on Mo(112)[END_REF]). Copyright line will be provided by the publisher

Considering the weak Si 3+ and Ti 3+ signals assigned to the suboxide cations, one can assume participation of asymmetric (O)3Si-Ti(O)3 bonds between the titania film and silica substrate. We notice in this connection that annealing at 300 °C promotes the formation of Ti-Si bonds in the Ti/O/Si system with a thin oxidized layer [START_REF] Butz | Chemical bonding and reactions at Ti/Si and Ti/oxygen/Si interfaces[END_REF] and plasma environment could activate the formation of Ti-Si bonds in our experiments. Furthermore, another solution can be proposed involving formation of (O)3Si-O-Ti(O)3 bonds. In fact, examination of Ti coordination number changes in (TiO2)x(SiO2)1-x layer has been undertaken by Gaultois and Grosvenor [START_REF] Gaultois | XANES and XPS investigations of (TiO2)x(SiO2)1-x: the contribution of final-state relaxation to shifts in absorption and binding energies[END_REF], who have shown that Ti2p and Si2p XPS band energies significantly decrease with Ti content increase due to nearest-neighbour and next-nearest-neighbour effects. In particular, the Si2p energy has decreased by 0.74 eV as x increased from 0 to 0.33. Therefore, somewhat large x-value in the interface region could explain the Si2p band position measured in our study. Although we cannot prefer one solution, both (O)3Si-Ti(O)3 and (O)3Si-O-Ti(O)3 covalent bonds assure strong adherence of the prepared nanocomposite coatings.

The degradation of cefixime was examined as a model reaction to compare photocatalytic activities of the prepared composite TiO2 films with different rutile/anatase ratios. The corresponding process followed first order reaction kinetics as shown in Fig. 4. Our experiments showed almost negligible photolytic decomposition of cefixime by UVA light on the timescale of four hours. In contrast, a strong disappearance of pollutant was observed in presence of the photocatalytic coatings. This effect was amplified in D2, D3 and D4 composite coatings containing rutile phase. The photocatalytic rates obtained from the linear fit of the first-order kinetic curves are shown in Fig. 5 as a function of rutile phase content in the anatase-rutile coatings. These results clearly show the activity enhancement with the increasing content of rutile phase. Furthermore, the activity gain in D2 (16 mol% rutile) compared to pure anatase TiO2 coating (D1) was 300 %. The following evolution was not so significant, and the overall gain in the best D4 (30 mol% rutile) was about 400 % over D1. More support to the similar architecture of the prepared films and Degussa P25 TiO2 nanopowders can be found in TRMC measurements, which reflect the behaviour of photoinduced charges. An important feature of the reference photocatalyst is an intense and long-lasted response of the photoexcited electrons in the conduction band [START_REF] Colbeau-Justin | Structural influence on charge-carrier lifetimes in TiO2 powders studied by microwave absorption[END_REF]. The TRMC signals of D1 (no O2) and D2 (2 ml/min O2) films and Degussa P25 TiO2 are shown in Fig. 6. In comparison with D2, the amplitude of D1 coating is 2.7 times lower suggesting a smaller number of long-lived charges, which is in agreement with its lower photocatalytic activity (Fig. 5). Except for the very beginning related to the photoinduced charges thermalisation, the both curves follow the power-law decay I=At -k , which has been previously related to the relaxation of a series of localised states observed via transient absorption spectroscopy [START_REF] Sobennikov | Femtosecond relaxation dynamics of photoinduced nanocrystalits TiO2 in aqua suspensions[END_REF] and later attributed to the liberation of charges from shallow traps (Ei) of the Urbach tail (EU) to the conduction band (Eg) with power constant k expressed as k=kBT/Ei [START_REF] Cheng | Photocatalytic nanoparticulate ZrxTi1-xO2 coatings with controlled homogeneity of elemental composition[END_REF], where Eg, EU and Ei are band gap, Urbach and trap states energies, kB is Boltzmann constant and T is ambient temperature. According to our measurements, the decay constants of the both materials were very similar and correspond to Eg-Ei  70 meV, which is twice larger than EU at room temperature [START_REF] Tang | Urbach tail of anatase TiO2[END_REF]. The similar amplitudes and decay time of the photoinduced charge kinetics (t>50 ns) suggest a similar energetic structure of the shallow states in both composite anatase-rutile materials. In contrast, evident differences in the thermalisation process kinetics (t<50 ns) was observed between films (Fig. 6a) and nanopowders (Fig. 6b), which understanding requires further studies.

Interesting that the rutile content of the prepared coatings was in the range reported for Degussa P25 TiO2 nanopowders: 2010 %. Furthermore, a net gain in the photocatalytic efficiency was observed in presence of rutile, in contrast to the opposite trend reported by Ohtani et al. [START_REF] Ohtani | What is Degussa (Evonik) P25? Crystalline composition analysis, reconstruction from isolated pure particles and photocatalytic activity test[END_REF], who has found no synergetic effect of the phases in contact. Further studies concerning photoexcited charges and reaction kinetics by using different pollutant molecules are therefore required in order to confirm the conclusion.

Conclusion

The nanocomposite TiO2 coatings composed of 10-nm crystallinities with anatase/rutile ratio varying between 0 and 30 mol% were prepared on glass substrates via sprayplasma deposition method and controlled via oxygen gas flow in plasma. The coatings formed strong covalent bonds (O)3Ti-Si(O)3 and (O)3Si-O-Ti(O)3 with silica substrate, suggesting their strong mechanical adherence. The degradation of cefixime as a model reaction in aqueous solutions under UVA light illumination was examined to compare photocatalytic activities of coatings. We showed that the coatings activity increased with an increase of rutile content by a factor of 3 to 4 compared to pure anatase TiO2 nanocoatings.

The use of immobilised nanostructured materials in large-scale applications is preferable since permits avoiding (or strongly suppressing) environmental pollution by nanoparticles. The prepared nanocomposite coatings can therefore be an attractive alternative to Degussa P25 powderlike photocatalyst, which rutile content is in the range of that in our materials.

Figure 1 :

 1 Figure 1: X-ray diffraction patterns (Cu-K radiation source) of TiO2 coatings D1 (a), D2 (b) and D3 (c). Peaks positions of anatase (COD -ID 9015929) and rutile (COD -ID 9015662) TiO2 polymorphs [27] are indicated below.

Figure 2 :

 2 Figure 2: HRTEM image of D4 coating (left) and analysis of interplane distances of anatase and rutile lattices (right).

2 Figure 3 :

 23 Figure 3: Ti2p, O1s and Si2p XPS spectra of TiO2 coatings on silica substrates (example of D1coating is given).

  photolysis D1 D2 D3 D4

Figure 4 :

 4 Figure 4: Photocatalytic kinetics of cefixime decomposition in aqueous solutions under UVA (=365 nm) light illumination.

Figure 5 :

 5 Photocatalytic rate, min -1

Figure 6 :

 6 Figure 6: TRMC signals of D1 and D2 coatings (a) and P25 Degussa TiO2 (b). The power least-squares fits with depicted power constants k are shown by solid lines.

Table 1 :

 1 Oxygen gas flow conditions and crystalline composition of coatings.

	coat-	oxygen	anatase	rutile crys-	rutile con-
	ing	gas flow	crystalline	talline size	tent (1)
		(mL/min)	size (nm)	(nm)	(mol%)
	D1	0	9.3	-	0
	D2	2	9.6	8.1	16
	D3	10	8.7	8.5	27
	D4	25	8.5	8.1	30

(1)

Complementary to anatase phase

Acknowledgments

Acknowledgments are due to the French Embassy in Tehran for a PhD grant supporting this work.

Declaration

The authors declare no conflict of interest.


Ibrahim Almuslimani
email: ibrahim.almuslimani@univ-rennes1.fr.

A fully adaptive explicit stabilized integrator for advection-diffusion-reaction problems

Keywords: advection-diffusion-reaction equations, explicit stabilized methods, Runge-Kutta Chebyshev methods, RKC, SK-ROCK, ARKC AMS subject classification (2010): 65L04, 65L20, 65M12

A novel second order family of explicit stabilized Runge-Kutta-Chebyshev methods for advection-diffusion-reaction equations is introduced. The new methods outperform existing schemes for relatively high Peclet number due to their favorable stability properties and explicitly available coefficients. The construction of the new schemes is based on stabilization using second kind Chebyshev polynomials first used in the construction of the stochastic integrator SK-ROCK. An adaptive algorithm to implement the new scheme is proposed. This algorithm is able to automatically select the suitable step size, number of stages, and damping parameter at each integration step. Numerical experiments that illustrate the efficiency of the new algorithm are presented.

Introduction

In this paper we use the idea of stabilization by combining first and second kind Chebyshev polynomials introduced in [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF] to derive explicit stabilized methods for advection-diffusion problems with, possibly, costly non-stiff reaction terms, ∂ t u(x, t) = ∇ • (D∇u(x, t)) -∇ • (vu(x, t)) + r(u(x, t)), (x, t) ∈ Ω × [0, T ], with initial and boundary conditions, where Ω ∈ R d , D is the matrix of diffusion coefficients, and v is the velocity vector. The function r represents non-stiff, but possibly costly, reaction terms. Note that in general, D and v may also depend on u leading to nonlinear diffusion and advection terms. In the linear one dimensional setting, the equation reduces to

∂ t u(x, t) = d∂ 2
x u(x, t) -a∂ x u(x, t) + r(u(x, t)), (x, t)

∈ Ω × [0, T ] (1)
where d and a are positive reals, and Ω is a real interval. The Peclet number is defined by a/d and is allowed here to be quite large. When discretizing the partial differential equation (PDE) (1) in space using centered finite difference for example, with mesh size ∆x, We obtain a system of ordinary differential equations (ODEs) of the form ẏ(t) = F D (y(t)) + F A (y(t)),

y(0) = y 0 ∈ R d , t ∈ [0, T ], (2)
where F D represents the diffusion term with eigenvalues of its Jacobian grow as 1/∆x 2 on the negative real axis, and F A represents the advection term (and possibly non-stiff reaction terms) with eigenvalues of its Jacobian are of size 1/∆x and located close to the imaginary axis and symmetric with respect to the origin. This means that the eigenvalues of the Jacobian of the obtained system are approximately located in an ellipse with the length of its minor axis proportional to the square root of the length of the major axis. Explicit stabilized Runge-Kutta-Chebyshev methods were originally introduced in the context of purely diffusive or diffusion dominated advection-diffusion problems (very small Peclet number) as a compromise between costly implicit methods and restrictive usual explicit schemes [START_REF] Abdulle | Chebyshev methods based on orthogonal polynomials[END_REF][START_REF] Abdulle | Fourth order Chebyshev methods with recurrence relation[END_REF][START_REF] Abdulle | Explicit Stabilized Runge-Kutta Methods[END_REF][START_REF] Abdulle | Second order chebyshev methods based on orthogonal polynomials[END_REF][START_REF] Van Der Houwen | On the internal stability of explicit, m-stage Runge-Kutta methods for large m-values[END_REF]. Due to their versatility, they were extended to many other types of problems such as advection-diffusion-reaction equations [START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF][START_REF] Tang | Improved runge-kutta-chebyshev methods[END_REF][START_REF] Verwer | An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations[END_REF][START_REF] Verwer | RKC time-stepping for advectiondiffusion-reaction problems[END_REF][START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF], stochastic differential equations (SDEs) [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF][START_REF] Abdulle | Explicit stabilized multirate method for stiff stochastic differential equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF][START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF], and optimal control problems [START_REF] Almuslimani | Explicit stabilized integrators for stiff optimal control problems[END_REF]. Other types of stabilized methods were studied in [START_REF] Torrilhon | Essentially optimal explicit Runge-Kutta methods with application to hyperbolic-parabolic equations[END_REF][START_REF] Jeltsch | Flexible stability domains for explicit Runge-Kutta methods[END_REF][START_REF] Ketcheson | Optimal stability polynomials for numerical integration of initial value problems[END_REF].

Typically, the stability domain of an explicit stabilized method contains a long narrow strip around the negative real axis. In the context of advection-diffusion problems, the authors of the article [START_REF] Verwer | RKC time-stepping for advectiondiffusion-reaction problems[END_REF] propose the usage of the RKC method with very large damping parameter to make the strip wider, which means that eigenvalues with slightly larger imaginary parts coming from the advection terms can be put in. This comes at the cost of a serious shortening of the strip, which means that, for a fixed time step, less eigenvalues with negative real parts can be put in. Later, a partitioned Runge-Kutta-Chebyshev method (PRKC) of order 2 was designed in [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF] based on the RKC method [START_REF] Sommeijer | RKC: an explicit solver for parabolic PDEs[END_REF] for the integration of ODEs that have a moderately stiff term (diffusion) and non-stiff terms (advection or costly reaction terms). PRKC has a limited stability for the advection term, and it shares with the standard RKC the same stability domain length over the negative real axis. In [START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF], the authors propose a partitioned implicit-explicit orthogonal Runge-Kutta method (called PIROCK) for the time integration of advection-diffusion-reaction problems with possibly severely stiff reaction terms and stiff stochastic terms. The diffusion terms are solved by the explicit, nearly optimal, second order orthogonal Chebyshev method (ROCK2). Applied to advection-diffusion problems, the method has order 2 of accuracy and can handle the large Peclet number regime but it needs very large damping that reduces a lot its stability domain length over the negative real axis. In addition, PIROCK relies on the ROCK2 method, for which no explicit formulas are available to compute the coefficients for a given stage number. Recently, the authors of [START_REF] Tang | Improved runge-kutta-chebyshev methods[END_REF] developed an improved version of the RKC method (called IMPRKC) for advectiondiffusion-reaction equations. Their idea is based on introducing an appropriate combination of RKC polynomials which leads to a significant increase of the width of the stability domain along the imaginary axis with almost no loss of its length along the real axis. This comes at the cost of a few additional function evaluations. The main drawback of the IMPRKC scheme is evaluating F A at each stage of the method. This rich literature shows that the domain of stabilized schemes is very active, and that the construction of an adaptively efficient explicit stabilized integrator for such important class of problems is quite challenging.

In [START_REF] Almuslimani | Explicit Stabilized Methods for Stiff Stochastic Differential Equations and Stiff Optimal Control Problems[END_REF]Sect. 3.7.2], the author of the present paper profited from the idea of stabilization using second kind Chebyshev polynomials, first introduced in [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF] in the context of SDEs, to construct a first order explicit stabilized method for advection-diffusion-reaction equations with optimal stability domain. What made this quite intuitive is the similarity between mean square stability for SDEs and the stability of the test equation for ODEs of the form (2) (see Sect. 4.1). In this work, we construct a second order integrator based on RKC that outperforms the existing methods in the literature. We propose a fully adaptive algorithm to implement the new second order method.

This paper is organized as follows: in Section 2, we give a fast revision on explicit stabilized methods. In Section 3 we present an optimal first order explicit stabilized method for advectiondiffusion-reaction equations that we constructed in the thesis [START_REF] Almuslimani | Explicit Stabilized Methods for Stiff Stochastic Differential Equations and Stiff Optimal Control Problems[END_REF] inspired by previous work on SDEs [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF]. In Section 4 we derive and analyze our new second order adaptive scheme in terms of stability and convergence, and we propose a local error estimator for automatic step size selection. Section 5 is dedicated to present and analyze some numerical experiments that illustrate the efficiency of the new schemes. Finally, we conclude in Section 6.

Preliminaries on explicit stabilized methods

This section is devoted to present useful standard materials.

In order to study the stability of a Runge-Kutta integrator applied to an ODE, the following approach is widely used [START_REF] Abdulle | Explicit Stabilized Runge-Kutta Methods[END_REF][START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF]. Consider the test ODE

ẏ(t) = λy(t), y(0) = y 0 , (3)
where λ ∈ C with negative real part. If we apply a Runge-Kutta method with step size h on (3), we get the relation y n+1 = R(hλ) n y 0 , where the rational function R(z) is called the stability function. Hence, the stability domain of the method is defined as

S := {z ∈ C; |R(z)| ≤ 1}.
In the particular case of explicit methods, R(z) is a polynomial which means that the stability domain is necessarily bounded. For example, the stability domain of the explicit Euler method is just a disk of radius 1 which explains the time step restriction it faces for stiff ODEs, while that of the implicit Euler method is the complementary of a disk of radius 1. This shows the advantage of implicit methods in terms of stability, however, for large dimensional problems and especially the nonlinear and ill-conditioned ones, implicit methods become very costly and difficult to implement.

Here appears the need to a compromise between classical explicit methods and implicit integrators. This compromise is "explicit stabilized methods" (see the survey [START_REF] Abdulle | Explicit Stabilized Runge-Kutta Methods[END_REF]). The idea is to construct explicit Runge-Kutta integrators with extended stability domain that grows quadratically with the number of stages s of the method along the negative real axis, and then allows to use large time steps typically for problems arising from diffusion dominant advection-diffusion-reaction PDEs for which the eigenvalues are close to the negative real axis and are very large in modulus.

Before proceeding, let us recall some useful facts on Chebyshev polynomials. The first kind Chebyshev polynomials are defined by

T 0 (x) = 1, T 1 (x) = x, T j (x) = 2xT j-1 (x) -T j-2 (x), j ≥ 2. ( 4
)
The second kind Chebyshev polynomials are defined by

U 0 (x) = 1, U 1 (x) = 2x, U j (x) = 2xU j-1 (x) -U j-2 (x), j ≥ 2. (5)
Moreover, the two kinds polynomials satisfy the following

U j-1 (x) = T j (x) j . (6)
The stabilization procedure is based on the above relations. The fact that both kinds share the same recurrence relation will be very useful in our analysis. Indeed, this allows to simultaneously derive the recurrence formulas of the methods, otherwise, the cost would be doubled. It was shown that for any explicit, consistent (order 1) Runge-Kutta method, the maximum stability domain length over the negative real axis is 2s 2 , where s is the number of stages of the method. The polynomial that achieves this optimal length is the shifted Chebyshev polynomial T s (1 + z/s 2 ). See, for example, [17, Chap. V, Th. 1.1]. For robustness reasons, a damping of this polynomial is introduced and the resulting scheme is recalled in the following subsection.

Optimal first order Chebyshev methods

Consider the ODE ẏ = f (y),

y(0) = y 0 , t ∈ [0, T ]. ( 7
)
Given y 0 , in order to compute y 1 ≈ y(h) using the optimal first order Chebyshev method applied to [START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF] with step size h, the following recurrence is applied

K 0 = y 0 , K 1 = K 0 + µ 1 hf (K 0 ), K j = µ i hf (K j-1 ) + ν i K j-1 + (1 -ν i )K j-2 , j = 2, . . . , s (8)
y 1 = K s ,
where ω 0 := 1 + η s 2 , ω 1 := Ts(ω0) T s (ω0) , and

µ 1 := ω 1 ω 0 , µ j := 2ω 1 T j-1 (ω 0 ) T j (ω 0 ) , ν j := 2ω 0 T j-1 (ω 0 ) T j (ω 0 ) , j = 2, . . . , s. (9)
The parameter η is called the damping parameter and it is necessary to avoid singularities in the stability domain which ensures the robustness of the method (See Figure 1). Typically for this method, η is fixed to 0.05. It can be easily verified, using the recurrence (4) and proceeding by induction, that applied to the test problem (3), the above method produces after one step

y 1 = K s = R s (hλ)y 0 with R s (z) = T s (ω 0 + ω 1 z) T s (ω 0 ) ,
for which the stability domain contains a narrow strip around the interval [-C η s 2 , 0] with

C η = 1 + ω 0 s 2 ω 1 2 -4/3η
is very close to 2 (the optimal value for order 1). For η = 0 the stability function is again T s (1 + z/s 2 ). The method has low memory requirements (only two stages have to be stored) and reasonable propagation of round-off errors even for large values of s needed in practice [START_REF] Van Der Houwen | On the internal stability of explicit, m-stage Runge-Kutta methods for large m-values[END_REF][START_REF] Verwer | Convergence properties of the runge-kuttachebyshev method[END_REF].

The fact that the length of the stability domain on the negative real axis enjoys a quadratic growth with respect to the number of stages s is crucial to the success of explicit stabilized Runge-Kutta methods.

Second order RKC methods

To design a second order method, we need the stability polynomial to satisfy *

R(z) = 1 + z + z 2 2 + O(z 3 ).
A correction to the first order shifted Chebyshev polynomials was introduced to construct a stabilized scheme of order 2 [START_REF] Bakker | Analytical aspects of a minimax problem[END_REF][START_REF] Van Der Houwen | On the internal stage runge-kutta methods for largemvalues[END_REF][START_REF] Verwer | Convergence properties of the runge-kuttachebyshev method[END_REF]. The obtained second order polynomial is the following

R s (z) = a s + b s T s (ω 0 + ω 2 z),
where,

a s = 1 -b s T s (ω 0 ), b s = T s (ω 0 ) (T s (ω 0 ) 2 ) , ω 0 = 1 + η s 2 , ω 2 = T s (ω 0 ) T s (ω 0 ) , η = 0.15. ( 10
)
For each s, |R s (z)| remains bounded by a s + b s = 1 -η/3 + O(η 2 ) for z in the stability interval (except for a small interval near the origin). The stability interval along the negative real axis is [-1+ω0 ω2 , 0] which is approximately [-0.65s 2 , 0], and covers about 80% of the optimal stability interval for second order stability polynomials, and the formula now for calculating s for a given time step h is

s := hλ max + 1.5 0.65 + 0.5 ,
where the brackets mean rounding to the nearest integer, and λ max is the spectral radius of the Jacobian of f that can be calculated at each step using power method for example. Using the recurrence relation of the Chebyshev polynomials, the RKC method as introduced in [START_REF] Verwer | Convergence properties of the runge-kuttachebyshev method[END_REF] is defined by

K 0 = y 0 , K 1 = K 0 + hb 1 ω 2 f (K 0 ), K j = µ j h(f (K j-1 ) -a j-1 f (K 0 )) + ν j K j-1 + κ j K j-2 + (1 -ν j -κ j )K 0 , y 1 = K s , (11)
where

µ j = 2b j ω 2 b j-1 , ν i = 2b j ω 0 b j-1 , κ j = - b j b j-2 , b j = T j (ω 0 ) T j (ω 0 ) 2 , a j = 1 -b j T j (ω 0 ), (12)
for j = 2, . . . , s. The parameters b 0 and b 1 are free (R 0 (z) is constant and only order 1 is possible for R 1 (z)) and the values b 0 = b 1 = b 2 are suggested in [START_REF] Sommeijer | A performance evaluation of a class of Runge-Kutta-Chebyshev methods for solving semidiscrete parabolic differential equations[END_REF].

Remark 2.1. For simplicity of the presentation, we consider the case of autonomous problems (with f independent of time) but we highlight that our approach also applies straightforwardly to non-autonomous problems ẏ(t) = f (t, y(t)). Indeed, a standard approach is to consider the augmented system with z(t) = t, i.e. dz dt = 1, z(0) = 0 and define ỹ(t) = (y(t), z(t)) T , see e.g. [14, Chap. III] for details. 3 Optimal first order scheme

Note that we have discussed the content of this section in the thesis [START_REF] Almuslimani | Explicit Stabilized Methods for Stiff Stochastic Differential Equations and Stiff Optimal Control Problems[END_REF], but we recall it here since it gives insight about the construction of the second order adaptive integrator in the next section.

Consider the linear test problem

ẏ = λy + iµy, y(0) = y 0 , (13)
where λ ∈ R -, µ ∈ R, and i = √ -1. Applying a Runge-Kutta method to the above equation, one gets an induction of the form y n+1 = R(p, q)y n , with p = hλ and q = hµ. We define the stability domain of a Runge-Kutta method applied to (13) by S = {(p, q) ∈ R 2 ; |R(p, q)| ≤ 1}.

Equation ( 13) can be seen as the test equation for linear SDEs with the iµ replacing the noise. Hence, inspired by SK-ROCK [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF], we consider the following stability polynomial

R(p, q) = A(p) + B(p)iq := T s (ω 0 + ω 1 p) T s (ω 0 ) + U s-1 (ω 0 + ω 1 p) U s-1 (ω 0 ) (1 + ω 1 2 p)iq, (14)
where T s and U s are the first and the second kind Chebyshev polynomials of degree s (the number of stages), and the coefficients ω 0 and ω 1 are the same as for the Chebyshev method [START_REF] Abdulle | Second order chebyshev methods based on orthogonal polynomials[END_REF]. The stability condition |R(p, q)| ≤ 1 is equivalent to A(p) 2 + B(p) 2 q 2 ≤ 1 which is exactly the mean square stability condition described in [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF]Sect. 2]. By [5, Theorem 3.2] and [5, Remark 3.6], for η > 0 and s ∈ N, |R(p, q)| ≤ 1 for all p ∈ [-2ω -1 1 , 0] and all q such that |q| ≤ √ -2p (See Figure 2). Remark 3.1. For SDEs, the condition |q| ≤ √ -2p guarantees the stability of the solution of the continuous problem, which means that numerical stability under this condition is sufficient. However, for our test equation [START_REF] Bakker | Analytical aspects of a minimax problem[END_REF], the exact solution is stable when λ ≤ 0 and µ ∈ R, hence even when |q| ≤ √ -2p is violated. Although stability under this condition is already a significant improvement comparing to standard explicit methods, we can use larger damping parameter η to increase the width of the stability region along the imaginary direction in the case of large advection.

The new order one method for space discretized advection-diffusion-reaction equations (2), is defined in the same way as the SK-ROCK method [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF], by replacing the noise term by the advection-reactions terms,

K 0 = y 0 K 1 = y 0 + µ 1 hF D (y 0 + ν 1 hF A (y 0 )) + κ 1 hF A (y 0 ) K j = µ j hF D (K j-1 ) + ν j K j-1 + κ j K j-2 , j = 2, . . . , s, y 1 = K s , (15)
where ν 1 = sω 1 /2, κ 1 = sω 1 /ω 0 and the rest of the coefficients are identical to those defined in [START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF]. Assuming enough regularity on F D and F A , the convergence proof is straightforward and based on [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF]Lemma 4.2]. The above scheme is also optimal in the sense that its stability region achieves the maximum possible length over the negative real axis for an explicit consistent scheme.

Note that the method requires only 1 evaluation of the advection-reaction terms F A per time step, while other existing methods for advection diffusion problems such as PIROCK [START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF] and PRKC [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF] require 3 and 4 evaluations of F A respectively per time step. It also outperform the mentioned methods by far in terms of stability. However, it has only order 1 of accuracy, which motivates the construction of a second order integrator that competes with existing schemes in terms of stability and convergence, taking advantage of the above analysis.

Note that we can adaptively increase the stability for advection by increasing η in the price of loosing some length on the negative real axis. This will be explained in details in the next section for the second order scheme, for which this adaptive increase of damping is a key feature.

New second order scheme

In this section, we introduce the new adaptive second order integrator based on the material presented in the first 3 sections and inspired by the SK-ROCK method introduced in [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF].

ARKC integrator

The ARKC integrator (for adaptive RKC) applied to (2) is defined as follows for s ≥ 2:

G = hF A (y 0 + h 2 F A (y 0 + ω 2 2 hF D (y 0 )) + h 2 F D (y 0 )) + hF D (y 0 + ω 2 -1 2 hF A (y 0 )) -hF D (y 0 ), K -1 = y 0 , K 0 = K -1 + ω 2 2 G, K 1 = K 0 + b 1 ω 2 hF D (y 0 )
+ αG, and for j = 2, . . . , s,

K j = µ j h(F D (K j-1 ) -F D (K 0 ) + (1 -a j-1 )F D (y 0 )) + ν j K j-1 + κ j K j-2 + (1 -ν j -κ j )K 0 , y 1 = K s , (16)
where α = 1 -ω2 2 b 1 sω 2 , and the other coefficients are identical to those defined in [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF] and [START_REF] Almuslimani | Explicit stabilized integrators for stiff optimal control problems[END_REF]. Note that for purely diffusive equations (F A ≡ 0), the method reduces to standard RKC [START_REF] Almuslimani | Explicit Stabilized Methods for Stiff Stochastic Differential Equations and Stiff Optimal Control Problems[END_REF].

Complexity The method requires s + 2 evaluations of F D and 3 evaluations of F A per time step. The standard RKC and the PRKC methods both require s evaluations of F D per time step, but RKC requires s evaluations of F A while PRKC requires only 4. The PIROCK merhod needs s + 2 + l evaluations of F D (l = 1 or 2) and 3 evaluations of F A per time step. Finally, the IMPRKC scheme requires s + ŝ evaluations of each F D and F A , where ŝ becomes quite large for very stiff problems. The advantage of our new scheme comes from having larger stability region that changes adaptively, which allows to use bigger time steps.

Construction In order to construct a second order scheme, we need our stability function to satisfy the following equality

R 2 (p, q) = 1 + p + iq + 1 2 (p + iq) 2 + O(((p + iq) 3 ) = 1 + p + iq + p 2 2 + ipq - q 2 2 + O((p + iq) 3 ). (17)
A natural approach would be to modify the last stage of the first order scheme to reach second order, but such naive modification will cause severe instability. Our idea is to start the stabilization using second kind Chebyshev polynomials from the beginning of the integration process. Hence, inspired by the previous section, we consider the following polynomial

R 2 (p, q) = A 2 (p) + B 2 (p)(iq - q 2
2 )

:= a s + b s T s (ω 0 + ω 2 p) + ω 2 2 + 1 - ω 2 2 U s-1 (ω 0 + ω 2 p) U s-1 (ω 0 ) (1 + ω 2 2 p)(iq - q 2 2 ), (18)
where all the coefficients are defined in [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF], from which we can derive the new adaptive second order ARKC method ( 16) for advection-diffusion-reaction problems, using the relations (4), ( 5), and [START_REF] Abdulle | Explicit stabilized multirate method for stiff stochastic differential equations[END_REF]. The first part of our stability function [START_REF] Jeltsch | Flexible stability domains for explicit Runge-Kutta methods[END_REF] corresponds to the stability polynomial of the standard second order RKC scheme, while the second part involves second kind Chebyshev polynomials (as in [START_REF] Hairer | Geometric numerical integration[END_REF]) in order to stabilize the advection part, and has the correct order. The construction of the method from the polynomial ( 18) is done by induction using the relations (4)- [START_REF] Abdulle | Explicit stabilized multirate method for stiff stochastic differential equations[END_REF]. See also Lemma 4.2 below. Unlike the case of standard ODEs where the term z 2 /2 is enough to have order two for nonlinear problems, some additional coupling conditions need to be satisfied here due to the partitioned nature of the scheme (see [START_REF] Hairer | Geometric numerical integration[END_REF]Sect. III.2]). These conditions result from the fact that the term ipq in [START_REF] Hundsdorfer | Numerical solution of time-dependent advection-diffusionreaction equations[END_REF] is in fact the sum of two terms: 1 2 ipq and 1 2 iqp, which are not necessarily equal for multidimensional and/or nonlinear problems. The quantity G is constructed carefully to get order 2 of convergence for general nonlinear problems. For linear equations, G = (1 + ω2

2 )(iq -q 2 2 ). Analogously to some existing second order methods [START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF][START_REF] Verwer | RKC time-stepping for advectiondiffusion-reaction problems[END_REF], our damping parameter is not fixed to a small value, but it is an increasing function of the stage number s. However, for relatively large Peclet numbers, no huge damping is needed, and the method still perform very well as will be shown in Section 5.

Stability analysis

Remark 4.1. The stability analysis of Runge-Kutta type methods is usually made on linear test problems, but we emphasize that such test equations only give insight on the stability of the stabilized method under study (here ARKC). This is because practical problems are usually nonlinear, and even for linear problems, the involved operators cannot in general be diagonalized simultaneously.

Lemma 4.2. The scheme (16) applied to the linear test problem (13) with step size h, produces the following recurrence

y n+1 = R s (hλ, hµ)y n ,
where the stability polynomial R s (p, q) is defined in [START_REF] Jeltsch | Flexible stability domains for explicit Runge-Kutta methods[END_REF]. Proof. By induction on j, it can be shown that, for every j ≥ 1, the internal stages satisfy

K j = R j (p, q)y n ,
where, R j (p, q) := a j + b j T j (ω 0 + ω 2 p)

+ ω 2 2 + 1 - ω 2 2 U j-1 (ω 0 + ω 2 p) b j sω 2 (1 + ω 2 2 p)(iq - q 2
2 ).

For j = s, we have that b s sω 2 = 1/U s-1 (ω 0 ) because s -1 T s (ω 0 ) = U s-1 (ω 0 ) and the proof is done.

In contrast to the first order scheme [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF] and the stochastic integrator SK-ROCK [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF], the damping parameter for the ARKC method is not fixed. Indeed, it is an increasing function of the Peclet number and the number of stages s. We provide numerical stability analysis that illustrates the nice features of the new scheme (see Figures 3 and4).

In what follows, we will respectively denote by δ s and α s the half width and the half height of the largest ellipse that can be put in the stability region of the corresponding method.

The method ( 16) is designed to handle quite large Peclet numbers, this means that we need the width of the stability region in the imaginary direction to be as large as possible. While the length over the negative real axis grows quadratically with s, we cannot achieve more than linear growth on the imaginary direction (this fact is given as an exercise in [16, Chap. IV]). In [START_REF] Verwer | RKC time-stepping for advectiondiffusion-reaction problems[END_REF], the RKC method is proposed with large damping value which reduces the length of the stability domain over the negative real axis to δ s = 0.34s 2 , and leads to a growth of only O( √ s) for the ellipse half-height α s . The PRKC method proposed in [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF] uses the standard small damping η = 0.15 that keeps δ s ≈ 0.65s 2 , but α s is fixed to 1.7 no matter how large is the number of stages s, which does not add much to the standard RKC method since this might be useful only in the small Peclet number regime. However, the main feature of the method is to reduce the number of evaluations of possible non stiff terms such as non-stiff advection or reaction terms. The PIROCK integrator [START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF], is notably better than the two mentioned methods RKC and PRKC. Two kinds of damping are proposed for PIROCK, the first lead to a nearly optimal length over the negative real axis δ s ≈ 0.81s 2 but α s is limited to 0.07696s + 1.878. For the second damping, α s = 0.5321s + 0.4996 but δ s reduces to 0.43s 2 . It shares as well the feature of reducing the evaluation of non diffusive terms. Nevertheless, PIROCK still have two drawbacks: the first one is the non-availability of explicit closed form formulas to compute its coefficients for a given number of stages s, the other is the limited choice of damping.

The stability region of the ARKC method changes adaptively with the spectrum of the Jacobian of the vector fields at each time step. This change is due to the damping parameter η which is not necessarily constant. In fact, the ellipse half width δ s still grows quadratically with s, while the ellipse half height α s grows linearly with s, in contrast to RKC and PRKC methods. For example, for the standard value of η = 2/13 used for RKC in [START_REF] Sommeijer | RKC: an explicit solver for parabolic PDEs[END_REF] and PRKC in [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF], δ s is still equal to 0.65s 2 and α s 0.17s (see Figure 3a), whereas for PRKC α s is fixed to 1.7 . In addition, increasing the value of the damping parameter η adds more space in the imaginary direction which is very favorable for the advection dominated problems. Figure 3c illustrates the stability region in the p-q plane for s = 20 and η = 3. We can see that δ s = 0.5s 2 and α s = 0.5s, while the best that PIROCK [START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF] could achieve for almost the same α s is δ s = 0.43s 2 which makes difference for large values of s. For η = 10 we have δ s = 0.35s 2 and α s = 0.9s as shown in Figure 3d, compare that with the case of standard RKC with infinite damping, considered in [START_REF] Verwer | RKC time-stepping for advectiondiffusion-reaction problems[END_REF], where δ s = 0.34s 2 while α s = O( √ s). The last method to compare with is IMPRKC introduced in [START_REF] Tang | Improved runge-kutta-chebyshev methods[END_REF] where the authors provide a plot for s = 50 stages and ŝ = 3 additional stages. The length over the negative real axis stays almost the same as RKC at 0.65s 2 , while the width along the imaginary axis is around 0.8s. As said before, the evaluation of the advection-reaction term at every step of IMPRKC is a disadvantage that limits its performance. For instance, for small Peclet numbers it needs time steps as much as the standard RKC scheme which makes it a bit more expensive because of the additional stages needed. In Section 5, our scheme ARKC is shown to perform better in all regimes.

Choice of damping

In fact, α s is a function of s and η, this introduces additional difficulty in the estimation of the value of α s . Therefore, we will introduce many choices of the range of Peclet number to simplify the implementation. Let ρ D and ρ A be the spectral radii of the Jacobians of F D and F A respectively. For ODEs coming from the discretization of the linear PDE (1), we have ρ D = 4d/∆x 2 and ρ A = a/∆x, then

ρ A √ ρ D = a 2 √ d .
We define the number

P e := a √ d , hence, ρ A / √ ρ D = 1 2 P e .
Similarly, we can show that q = 1 2 P e √ -hp. We will use this new number P e to adapt our choice of damping according to the parameters of the problem. Notice that this is just another measure of the advection dominance that is proportional to the Peclet number and more convenient to use in our algorithm.

For P e ≤ 1/10, i.e, ρ A / √ ρ D ≤ 1/20, we fix η = 0.15 up to s = 200, and for s between 200 and 500 we set η = 0.6 (we do not allow s to be more than 500). For 1/10 < P e ≤ 1/2 which means that 1/20 < ρ A / √ ρ D ≤ 1/4, we consider the following choices for s and η:

2 ≤ s ≤ 30 31 ≤ s ≤ 60 61 ≤ s ≤ 110 111 ≤ s ≤ 160 η = 0.2 η = 0.45 η = 1 η = 1.5 161 ≤ s ≤ 260 261 ≤ s ≤ 360 361 ≤ s ≤ 500 η = 2.4 η = 3 η = 4 Table 1: 1/10 < P e ≤ 1/2, 1/20 < ρ A / √ ρ D ≤ 1/4.
Tables for other values of P e are given in appendix A. Infinitely many choices could be done, but we will present and use only a few choices since they are enough for the method to perform very well. The methodology to compute the above values is easy, it is enough to plot R s (p, q) for the corresponding choice of q and to vary η in a way that we stay stable for the given value of s (see Figure 4).

Convergence analysis

In this section we will prove the second order of convergence of the scheme ( 16) when applied to ODEs of the form (2) arising from the discretization of advection-diffusion-reaction problems.

Theorem 4.3. Let T > 0 and consider the system of ODEs (2) on the time interval [0, T ], where F D and F A are of class C 2 and are Lipschitz continuous. Suppose in addition that the first and second derivatives of F D and F A are bounded. Let N ∈ N, h = T /N , and t n = nh, n = 1, . . . , N , and consider the method (16) applied to (2) with step size h, such that the number of stages s and the damping parameter η are chosen appropriately to guarantee stability. Then, we have for all n = 1 . . . , N ,

y(t n ) -y n ≤ Ch 2 , ( 19
)
where C is independent of h and n. In other words, the method converges with order 2.

Proof. Let us prove first that the local error (the error after one step) satisfies

y(h) -y 1 = O(h 3 ). (20)
• The proof for γ j 3 is very similar to that of γ j 2 . Indeed, we can prove for each j = 0, . . . , s the equality γ j 3 = 2ω 2 2 b j Q j (ω 0 ), where

Q 0 (x) = 0, Q 1 (x) = 0, Q j (x) = T j-1 + 2xQ j-1 (x) -Q j-2 (x), j ≥ 2.
Using the relation (4), we observe that for all j ≥ 0,

Q j (x) = T j (x) 4 , which implies that γ s 2 = 2ω 2 2 b s T s (ω 0 )/4 = 1.
Thus, ( 20) is proved, and using regularity assumptions made on the vector fields, [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF]Theorem 3.6] implies the global convergence estimate [START_REF] Ketcheson | Optimal stability polynomials for numerical integration of initial value problems[END_REF].

Variable step size control and the fully adaptive algorithm

We introduce the following local error estimator that allows us to adaptively select the time step size in order to reach a given accuracy,

Est n+1 = C(12(y n -y n+1 ) + 6h(F D (y n ) + F A (y n ) + F D (y n+1 ) + F A (y n+1 ))),
where

C = 1/6 -c 2 + (1/2 -c 1 )ζ -1/6ζ, with ζ = 0 if F A ≡ 0 and ζ = 1 otherwise, and
c 1 = ω 2 2 1 - ω 2 2 1 + ω 2 U s-1 (ω 0 ) U s-1 (ω 0 ) , c 2 = sb s U s-1 (ω 0 ) ω 3 2 6 .
To get an intuition about the above coefficients, compare the third order term in the exact polynomial which is (p + iq) 3 /6 = p 3 /6 -iq 3 /6 + iqp 2 /2 -pq 2 /2 and the third order term of stability polynomial [START_REF] Hairer | Geometric numerical integration[END_REF] that is equal to c 2 p 3 + c 1 iqp 2 -pq 2 /2. Note that the above estimator is inspired by the one considered for RKC in the paper [START_REF] Sommeijer | RKC: an explicit solver for parabolic PDEs[END_REF], and they coincide for ζ = 0. In contrast to the estimators introduced for PRKC [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF] and PIROCK [START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF], we do not consider two separate estimators for F D and F A , since our method is defined in a different way. Indeed, for improved stabilization, the advection-reaction terms are computed at the beginning and not separately at the end, which makes the method more similar to RKC. We adopt the standard step size selection strategy proposed in [START_REF] Sommeijer | RKC: an explicit solver for parabolic PDEs[END_REF] for RKC (see also [15, page 167]). Now, we are ready to present our fully adaptive algorithm, Algorithm 4.4 (y 0 → y 1 ). Given a time step size h and an initial value y 0 :

• Calculate ρ A and ρ D at the current value of the solution.

• Calculate ρ A / √ ρ D and choose the corresponding table among Tables 1, 3,

• Search the minimum s (and the corresponding η) in the chosen table such that 1+ω0 ω1 > hρ D . • Generate the coefficients (10) and (12) and apply the recurrence (16) to calculate y 1 .

• Update the step size according to the automatic step size selection procedure and repeat until reaching the final time.

Remark 4.5. To ensure stability, it is enough to choose s such that 1+ω0 ω1 > hρ D , because the relation between η and the corresponding range for s in each table is built to ensure that, once -hρ D lies inside the stability domain, the whole ellipse containing the eigenvalues for the given spectral radii fits inside. Remark 4.6. In the case where the eigenvalues have nonzero imaginary part and very small real part, this is close to purely advective regime, which is out of the scope of the paper. However, this case can be treated in two different ways, either by using the adaptive algorithm without any modification and then the error estimator will choose a very small time step, or by modifying the algorithm to integrate such systems (extremely large Peclet number) with RK3 or RK4 (for which the stability domain includes a part of the imaginary axis).

For the calculation of the spectral radii ρ D and ρ A , we use the Matlab function "eig". In other programming languages, one can use nonlinear power method as in [START_REF] Abdulle | ROCK2 and ROCK4: software for stiff differential equations (discretized parabolic problems)[END_REF] for example. However, the cost of such methods is a different issue and is beyond the scope of the present paper. The code of the ARKC integrator as well as the drivers that reproduce the numerical experiments will be made publicly available on the page: https://sites.google.com/view/ ibrahim-almuslimani.

Numerical experiments

We will compare ARKC with IMPRKC, PIROCK, and PRKC as they were shown to outperform the other stabilized methods for advection-diffusion-reaction problems where the reaction term is not stiff.

Linear 1D advection-diffusion problem

An excellent example to compare the performance of ARKC with IMPRKC, PIROCK, and PRKC is the following 1-dimensional advection-diffusion equation with periodic boundary conditions

∂ t u + a∂ x u = ∂ 2
x u, u(x, 0) = sin(2πx), u(0, t) = u(1, t), [START_REF] Tang | Improved runge-kutta-chebyshev methods[END_REF] where a is a positive constant (P e = a), x ∈ [0, 1], and t ≥ 0. We discretize the space interval [0, 1] to a uniform grid {x k } N k=0 with x k = kh, and h = 1/N . We use second order central differences for the advection and the diffusion terms. We denote by u k (t) the approximation of u(x k , t), and the periodic boundary conditions propose that u 0 (t) = u N (t). The eigenvalues of the obtained matrix are

λ k = 2 h 2 (cos(2kπh) -1) - ia h sin(2kπh), k = 1, . . . , N,
and are located in an ellipse in the left half-plane C -, which makes the problem typical for the comparison of the three schemes.

For the numerical experiments, we take N = 150, t ∈ [0, 1/2], and we fix Atol = Rtol = tol. The initial step is fixed to 10 -3 . The number of rejected steps is always very small and thus neglected. We can clearly see in Table 2 that PRKC can compete with our scheme ARKC only in the very small Peclet number regime, while For moderate and large Peclet number, ARKC is notably better. On the other hand, the cost of ARKC is close to that of PIROCK for moderate P e with small advantage for the latter when using large tolerance 10 -2 . However, ARKC becomes cheaper and more accurate at the same time for large values of P e (in our experiment for P e = 10 and 12). For small tolerance 10 -5 ARKC outperforms PIROCK in all Peclet number regimes. This is expected because the damping (and so the vertical with of the stability region) is fixed for PIROCK (when F A = 0, PIROCK has δ s = 0.43s 2 and α s 0.53s + 0.5). The adaptivity of ARKC allows to continuously increase the vertical width of its stability region which lets it to be more flexible with respect to the change in P e . One should not forget that the explicit availability [START_REF] Tang | Improved runge-kutta-chebyshev methods[END_REF]. The numbers on the left correspond to tol = 10 -2 , while those on the right correspond to tol = 10 -5 .

of ARKC coefficients helps a lot in making the scheme adaptive with respect to the change in P e , the feature that is missing in PIROCK.

In the small Peclet number regime, IMPRKC performs very similar to standard RKC method, which is expected, and thus it is outperformed by ARKC and the other two schemes. For large Peclet numbers, IMPRKC outperforms PRKC and its cost is close to PIROCK, while ARKC performs much better due to the very small number of evaluations of the advection term compared to that of IMPRKC.

The smaller number of steps needed in ARKC compared to the three other methods, leads to a significantly lower number of F A evaluations no matter how big is the number of stages.

Burgers equation with a nonlinear reaction term

As an example of a PDE with variable Peclet number, we consider the following Burgers equation with a nonlinear reaction term and periodic boundary conditions

∂ t u + 10u∂ x u = ∂ 2 x u + sin(u 2 ), (x, t) ∈ [0, 1] × [0, 1/2], u(x, 0) = 1 + sin(2πx), u(0, t) = u(1, t) ∀ t ∈ [0, 1/2]. (23)
We discretize the above equation in space using second order central differences with ∆x = 10 -2 , into N + 1 grid points with N = 100. We calculate a reference solution using the Radau IIA method of order 5 [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF]. In Figure 5, we plot the Peclet number of equation [START_REF] Torrilhon | Essentially optimal explicit Runge-Kutta methods with application to hyperbolic-parabolic equations[END_REF] as a function of time. we see that it is variable and of quite large magnitude. In Figure 7, we compare the number of functions evaluations needed to obtain a given accuracy for the solution of the Burgers equation (23) using 3 different numerical methods: ARKC, IMPRKC, and PIROCK. The results are obtained for tol = 10 -r , r = 1, . . . , 6. The advantage of our scheme ARKC ( 16) is very clear. For IMPRKC, the number of F D evaluations is reasonable, while that of F A evaluations is very large compared to ARKC and PIROCK. For PIROCK, the fixed damping increases the number of steps and of functions evaluations with respect to ARKC. The flexibility of ARKC gives it remarkable advantage over the other schemes. We can also see that even when the number of F D evaluations for ARKC is close to that of the other two schemes, the number of evaluations of F A containing the nonlinear advection and reaction terms stays much lower, that is because of the low number of time steps needed.

Conclusion

In this paper, we have constructed a fully adaptive second order explicit stabilized Runge-Kutta-Chebyshev time integrator for advection-diffusion-reaction PDEs, called ARKC. The new scheme is implemented using an algorithm that is able to adaptively choose the step size, the number of stages and the damping parameter of the method according to the Peclet number. The new scheme is shown to outperform existing methods in the literature for the same type of problems. This high performance is a result of the full adaptivity of Algorithm 4.4, in particular, the adaptive damping that allows significant control of the form of the stability region as a function of the Peclet number.

A Damping and number of stages for some choices of Peclet number One can get more tables and increase the adaptivity of the algorithm with respect to damping. This will for sure increase the performance of the method. However, the method performs already very well with the tables we provided.

s

Figure 1 :

 1 Figure 1: Stability domains of the Chebyshev method (8) for s = 7 and different damping values η = 0, 0.05, 3.98.

Figure 2 :

 2 Figure 2: Stability domain of the new optimal first order method (15) in the p -q plane for s = 10 and η = 0.05. The dashed lines correspond to ± √ -2p.

  (a) η = 0.15, δs = 0.65s 2 , αs = 0.17s. (b) η = 1.5, δs = 0.56s 2 , αs = 0.35s. (c) η = 3, δs = 0.5s 2 , αs = 0.5s. (d) η = 10, δs = 0.35s 2 , αs = 0.9s.

Figure 3 :

 3 Figure 3: Illustration of the stability regions of the ARKC method (16) in the p -q plane for s = 20 and different values of the damping parameter η. For each stability region, we plot the largest possible ellipse that fits inside.

5 (b) η = 1 , 5 (c) η = 5 ,

 5155 η = 0.15, δs = 0.65s 2 . δs = 0.59s 2 . δs = 0.4345s 2 .

Figure 4 :

 4 Figure 4: Stability polynomial (18) for s = 200 (left) and s = 30 (middle and right), and different values of η and q.

Figure 5 :

 5 Figure 5: Change in Peclet number with respect to time in problem (23).

Figure 6 :

 6 Figure 6: Solution of equation (23) obtained using ARKC method (16) at different time moments.

Figure 6

 6 Figure6shows the solution u(x, t) of equation (23) obtained using ARKC method (16) at different time moments.

Figure 7 :

 7 Figure 7: Comparison between ARKC, IMPRKC, and PIROCK in terms of cost for problem (23).

Table 2 :

 2 Comparison of ARKC, PRKC, and PIROCK for the linear advection-diffusion problem

Table 3 :

 3 1/2 < P e ≤ 1, 1/4 < ρ A / √ ρ D ≤ 1/2.

		< 11	< 21	< 31	< 41	< 51	< 61	< 71	< 81	< 91 < 101
	η	0.15	0.6	1	1.4	1.7	2.1	2.4	2.7	3	3.3
	s < 121 < 141 < 161 < 181 < 201 < 251 < 301 < 401 < 501
	η	3.7	4.1	4.5	4.9	5.3	6	6.6	7.7	8.8
		s < 11 < 21	< 31	< 41	< 51	< 61	< 71	< 81
		η	0.7	1.5	2.3	2.9	3.5	4	4.5	4.9
		s < 91 < 101 < 141 < 181 < 251 < 301 < 401 < 501
		η	5.2	5.5	6.7	7.7	8.8	9.8	11	12

Table 4 :

 4 1 < P e ≤ 3/2, 1/2 < ρ A / √ ρ D ≤ 3/4.

	s < 11 < 21 < 31 51 < 71 < 111 < 151 < 311 501
	η	1	2.5	3.5	4.8	6	7.8	9	12.5	15

Table 5 :

 5 3/2 < P e ≤ 2, 3/4 < ρ A / √ ρ D ≤ 1.

	s < 11 < 21 < 31 51 < 71 < 111 < 151 < 311 501
	η	2	3.8	5	6.8	8	10.4	12		16	19
			Table 6: 2 < P e ≤ 2 √	2, 1 < ρ A /	√ ρ D ≤	√	2.
		s < 11 < 31 < 71 < 151 < 311 < 501
		η	4	9	13.5		18		23		27

Table 7 :

 7 P e > 2 √ 2, ρ A / √ ρ D > √ 2.

Acknowledgements The author is grateful to Gilles Vilmart for helpful discussions and comments.

Throughout the proof, F i (y 0 ) will be simply denoted by F i , and F i (y 0 )F j (y 0 ) will be denoted by F i F j with i, j ∈ {D, A}.

Using Taylor expansion, we can easily see that

Now, let us suppose that

The first two stages of the method satisfy:

. By performing a Taylor expansion of the stages K j defined in ( 16) and replacing K j-1 and K j-2 by the expansion defined in [START_REF] Sommeijer | A performance evaluation of a class of Runge-Kutta-Chebyshev methods for solving semidiscrete parabolic differential equations[END_REF], and finally identifying the coefficients, we get the following relations

In order to prove [START_REF] Sommeijer | RKC: an explicit solver for parabolic PDEs[END_REF], it is sufficient to show that γ s 1 = γ s 2 = 1 and γ s 3 = γ s 4 = γ s 5 = 1 2 . Obviously, γ j 6 = γ j 5 for all j ≥ 0. We will provide proofs for the first 3 coefficients, the other two can be done using similar arguments.

• For γ j 1 we have:

, thus γ j 1 are the internal stages of the RKC method [START_REF] Almuslimani | Explicit Stabilized Methods for Stiff Stochastic Differential Equations and Stiff Optimal Control Problems[END_REF] applied to the problem ẏ = 1, y(0) = 0 with step size h = 1. This means that for all 0 ≤ j ≤ s, γ j 1 = c j , where c j is the j th node of the RKC method. Therefore, γ s 1 = c s = 1. Moreover, according to [START_REF] Sommeijer | RKC: an explicit solver for parabolic PDEs[END_REF]Sect. 2], for all j = 2 . . . , s we have

• It can be proved that γ j 2 , j = 2, . . . , s are given by γ j 2 = ω2 2 + α bj b1 P j (ω 0 ), where P j (x) are polynomials that satisfy the following two term recurrence relation P 0 (x) = 0, P 1 (x) = 1, P j (x) = 2xP j-1 (x) -P j-2 (x), j ≥ 2.

Comparing with the relation [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF], it can be easily seen that for all j ≥ 1, P j (x) = U j-1 (x) = T j (x)/j. Hence, γ s 2 = ω2 2 + α bs sb1 T s (ω 0 ) = ω2 2 + (1 -ω2 2 )ω 2 b s T s (ω 0 ) = 1.


Marine Soler

Celine Barchasz

Vasily Tarnopolskiy

Frédéric Le Cras

Frédéric Le

Evaluation of the Chemical Stability of Sulfide Solid Electrolytes with Lithium Metal Using In-Situ, Operando and Post-Mortem Analyses

Evaluation of the Chemical Stability of Sulfide Solid Electrolytes with Lithium Metal Using In-Situ, Operando and Post-Mortem Analyses

The high demand for portable energy storage is driving the development of new, more efficient and safer portable energy storage systems. A lithium-sulfur battery is based on a lithium (Li) metal anode and an elemental sulfur cathode. It is a very promising system thanks to the very high theoretical energy density, while sulfur is also an inexpensive, abundant, non-critical and non-toxic raw material. On the other hand, lithium-sulfur batteries with conventional liquid electrolytes suffer from several problems preventing them from being commercialized today. The development of so-called "all-solid-state" lithium-sulfur batteries with the use of a highly conductive solid electrolyte (SE) could resolve some of these difficulties [START_REF] Wang | Research Progress of the Solid State Lithium-Sulfur Batteries[END_REF] . In particular, solid sulfide electrolytes (SSEs) have received increasing interest in recent years due to their promising ionic conductivities around 10 -2 -10 -3 S/cm [START_REF] Busche | In Situ Monitoring of Fast Li-Ion Conductor Li 7 P 3 S 11 Crystallization Inside a Hot-Press Setup[END_REF] and their ductility. However, while it solves some of the problems associated with liquid electrolytes, the all-solid-state system also presents a number of specific challenges.

One major challenge of the all-solid-state systems is the study and the control of the interfacial issues, in particular of the SE/Li metal anode interface. Indeed, stability issues of the sulfide-based electrolytes with Li metal may be encountered, as being thermodynamically unstable against lithium metal, as well as the possible growth of lithium dendrites through the solid electrolyte separator [START_REF] Wu | Interfacial Behaviours between Lithium Ion Conductors and Electrode Materials in Various Battery Systems[END_REF][START_REF] Yue | Progress of the Interface Design in All-Solid-State Li-S Batteries[END_REF] . Since both of these phenomena can severely impact the performance and safety of the battery (poor electrochemical stability of the SE against Li metal leads to decomposition of the SE thus increasing the impedance; growth of lithium dendrites leads to short-circuit of the battery), they must be understood and suppressed.

In this work, we study two sulfide electrolytes: argyrodite-type crystalline Li6PS5Cl and glass-ceramic Li7P3S11. In order to assess and compare their stability against Li metal, Electrochemical Impedance Spectroscopy measurements were carried out over time on symmetric Li/SSE/Li cells. While the impedance is increasing constantly in the case of Li7P3S11, it remains relatively stable for argyroditetype SE. In addition, the nature and the morphology of the decomposition products of the two electrolytes in contact with Li metal have been characterized post mortem by XPS, ToF-SIMS and SEM analyses, as well as in-situ and operando using synchrotron radiation to perform X-Ray Diffraction Computed Tomography and Absorption Tomography. We demonstrate that in the case of argyrodite, the decomposition products limit the further reaction of electrolyte with Li metal, whereas they promote the propagation of the decomposition reaction for Li7P3S11.