Gregory Di

Guillaume Prenat

Bernard Dieny
email: bernard.dieny@cea.fr

InMRAM: Introductory course on Magnetic Random Access Memories for microelectronics students and engIneers

Magnetic Random Access Memories (MRAM) interest is growing fast in the microelectronics industry. Commercial MRAM products already exist and all major industrial players have launched large R&D efforts to bring Spin Transfer Torque MRAM to production. The principle goal is to replace part of the memory hierarchy (in particular DRAM) by MRAM below the 20nm node. Interest in bringing non-volatile elements closer to the logic is also growing, in order to improve the performance of electronic circuits by increasing the bandwidth between logic and memory and reducing the power consumption. Despite this nsmg interest, very few microelectronics engineers have a background in magnetism and even more in spinelectronics so that it is very difficult for them to get into this emerging field. Magnetism and microelectronics communities have worked separately so far. We are convinced it is time to develop more relationships between these two communities through joint symposia, workshops, conferences, as well as to introduce some courses in magnetism and spin electronics in the education of young microelectronics engineers. In this context, SPINTEC, a French lab associated to CEA/CNRS/Grenoble University, started organizing last year an annual Introductory Course on MRAM (InMRAM) for students, engineers and researchers having no or little background in magnetism. The second edition will take place in Grenoble from 2nd to 4th July 2014.

I.

INTRODUCTION

Spinelectronics is a merging between magnetism and electronics. [t started in 1988 with the discovery of Giant Magnetoresistance (GMR) [START_REF] Baibich | Giant Magnetoresistance of (OOl)Fe/(OOl)Cr Magnetic Superiattices[END_REF]. The first application of GMR has been as very sensitive magnetoresistive sensors called spin In a second step, as the switching speed of MRAM will increase for instance thanks to precessional srr switching, MTJ may be used even closer to the logic units. In a third step, new logic-in-memory architectures may be introduced wherein the memory is distributed with a much finer granularity within logic circuits.

Although most of the interest in MRAM is now focused on STT-MRAM, other companies in particular Crocus

Technology are still developing MRAM based on field induced magnetization switching but in combination with thermal assistance. The idea is to take advantage of the Joule heating produced in the MTJ around the tuunel barrier when a current is flowing through it. At current densities of the order of a few can be used to enable the magnetization switching with a weak magnetic field [13]. Thanks to this combination of field and heating, the field pulses can be shared between various bits (for instance writing all "0" in a 64 bit word with a single pulse of magnetic field, then all "I" with a second pulse), thereby strongly reducing the power consumption associated with field writing. This approach has very nice perspectives of applications in the field of security of data, but also for Content Addressable Memories (Fig. 3) which are used in search engines, routers ... [I.

HOW TO FAVOR LINKS BETWEEN MAGNETISM AND MICROELECTRONICS COMMUNITIES?

[t is interesting to note that the places where the R&D in spinelectronics and MRAM has grown the fastest and the strongest are in places where expertise in both magnetism and microelectronics were present in close proximity. This is true for instance at Tohoku University in Japan or in Grenoble France. This is also true in companies which are producing or have produced both electronic devices and magnetic hard disk drives (Toshiba, Samsung, IBM, Hitachi ... ).

At the international level, it is important to organize joint events gathering magnetism and microelectronics experts.

  valves [2] which have been used as read heads in Hard Disk Drives between [996 and 2004. Two other major breakthroughs have further stimulated this field. The first one was the observation of large tunnel magnetoresistance (TMR) at room temperature in magnetic tunnel junctions (MTJ), first in amorphous Alumina based MTJ (TMR amplitude in the range 40-70%) [3, 4] and more recently in crystalline MgO based MTJ (TMR amplitude in the range [00-600%) [5, 6]. MTJ s consist of two magnetic metallic layers separated by a thin tunneling barrier. When a bias voltage is applied across the stack, the electron tunneling probability depends on the relative orientation of the magnetization in the two adjacent magnetic layers. [n most common cases, the tunneling probability is 978-1-4799-4016-5/14/$31.00 ©2014 IEEE higher when the magnetizations are parallel (yielding a minimum resistance) and lower when they are antiparallel (yielding a higher resistance). Somehow, the magnetic layers play the role of polarizer and analyzer for the spin of the electrons. MTJs offer very nice possibilities of integration with CMOS components for several reasons. Their resistance can be tuned by adjusting the thickness of the tunnel barrier to be of the order of a few kQ i.e. compatible with the typical resistance of MOSFET. Furthermore, they can be grown above the CMOS technological levels and electrically connected to the CMOS components through Cu via. [n this integration, there is no contamination of the CMOS components. The MTJs only require annealing at temperature of the order of 300°C compatible with the CMOS process. The MTJ can be used as variable resistance controlled by magnetic field or by current/voltage thanks to the spin transfer torque phenomenon (see further). Soon after this discovery of large TMR at room temperature, it was realized that MTJ could be used to make arrays of memory elements called MRAM (Magnetic Random Access Memories). A significant R&D effort then started on MRAM written by pulses of magnetic field, mostly led by Freescale/Everspin. This company successfully launched a 4Mbit MRAM product in 2006 which is used in microcontrollers, automotive applications or space applications (hardness to radiations). However, this technology is not scalable beyond the 65nm technology node because of the large current required to generate the pulses of magnetic field during write. This issue has been circumvented thanks to a second major discovery: the spin transfer torque (STT) phenomenon [7, 8]. This is the action that a spin polarized current exerts on the magnetization of a magnetic nanostructure due to the exchange interaction between the spin of the conduction electrons and the spin of the electrons responsible for the local magnetization. If the current density is large enough, this phenomenon can be used to switch the magnetization of a magnetic nanostructure. [t is particularly useful in the context of MRAM to switch the magnetization of the storage layer. This phenomenon has been a major discovery because it provides a new way to manipulate a magnetization with much less energy and in a much more localized way than with a magnetic field [9]. STT induced magnetization switching has been first observed in fully metallic magnetic pillars [9] but also, after some optimization, in MTJs [10]. In STT, the magnetization switching is determined by the current density flowing through the MTJ. As a result, the current required to write in STT MRAM scales as the MTJ area down to very small values of the order of lOI-lA where it reaches a bottom limit set by the memory retention. This is very favorable for the downsize scalability of MRAM based on STT switching. This property renewed the interest for MRAM since the possibility for this technology to be scalable to and beyond the 20nm node then became clearer. Now most of the major companies in microelectronics have launched significant R&D efforts in STT MRAM (Samsung, Hynix/Toshiba, IBM, TDK, NEC, Qualcomm, TSMC ... ) and a number of start-up companies have emerged in this field Transfer Technology, Avalanche ... ). The interest for MRAM has been strongly rising over the past few years as illustrated in Figure 1 which represents the number of papers and patents published per year containing the word "MRAM". These data are derived from the scirus database (www.scirus.com). MRAM have now a clear place in the ITRS. The ITRS ERD and ERM working groups have also published a report in July 2010 in which "STT MRAM and Redox RAM were identified as emerging memory technologies recommended for accelerated R&D leading to scaling and commercialization of non-volatile RAM to and beyond the 16nm generation" [11]. The goal is here to replace parts of the memory hierarchy by non-volatile MRAM as illustrated in Fig.2. This would have lot of advantages in terms of reducing the power consumption (both static and dynamic), reducing the footprint on the wafer, speeding up the communication between logic and memory functions, introducing ultrafast reprogrammability capabilities in electronic circuits. In a further step, the conventional memory hierarchy could even be revisited with new logic-in memory architecture wherein the memory would be distributed with a much fmer granularity within the logic circuits as illustrated in Figure 2. These architectures would allow a Normally-Off/Instant-On electronics with much reduced power consumption and faster bandwidth between logic and memory [12].

I

  Fig. I. Number of publications and patents per year containing the word "MRAM

  Fig. 3. concept of Magnetic Logic Unit (MLU) as developed by Crocus Technology (http://www.crocustechnology.comlpdf/Crocus MIP White Paper v6.pdf) . A MLU is a MT J comprising a soft reference layer in contrast to conventional MTJs which have a pinned reference layer. MLUs intrinsically combine a storage function together with a comparison (XOR) function. This is called the Match-In-Place function. MLUs are very useful for a number of applications in particular in NAND type of chains used for Content Addressable Memory (CAM) applications. The purpose of a CAM is to find at which address a particular word is located. The example above illustrates how to find the written word "1010" in a set of 4 addresses. By sending pulses of current corresponding to "1010" in the field generating lines located below the MTJ, the magnetization in all soft layers along each line get oriented along the corresponding field direction. Only the chain containing the searched word has the minimum resistance since all the MLU elements along this chain are in parallel magnetic configuration i.e. are in their minimum resistance state (3, d chain from the top in this example).

o

  Several initiatives have been taken along this line in the past couple of years. For instance, a one-day workshop was organized in May 20 I 0 at Hsinchu Taiwan, the eve of the start date of the VLS[ and [ntermag 2010 conferences which happened to take place the same week in Taiwan. This has been an excellent opportunity for comparing FLASH, ReRAM and MRAM technologies and getting reciprocally educated. Conferences such as IEDM, [MW, VLS[ have clearly expressed their willingness to include more MRAM contributions in their technical program. So far, most of the MRAM reports are presented in magnetism conferences (Intermag and MMM). We should encourage the teams working on MRAM to also present their reports in microelectronics conferences. The IEEE Magnetic Society is very supportive of this endeavor. [t has for instance sponsored a half-day short-course on MRAM at IMW2012. Reciprocally, it is important for magnetism experts working in MRAM and hybrid CMOSlMagnetic technology, who are most of the time physicists, to [earn more about microelectronics so as to better orient their R&D efforts. III. [NMRAM: [N TRODUCTORY COURSE ON MRAM FOR STUDENTS AND ENGINEERS IN MICROELECTRONICS Magnetism is playing growing role in industry in general. Traditionally, the most important areas of use of magnetism has been in hard disk drives industry, in permanent magnets used for instance in electric motors, soft magnets used in inductors ... in magnetic sensors implemented in a large variety of applications (automotive industry, robotics, electronic compass, ... ). Other areas of applications have more recently emerged related to biotechnology (magnetic particles for diagnostics or treatments, medical imaging ... ), or energy harvesting. Besides, MRAM and more generally hybrid electronic circuits combining CMOS and magnetic technologies is becoming a major direction of R&D at the merge between magnetism and microelectronics. Its potential of industrial development is huge. In this regard, the worldwide number of engineers educated in magnetism and spinelectronics looks quite small. [t seems time to introduce courses in magnetism and spinelectronics in electrical engineering education. This would give to the young engineers in this field the basis for understanding the materials, the phenomena, the various technical issues, related to MRAM technologies. To contribute in this direction, SPINTEC, which is a French lab affiliated to CEA (Commissariat a l'Energie Atomique et aux Energies Alternatives), CNRS (Centre National de la Recherche Scientifique) and Grenoble University, organizes an annual Introductory Course on MRAM (InMRAM) whose first edition took place in 2013. The second edition will take place from July 2il to 4 th 2014 in Grenoble, France. This summer school addresses several topics Combining the best of CMOS and magnetic worlds. e.g., magnetic full adder, content addressable memory, M-FPGA, ... o Logic-in-memory architecture o Normally-OFFlInstant-ON computing. o Low power electronics based on spintronics o Benefits from radiation hardness [nMRAM speakers are experts coming from [nstitut Neel, CNRS/Thales, CNRS, CEA, Yole and University of Montpellier 2. More details are given at www.inMRAM.com The [nMRAM first edition feedbacks from the 1 [2 attendees were extremely positive. The attendees came from universities, research laboratories and companies as illustrated in Figure 4. This figure also shows that about 2/3 of the attendees were originating from the microelectronic community with no magnetism background and 113 from the magnetism community. Concerning the geographical origin of the attendees, the latter was quite diverse. Approximately half of the audience was from France and the other half from 20 different countries, as illustrated in Figure 5.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Repartition of attendees in terms of (a) professional origin and (b) scientific/technical background.

A/cm2, the MTJ temperature can increase by 200°C in a few ns. By properly designing the MTJ stack, this temperature rise

ACKNOWLEDGMENT

Liesl Folks, Bruce Terris, Ron Goldfarb and Takao Suzuki from the IEEE Magnetic Society are acknowledged for their continuous support in this initiative to foster more relationships between magnetism and microelectronics conferences. This effort is also supported by the European Commission through the ERC Adv grant HYMAGINE (ERe 246942).


Emanuele Congiu

Marco Minissale

Saoud Baouche

Audrey Moudens

Stephanie Cazaux

Giulio Manicò

Valerio Pirronello

Francois Dulieu

Henda Chaabouni

François Dulieu

Efficient diffusive mechanisms of O atoms at very low temperatures on surfaces of astrophysical interest

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the cold regions of the Universe, where temperatures are lower than 10 K and densities are far weaker than those attainable on earth, a rich chemistry is initiated on the surfaces of minuscule dust particles. [1][2][3][4] The species weakly bound to the surface, most of all, play a central role in the evolution of the pristine chemistry governed by the diffusion of reactive species. 5 In space, thermal atom-addition induced chemistry occurs mostly at low temperatures ( 10 K) newly formed species are protected from radiation to a great extent by dust particles. These regions are parts of collapsing envelopes that feed young stellar objects and that provide the original material from which comets and ultimately planets are made. 6 Hydrogenation of interstellar ices can induce the formation of new species in the solid phase and, therefore, it has been the topic of recent laboratory-based studies. [7][8][9][10] The efficient surface formation of the bulk of interstellar ices, i.e., water, methanol, formaldehyde, and formic acid has been demonstrated through H-atom additions of CO-and/or O 2 -ices under interstellar relevant conditions. [11][12][13] In particular, the formation of the most important and abundant ice of all, amorphous solid water (ASW), proved to be the principal product of all the possible chemical pathways involving O and H atoms/molecules or the hydroxyl radical (O 3 , O 2 , O, OH, H and H 2 ). [14][15][16][17][18][19][20] So far, diffusion has only been explored experimentally for H atoms. 21,22 Nevertheless, if we assume that only H atoms are mobile and can diffuse over the dust grain surface, it becomes difficult to meet the observational evidence for CO 2 -the second most abundant condensed species in grain mantles -as well as for the rich molecular diversity for the general interstellar medium. CO 2 is believed to form in the solid phase via several energetic [23][24][25][26] and non-energetic [23][24][25][26][27][28][29][30][31][32] mechanisms, i.e., CO+OH and CO+O reactions at 10 K. The reaction CO+OH leading to CO 2 , however, occurs in competition with H 2 O formation via the H+OH (barrierless) pathway, 30 which is a much faster route whenever H is the only mobile species. For this reason, the CO+OH reaction alone cannot account for the CO 2 ice observed in quiescent regions. 33,34 Furthermore, it seems that a sort of depth (i.e., age) segregation exists between the three most abundant ices (H 2 O, CO 2 , CO): 35,36 water tends to be concentrated in the layers forming the inner (and older) part of the mantles, while CO 2 and CO abundances increase in the outer (and more recent) layers. CO accretes onto icy grain surfaces at a rate proportional to the gas density, 37,38

Experimental

The FORmation of MOLecules in the InterStellar Medium (FORMOLISM) experimental set-up has been developed with the purpose of studying the reaction and interaction of atoms and molecules on surfaces simulating dust grains under interstellar conditions (the relevance of substrate, low density, and very low temperatures 10 K). ∼ 42 FORMOLISM is composed of an ultrahigh vacuum chamber with a base pressure of 10 ∼ -10 mbar, a rotatable quadrupole mass spectrometer (QMS) and an oxygen-free high-conductivity copper sample holder. The sample holder is attached to the cold finger of a closed-cycle He cryostat and can be cooled to 6 K. The temperature is measured with a calibrated silicon diode clamped on the sample holder and controlled by a Lakeshore 334 controller to ±0.2 K with an accuracy of ±1 K in the 8-400 K range. The apparatus is also equipped with a reflection-absorption infrared spectroscopy (RAIRS) facility used to probe the deposited or produced species in situ. 7 Reactants are introduced into the vacuum chamber via two separated triply differentially pumped beam lines aimed at the cold surface. molecules exiting the source are in their ground state 3 P and X 3 Σ g -, respectively. The beam flux was calibrated by using temperature-programmed desorption (TPD) to determine the O 2 exposure time to saturate the O 2 monolayer.

Five surfaces were investigated in this study: porous ASW (H2O (p) ), non-porous ASW (H2O (np) ), crystalline ice (H2O (c) ), amorphous olivine-type silicate (SiO x ), and an oxidised slab of highly oriented pyrolytic graphite (HOPG). The SiO x , HOPG, and other carbon-based surfaces, mimic bare dust grains in molecular clouds and have been previously used in investigations of heterogeneous chemistry on dust grain analogues. The silicate sample is amorphous in nature, as evidenced by infrared spectroscopic studies, 43 while TPD experiments reveal the surface to be non-porous on the molecular scale. 44 The HOPG surface used in the experiments is a ZYA-grade HOPG sample, which had been previuosly exposed to an O-atom beam (oxidised) to avoid surface changes during the experimental sequences. For the water substrates, ice films were grown on top of the silicate surface by spraying water vapour from a microchannel array doser located 2 cm in front of the surface. 45 To form H 2 O (c) , the surface was held at 120 K during the deposition, then flash heated at 50 K min -1 to 145 K. For each type of ice surface, the temperature was held constant until the background pressure in the chamber stabilised, before cooling it back down in the range 6 -25 K, at which temperature O atoms were dosed onto the respective surfaces.

Results

In Fig. 1 we show the type of raw experimental data we obtain: O 2 and O 3 mass spectra (TPD performed at a heating rate of 10 K min -1 ) are recorded after deposition of a O+O 2 dose on graphite at a given surface temperature (Fig. 1, left panel). O 2 desorbs between 27 and 50 K, while ozone desorption is observed between 55 and 75 K (directly at mass 48 a.m.u., or via the O 2 + fragments at mass 32 a.m.u.). O (16 a.m.u.) desorption was never observed. The area of the peaks (proportional to the amount of the species formed on the surface) changes depending on the coverage (O+O 2 dose) and on the surface temperature. at several O+O 2 initial doses. In two previous papers 39,40 dealing with the diffusion of oxygen atoms on water ice and amorphous silicate, we showed the complete series of TPDs performed at several O+O 2 coverages and at various surface temperatures. From the resulting O 3 and O 2 yields, we inferred that the O 3 /O 2 ratio increases with initial coverage, as an incoming O-atom is more likely to find O 2 molecules at higher coverages, and with surface temperature because at higher temperatures the mobility of O atoms is favoured and ozone formation is more efficient.

If we focus on the diffusion dependency on T s , and at coverages less than 0.5 ML -more suitable to the astrophysical context -the LH mechanism can be fairly considered the main process governing the surface reactions involving oxygen atoms and molecules. One may argue, however, that our experiments -based on linear thermal ramps starting from deposition temperature T s -can affect the diffusion of O atoms and the actual O 2 and O 3 yields. To rule this possibility out, we planned an experiment in which TPD and RAIRS results could be compared in order to confirm or discard the possible role of the heating ramp in O diffusion. In Fig. 2 we present the amount of ozone formed after depositing a given dose of O+O 2 on graphite. The experiment was performed using two different procedures. First, after a 0.3 ML of O+O 2 was dosed at various temperatures, a RAIR spectra was recorded then a TPD was started for each deposition. Blue squares in Fig. 2 indicate the TPD ozone yields after deposition of 0.3 ML of O+O 2 on oxidised HOPG kept at 6.5, 10, 15, 20 and 25 K. Green squares were obtained by integration of the band area at 1043.5 cm -1 of ozone, recorded after each deposition. Second, after depositing 0.3 ML of O+O 2 on graphite at 6.5 K, we monitored the evolution of the ozone IR band intensity with temperature. By this second method we could estimate the formation of ozone during the thermal ramp. The result of this experiment is shown by the black circles in Fig. 2. They give the ozone IR band area as a function of temperature after a single 0.3-ML dose of O+O 2 at 6.5 K. For temperatures greater than 40 K, the O 3 band does not increase because O 2 desorbs and the O+O 2 reaction can no longer take place. The data presented in Fig. 2, in summary, show the yield of ozone formed following a given dose of O+O 2 deposited at different T s between 6 and 25 K (T s greater than 25 K would make O 2 mobile as well, hence add another degree of complexity to this study). It is clear that ozone formation efficiency grows fast with the deposition temperature (see squares in Fig. 2), while the contribution to ozone provided by the thermal ramp -which is likely to induce diffusion of the residual O atoms -is very small (see black circles in Fig. 2). In addition, if we consider the IR data at 15 K, the O 3 yield obtained after O+O 2 deposition at 15 K (green square) is much higher than the ozone yield after deposition at 6.5 K and heated to 15 K (black circle). This fact confirms that all the chemistry has occurred at the deposition temperature; and if it is true at Ts = 6.5 K, then it is the case at any other temperature higher than 6.5 K.

Experimental data are then inserted into a model composed by a series of rate equations used to simulate the O 2 and O 3 formation yields according to coverage and surface temperature. The model includes both LH and ER mechanisms, and it allows reactions to occur during the deposition phase, as well as during the heating phase (TPD). A complete account of our model is given in Minissale et al 2014. 40 Here we will focus on the diffusion rates k and the different methods by which they are calculated. We already alluded to the fact that reactions mostly occur during the exposure phase. The diffusion of atoms during the heating phase is small because not more than a few percent of the deposited O atoms remain available on the surface in the low coverage regime. The effect of a possible diffusion during the TPD lies within the error bars of the experimental data, and can be neglected. For this reason, in what follows, we will address only the diffusion constants at a fixed temperature for each one of the substrate investigated.

The diffusion coefficients k include two components due to quantum tunnelling and thermal motion:

46 k = k qt + k tm
In our model, k can be treated as a free numerical parameter during the deposition phase at constant temperature, owing to the fact that the evolution of the coverage with time is known and provides a strong constraint. Therefore, resulting k values are a set of constants giving the diffusion rate at given temperatures, although no information can be inferred about the nature of the diffusive process. In

The diffusion coefficients given by eqn (1) and eqn (2) provide the diffusion probability of exploring a fraction of 1 ML per unit time, and can be converted into the usual units cm 2 s -1 (Fig. 3, Fig. 4, and Fig. 6) by assuming that there are 10 15 sites cm -2 .

Fig. 4 displays a fit of the experimental values obtained on amorphous silicate according to the empirical law given in Eq. (1). k 0 (s -1 ) can be considered the minimum value of k, or the value k must tend to near T = 0 K. α is a free parameter with values between 0 ad 1 K -1 s -1 , and it accounts for diffusion efficiency differences between the various substrates. α is one for water ice while is about 0.1 for graphite and silicate. The dependency on the surface temperature is governed by the factor (T/10) β , the exponent β can have a value between 3 and 4, with variations due to the surface nature, although the best fits give typical values of β ~ 3.5.

Case b), the classical Arrhenius law used to model the diffusion coefficients k is

k Arr (T) = ν 0 exp[-E diff (T)/T]. ( 2
)
E diff is the diffusion barrier expressed in kelvins (eV/k b ) and ν 0 (= 10 12 s -1 ), the pre-exponential factor, can be seen as a trial frequency for attempting a new event. In Fig. 5 we present a fit of the diffusion coefficients k on non-porous ASW obtained by using the Arrhenius law. Fig. 5 actually displays the activation energies for diffusion (E diff ) as a function of temperature. In fact, according to Eq. ( 2), a suited set of E diff can be used to derive one diffusion coefficient for each temperature. It is thus possible to link each of these diffusion coefficients to an Arrhenius behaviour, and find one energy barrier at each temperature as shown in Fig. 5 (see also dashed lines in Fig. 3). It should be noted, however, that an Arrenius-law form in which E diff is fixed (indipendent of T), or where a distribution of E diff is given, is not able to fit the data. This is why we discarded the Arrhenius-type behaviour of k as it made no physical sense to us. In fact, a systematic increase of the Arrhenius barrier with temperature seemed to us an ad hoc solution. Also, this implies that at low temperatures (~ 6 K) diffusion occurs through low diffusion barriers (e.g., E diff = 170 K). If such low barriers actually exist, they represent fast connections between adsorption sites. Why then would these low energy barriers vanish at higher temperatures? To put it in other terms, why and how atoms would diffuse through slow pathways (high diffusion barriers) at high temperatures (~ 20 K), if faster pathways exist? We consider this unlikely and not physically reasonable.

In Fig. 6 we show a comparison between the classical behaviour (described by an Arrhenius-type law) of the diffusion coefficients as reported by Karssmeijer et al 47 for CO molecules on hexagonal water ice, and the trend that we find experimentally for O atoms on amorphous silicate and crystalline water ice. It is clear that our experimental values do not follow an Arrhenius behaviour, suggesting that a classical description is incomplete to explain the experimental data (squares and triangles in Fig. 6). In fact, in a pure thermal diffusion the slope is very different, and if we fit the data by using a classical Arrhenius law, we find values of ν 0 and E diff not physically acceptable.

Therefore, a quantum mechanical approach ought to be used to account for the deviations from the classical trend. Our results on oxygen atoms are consistent with a tunnelling-dominated diffusion found for H atoms on H 2 O (np) by Senevirathne et al 48 in the 6 -13 K temperature range (the slopes of H-and O-diffusion constant behaviours are similar). They also found that diffusion of H atoms is enhanced around 13 K, as occurs to O atoms around 22 K, just where classical thermal motion begins to predominate over quantum processes. 39,48 Hama et al 49 found that this temperature border between quantum and classical diffusion of H atoms is likely to be at T s < 10 K. We found that at very low temperatures the diffusion of O atoms is better simulated by quantum tunneling through a square barrier. 39,[START_REF] Messiah | Quantum Mechanics[END_REF] The physical parameters we use to describe such a quantum jump are the width a and the height E a of the barrier. The choice of a square barrier, the simplest shape of a potential, was made on purpose to show that the right trend is obtained if one uses quantum-tunneling diffusion, not because we believed that a square barrier was the right one. We believe that any other more realistic potential shape we could use, would not fundamentally change the results, and it would still be unrealistic given the complexity of the distribution of diffusion barriers. We did not try to obtain the best fit of our data, but tried to show that the right trend is obtained if one uses quantum-tunneling diffusion (see solid lines in Fig. 3). Hence, we chose to model the quantum diffusion with two parameters which have a simple physical meaning, although they correspond to macroscopic values that come from the interplay of many microscopic different situations.

The values of k and of all the parameters used to fit the diffusion coefficients on each substrate, using the three methods, are listed in Table 1.

The diffusion coefficients of O atoms calculated on water ices are one order of magnitude greater than those found on amorphous silicate and oxidised HOPG, namely O diffusive mechanism is more efficient on icy grains. Also, as opposed to the case of H atoms, there is no difference between the efficiency of O mobility on the three types of water ices investigated (H 2 O (p) , H 2 O (np) , and H 2 O (c) ). In the light of our experimental results, we can only observe and simply report this finding. In fact, to deal with atoms makes it very difficult to derive key parameters such as the energy barrier for diffusion, or even the energy barrier for desorption, hence no pertinent assumption can be made to explain these findings from a physicochemical point of view. However, to give a physical explanation of our results is beyond the scope of this paper, since we believe that quantum calculations and simulations will be necessary to thoroughly describe O diffusion mechanisms at low temperatures.

Astrophysical implications

As far as the diffusion of O atoms is concerned, it turns out that, whenever a diffusive process exists, this has an impact on the chemistry occurring at the surface of dust grains. In fact, either the formation of some species may be enhanced or at least the relative abundances of the final products is affected if O diffusion is efficient. An important example of how O-atom mobility can module the abundances of key species of ices in the ISM is the case of the H 2 O/CO 2 ratio.

In dense quiescent molecular clouds, hydrogen atoms have always been thought to be the only mobile species on the surface of icy grains. Most of the molecular variety observed in interstellar ices has long been considered the outcome of H-atom addition reactions involving O, O 2 , O 3 , CO, N, and NO. Water formation, for example, is the final and most stable species of all the chemical network between H and O, O 2 and O 3 , which justifies its role as the most abundant ice in the Universe. If the reactive partner of H is CO, then CH 3 OH is obtained via a series of successive hydrogenations. On the other hand, if H is the only mobile species able to scan the entire surface of the grain, 21 it may be difficult to explain the abundance of CO 2 , the second most abundant condensed species. CO 2 can also be formed via energetic processes by irradiating ice mixtures of H 2 O and CO with UV photons or ions. In the dense core of molecular clouds, however, these processes may not apply, and CO 2 can only be formed via non-energetic mechanisms, i.e., the reactions CO+OH and CO+O. If these chemical routes leading to CO 2 involved only species not mobile at 10 K, then CO 2 formation would be greatly hindered by the rate of accretion and the high mobility of H atoms, able to reach CO, OH, and O long before these species can meet to form carbon dioxyde. Our present and previous works introduce strong arguments to believe that O atoms too are mobile at very low temperatures. This implies that the formation rate of CO 2 in dense clouds is governed by a balance between the accretion rate of H atoms and the diffusion rate of O atoms on the surface of dust grains. The cartoon in Fig. 7 shows that when the accretion rate of H atoms is dominant, H 2 O and CH 3 OH are for the most part the final products; when the diffusion rate of O atoms prevails, formation of CO 2 (and O 3 ) is favoured.

With this in mind, we made some calculations to show the evolution of the relative abundances of H atoms and O atoms on the surface of dust grains and -assuming that both species are mobile at low temperaures -how this balance can affect the chemistry within interstellar clouds of various densities. In fact, different environments are characterized by different densities, the abundances of species in the gas phase change and this entails a change in the accretion time-scales of particles on dust grains. In diffuse clouds, hydrogen is mainly in its atomic form and is by far the most abundant atomic species. In dark clouds, hydrogen is mainly in its molecular form, so H atoms become a rather rare reactant with [H]/[H 2 ] 10 ∼ -3 (see, e.g., Li & Goldsmith). [START_REF] Li | [END_REF] The number density of H atoms is mostly governed by the destruction of H 2 due to cosmic rays. This value, almost indipendently of the density of the cloud, is of the order of 1 H cm -3 . On the other hand, the [O]/[H 2 ] ratio remains approximately constant (10 -4 ), thus the number of atomic O, unlike H, is proportional to the density of the cloud (see, for example, Table 1 of Caselli et al). 52 For a cloud with number density of 10 In Fig. 8 we show the time interval between two impacts of particles of the same species (H or O) on a single dust grain, as a function of the density n of the cloud. The time intervals between two arrivals are derived from the actual particle flux of a given species. The interstellar flux of species accreting on dust grains can be calculated as follows:

Ф x = 1/4 n x v x ( 3
)
where n x is the density of species x in gas phase and v x = (8k b T/πm x ) 0.5 is its mean velocity. Ф x is thus expressed in particles cm -2 s -1 . For our calculation, we can approximate the dust grains to spheres with typical radius r = 0.1 μm, with accessible surface area A = 4 πr 2 . The time interval between the impacts of two particles then is

t = (n x v x A/4) -1 . (4)
In Fig. 8, the grey solid line rappresents the time interval between the impact of two hydrogen atoms, calculated by assuming a constant density of H atoms n H = 2.3 cm -3 (from Li & Goldsmith). [START_REF] Li | [END_REF] The density of O atoms is proportional to the density of the clouds n, namely, n O = 5•10 -4 n. The time interval between the arrival of two O atoms is displayed as a red solid line, which clearly shows that arrivals of oxygen atoms become more frequent (shorter time between two impacts) with the density of the cloud. The grey and red lines cross around a density n of 10 4 cm -3 . This suggests that for cloud densities of ~ 10 4 cm -3 the accretion rates of H and O are comparable and, given that both species can diffuse, oxidation reactions on grains may play a role, although H-atom additions are still dominant owing to the higher mobility of H. In Fig. 8 we also indicate the mean time O atoms need for completing a scan of all the adorption sites on the surface of one typical grain used above, with radius = 0.1 μm and 10 6 absorption sites (10 15 sites cm -2 ).

Mean times needed for a complete scan of the grain surface were calculated for a surface temperature of 10 K by using the diffusion constants k of O atoms on each substrate presented in this work, taking into account that k = 10 -15 cm 2 s -1 corresponds to one jump per unit time. As to H atoms, the mean time for scanning the entire surface of water ice was derived by the energy barrier for diffusion of 255 K (at 10 K)

given by Matar et al. 21 Again, in Fig. 8, it is interesting to observe the intersection occurring at n ~ 10 5 cm -3 between the red line and the band giving the mean time H atoms employ to scan the whole surface of the grain. This implies that at cloud densities of ~ 10 5 cm -3 or greater the diffusion and accretion rate of H atoms are smaller than the accretion rate of O atoms. Therefore, in very dense clouds oxygen atoms may become the dominant reaction partner able to react with CO and produce CO 2 , as well as with H and produce OH. Since H atoms are rare in this environments, OH will not be readily transformed into water via hydrogenation, and also the hydroxil radical is likely to react with the abundant CO molecules to form CO 2 .

Fig. 1 ,

 1 Fig. 1, right panel: RAIR spectrum of the ν 3 -ozone band at 1043.5 cm -1 recorded at 6.5 K after deposition of 0.3 ML of O+O 2 on graphite at the same temperature. O 3 and O 2 yields are calculated after depositions of O+O 2 performed by a) varying the dose (coverage) in the range 0.05 -1 ML and keeping the surface temperature (T s ) fixed, and b) dosing with a fixed amount and varying the deposition temperature between 6 and 30 K. The acquisition of two sets of data is motivated by the nature of the mechanisms through which surface reactions leading to O 2 and O 3 proceed. The two main mechanisms invoked for surface catalysis are the Langmuir-Hinshelwood (LH) mechanism, in which O 2 and O 3 form in an O-atom diffusion process, and the Eley-Rideal (ER) mechanism, or 'prompt' mechanism, in which an impinging O-atom reacts directly with an adsorbed O or O 2 . The LH mechanism is highly dependent on the surface temperature as it affects the mobility of species on the surface. This justifies depositions of O 2 and O 3 at different surface temperatures. The ER mechanism, on the other hand, is indipendent of T s , becomes more efficient with the increase in surface coverage, and therefore we also investigated the O 2 and O 3 formation

Fig. 3

 3 Fig.3the diffusion constants k that we obtained for various substrate compositions are plotted as a function of temperature. An important finding of this comparative study is that diffusion coefficients on water ices (regardless of its morphology) are about one order of magnitude greater than those on silicate and graphite.Diffusion coefficients k vs T s can be displayed in several ways according to the law used to describe them, namely, k may have a) an empirical law built for fitting the experimental values, b) an Arrhenius-law form, with an activation energy E diff free to vary, or c) a quantum-tunnelling form with a width and height of the barrier. A detailed analytical or numerical solution of the dependence of k with T s can help have some insight into the physical nature of the diffusive process at play.

Fig. 1

 1 Fig. 1Left panel: O2 and O3 TPD traces obtained after deposition of 0.3 ML of O+O2 on oxidised HOPG held at 6.5 K. Right panel: RAIR spectrum recorded after deposition of 0.3 ML of O+O2 on oxidised HOPG at 6.5 K; the absorption band at 1043.5 cm -1 is due to the ν3 asymmetric stretching mode of O3.

Fig. 2

 2 Fig. 2 Ozone yields derived from TPD peak areas (blue squares) and RAIRS ozone ν3-band area (green squares) after deposition of 0.3 ML of O+O2 on oxidised HOPG at 6.5, 10, 15, 20, and 25 K. Black circles represent RAIRS ozone ν3-band areas obtained by depositing a unique 0.3-ML dose of O+O2 on oxidised HOPG held at 6.5 K then by recording a RAIR spectrum at different surface temperatures (6.5, 15, 25, and 40 K).The dashed and the solid red lines are fits of the experimental values and serve as a guide to the eye. All ozone yields were normalised to the O 3 yield obtained from oxidised HOPG after deposition of 0.3 ML O+O2 at 6.5 K.

Fig. 3

 3 Fig. 3 Diffusion constants k of O atoms obtained on H2O(p) (open pink circles), on H2O(np) (black circles), H2O(c) (blue triangles), amorphous silicate (red squares), and oxidised HOPG (green stars), plotted as a function of surface temperature. Dashed lines represent a series of Arrheniustype laws generated by using five values of Ediff (from 150 to 700 K). The two solid lines are best fits of the experimental values obtained through the quantum-tunneling diffusion law for O atoms on H2O(np) (black solid line) and amorphous silicate (red solid line); seeTable 1 for best fit values of a (barrier width) and Ea (barrier height).

Fig. 4

 4 Fig. 4Black squares represent diffusion constants of O atoms on amorphous silicate as a function of temperature. The red solid line is a best fit of diffusion constants vs temperature obtained by using the empirical law given in Eq. (1); see Table1for best fit values of k0, α, and β.

Fig. 5

 5 Fig. 5Energy barrier for diffusion on H2O(np) as a funcion of surface temperature in the case diffusion constants are derived from the Arrhenius-type law given in Eq. (2). The red solid line represents a linear fit of Ediff (T). A single value of Ediff cannot be a solution satisfying the whole set of diffusion constants observed in the 6 -25 K temperature range; see Table1for the interval of E diff values needed to obtain kArr between 6 and 25 K.

Fig. 6

 6 Fig. 6 Comparison between diffusion behaviours of CO molecules, H atoms, and O atoms. Blue squares and green triangle represent O-atom diffusion constants as a function of surface temperature on amorphous silicate and crystalline water ice, respectively (this work). The black dashed line represent the thermal diffusion (Arrhenius behaviour) of CO molecules on hexagonal water ice found in Karssemeijer et al. 47 The red solid line and the red dashed line display the H-atom tunnelling (6 -13 K) and the H thermal diffusion (Ts > 13 K), respectively, obtained by Senevirathne et al 48 on compact amorphous water ice. The difference between the slopes of CO and O behaviours, and the similarity between the slopes of H tunnelling and O data, corroborates the conclusion that O atoms diffuse via quantum tunnelling in the surface temperature domain between 6 and 22 K.

Fig. 7

 7 Fig. 7 Schematic view of the two possible scenarios concerning the H2O/CO2 balance on icy grains in dense molecular clouds. If H atoms (black dots) are the only mobile species, H-addition mechanisms are dominant and the formation of water (H + OH/O/O 2/O3 → H2O) and of other H-saturated species (e.g., CO + 3H → CH3OH) is favoured. On the other hand, if also O atoms (red circles) can diffuse at very low temperatures, the formation of CO2 in dense clouds may proceed via non-energetic reactions (CO+O and CO+OH) as well, making possible that CO 2 is the second most abundant ice in the ISM.

Fig. 8

 8 Fig. 8 Time intervals between two impacts of H and O, and times employed to scan a whole grain (shaded horizontal bands) on various surfaces of interest, are plotted as a function of the density of the cloud. The time interval between two arrivals of H is constant as the density of H atoms remains rather constant regardless of the density of the medium, while O atoms abundance grows with cloud density.

  , i.e., in the innermost part of the clouds where ∼

a LERMA Université de Cergy-Pontoise and Observatoire de Paris, ENS, UPMC, UMR 8112 du CNRS 5, mail Gay Lussac, 95000 Cergy Pontoise cedex, France; E-mail: emanuele.congiu@u-cergy.fr b Kapteyn Astronomical Institute, PO box 800, 9700AV Groningen, The Netherlands; c Dipartimento di Fisica e Astronomia, Università di Catania, via Santa Sofia 64, 95123 Catania, Sicily, Italy.

  almost certainly because dust grains are not cold enough in low density regions (e.g., diffuse clouds) so that CO residence time on the surface is extremely short. Conversely, why H 2 O (together with other fully hydrogenated species like CH 3 OH and NH 3 )35 is mostly found in the inner part (near to the silicate/carbonaceous core), and CO 2 in the outer layers, is not fully understood. This could actually be explained by assuming the presence of efficient diffusive processes of O atoms on the icy grains, as recent laboratory works by our group demonstrated on realistic dust grain analogues, such as water ice and silicate at very low temperatures.39,40 Also, this scenario is consistent with the CO 2 formation timeline proposed in Noble et al,41 that they tested through concurrent observations of CO, CO 2 , andH 2 O ices.Here we present new results of diffusion constants of O atoms calculated on graphite, and compare them to the diffusion constants previously calculated on amorphous silicate and three different surfaces of water ice: porous, compact amorphous, and crystalline water ice. In this comparative study, the measured

diffusion constants on each surface are given as a function of surface temperature. These data are modelled to determine a simple analytical expression that accurately reproduces the diffusion constants of O atoms with surface temperature. Also, by changing the model to incorporate both a classical (Arrhenius-type) law and a quantum tunnelling description, we are able to address some key physical questions, namely what diffusive process is at play in the 6 -25 K temperature range.

  4 cm-3 , the [H]/[O] ratio is 1/0.75, while for a denser cloud with a density of 10 ∼ 5 cm -3 , the [H]/[O] ratio is 1/7. Therefore, for very dense clouds, O is the most abundant species in atomic form, can ∼ accrete on grains and, provided O atoms are mobile, subsequently react with other species before these get fully saturated by H-additions. Accretion rates of H atoms and diffusion coefficients of O are then the key factors to be compared in order to determine at what density of the medium oxidation reactions become comparable to H-atom additions.

Table 1

 1 Best fit parameters of the three methods used to model the diffusion constants for O diffusion on five different grain surface analogues.

		Quantum tunnelling	Arrhenius law	Empirical law
		a (Å)	Ea (K)	Ediff [6<Ts<25] (K)	k0/10 -15	α	β
	Porous ASW	0.69±0.10	530±70	170 < Ediff < 600	1.30	1	3
	Non-porous ASW	0.70±0.05	520±60	170 < Ediff < 600	1.21	1	3
	Crystalline water ice	0.69±0.05	500±50	170 < Ediff < 600	1.42	1	3
	Amorphous silicate	0.67±0.10	720±70	290 < Ediff < 740	0.15	0.1	4
	Oxidised HOPG	0.67±0.10	740±60	290 < Ediff < 740	0.1	0.1	4

Table 1

 1

for best fit values of k0, α, and β.

Acknowledgements

The LERMA-LAMAp team in Cergy acknowledges the support of the national PCMI programme founded by CNRS, and the Conseil Regional d'Ile de France through SESAME programmes (contract I-07597R). MM acknowledges financial support by LASSIE, a European FP7 ITN Communitys Seventh Framework Programme under Grant Agreement No. 238258. MM also thanks Prof. Tomellini (Università di Roma Tor Vergata) for fruitful discussions.


Radouan Boukharfane

Pedro José

Martínez Ferrer

Arnaud Mura
email: arnaud.mura@ensma.fr

Vincent Giovangigli

On the role of bulk viscosity in compressible reactive shear layer developments

Keywords:

Despite 150 years of research after the reference work of Stokes, it should be acknowledged that some confusion still remains in the literature regarding the importance of bulk viscosity effects in flows of both academic and practical interests. On the one hand, it can be readily shown that the neglection of bulk viscosity (i.e., κ = 0) is strictly exact for mono-atomic gases. The corresponding bulk viscosity effects are also unlikely to alter the flowfield dynamics provided that the ratio of the shear viscosity µ to the bulk viscosity κ remains sufficiently large. On the other hand, for polyatomic gases, the scattered available experimental and numerical data show that it is certainly not zero and actually often far from negligible. Therefore, since the ratio κ/µ can display significant variations and may reach very large values (it can exceed thirty for dihydrogen), it remains unclear to what extent the neglection of κ holds.

The purpose of the present study is thus to analyze the mechanisms through which bulk viscosity and associated processes may alter a canonical turbulent flow. In this context, we perform direct numerical simulations (DNS) of spatially-developing compressible non-reactive and reactive hydrogen-air shear layers interacting with an oblique shock wave. The corresponding flowfield is of special interest for various reactive high-speed flow applications, e.g., scramjets. The corresponding computations either neglect the influence of bulk viscosity (κ = 0) or take it into consideration by evaluating its value using the EGlib library. The qualitative inspection of the results obtained for two-dimensional cases in either the presence or the absence of bulk viscosity effects shows that the local and instantaneous structure of the mixing layer may be deeply altered when taking bulk viscosity into account. This contrasts with some mean statistical quantities, e.g., the vorticity thickness growth rate, which do not exhibit any significant sensitivity to the bulk viscosity. Enstrophy, Reynolds stress components, and turbulent kinetic energy (TKE) budgets are then evaluated from three-dimensional reactive simulations. Slight modifications are put into evidence on the energy transfer and dissipation contributions. From the obtained results, one may expect that refined large-eddy simulations (LES) may be rather sensitive to the consideration of bulk viscosity, while Reynolds-averaged Navier-Stokes (RANS) simulations, which are based on statistical averages, are not.

Introduction

The bulk (or volume) viscosity κ, which can be related to the second (or dilatational) viscosity coefficient λ, is associated to the vibrational and rotational energy of the molecules. From the macroscopic viewpoint, it characterizes the resistance to dilatation of an infinitesimal bulk element at constant shape [START_REF] Balucani | Dynamics of the Liquid State[END_REF]. It is strictly zero only for dilute monoatomic gases and this theoretical result is often used to discard it, regardless of the nature or internal structure of the fluid as well as the flowfield conditions. However, acoustic absorption measurements performed at room temperature have shown that the ratio of the volume to the shear viscosity κ/µ may be up to thirty for hydrogen at room temperature [START_REF] Cramer | Numerical estimates for the bulk viscosity of ideal gases[END_REF], and recent analyses of reactive multicomponent high-speed flows have confirmed that it is not justified to neglect it, except for the sake of simplicity [START_REF] Billet | Impact of volume viscosity on a shock-hydrogen-bubble interaction[END_REF]. The dilatational viscosity is important in describing sound attenuation in gaseous media, and the absorption of sound energy into the fluid depends itself on the sound frequency, i.e., the rate of fluid expansion and compression. For polyatomic gases, the available measurements of κ, which remains quite seldom due to the complexity of its determination, show that it is certainly not zero and actually far from negligible. It is also noteworthy that theoretical analyses do show that κ/η is at least of the order of unity. Therefore, since the ratio κ/µ can display significant variations and may reach very large values, it is unclear to what extent the Stokes hypothesis (i.e., λ = -2µ/3 or κ = 0) may be used for compressible and turbulent flows of gases featuring a ratio κ/µ greater than unity.

In either an expansion or a contraction of the gas mixture, the work done by the pressure modifies immediately the translational energy of the molecules, while a certain time-lag is needed for the translational and internal energy to re-equilibrate through inelastic collisions [START_REF] Chapman | The Mathematical Theory of Non-uniform Gases[END_REF]. This can be described through a system of two coupled partial differential equations written for the internal and translational temperatures, with a pressure-dilatation term that acts as a source term in the translational temperature budget. The volume (or bulk) viscosity is associated to this relaxation phenomenon and it is evaluated from this internal energy relaxation time-lag. The evaluation of this property for a mixture of polyatomic gases is far from being an easy task since the kinetic theory of gases does not yield an explicit expression for this transport coefficient, but instead linear systems that must be solved [START_REF] Ern | Volume viscosity of dilute polyatomic gas mixtures[END_REF]. The corresponding systems are derived from polynomial expansions of the species' perturbed distribution functions. The bulk viscosity is obtained here using the library EGlib developed by Ern and Giovangigli [START_REF] Ern | Fast and accurate multicomponent transport property evaluation[END_REF][START_REF] Ern | EGlib users manual[END_REF]. It is evaluated as a linear combination of the pure species volume viscosities, which require the evaluation of various collision integrals [START_REF] Ern | Volume viscosity of dilute polyatomic gas mixtures[END_REF].

The impact of bulk viscosity effects has been previously analyzed in several situations including shock-hydrogen bubble interactions [START_REF] Billet | Impact of volume viscosity on a shock-hydrogen-bubble interaction[END_REF], turbulent flames [START_REF] Fru | Impact of volume viscosity on the structure of turbulent premixed flames in the thin reaction zone regime[END_REF], compressible boundary layers [START_REF] Cramer | Effect of large bulk viscosity on large-Reynoldsnumber flows[END_REF], shock-boundary layer interaction [START_REF] Bahmani | Suppression of shock-induced separation in fluids having large bulk viscosities[END_REF], and planar shockwave [START_REF] Chikitkin | Effect of bulk viscosity in supersonic flow past spacecraft[END_REF]. All these studies confirm that the bulk viscosity effects may be significant. The purpose of the present work is to assess its influence in regard to both the instantaneous and statistical features of canonical compressible turbulent multicomponent flows. Using direct numerical simulation (DNS), we investigate the impact of the bulk viscosity coefficient κ on the spatial development of reactive and non-reactive compressible mixing layers interacting with an oblique shock wave. Such a canonical flowfield is typical of the shock-mixing layer interactions that take place in compressible flows of practical interest [START_REF] Huete | Diffusion-flame ignition by shock-wave impingement on a hydrogen/air supersonic mixing layer[END_REF]. For instance, supersonic jets at high nozzle-pressure ratio (NPR) give rise to complex cellular structures, where shocks and expansions waves interact with the turbulent outer shear layer [START_REF] Buttay | Analysis of small-scale scalar mixing processes in highly under-expanded jets[END_REF]. It is also relevant to scramjet intakes and combustors, where shock waves interact with the shear layers issued from the injection systems. On the one hand, it is clear that the occurrence of shock waves in supersonic combustors induces pressure losses that cannot be avoided but, on the other hand, the resulting shock interactions with mixing layers contribute to scalar dissipation (i.e., mixing) rates enhancement [START_REF] Boukharfane | Evolution of scalar and velocity dynamics in planar shock-turbulence interaction[END_REF], and may favor combustion stabilization in high-speed flows.

The present manuscript is organized as follows: the mathematical model is presented in the next section (i.e., Section 2), which also includes a short description of the numerical methods. The details of the computational setup are subsequently provided in Section 3. Section 4 gathers all the results issued from (i) twodimensional numerical simulations of both inert (4.1) and reactive (4.2) cases, and (ii) the three-dimensional case, which is analyzed in 4.3. Finally, some concluding remarks and perspectives for future works are presented in Section 5.

Mathematical description and computational model

In this work, the in-house massively parallel DNS solver CREAMS is used. It solves the unsteady, three-dimensional set of compressible Navier-Stokes equations for multicomponent reactive mixtures [START_REF] Ferrer | A detailed verification procedure for compressible reactive multicomponent Navier-Stokes solvers[END_REF]:

∂ t (ρ) + ∇ • (ρu) = 0, (1a) ∂ t (ρu) + ∇ • (ρu ⊗ u) = ∇ • σ, (1b) ∂ t (ρE t ) + ∇ • (ρuE t ) = ∇ • (σ • u -J ) , (1c) ∂ t (ρY α ) + ∇ • (ρuY α ) = -∇ • (ρV α Y α ) + ρ ωα , (1d)
where t denotes the time, ∇ is the spatial derivative operator, u is the flow velocity, ρ is the density, E t = e+u•u/2 is the total specific energy (obtained as the sum of the internal specific energy, e, and kinetic energy), Y α is the mass fraction of chemical species α (with α ∈ S = {1, . . . , N sp }), V α is the diffusion velocity of species α, J is the heat flux vector and ωα represents the chemical production rate of species α. The integer N sp denotes the number of chemical species.

The above set of conservation equations (1) requires to be completed by constitutive laws. In this respect, the ideal gas mixture equation of state (EoS), P = ρRT /W with R the universal gas constant, is used to relate the pressure P to the temperature T . In this expression, the quantity W denotes the molar weight of the multicomponent mixture, which is obtained as the sum of the molecular mass of each individual species

W -1 = ∑ Nsp α=1 Y α /W α .
Within the framework of the kinetic theory of dilute polyatomic gas mixtures, the molecular diffusion velocity vector V α , α ∈ S , heat flux vector J , and second-order stress tensor σ are expressed as follows:

ρV α Y α = - ∑ β∈S ρY α D α,β ( d β + χ β X β ∇ ( ln T ) ) , (2a)
J = ∑ α∈S ρV α Y α ( h α + RT χ α W α ) -λ T ∇T , (2b)
σ = -PI + τ = -PI + µ (∇u + ∇u ⊺ ) + λ (∇ • u) I , (2c)
where D α,β are the multicomponent diffusion coefficients, d α , α ∈ S , the species diffusion driving forces, χ α the rescaled thermal diffusion ratios, X α the species mole fractions, h α the enthalpy per unit mass of the αth species, and λ T the thermal conductivity. The diffusion driving force d α of the αth species is given by

d α = ∇X α + ( X α -Y α ) ∇ ( ln P ).
The quantity µ denotes the shear viscosity and λ denotes the second (or dilatation) viscosity coefficient. The bulk viscosity coefficient κ appears explicitly in the expression of the viscous stress tensor τ. A relationship between the bulk viscosity κ and viscosity coefficients µ and λ can be deduced from the expression of the total pressure, which can be evaluated as the component of the spherical tensor based on the trace of the total stress tensor σ:

- tr(σ ) 3 = - i=3 ∑ i=1 σ ii 3 = P - ( λ + 2 3 µ ) ∇ • u = P -κ∇ • u (3)
The second term in the right-hand-side of the above expression is the dilational contribution, which defines the bulk viscosity as κ = λ + 2µ/3. As mentioned above, the Stokes' hypothesis, stating that λ = -2µ/3 (and hence κ = 0), is often retained as a simplifying assumption. Many efforts have been devoted to the derivation of relationships between the bulk viscosity and fundamental fluid properties [START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF][START_REF] Prangsma | Ultrasonic determination of the volume viscosity of N 2 , CO, CH 4 , and CD 2 between 77 and 300 K[END_REF]. If we consider a single polyatomic gas with a unique internal energy mode, the internal energy relaxation time τ int can be related to the bulk viscosity [4,18]:

κ = ( P R/c 2 v ) • c int τ int , (4)
when the relaxation time is smaller than fluid characteristic times.

In the above expression, c int denotes the internal heat capacity and c v the specific heat at constant volume. When there are several internal energy modes and/or several species present in the mixture, the above simple expression is replaced by the solution to a linear system [START_REF] Ern | Multicomponent Transport Algorithms[END_REF]. Within the Monchick and Mason approximation [START_REF] Monchick | Transport properties of polar gases[END_REF], neglecting complex collisions characterized by more than one quantum jump, the reduced system is diagonal and yields κ [START_REF] Billet | Impact of volume viscosity on a shock-hydrogen-bubble interaction[END_REF]:

κ = ( P R/c 2 v ) • ∑ k∈P X k c int k τ int k , (5)
where the integer P = 1, . . . , n p is the polyatomic species indexing set. The average relaxation time for internal energy of the kth species τ int k is then expressed as:

c int k /τ int k = ∑ l∈N c l k /τ l k , (6)
where τ l k denotes the average relaxation time of internal energy mode l for the kth species, and N is the internal energy mode indexing set.

The CREAMS solver is coupled with the EGlib library to estimate transport coefficients from the kinetic theory of gases [START_REF] Ern | EGlib: A general-purpose fortran library for multicomponent transport property evaluation[END_REF]. In this library, the optimized subroutines EGSKm are used to evaluate the bulk viscosity. The integer m ∈ 2, 6 associated to the subroutine name refers to retained level of approximation. The higher the value of m, the more expensive the algorithm but also the more accurate the bulk viscosity expression. Following the work of Billet et al. [START_REF] Billet | Impact of volume viscosity on a shock-hydrogen-bubble interaction[END_REF], the value m = 3 is retained for the purpose of the present study. The shear viscosity and diffusion velocities are evaluated with the routines EGFE3 and EGFYV, respectively. EGFLCT3 is used to determine the thermal conductivity λ T and rescaled thermal diffusion ratios χ α . In this respect, some additional computations showed that the extra time needed for m = 4 and m = 5 is approximately 30% in comparison with the one needed for m = 2 and m = 3, while an additional time larger than 250% is required for m = 6 compared to m = 2 and m = 3.

The above system (1) is discretized on a Cartesian grid. A seventh-order accurate WENO scheme is used to approximate inviscid fluxes, while an eighth-order accurate centered difference scheme is retained to approximate viscous and diffusive contributions. Time integration is performed with a third-order accurate TVD Runge-Kutta scheme. The stiffness associated to the wide range of time scales involved in the description of the chemical system is addressed using the Sundials CVODE solver [START_REF] Hindmarsh | Sundials: Suite of nonlinear and differential/algebraic equation solvers[END_REF]. A standard splitting operator technique, similar to the one previously retained in Ref. [START_REF] Ziegler | An adaptive high-order hybrid scheme for compressible, viscous flows with detailed chemistry[END_REF], is used. A detailed verification of the solver may be found in Refs. [START_REF] Ferrer | A detailed verification procedure for compressible reactive multicomponent Navier-Stokes solvers[END_REF][START_REF] Boukharfane | A combined ghost-pointforcing/direct-forcing immersed boundary method (IBM) for compressible flow simulations[END_REF].

Problem statement and computational setup

We study the interaction of an oblique shock with a spatiallydeveloping shear layer. The upper stream corresponds to the fuel inlet, i.e., a mixture containing hydrogen, and the bottom inlet stream to vitiated air. Both two-and three-dimensional computations are performed. Fig. 1 provides a typical sketch of the corresponding computational geometry and Table 1 gathers the values of the main parameters relevant to the present numerical simulation. The flow initialization is similar to the one retained in Ref. [START_REF] Ferrer | A detailed verification procedure for compressible reactive multicomponent Navier-Stokes solvers[END_REF]. Assuming equal free-stream specific heat capacity ratios, the convective Mach number may be evaluated from M c = (U 1 -U 2 )/(a 1 + a 2 ), where a 1 and a 2 denotes the sonic speeds of streams 1 (oxidizer inlet stream) and 2 (fuel inlet stream) respectively. For the present set of computations, it is equal to M c = 0.48. The values of the velocity at the inlets are U 1 = 1634.0 m/s at bottom (oxidizer stream) and U 2 = 973.0 m/s at top (fuel stream).

The mixing layer flow is impinged by an oblique shock wave that is issued from the oxidizer inlet stream (1) at the bottom boundary. The oblique shock wave angle is β = 33 • , see Fig. 1. The geometrical parameters relevant to the present set of numerical simulations are provided in Table 2. The quantities L x 1 , L x 2 , and L x 3 denote the computational domain lengths in each direction normalized by the initial vorticity thickness δ ω,0 , while

N x 1 , N x 2 ,
and N x 3 are the corresponding numbers of grid points. In the two-dimensional computations, only the x 1 -and x 2 -directions are considered.

The flow is initialized with a hyperbolic tangent profile for the streamwise velocity component, while the other velocity components are set to zero. Species mass fractions and density are also set according to the following general expression:

ϕ(x 1 , x 2 , x 3 ) = ϕ 1 + ϕ 2 2 + ϕ 1 -ϕ 2 2 tanh ( 2x 2
δ ω,0

) ,

where ϕ denotes any of the flow variables mentioned above (i.e., species mass fraction or streamwise velocity component). The value of the Reynolds number Re δω , based on the initial vorticity thickness and inlet velocity difference ∆U = U 1 -U 2 is Re δω = 640.

One of the fundamental statistical quantities that characterizes the mixing layer development is its normalized growth rate [START_REF] Ramshaw | Simple model for mixing at accelerated fluid interfaces with shear and compression[END_REF].

Although the definition of this growth rate is not unique, one standard expression relies on the vorticity thickness definition [START_REF] Ferrer | Compressibility and heat release effects in high-speed reactive mixing layers. Growth rates and turbulence characteristics[END_REF][START_REF] Pantano | Mixing of a conserved scalar in a turbulent reacting shear layer[END_REF]:

δ ω (x 1 ) = U 1 -U 2 ∂ũ 1 /∂x 2 | max . ( 8
)
Dirichlet boundary conditions are applied at the two supersonic inlets, perfectly non-reflecting boundary conditions are set at the outflow, and periodic boundary conditions are settled along the x 3direction. A slip boundary condition is imposed at the top, while the bottom boundary condition is set by using Rankine-Hugoniot relations, generalized for a multicomponent mixture [START_REF] Mitchell | General-purpose computer code for predicting chemical-kinetic behavior behind incident and reflected shocks[END_REF]. In order to trigger flow transition, a slight white noise fluctuation is superimposed to the transverse velocity component along the line slight overestimate compared to the case where the effects of κ are considered. The region that extends from x 1 /δ ω,0 = 150.0 until the interaction with the reflected shock is characterized by a significant change of behavior and the values obtained with κ ̸ = 0 are larger than those obtained with κ = 0. This region is characterized by strong pressure wave reflection from the upper limit of the computational domain. After the interaction with the reflected shock wave, the maxima of the TKE obtained with or without taking into account the bulk viscosity effects tend to become similar as the end of the computational domain is approached. It can be concluded that, in the absence of the second shock wave interaction and associated parasitic pressure waves issued from the top of the computational domain, the TKE levels would be slightly underestimated if the effects of κ were not taken into account. In an attempt to better understand the behavior of the TKE, the analysis of the main terms involved in its transport equation is now carried out. The transport equation for the turbulent kinetic energy K is given by

(x 1 , x 2 ) = (4δ ω,0 , 0).
∂ t (ρK) + ∇ • (ρũK) = P + ε + T + Π + Σ (12)
In this equation, P is the production term, ε is the dissipation term, T denotes the turbulent transport term, Π is the pressure-strain term, and finally Σ the mass flux term. The budget ( 12) is deduced from the transport equation of the Reynolds tensor components: with

∂(ρR ij ) ∂t + ∂(ρũ k R ij ) ∂x k = P ij + ε ij + T ij + Π ij + Σ ij , (13)
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ P ij = -ρ ( R ik ∂ũ j ∂x k + R jk ∂ũ i ∂x k ) , (a)
ε ij = -τ ′ ik ∂ u ′′ j ∂x k -τ ′ jk ∂ u ′′ i ∂x k , (b)
T ij = - ∂ ∂x k ( ρu ′′ i u ′′ j u ′′ k + P ′ u ′′ i δ jk + P ′ u ′′ j δ ik -τ ′ jk u ′′ i -τ ′ ik u ′′ j ) , (c)
Π ij = P ′ ∂ u ′′ i ∂x j + P ′ ∂ u ′′ j ∂x i , (d)
Σ ij = ( u ′′ i ∂ τ jk ∂x k + u ′′ j ∂ τ ik ∂x k ) - ( u ′′ i ∂ P ∂x j + u ′′ j ∂ P ∂x i ) (e) (14)
In the above set of equations, it has been chosen to split the pressure into a mean and a fluctuating part but it should be recognized that there exist other ways to handle pressure terms in the second-order moment transport equations. Here, we are following the same procedure as the one previously considered by Pantano et Sarkar [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF]. However, it is acknowledged that, for other conditions related, for instance, to transport modeling in turbulent premixed flames in the flamelet regime, keeping the instantaneous pressure can be a better choice since local flamelet relationships may provide relevant closures for the corresponding effects, see for instance Bray et al. [START_REF] Bray | Premixed flames in stagnating turbulence. Part IV: a new theory for Reynolds stresses and Reynolds fluxes applied to impinging flows[END_REF] or Robin et al. [START_REF] V. Robin | A second-order model for turbulent reactive flows with variable equivalence ratio[END_REF].

The analysis of the main terms involved in the TKE transport equation is carried out at two distinct locations to infer the impact of the volume viscosity. Figure 20 shows that the most important contributions are associated to the production and dissipation terms. Their amplitude is found to be slightly smaller when κ is not taken into account. The turbulent transport term is positive at the periphery of the mixing layer while it tends to be negative within the mixing layer. This quantity, which is larger in the case featuring κ ̸ = 0, removes energy from regions characterized by large fluctuations levels to transfer it in regions characterized by lower levels of TKE. Figure 20 also shows that the contributions due to pressure-strain and mass flux terms remain negligible compared to the others, for both cases.

Figure 21 reports the distribution of the Reynolds stress components together with the variance of the passive scalar and the scalar flux components for both cases κ ̸ = 0 and κ = 0. The three streamwise positions under consideration are representative of the variations observed on the TKE profile reported in Fig. 19. The profiles of the Reynolds stress components show that the maxima of its diagonal components follow the trends reported in Fig. 19. Figure 21(f), which displays the longitudinal evolution of the scalar flux component ũ′′ 1 ξ ′′ /(u 1,RMS ξ RMS ), reveals that the maximum value of the correlation between the longitudinal velocity fluctuation and the scalar fluctuation is slightly underestimated when the effects of bulk viscosity are not considered.

Figure 22 reports the variance of the mass fractions of chemical species present in the mixture. The hydrogen, which is characterized by the highest ratio κ/µ is the one that displays the largest differences (up to approximately ten percent) between the two cases, i.e., κ = 0 and κ ̸ = 0. The differences observed at the three locations concern both the shape and maximum levels, which depend on the species under consideration. Indeed, it is found that the distribution of the profiles for all chemical species is slightly wider -indicating that the fluid is incorporated more markedlywhen the effects of bulk viscosity are not taken into account, which leads to a reduction of fluctuations around the averaged value. A similar effect is observed when the convective Mach number values are increased [START_REF] Ferrer | Compressibility and heat release effects in high-speed reactive mixing layers. Growth rates and turbulence characteristics[END_REF][START_REF] Mahle | Direct and Large-eddy Simulation of Inert and Reacting Compressible Turbulent Shear Layers[END_REF].

Summary and conclusions

In the present manuscript, two-and three-dimensional numerical simulations of spatially-developing compressible mixing layers impacted by an oblique shock wave are conducted for a convective Mach number M c = 0.48. The emphasis is placed on the possible influence of the bulk viscosity on the mixing processes. Thus, a mixture of hydrogen and air is considered in conditions that are representative of experimental benchmarks relevant to high-speed flow combustion. In a first step of the analysis, twodimensional computations of inert and reactive mixing layers are performed. A significant impact of the bulk viscosity is observed on the instantaneous flowfields while averaged quantities do not exhibit any remarkable modification. It is also worth noting that the reactive cases only display slight differences with respect to inert cases: this is especially true for the longitudinal evolutions of the vorticity thickness and turbulent kinetic energy. Threedimensional simulations of inert mixing layers are subsequently conducted. The influence of the bulk viscosity is more visible in these three-dimensional cases: it tends to reduce the mixing layer growth rate compared to the case where it is not taken into account. The comparison is also performed in terms of higherorder statistical moments. This last part of the analysis shows that the bulk viscosity effects tend to amplify the velocity gradients at the boundaries of the mixing layer, and consequently favor the return to equilibrium. From the above synthesis of the obtained results, one may expect that refined large-eddy simulations (LES) may be rather sensitive to the consideration of bulk viscosity, while Reynolds-averaged Navier-Stokes (RANS) simulations, which are based on statistical averages, are not. One possible perspective of the present work concerns the filtering of the present dataset, which may provide further insights so as to assess this conclusion. Finally, from the present set of results, it is recommended to take the bulk viscosity effects into account especially when highlyresolved large-eddy simulations (LES) are considered.

Fig. 21 .

 21 Fig. 21. Profiles of the Reynolds stress tensor components normalized by ∆U together with the mixture fraction variance ξ ′′ ξ ′′ and longitudinal component of the scalar flux ũ′′ 1 ξ ′′ at three abscissae, same symbols as those retained in Fig.17.

17

 17

  .

Fig. 22 .

 22 Fig. 22. Profiles of the variances of chemical species mass fractions at three abscissae, same symbols as those retained in Fig. 17.

  Fig. 22. Profiles of the variances of chemical species mass fractions at three abscissae, same symbols as those retained in Fig. 17.

Table 1

 1 Parameters of the shock-mixing layer interaction case.

			Fuel		Oxidizer		Bottom
	P (Pa)		94 232.25		94 232.25		129 951.0
	T (K)		545.0		1475.0		1582.6
	u 1 (m/s)		973.0		1634.0		1526.3
	u 2 (m/s)		0.0		0.0		165.7
	u 3 (m/s)		0.0		0.0		0.0
	Mach (-)		1.6		2.12		1.93
	Y H 2 (-)		0.05		0.0		0.0
	Y O 2 (-)		0.0		0.278		0.278
	Y H 2 O (-)		0.0		0.17		0.17
	Y H (-)		0.0		5.60 • 10 -7		5.60 • 10 -7
	Y O (-)		0.0		1.55 • 10 -4		1.55 • 10 -4
	Y OH (-)		0.0		1.83 • 10 -3		1.83 • 10 -3
	Y HO 2 (-)		0.0		5.10 × 10 -6		5.10 × 10 -6
	Y H 2 O 2 (-)		0.0		2.50 × 10 -7		2.50 × 10 -7
	Y N 2 (-)		0.95		0.55		0.55
	Table 2
	Computational mesh description.
	L x 1	L x 2	L x 3	N x 1	N x 2	N x 3	δ ω,0 (m)
	280	130	15	1640	750	180	1.44 × 10 -4

  The value of the CFL number is set to 0.75.

Reactive flow simulations are conducted with the detailed mechanism of O'Conaire et al.

[START_REF] Conaire | A comprehensive modeling study of hydrogen oxidation[END_REF]

. It consists of nine chemical species (H 2 , O 2 , H 2 O, H, O, OH, HO 2 , H 2 O 2 , and N 2 ) and 21 elementary reaction steps. The concentrations of these species at the inlet have been determined from equilibrium conditions so as to reach favorable self-ignition conditions within the extension of the computational domain.

Throughout this manuscript, the Reynolds and Favre averages of any quantity ϕ are denoted by ϕ and φ, while the corresponding

Acknowledgments

The computations were performed using the High Performance Computing resources from the mésocentre de calcul poitevin and from genci under allocations x20142a0912 and x20142b7251. The first author benefited from interesting discussions with Aimad Erraiy. This work has been presented at the tenth International Conference on Computational Fluid Dynamics (Barcelona, July 2018).