Mattia Farronato

Johann Lüder

Danilo Longo

Hervé Cruguel

Marcel Bouvet

Barbara Brena

Nadine Witkowski
email: nadine.witkowski@sorbonne-universite.fr

High Tolerance of Double-Decker Phthalocyanine Towards Molecular Oxygen

Because organic electronics suffer from degradation inducing oxidation processes, oxygen tolerant organic molecules could solve this issue and be integrated to improve the stability of devices during operation. In this work, we investigate how lutetium double-decker phthalocyanine (LuPc 2 ) reacts towards molecular oxygen and we report microscopic details of its interaction with LuPc 2 film by combining X-ray Photoemission Spectroscopy, Near Edge X-ray Absorption Fine Structure Spectroscopy and Density Functional Theory. Surprisingly, LuPc 2 molecules are found to weakly physisorb below 120 K and appear rather inert to molecular oxygen at more elevated temperatures. We are able to draw a microscopic picture at low temperature in which oxygen molecules stick on top of the pyrrolic carbon of LuPc 2 . Our work sheds light on a class of semiconducting molecules, namely double-decker phthalocyanines, that present a high tolerance towards molecular oxygen, opening promising perspectives for the design of stable materials to be applied in the next generation of organic based electronic devices operating under ambient conditions.

Introduction

Recently, devices based on organic electronics have attracted a lot of attention both on the market and in the research laboratories. [START_REF] Kaltenbrunner | Ultrathin and Lightweight Organic Solar Cells with High Flexibility[END_REF][START_REF] Jean | In Situ Vapor-Deposited Parylene Substrates for Ultra-Thin, Lightweight Organic Solar Cells[END_REF][START_REF] Nakanotani | High-Efficiency Organic Light-Emitting Diodes With Fluorescent Emitters[END_REF] Their advantages versus the inorganic silicon based devices are multiples: they are cheap, flexible, and they offer almost infinite tayloring possibilities making these devices ideal for large scale implementation in urban areas.

Among the organic materials used both for research and for devices application, one of the most interesting is represented by phthalocyanines (Pcs) and their derivatives. They are thermally stable, and present interesting electronic and optical properties. For instance, copper phthalocyanine, known as p-type semiconductor, has been extensively studied in organic light-emitting diode (OLED) or in organic photovoltaic cells. [START_REF] Bof Bufon | Determination of the Charge Transport Mechanisms in Ultrathin Copper Phthalocyanine Vertical Heterojunctions[END_REF][START_REF] Zhao | Emission Spectral Stability Modification of Tandem Organic Light-Emitting Diodes through Controlling Charge-Carrier Migration and Outcoupling Efficiency at Intermediate/Emitting Unit Interface[END_REF][START_REF] De Arco | Continuous, Highly Flexible, and Transparent Graphene Films by Chemical Vapor Deposition for Organic Photovoltaics[END_REF] One of the major drawbacks that hinders phthalocyanine based materials from taking over the market, is their reduced lifetime originating from photobleaching, exposure to elevated temperature (in the range of 70-100 • C) or exposure to reactive gases or radiation. In particular the exposure towards atmospheric gases can cause many different problems as reactivity, [START_REF] Seemann | Reversible and Irreversible Degradation of Organic Solar Cell Performance by Oxygen[END_REF][START_REF] Volonakis | Impurity-related Degradation in a Prototype Organic Photovoltaic Material: A First-principles Study[END_REF] and even polymerization of the organic material in presence of oxygen at operating temperatures. [START_REF] Lessmann | Aging of Flat Heterojunction Zinc Phthalocyanine/Fullerene C-60 Organic Solar Cells[END_REF] In the case of phthalocyanine few studies have been published on its reactivity towards gases as O 2 on FePc, [START_REF] Sedona | Tuning the Catalytic Activity of Ag(110)-supported Fe Phthalocyanine In the Oxygen Reduction Reaction[END_REF] NO 2 on CuPc [START_REF] Lozzi | Photoemission and Theoretical Investigations on NO 2 Doping of Copper Phthalocyanine Thin Films[END_REF] and NO on CoPc, [START_REF] Flechtner | NO-Induced Reversible Switching of the Electronic Interaction between a Porphyrin-Coordinated Cobalt Ion and a Silver Surfacee[END_REF] all evidencing that the reaction takes place on the central metallic cation, with the macrocycle only able to increase or decrease this reactivity by modifing the molecular orbitals. [START_REF] Flyagina | DFT Study of the Oxygen Reduction Reaction on Iron, Cobalt and Manganese Macrocycle Active Sites[END_REF] A study on the reactivity of metal free naphthalocyanine towards oxygen was presented by Ottaviano et al. [START_REF] Ottaviano | Interaction of Naphthalocyanine with Oxygen and with Si(111) 7x7: An in-situ X-ray Photoelectron Spectroscopy Study[END_REF] who found that the reaction takes place on the azabridge nitrogen atoms which are saturated by H in metal-free phthalocyanine and derivatives; the two hydrogens are dislodged by the oxygen molecules binding to the phthalocyanine. Also studies on the interaction between single decker phthalocyanine and transition metal oxide surface were performed. For example, the investigation from Glaser et al. [START_REF] Glaser | Electronic Structure at Transition Metal Phthalocyanine-Transition Metal Oxide Interfaces: Cobalt Phthalocyanine on Epitaxial MnO Films[END_REF] showed a charge transfer from the MnO substrate to the CoPc overlayer, and a subsequent shift of the Co 2p towards lower Binding Energy (BE) due to the rehybridization of the Co and MnO orbitals, while no modifications were seen on the N and C 1s spectra. Liu and coworkers [START_REF] Liu | A Direct Fe-O Coordination at the FePc/MoOx Interface Investigated by XPS and NEXAFS Spectroscopies[END_REF] showed instead that deposition of FePc over MoO x led to oxidation of the metallic centre. All these works evidenced that the crucial location for the molecular reactivity lies on the central cation (or the N atoms saturated by hydrogens in the case of metal free Pcs). The central metallic cation has a dominant role in the charge redistribution during the charge transfer process [START_REF] Chen | Charge Transport in Stacking Metal and Metal-Free Phthalocyanine Iodides. Effects of Packing, Dopants, External Electric Field, Central Metals, Core Modification, and Substitutions[END_REF] and therefore influences the charge carriers mobility.

To improve the stability of phthalocyanine based organic devices, a class of molecules that would naturally prevent the above mentioned reaction mechanism on the cation, is represented by double-decker phthalocyanine. The geometry of LuPc 2 is shown in Figure 1, the coordination sphere of the central Lu 3+ cation is saturated by the eight pyrrolic nitrogen atoms. This saturation should prevent the central ion from reacting, making of this class of phthalocyanine a promising alternative to the widely used single decker ones. Only few studies were published on the reactivity of double-decker phthalocyanines [START_REF] Passard | Gaseous Oxidation and Compensating Reduction of Lutetium Bis-phthalocyanine and Lutetium Phthalo-naphthalocyanine Films[END_REF][START_REF] Yamana | Effects of Heat Treatment on Oxidative Gas Adsorption for Lead Naphthalocyanine Thin Films[END_REF][START_REF] De Saja | Sensors Based on Double-Decker Rare Earth Phthalocyanines[END_REF] and mostly from a device point of view. For example the study of the conductivity before and after the exposition to highly oxidizing gases showed at first an increase of the conductivity due to the partial oxidation of the film [START_REF] Passard | Gaseous Oxidation and Compensating Reduction of Lutetium Bis-phthalocyanine and Lutetium Phthalo-naphthalocyanine Films[END_REF] which produced an increased hopping possibility caused by the increased number of accepting molecules. A further oxidation led then to a decrease in conductivity, due to the increased fraction of oxidized molecules, and therefore a reduced number of accepting sites.

So far very little is known about the microcopic molecular processes leading to such reactivity; in the present study, we propose to adress this issue by investigating the reactivity of LuPc 2 thin films towards molecular oxygen at low temperature. Surface spectroscopies as X-ray Photoelectron Spectroscopy (XPS) and Near Edge X-ray Absorption Fine Structure Spectroscopy (NEXAFS) are the perfect tools to reveal the intimate mechanisms, as they guarantee great surface sensitivity (i.e. only the top few layers are probed), chemical selectivity and site specificity as well as information about the charge redistribution after the reaction. By combining experimental techniques and advanced Density Functional Theory (DFT) simulations, we were able to define precisely the location and energy adsorption of O 2 molecules on LuPc 2 .

Experimental and/or Theoretical Methods

The experiment was performed at the beamline D1011 in MAX IV Laboratory in Lund at the front station. [START_REF] Nyholm | A Soft-x-ray Monochromator for the Max Synchrotron Radiation Facility[END_REF] The D1011 beamline is a bending magnet beamline which allows a flux of 10 10 -10 11 linearly polarized (85%) photons•s -1 . The endstation is composed of two chambers: a preparation chamber and an analysis chamber. In both chambers the base pressure was in the low 10 -10 mbar range. The preparation chamber was equipped with an ion gun, an e-beam heating to prepare the sample, and Low Energy Electron Diffraction (LEED). The analysis chamber was equipped with electron spectrometer (VG Scienta SES200 concentric hemispheric analyser), electron detector with partial and total electron yield for NEXAFS measurement. Au(111) single crystal was purchased by Surface Preparation Laboratory (SPL). It was cleaned by repeated cycles of sputtering and annealing until no N, C signal was seen in the NEXAFS or in the XPS spectrum. The last annealing was performed at 873 K and slowly cooled to obtain the ( √ 3 × 22) herringbone reconstruction, [START_REF] Barth | Scanning Tunnelling Microscopy Observations on the Reconstructed Au(111) Surface: Atomic Structure, Long-range Superstructure, Rotational Domains, and Surface Defects[END_REF] which was checked by LEED.

LuPc 2 molecules were synthesized as reported in literature. [START_REF] Bouvet | Electrical-properties of Rare-earth Bisphthalocyanine and Bisnaphthalocyanine Complexes[END_REF][START_REF] Clarisse | Synthesis and Characterization of some Lanthanide Phthalocyanines[END_REF] Before the evaporation, the crucible and the molecules were carefully outgassed at 250 • C for several hours to remove the residual H 2 Pc, until the base pressure was recovered. The organic layer was prepared by thermal evaporation of the molecules from a homemade Ta pocket heated by direct current.

Evaporator temperature was kept around 573 K and was measured with a thermocouple coiled around the crucible. During the deposition, the sample was kept at room temperature and the total pressure was kept in the 10 -8 mbar regime for about 30 minutes. The layer was thick enough so that no Au 4f peak was seen with the XPS ensuring a covering of at least 3 nm. The XPS data were collected with pass energy of 100 eV, and the overall resolution was about 0.1 eV. The spectra were recorded with Photon Energy (PE) of 500 eV and 450 eV for N 1s and C 1s spectra respectively. Electrons were collected in normal emission and the BE scale was calibrated on the main C 1s line, which was found at the same BE as previously published, [START_REF] Bidermane | Experimental and Theoretical Study of Electronic Structure of Lutetium Bi-phthalocyanine[END_REF][START_REF] Toader | Initial Growth of Lutetium(III) Bis-phthalocyanine on Ag(111) Surface[END_REF] or on Lu 4f which both did not show any shift after oxygen exposure. A mixed Shirley and linear background was subtracted from XPS data, and the data were analysed with IgorPro software and the SPANCF procedure. [START_REF] Kukk | Spectral Analysis by Curve Fitting Macro Package SPANCF[END_REF] Voigt curves with fixed Lorentzian FWMH of 0.1 eV and 0.132 eV were used for C 1s and N 1s respectively while the width of the Gaussian curve was left free to vary in the fitting procedure. The intensities have been normalized to unity on the main feature height.

NEXAFS spectra were taken in partial electron yield mode using a channel plate with a retarding voltage of 150 eV and an energy step of 50 meV for the π resonances and 100 meV

for the σ resonances. The resolution E/∆E exceeded 6000 for the N and O K-edge thresholds, [START_REF] Vinogradov | Low-lying Unoccupied Electronic States in 3d Transition-metal Fluorides Probed by NEXAFS at the F 1s Threshold[END_REF] equivalent to a resolution of around 50 meV for N K-edge, and 80 meV for O K-edge.

Two polarizations were investigated: in-plane polarization (IPP) for which the polarization E vector of the radiation was nearly parallel to the surface (misalignment of about 4 • ) and out-of-plane polarization (OPP). The photon energy was calibrated at the beginning and at the end of the first N K-edge spectrum with C 1s first and second order, and successive spectra were calibrated on the first π * resonance. The spectra were first divided by the references on clean Au substrate taken at the beginning of the experiment, and then set to 0 below threshold. [START_REF] Stöhr | NEXAFS Spectroscopy[END_REF] Oxygen with purity >99.995%, purchased by Air Liquid, was introduced through a leak valve in the preparation chamber. During the dose, the pressure was kept in the low 10 -5 mbar.

Density Functional Theory (DFT) was employed to gain insights into the interaction between the O 2 molecule and LuPc 2 . The computations, i.e. structure optimization and XPS simulations, were performed on single molecules and for the two molecules, i.e. O 2 adsorbed on LuPc 2 , interacting with each other. For this purpose, the VASP package [START_REF] Kresse | Ab initio Molecular Dynamics for Liquid Metals[END_REF][START_REF] Kresse | Efficient Iterative Schemes for ab initio Total-energy Calculations using a Plane-wave Basis Set[END_REF][START_REF] Kresse | Efficiency of ab-initio Total Energy Calculations for Metals and Semiconductors using a Plane-wave Basis Set[END_REF] was used with the implemented projector augmented wave method (PAW), [START_REF] Kresse | From Ultrasoft Pseudopotentials to the Projector Augmentedwave Method[END_REF] the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF][START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] and the Tkatchenko-Scheffler (TS) dispersion correction, [START_REF] Tkatchenko | Accurate Molecular van der Waals Interactions from Ground-State Electron Density and Free-atom Reference Data[END_REF] which is required in DFT to accurately describe weakly interacting systems, e.g. such that are physisorbed. To minimize cross talk between periodic images, the simulations cells contained 10 Å of vacuum in each direction which resulted in box sizes of at least, e.g. 25×25×16 Å3 for a single LuPc 2 . The interaction of LuPc 2 with one O 2 molecule was investigated in a two steps optimization approach and with several initial geometries. The initial geometries were pseudo-randomly generated with distances between O 2 and LuPc 2 ranging from 2.4 to 3.5 Å. Due to symmetry of the LuPc 2 structure, the initial O 2 positions were located at only one of the isoindene rings. In total, 90 initial geometries were generated with O 2 at different sites and with different orientations. Then, preoptimizations were carried out for these structures with an energy cutoff of 300 eV and a low force convergence of 0.03 eV/ Å. In step two, the six most stable structures were selected for continuation of their geometry optimization being carried out with an energy cutoff of 500 eV until all forces on the atoms were less than 0.02 eV/ Å. The isolated LuPc 2 and the isolated O 2 molecule were also optimized with the latter convergence criteria. For all calculations, the Brillouin zone was sampled at the Γ point only.

E ad = E LuP c 2 O 2 -(E O 2 + E LuP c 2 ).
N 1s and C 1s XPS simulations were performed with the modified PAW by Köhler [START_REF] Köhler | Density Functional Study of CO on Rh(111)[END_REF] in VASP, which allows to relax the otherwise frozen core states. We chose to compute energy shifts in XPS via a ground state approximation for which we estimate the 1s ionization energy of a state being equal to its computed eigenvalue. In this way, computed 1s XPS signal is, however, a relative measure and only comparable to the same element within the same computational setup while absolute excitation energies cannot be obtained. Nonetheless, relative shifts in photoemission spectra, for example in C 1s, can provide information about chemical environment of atoms and changes of it which makes this approach suitable to study the interaction between LuPc 2 with O 2 and it allows to interpret possible changes in experimental data. We distinguish between two different shifts: i) a shift (∆) due to the chemical bonds within a molecule, e.g., the well-known shift of about 1.3 eV between pyrrolic and benzene C 1s signal (which is well-known for many phthalocyanines 25 from experimental measurements) and ii) a shift (δ) due to adsorption of other molecules like O 2 which is measured at the same atom before and after O 2 adsorption. We also introduce the XP range, that corresponds to the maximal difference computed for the energy position of C 1s or N 1s core levels.

Results and discussion

XPS results

Figure 2 presents the experimental XPS spectra of C 1s and N 1s (green dots for spectra (C1) and (N1)) for the as deposited LuPc 2 film, after a dose 3600 L O 2 (spectra (C2) and (N2))

and 90000 L O 2 (spectra (C3) and (N3)) at 90 K, respectively. The N 1s spectrum for the as deposited sample shows a main peak at 398 eV with a single shoulder at higher BE around 399 eV. The C 1s spectrum instead shows a main peak at 284.5 eV with two shoulders at higher BE. In both spectra, a shoulder appears at higher BE next to the main peak after the first oxygen dose of 3600 L (1L = 1 × 10 -6 Torr×1s) (spectra (C2) and (N2)), that increases further after the second dose of 90000 L (spectra (C3) and (N3)). The measurements at 90 K were followed by a warming up to RT that revealed a complete restoration of the XPS spectra for both carbon and nitrogen species (see Supporting Information) and also NEXAFS at the nitrogen K-edge. The original NEXAFS was already retrieved for a temperature around 120 K (not shown), evidencing the weak bonding between the molecular oxygen and the LuPc 2 molecules which affects both C 1s and N 1s. To go further in the analysis of the XPS, spectral components for the N 1s and C 1s peaks have been extracted and are displayed in

C 1s hν=450eV

N 1s hν=500eV

Binding Energy (eV) For the N 1s level in the pristine LuPc 2 film (Figure 2 spectrum (N1)), the main peak is composed of both the azabridge and the pyrrolic nitrogen atoms contributions, which are too close (30 meV) to be resolved, [START_REF] Bidermane | When the Grafting of Double Decker Phthalocyanines on Si(100)-2 x 1 Partly Affects the Molecular Electronic Structure[END_REF] and their associated shake-up at 399. For the C 1s level in the pristine film (Figure 2 spectrum (C1)) the main peak is mostly due to contributions from benzenic carbon found at 284.25 eV. A slight asymmetry is used to take into account for the contribution arising from the two slightly different carbon atoms:

the four linked directly to pyrrolic carbons and the four bound to other benzenic carbon as discussed in ref. [START_REF] Brena | Equivalent Core-hole Time-dependent Density Functional Theory Calculations of Carbon 1s Shake-up States of Phthalocyanine[END_REF] The shoulder at higher BE is composed by contributions from pyrrolic carbons at 285.40 eV and shake-up of the benzenic carbon 39 at 286.30 eV. The last shoulder

DFT results

NEXAFS results

In Figure 4 Bidermane et al. [START_REF] Bidermane | Experimental and Theoretical Study of Electronic Structure of Lutetium Bi-phthalocyanine[END_REF] The NEXAFS spectra show a clear dichroism, with the π * (σ * ) resonances dominating in the out of plane (inplane) polarization for the pristine LuPc 2 film.

As the NEXAFS cross section is strongly dependent on the mutual orientation of the light polarization and the final orbital, [START_REF] Stöhr | NEXAFS Spectroscopy[END_REF] this is a clear indication that the molecules lie flat on the substrate for this thickness (>3 nm) . The dichroism is conserved after the 90000 L O 2 dose, confirming that the molecules still lie flat after the dose and the almost unchanged intensity ratio between the σ * and π * resonances before and after the deposition is due to a limited rehybridization of the molecular orbitals after interaction with O 2 . After the oxygen dose, a new peak is visible at 399 eV in the IPP geometry and also the intensity of the π resonance at 400 eV, corresponding to the pyrrolic atom, becomes stronger, as clearly seen in Figure 4 (B). The origin of the new feature will be discussed in the next section. Note that the possibility of a reaction with condensed water, atomic oxygen or O 3 (created by the x-ray irradiation) was ruled out thanks to NEXAFS measurement on the oxygen K-edge (see Supporting Information).

From the NEXAFS data, it is possible to evidence that the LuPc 2 molecules lie flat on the surface before and after the oxygen dose, and that a perturbation of the electronic structure is induced by the oxygen adsorption at low temperature.

Discussion

The first adsorption mechanism evidenced from the experiment, relies in the weakness of the LuPc 2 -O 2 bond, shown by the low desorption temperature of around 120 K and the complete restoration of the photoemission spectra after desorption. This behavior points towards a physisorption rather than a chemisorption of the oxygen molecule on the LuPc 2 film. This hypothesis is corroborated by the analysis of the binding energy of the XPS additional peaks.

No literature could be found on this specific system, but the N+O feature at 399.65 eV in the N 1s XPS spectrum can not be related to N atoms directly bound to oxygen which are usually found for BE around 401.3 eV. [START_REF] Bouvet | Detection and Titration of Ozone using Metallophthalocyanine Based Field Effect Transistors[END_REF][START_REF] Dash | An XPS Investigation on a Series of Schiff-base Dioxime Ligands and Cobalt Complexes[END_REF][START_REF] Hueso | XPS Investigation of the Reaction of Carbon with NO, O-2, N-2 and H2O Plasmas[END_REF] The DFT results give further evidence of weak interaction, i.e. physisorption, between O 2 and LuPc 2 with adsorption energies of 270 meV or less and adsorption distances of at least 3.1 Å.

The preferential adsorption sites are deduced from the combined analysis of XPS core levels and DFT simulations. There are no studies reporting on similar systems but analogies can be drawn based on N 1s and C 1s binding energies. Whereas BE of about 401 eV is found for N atoms chemically bound to oxygen, the BE falls around 399 eV [START_REF] Clegg | Hydrogen-bonding, Self-assembled Monolayers: Ordered Molecular Films for Study of Through-peptide Electron Transfer[END_REF][START_REF] Baumgarten | Synthesis and Characterization of a New Platinum Supported Catalyst Based on Poly-(acrylamideco-3-(acryloylamino)propyltrimethylammoniumchloride) as Carrier[END_REF][START_REF] Lee | X-ray Photoelectron Spectra and Electronic Structure of some Diamine Compounds[END_REF] for N atoms in the vicinity of carbonyl groups, suggesting such kind of proximity in our system. In the case of the C 1s spectrum the additional peak can be related to a C-O bond in polymers, in small molecules or on graphite [START_REF] Bouvet | Detection and Titration of Ozone using Metallophthalocyanine Based Field Effect Transistors[END_REF][START_REF] Hueso | XPS Investigation of the Reaction of Carbon with NO, O-2, N-2 and H2O Plasmas[END_REF][START_REF] Grunze | Adhesion of Vapor Phase Deposited Ultra-thin Polyimide Films on Polycrystalline Silver[END_REF][START_REF] Marchant | Degradation of a Poly(ether Urethane Urea) Elastomer: Infra-red and XPS Studies[END_REF][START_REF] Steiner | Ultrathin Layers of Low and High-molecular-weight Imides on Gold and Copper[END_REF][START_REF] Gervais | Surface Analysis of Lipopeptides using X-ray Photoelectron Spectroscopy: I. Lipopeptides with Polysarcosine Peptidic Chains[END_REF][START_REF] Sundberg | On the XPS Spectra of some Carbonyl Complexes and the Effective Charge of the Metal and Ligand Atoms[END_REF] respectively. The tilt of molecules in the stacking with respect to the surface could also play a role in the reaction making the configuration (f) accessible to oxygen molecules for an edge-on orientation of the molecules. Unfortunately, this configuration is not expected to lead to an additional feature in the XPS spectrum and more appropriate measurements such as transport would be more relevant to evidence the influence of the stacking on the reactivity.

Therefore, it can be concluded that at least one of these configuration is present in the experiment and that O 2 is adsorbed in the vicinity of the pyrrole ring.

A detailed understanding of the electronic structure derived from the XPS and NEXAFS measurements can be drawn. Looking more carefully at the zoom of the IPP NEXAFS in the Figure 4 for the spectra collected before and after the exposition to oxygen. Instead, as the NEXAFS peaks intensity is strongly dependent on the occupancy of the final state orbitals and their energy density, this is probably an evidence of charge withdrawal 51 from the molecule to the oxygen. The effect is clearly seen in the IPP but not in the OPP (see Figure 4 The oxygen adsorption lifts the degeneracy between the 8 equal pyrrolic nitrogens, as shown by the additional peak in the XPS spectrum, which is relative to a chemically inequivalent N atom, giving rise to an additional transition in the NEXAFS spectrum. This effect was already shown by Solomon and coworkers on a similar system. [START_REF] Solomon | Orientation and Absolute Coverage of Benzene, Aniline, and Phenol on Ag(110) Determined by NEXAFS and XPS[END_REF] The lower energy (PE) of the additional peak in the NEXAFS spectrum with respect to the energy (BE) of the N+O component in the XPS spectrum might be due to the better screening of the final state in the NEXAFS process which retains the electron around the molecule, allowing it to still screen the core hole. [START_REF] Stöhr | NEXAFS Spectroscopy[END_REF] Surprisingly, the physisorption of molecular oxygen induces a reversible charge transfer from the nitrogen atom of the pyrrolic group to the oxygen. Such intriging behavior would require further calculations to fully understand the underlying mechanisms.

Conclusion

Thanks to a combined theoretical and experimental investigation, we are able to draw a microscopic picture for the adsorption of molecular oxygen on LuPc 2 film. We demonstrate that double-decker phthalocyanines of lutetium are poorly reactive to molecular oxygen and appears as a high tolerant organic molecules leaving the Lu ion unaffected by the process.

The experimental XPS data supported by computed XPS carried out for the first time on individual molecules, confirm that oxygen molecules weakly physisorb nearby the pyrrolic nitrogen of the ligand plateau.

The high tolerance of double-decker phthalocyanine towards molecular oxygen can open promising perspectives for the design of novel materials to be applied in the next generation of organic based electronic devices operating under ambient conditions.

Figure 1 :

 1 Figure 1: DFT optimized structure of an isolated LuPc 2 molecule.

Figure 2 .

 2 Figure 2. Fitting results are detailed in Tab. 1 and Tab. 2 for the N 1s and C 1s levels respectively and are similar to already published data[START_REF] Bidermane | Experimental and Theoretical Study of Electronic Structure of Lutetium Bi-phthalocyanine[END_REF][START_REF] Toader | Initial Growth of Lutetium(III) Bis-phthalocyanine on Ag(111) Surface[END_REF] on pristine LuPc 2 films, except for the resolution which is inferior in our case.

Figure 2 :

 2 Figure 2: XPS spectra of C 1s (left) and N 1s (right) for increasing oxygen exposure at 90 K on a LuPc 2 film.

  10 eV. No asymmetry is taken into account while fitting the N 1s peak as the two peaks are supposed to have the same intensity (i.e. same number of atoms contributing to the peak). After the first oxygen exposition (Figure 2 spectrum (N2)), a third component is evidenced, at BE of 399.65 eV, referred to as N+O. After the second exposition (Figure 2 (N3)), an increase of this N+O component is clearly observed together with the appearance of a second additional component at 401.25 eV, referred to as N+O sh-up .

Figure 3

 3 Figure 3 shows the six most stable DFT optimized structures of one O 2 molecule adsorbed on one LuPc 2 labeled from a) to f). In Tab. 3, the adsorption energies of the structures shown in Figure 3 are given together with their O 2 binding distances. For these structures, the computed adsorption energies are found between -270 and -200 meV, the distances between O 2 and the Pc rings are in the range of 3.1 to 3.7 Å and the distances between Lu ion

Figure 3 :Table 3 : 3 .

 333 Figure 3: Optimized adsorption structures, superimposed onto each other, a) to f) of O 2 on LuPc 2 in i) top view and ii) side view.

Figure 4 :

 4 Photon Energy (eV) 395 420 415 410 405 400

Figure 5 :

 5 Figure 5: Energy diagram of the allowed NEXAFS transitions for the IPP geometry before (left) and after (right) oxygen exposure derived from XPS and NEXAFS.

  (A)), and it has to be remembered that in this polarization the two π * resonances can be attributed to different N atoms: the first one is mostly due to N azabridge →LUMO, while the second one is mostly due to N pyrrolic →LUMO.[START_REF] Bidermane | Experimental and Theoretical Study of Electronic Structure of Lutetium Bi-phthalocyanine[END_REF] This finding is reinforced by the intensity and the BE of the additional component in the N 1s spectra (Figure2spectrum (N3)), which suggests a strong charge withdrawal from the N atoms due to the spatial proximity. Both XPS and NEXAFS strongly suggest a charge withdrawal from the pyrrolic nitrogen, where most of the LUMO is localized. The nature of the additional component that appears at 398.9 eV can be understood by a combined interpretation of the XPS and NEXAFS spectra and is depicted by the diagram in Figure5. It is assigned to the same transition N 1s → π * as the other two resonances, but from the nitrogen atom close to the oxygen adsorption site (illustrated by violet arrow and responsible for the N+O component in the XPS spectrum).

Table 1 :

 1 Peak fitting parameters of the N 1s spectra including BE, Area (normalized on the height of the first peak). Sh-up indice refers to shake-up satellite.

	Peak	Position (eV) Area
		Clean
	N 1s	398.15	0.95
	N 1s sh-up	399.10	0.30
	After 3600 L of O 2
	N 1s	398.15	0.95
	N 1s sh-up	399.10	0.30
	N+O	399.65	0.19
	After 90000 L of O 2
	N 1s	398.15	0.90
	N 1s sh-up	399.10	0.35
	N+O	399.55	0.30
	N+O sh-up	401.25	0.05

Table 2 :

 2 Peak fitting parameters of the C 1s spectra including BE and Area. Sh -up indice refers to shake-up satellite.

	Peak	Position (eV) Area
	As deposited
	C benz	284.25	1.35
	C benz/sh-up	286.30	0.10
	C pyrr	285.40	0.3
	C pyr/sh-up	287.15	0.10
	After 3600 L of O 2
	C benz	284.25	1.35
	C benz/sh-up	286.30	0.10
	C pyrr	285.40	0.35
	C pyr/sh-up	287.15	0.10
	C+O	286.55	0.05
	After 90000 L of O 2
	C benz	284.25	1.45
	C benz/sh-up	286.30	0.10
	C pyrr	285.35	0.35
	C pyr/sh-up	287.15	0.10
	C+O	286.50	0.15
	C+O sh-up	288.10	0.05

is due to the shake-up of the pyrrolic carbon main peak, found at 287.15 eV. After the first oxygen exposition (Figure

2

spectrum (C2)) a fifth component is present, at BE of 286.55 eV, called C+O. After the 90000 L O 2 dose, the C+O peak increases in intensity, and also a sixth peak is appearing at 288.10 eV, referred to as C+O sh-up . No changes could be observed in the Lu 3d and 4f peaks. The XPS results reveal the presence of additional features in the N 1s and C 1s core levels whose intensities increase with the oxygen dose at 90 K.

Table 4 :

 4 XP ranges for N 1s (XP N ) and C 1s (XP C ) levels of the single LuPc 2 and O 2 optimized structures according to Figure3. Large values leading to XPS shifts are highlighted.

		XP N (eV) XP C (eV)
	LuPc 2	0.044	1.759
	a)	0.025	1.755
	b)	0.033	1.769
	c)	0.378	2.584
	d)	0.044	1.789
	e)	0.466	2.090
	f)	0.026	1.772

  but again, no literature could be found on this specific or closely related systems. Therefore, the analysis of the BE in the XPS data suggest the adsorption of the oxygen molecule close to carbon atoms, and most likely the pyrrolic carbon atoms, since N 1s is largely affected. This finding is supported by our DFT calculations for which the two adsorption geometries presenting the larger changes in BE are obtained for oxygen molecules adsorbed near or on a pyrrole ring. However, our DFT results do not address the attachment of more than one O 2 molecule at one LuPc 2 which may cause slight changes in the adsorption energies. Because the XPS measurements detect large shifts in both N 1s and C 1s resonances due to O 2 adsorption, which was not seen in the simulated XPS for the most stable adsorption site computed (a), it can be assumed that other effects (e.g. vibrations) prevent the structure a) from being the most stable structure or/and that several O 2 molecules can adsorb on one LuPc 2 simultaneously, molecules in a) configuration leaving the experimental spectra unaffected. Then, other structures must be present on the LuPc 2 to account for experimental findings. In this regard, structures c) and e) give qualitative agreement between computed and experimental XPS changes. Moreover, the adsorption energies of c) and e) are only slightly larger than for a), i.e. 50 and 60 meV,

Acknowledgments

We are grateful to D1011 beamline staff for providing support and assistance during the beamtime. We acknowledge the financial support from the Knut and Alice Wallenberg foundation, the European Union Seventh Framework Programme under grant agreement n • 607232 THINFACE, the Swedish Research Council (Grant 2014-3776) and the Swedish National Infrastructure for Computing (SNIC).

Supporting Information

XPS spectra of C1s and N1s for increasing dose of oxygen and after annealing at room temperature.

NEXAFS at the O K edge.


Signe Stéphanie

Farzam Mamba

Michael Zoueshtiagh

Baudoin

S Signe Mamba

F Zoueshtiagh

M Baudoin

Pressure-driven dynamics of liquid plugs in rectangular microchannels: influence of the transition between quasi-static and dynamic film deposition regimes

Keywords: Liquid plug, rectangular channel, Hele-Shaw geometry, wet fraction, periodic forcing

come

Introduction

Two-phase gas-liquid flow in microfluidic devices is a hydrodynamic problem with practical applications in a variety of engineered systems including flows in microreactors [START_REF] Song | A microfluidic system for controlling reaction networks in time[END_REF][START_REF] Gunther | Transport and reaction in microscale segmented gas-liquid flow[END_REF][START_REF] Assmann | Extraction in microreactors: intensification by adding an inert gas phase[END_REF][START_REF] Sobieszuk | Hydrodynamics and mass transfer in gas-liquid flows in microreactors[END_REF], enhanced oil recovery [START_REF] Havre | Taming slug flow in pipelines[END_REF], flow in porous media [START_REF] Lenormand | Mechanisms of the displacement of one fluid by another in a network of capillary ducts[END_REF][START_REF] Dias | Network models for two-phase flow in porous media part 1. immiscible microdisplacement of non-wetting fluids[END_REF][START_REF] Stark | The motion of long bubbles in a network of tubes[END_REF], film coating and biomechanical systems as pulmonary flows [START_REF] Kamm | Is airway closure caused by a liquid film instability?[END_REF][START_REF] Heil | The mechanics of airway closure[END_REF][START_REF] Grotberg | Respiratory fluid mechanics[END_REF]. The significant scope of the topic motivated early studies on the dynamics of these interfacial flows, first in cylindrical tubes [START_REF] Fairbrother | Studies in electro-endosmosis[END_REF][START_REF] Taylor | Deposition of a viscous fluid on the wall of a tube[END_REF][START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] and soon after in polygonal channels [START_REF] Saffman | The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid[END_REF][START_REF] Jensen | Effect of gravity on the saffman-taylor meniscus: Theory and experiment[END_REF][START_REF] Ratulowski | Transport of gas bubbles in capillaries[END_REF]Wong et al., 1995b;[START_REF] Kolb | Coating the inside of a capillary of square cross section[END_REF][START_REF] Thulasidas | Bubble-train flow in capillaries of circular and square cross section[END_REF][START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF][START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF][START_REF] Han | Measurement of liquid film thickness in micro square channel[END_REF][START_REF] Han | Measurement of liquid film thickness in a micro parallel channel interferometer and laser focus displacement meter[END_REF]. In the last decades, the interest in segmented gas-liquid flows in polygonal, and in particular, rectangular microchannels has further grown with the development of soft lithography techniques in microfluidics [START_REF] Duffy | Rapid prototyping of microfluidic systems in poly (dimethylsiloxane)[END_REF][START_REF] Quake | From micro-to nanofabrication with soft materials[END_REF], which enable simple design of complex microchannels with rectangular cross-sections (Anderson et al., 2000).

A finite volume of liquid (liquid plug or slug) that is displaced by an air finger at constant flow rate or pressure head (Taylor flow) in a tube, leaves on the walls a trailing liquid film. Its thickness can be quantified by the so-called wet fraction, which is the proportion of the tube section occupied by the liquid film. In cylindrical capillary tubes, this parameter increases monotonically with the dimensionless velocity of the meniscus (the so-called capillary number Ca) with a Ca 2/3 law [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] at low capillary number. This law can be further extended to larger capillary number as demonstrated by [START_REF] Aussillous | Quick deposition of a fluid on the wall of a tube[END_REF]. In polygonal microchannels however, a transition occurs at a critical capillary number between two radically different regimes: Under this critical parameter, the shape of the meniscus remains close to the static one and most of the fluid deposition occurs in the corners of the tube. Indeed, the static shape of the meniscus cannot follow the singular shape of the sharp edges, which leads to significant deposition in the tube corners, even at vanishingly small capillary number. Thus, the wet fraction slightly relies on the plug dynamics and the relative variation of the wet fraction with the capillary number remains weak. (Wong et al., 1995a,b). Above this critical number, the fluid deposition resulting from the deformation of the rear meniscus induced by the flow overcomes the static one. In this case, the wet fraction becomes again strongly dependent over the capillary number, similarly to what is observed in cylindrical channels. Nevertheless, this process in polygonal channels also depends on the tube geometry [START_REF] Jensen | Effect of gravity on the saffman-taylor meniscus: Theory and experiment[END_REF]Wong et al., 1995a,b;[START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF][START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF]. [START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF][START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF] showed that it is possible to extend the laws introduced for cylindrical tubes to rectangular tubes, providing the introduction of an aspect-ratio-dependent capillary number .

This liquid film deposition process induces a dramatic acceleration of a liquid plug when it is pushed at constant pressure head [START_REF] Baudoin | Airway reopening through catastrophic events in a hierarchical network[END_REF]. Indeed the diminution of the plug size leads to a reduction of the viscous resistance of the plug to motion, itself leading to an acceleration of the plug and thus more fluid deposition. More recently, it has been shown experimentally by [START_REF] Magniez | Dynamics of liquid plugs in prewetted capillary tubes: from acceleration and rupture to deceleration and airway obstruction[END_REF] that the inverse behaviour (progressive slow down and growth of the liquid plug) might also be observed in prewetted capillary tubes depending on the value of the driving pressure and the thickness of the prewetting film. The acceleration and rupture of a liquid plug has also been evidenced in complex tree geometries [START_REF] Baudoin | Airway reopening through catastrophic events in a hierarchical network[END_REF][START_REF] Song | The air-liquid flow in a microfluidic airway tree[END_REF] and for cyclic pressure forcing [START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF]. In the latter case however, both the diminution of the viscous resistance and interfacial resistance (due to lubrication effects) at each cycle contribute to the plug acceleration and breaking.

Nevertheless, all the aforementioned studies were conducted in cylindrical tubes or at capillary numbers well above the critical capillary number. In this paper, we study experimentally and theoretically the influence of the transition between quasi-static and dynamic fluid deposition process on the dynamics of liquid plugs in rectangular channels pushed either with a unidirectional or a cyclic pressure forcing. In both cases, it is shown that the transition between these two regimes leads to a dramatic acceleration of the plug eventually leading to its rupture. For cyclic forcing, it is shown that under a critical length the plug dynamics is unstable and leads to the plug rupture while above it is stable and periodic; the experimental results are recovered from a reduced dimension model, inspired from previous theoretical developments by [START_REF] Baudoin | Airway reopening through catastrophic events in a hierarchical network[END_REF], [START_REF] Magniez | Dynamics of liquid plugs in prewetted capillary tubes: from acceleration and rupture to deceleration and airway obstruction[END_REF] and [START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF] adapted here to take into account (i) the modifications of the laws in the rectangular geometry, (ii) lubrication effects resulting from the back and forth motion of the liquid plug on a prewetted tube, and (iii) the transition between quasi-static and dynamic film deposition. The first and second sections provide the experimental and model details. The third and fourth sections explore respectively the response of liquid plugs to unidirectional and periodic cyclic pressure forcings and compare the observed dynamics to results in cylindrical tubes.

Method

The experimental setup is represented on Fig. 1. The experiments are conducted in rectangular polydimethylsiloxane (PDMS) microfluidic channels of high aspect ratio obtained by standard photolithography techniques: A mold is etched by depositing a layer of photoresist resin (Microchem, SU8-2035) on a silicon wafer. This layer is spin-coated and patterned by standard photolithography. The spin-coating speed combined with the choice of the photoresist set the height h = 45 ± 2µm of the microfluidic channels, while the width w = 785µm is controlled by the design of the patterned masks which is used during the UV exposure. These values of the channel width and height were measured afterwards with a profilometer (Dektak XTL). After exposure, the film is developed in an organic solvent solution (SU-8 developer) to yield the negative of the channel design. This SU8 mold was used to pour PDMS (Dow Corming, Sylgard 184) whose polymerisation was obtained by curing it at 100 • C. The microfluidic channel is then cut out and bonded on a glass microscope slide by passing the two surfaces in an oxygen plasma. The microscope slides are covered by a thin PDMS membrane in order to guarantee identical boundary conditions to the four channel walls. The total length of the channel is 6 cm.

Then, perfluorodecalin (PFD) liquid plugs are created in the channel by pushing alternatively some liquid and some air at a Y-junction with a syringe pump and a MFCS Fluigent pressure controller respectively, connected to both entrances of the microfluidic device. Perfluorodecalin was used for its hydrophilic properties with PDMS (static contact angle θ s = 23 ± 1 • ) and since it does not swell PDMS [START_REF] Lee | Solvent compatibility of poly (dimethylsiloxane)-based microfluidic devices[END_REF]. This fluorocarbon has a dynamic viscosity µ = 5.1 × 10 -3 Pa.s, surface tension σ = 19.3 × 10 -3 N/m and density ρ = 1.9 × 10 3 kg/m 3 . Then, air is blown in the channel at low pressure to bring the liquid plug to the center of the microfluidic channel and stopped manually when the target position is reached. Finally, the plug motion is forced with either a unidirectional or cyclic periodic pressure forcing with the MFCS programmable pressure controller. For cyclic forcing one entrance and the exit of the channel are connected to two channels of the MFCS pressure controller.

Then an overpressure (compared to atmospheric pressure) is applied alternatively to each end of the channel while the other is set to atmospheric pressure. The resulting shapes of the pressure forcing measured with an internal pressure sensors in these two cases are represented on Fig. 2. For cyclic forcing the period was fixed to 2T = 4 s or 2T = 6 s, with T the duration of a half cycle. The measured unidirectional pressure driving can be approximated by the following analytical expression based on Gompertz functions:

P t = 900e -3e -10t , (1)
while the cyclic forcing can be approximated by the expression : Pc = 1200e -2.5e -3(t-nT ) Pa and P d = 1200e -1.2(t-nT ) e -0.02e -1.2(t-nT ) Pa (red dotted line) and T = 2.12 is the half period.

∆P t = 1200e -6e -3t Pa for t ∈ [0; T ] (2) ∆P t = (-1) n (P c -P d ) for t ∈ [nT ; (n + 1)T ] (3)
with P c = 1200e -2.5e -3(t-nT ) Pa (4)

and P d = 1200e -1.2(t-nT ) e -0.02e -1.2(t-nT ) Pa (5)

Initially the microfluidic channel is dry. Hence, the liquid plug moves on a dry portion of the channel for unidirectional forcing. Nevertheless, the motion of the liquid plug leaves a trailing liquid film on the walls and in the corners of the channel. Thus, for cyclic forcing the liquid plug moves on a prewetted channel after the first cycle as long as it is moves on a portion of the channel already visited by the liquid plug in the previous back and forth motions.

Experiments are recorded with a Photron SA3 high speed camera mounted on a Z16 Leica Microscope. The resolution of the camera used in the experiments is 1024 × 64 pixels, the acquisition frame rate 125 images/s and the shutter time 1/3000 s. The image analysis is then performed using ImageJ software and Matlab. The plug evolution is characterised by monitoring the positions of the rear meniscus x r and front meniscus x f (see Fig. 1), and deducing the evolution of the plug length L p (t) = x f (t) -x r (t) and the speed of the rear meniscus of the plug U r = dx r /dt.

Model of a plug flow in dry and prewetted rectangular microfluidic channels

The model derived in this paper to describe the dynamics of liquid plugs under pressure forcing combines previous theoretical developments by [START_REF] Baudoin | Airway reopening through catastrophic events in a hierarchical network[END_REF], [START_REF] Magniez | Dynamics of liquid plugs in prewetted capillary tubes: from acceleration and rupture to deceleration and airway obstruction[END_REF] and [START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF], and integrates additional elements to include the transition between quasi-static and dynamic liquid film deposition.

Dimensional analysis and characterisation of the regime

In this problem, we consider a single liquid plug of initial length L 0 set into motion in a rectangular microfluidic channel under the unidirectional or periodic forcings represented on Fig. 2.

Dimensionless number Formula Estimation

τ c /τ exp l c /U 1.8 × 10 -2 τ d /τ exp ρl c 2 /µ 1.3 × 10 -2 Re τ d /τ c 0.7 W e ρU 2 l c /σ 1.9 × 10 -3 Ca µU /σ 2.6 × 10 -3 Bo ρgh 2 /σ 2 × 10 -3
Table 1. Values of the key dimensionless parameters associated with the mean characteristic velocity U = 1cm/s.

The characteristic parameters in this problem are the width of the microfluidic channel w, its height h, the viscosity of the liquid plug µ, the surface tension σ, the speed of the liquid plug U and the characteristic time associated with the plug evolution in the experiments τ exp . For cyclic forcing, this time is simply the half period of the signal τ exp = T , while for unidirectional forcing, it is the time required for the plug to rupture. In the following, the geometry will be characterised by the aspect ratio α = w/h and the characteristic length scale l c = √ wh. From these parameters, we can construct the characteristic convection time τ c = l c /U , and the characteristic viscous diffusion time τ d = ρl c 2 /µ. Then we can characterise the flow regime by introducing the following dimensionless numbers: the Reynolds number (Re = τ d /τ c ) which compares convection to viscous diffusion, the Weber number (W e = ρU 2 l c /σ) which compares inertia to surface tension, the capillary number (Ca = µU/σ) which compares viscous effects to surface tension effects, the Bond number which compare gravity to surface tension (Bo = ρgh 2 /σ), and finally the ratio of the experimental characteristic time τ exp to the convective and diffusion times τ exp /τ c and τ exp /τ d .

In the experiments the average velocity is typically U mean = 1 cm/s and the maximal velocity U max = 4.5 cm/s. The time required for the plug to rupture varies between 0.5 s and 5 s for unidirectional forcing and the half period of periodic forcing is T = 2 s or T = 3 s. These values enable the estimation of the dimensionless numbers introduced previously and summarised in Table 1 (for these estimations, we take τ exp = 1 s and U = 1 cm/s).

These values of the dimensionless numbers indicate that, globally, surface tension effects are dominant over viscous effects, themselves being dominant over inertial effects. This means that away from the walls, the shape of the meniscus is mainly dictated by the minimisation of the interfacial capillary energy. This undeformed part of the meniscus is called the static meniscus. Nevertheless, as demonstrated first by Bretherton [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] in cylindrical tubes, viscous effects still play an important role close to the walls due to the incompatibility between the no-slip condition (null velocity at the walls) and the homogeneous motion of an undeformed meniscus. This singularity leads to the existence of large shear stresses close to the walls and, hence, the deposition of a thin trailing film behind the rear meniscus and a change in the apparent contact angle for the front meniscus. Thus viscous effects still play an important role close to the walls despite the low value of the capillary number. Then, the relative importance of unsteady effects can be estimated from the ratios τ c /τ exp and τ d /τ exp . Since, these two ratios are small, the unsteady terms can be neglected in Navier-Stokes equations even if the plug evolves over time: the flow is quasistatic. Finally, the small value of the Bond number indicates that gravity effects can be neglected. Based on this analysis, we expect the liquid plugs dynamics in the present problem to be considered as a quasi-static visco-capillary flow governed by steady Stokes equation.

Model of the plug dynamics.

The model describing the pressure drop in the microfluidic channel is obtained by equalizing the driving pressure head ∆P t (Fig. 2) to the sum of the pressure drop resulting from viscous dissipation in the bulk of the liquid plug ∆P bulk visc , the pressure drops at the front and rear meniscus of the plug ∆P men f ront , ∆P men rear and the pressure drop ∆P bubble inside the air. Simple estimation of these pressure drops show that this latest contribution can be ignored compared to the other ones (see e.g. Kreutzer et al. (2005b)). Thus, the steady state balance of pressure across the liquid plug becomes:

∆P t = ∆P bulk visc + ∆P men f ront + ∆P men rear (6)

Viscous pressure drop

The pressure drop resulting from a laminar flow of a fluid in a rectangular geometry is given by [START_REF] White | Fluid mechanics[END_REF]:

∆P bulk visc = aµQL wh 3 (7) with a = 12 1- 192 π 5 α tanh πα 2 -1
, L the portion of the tube considered and Q the flow rate. In the limit α 1 considered here (α = 17.5 in the experiments) and for a liquid plug of length L p , this expression becomes:

∆P bulk visc = 12µQL p wh 3 (8)
where Q = U r S r is the flow rate, S r is the the cross sectional area open to air behind the liquid plug and U r the speed of the rear meniscus. This expression relies on two assumptions: (i) it assumes that the pressure drop inside the plug follows a Hagen-Poiseuille law despite the finite size of the plug and the recirculation occurring close to the menisci and (ii) it assumes the same speed for the front and rear meniscus. The validity of the first approximation has been tested numerically with the OpenFoam Volume of Fluid code [START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF] in 2D. These 2D simulations show that equation ( 8) is an excellent approximation of the viscous pressure drop (error < 4.5%) as long as the length of the plug remains larger than the height of the channel. For smaller plugs the discrepancy increases progressively but, in this case, the pressure drops at the menisci strongly dominate over bulk viscous pressure drop leading to minor effects of the error on the overall plug dynamics. The second approximation amounts to neglect the evolution speed of the plug length dL p /dt compared to the translational speed of rear meniscus dx r /dt, since dx f /dt = dx r /dt+dL p /dt. Experimental measurements of the speed of the front and rear meniscus show that this approximation holds within a few percent of accuracy. In the remaining part of the manuscript, we will therefore neglect the difference between the front and rear menisci in the estimation of the pressure drops, and the capillary number Ca will therefore be constructed on the rear meniscus velocity: Ca = µU r /σ.

Front meniscus pressure drop

When a liquid plug is at rest (no pressure head) in a rectangular tube, its front and rear menisci adopt complex shapes minimising the interfacial energy. This minimisation problem can be solved with the method of Lagrange multipliers. The solution is a constant mean curvature (CMC) surface verifying the wetting conditions at the walls. In a rectangular geometry, this shape is rather complex (Wong et al., 1995a) and there is no analytical expression of these surfaces geometry. Nevertheless, what matters when we consider the motion of a liquid plug is not the static shape of the meniscus but the departure from this static shape when the plug moves. Indeed the Laplace pressure jumps resulting from the curvatures of the front and rear meniscus at rest compensate one another leading to a zero contribution. Therefore we will only consider in the following the dynamic pressure jumps at the meniscus, that is to say the Young-Laplace pressure jumps when the plug is moving minus their value at rest.

The computation of the dynamic pressure jumps at the front meniscus can be greatly simplified in a rectangular geometry with high aspect ratio α 1. In this case, the principal curvature in one direction κ h ≈ 2 cos(θ a d )/h is strongly dominant over the curvature in the other direction κ w ≈ 2 cos(θ a d )/w, where θ a d is the advancing dynamic apparent contact angle. Based on Young-Laplace equation, the dynamic pressure jump at the front meniscus can thus be estimated from the formula:

∆P men f ront ≈ -σ(κ h -κ s h ), with κ h ≈ 2 cos(θ a d )/h (9)
and κ s h ≈ 2/h the principal curvature in the vertical direction at rest. In the limit of low capillary numbers, asymptotic expansion leads to:

cosθ a d ∼ (1 -θ a d 2 /2),
and the dynamic pressure drop becomes:

∆P men f ront = 2σ h θ a d 2 2 (10)
Of course, this expression is an approximation since (i) it neglects the horizontal curvature compared to the vertical one and (ii) it neglects the thickness of the prewetting film (if the plug is moving on a prewetted capillary tube). This expression is therefore only valid for high aspects ratios h/w, low capillary numbers and thin prewetting films. The next step is to determine the value of the dynamic apparent contact angle θ a d as a function of the capillary number. On a dry substrate, the dynamic contact angle can be estimated from Hoffman-Tanner's law:

θ a d = ECa 1/3 (11)
with E a constant of order (4 -5) for a cylindrical dry tube as reported by [START_REF] Hoffman | A study of the advancing interface. i. interface shape in liquid-gas systems[END_REF] and [START_REF] Tanner | The spreading of silicone oil drops on horizontal surfaces[END_REF]. [START_REF] Bico | Falling slugs[END_REF] 2013) measured a value of E = 4.9 for a perfluorodecalin liquid plug moving in PDMS rectangular capillary tubes with high aspect ratio. This value is in good agreement with the experiments performed in this paper.

A theoretical expression of the dynamic contact angle on a prewetted surface was proposed theoretically by [START_REF] Chebbi | Deformation of advancing gas-liquid interfaces in capillary tubes[END_REF] and validated experimentally by [START_REF] Magniez | Dynamics of liquid plugs in prewetted capillary tubes: from acceleration and rupture to deceleration and airway obstruction[END_REF] for the motion of a liquid plug in a prewetted cylindrical capillary tube:

tanθ a d = (3Ca) 1/3 f ((3Ca) -2/3 cosθ a d h f /R) (12)
with h f /R the thickness of the liquid film ahead of the liquid plug, f (y) = 3 j=0 b n [log 10 y] n and the coefficients b n are tabulated in [START_REF] Chebbi | Deformation of advancing gas-liquid interfaces in capillary tubes[END_REF]. This expression enables to integrate a lubrication effect induced by the presence of a prewetting film, which facilitates the displacement of the front meniscus and thus reduces the pressure jump. To the best of our knowledge, no rigorous derivation of an analytical formula exists for rectangular geometries. Nevertheless, a similar expression as equation ( 12) is expected in rectangular geometries with high aspect ratios α 1. Indeed, Chebbi (similarly to Bretherton) derived the above theoretical expression in the approximation of thin prewetting liquid film h f compared to the radius of the tube (h f /R 1). In this approximation, the radial curvature of the tube is locally neglected and the problem solved is identical to a 2D planar problem.

In the rectangular configuration, an estimation of the thickness of the prewetting film in the vertical direction is nevertheless missing. To adapt this formula to rectangular geometries, Baudoin et al. [START_REF] Baudoin | Airway reopening through catastrophic events in a hierarchical network[END_REF] proposed to estimate the relative thickness of the prewetting film by the formula 1 -Sf with Sf = S f /(h w) the dimensionless cross sectional area open to air in front of the plug. For rectangular channel with large aspect ratios, it is expected that this expression slightly overestimates lubrication effects since the prewetting film is always thicker on the lateral walls than in the central ones.

Finally, as demonstrated by Signe-Mamba et al. ( 2018), a good approximation of the implicit formula (12) in the limit of low capillary numbers is:

θ a d = F Ca -1/3 (13) with F = 3 1/3 b 0 + b 1 log 10 (A) + b 2 [log 10 (A)] 2 + b 3 [log 10 (A)] 3 and A = (3Ca) -2/3
1 -Sf . This expression shares some similarities with Hoffman-Tanner's law but this time the coefficient F depends on the capillary number, underlining the lubrication effect.

Rear meniscus pressure drop

The dynamical pressure drop at the rear meniscus was calculated theoretically by [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] in cylindrical geometries, Wong et al. (1995a,b) for polygonal channels, and later on numerically by [START_REF] Hazel | The steady propagation of a semi-infinite bubble into a tube of elliptical or rectangular cross-section[END_REF] for rectangular microfluidic channels of different aspect ratios at finite capillary numbers. Based on the results of [START_REF] Hazel | The steady propagation of a semi-infinite bubble into a tube of elliptical or rectangular cross-section[END_REF], it is possible to infer the following formula for the dynamic pressure jump at the rear meniscus [START_REF] Baudoin | Airway reopening through catastrophic events in a hierarchical network[END_REF]:

∆P men rear = 2σ h Df (α)Ca 2/3 (14)
with f (α) = (0.52 + 0.48/α) and D = 4.1, a constant obtained by [START_REF] Baudoin | Airway reopening through catastrophic events in a hierarchical network[END_REF] from least square fit of the numerical data points of Fig. 8 of [START_REF] Hazel | The steady propagation of a semi-infinite bubble into a tube of elliptical or rectangular cross-section[END_REF] in the range Ca ∈ [10 -3 , 0.3], with Broyden-Fletcher-Goldfarb-Shanno (BFGS) minimization algorithm.

Total pressure drop

If we combine equations ( 8), ( 10), ( 11), ( 13) and ( 14), the total pressure drop across the liquid plug becomes:

Dry tube: ∆P t = 12σS r L p wh 3 Ca + σ h E 2 + 2Df (α) Ca 2/3 (15) Wet tube: ∆P t = 12σS r L p wh 3 Ca + σ h F 2 + 2Df (α) Ca 2/3 with E = 4.9, F = 3 1/3 b 0 +b 1 log 10 (A)+b 2 [log 10 (A)] 2 +b 3 [log 10 (A)]
3 and A = (3Ca) -2/3 1 -Sf . These equations can be written under dimensionless form by introducing the characteristic length l c = √ wh and the characteristic pressure variation σ/h:

Dry tube: ∆ Pt = 12 √ α Sr Lp Ca + E 2 + 2Df (α) Ca 2/3 (16) Wet tube: ∆ Pt = 12 √ α Sr Lp Ca + F 2 + 2Df (α) Ca 2/3
where the tildes indicate dimensionless functions. To achieve a closed set of equations, two equations are missing: one determining the evolution of the plug length Lp and one determining the fluid deposition process on the walls and consequently Sr and Sf .

Evolution of the plug length

The first equation is simply obtained from a mass balance between the amount of liquid that the plug collects and loses. Let S 0 = wh be the cross section of channel, and S r and S f the sections of the tube open to air behind and in front of the plug respectively. The mass balance becomes:

S 0 dL p = (S 0 -S f )dx f -(S 0 -S r )dx r ( 17
)
with dx r and dL p the displacement of the rear interface and the variation of the plug length during an infinitesimal time step dt and dx f = dL p +dx r . Therefore, the equation giving the evolution of the length of the liquid plug takes the form:

dL p dt = S r S f -1 U = σ µ S r S f -1 Ca (18)
where S r depends only on the dimensionless speed of the rear meniscus S r = S r (Ca) and S f depends on the history of the liquid deposition on the walls and the position of the front meniscus S f = S f (x f , t). If we introduce the characteristic length l c and the viscocapillary time scale τ = µl c /σ, we obtain the dimensionless equation:

d Lp d t = Sr (Ca) Sf ( xr , t) -1 Ca (19)
When the plug moves on a dry surface, then Sf = 1. Otherwise the value of Sf is either inferred from the initial condition if the plug moves on a prewetted tube or from a memory of the liquid deposition by the liquid plug when the plug undergoes cyclic motion (see [START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF] for cylindrical tubes). In many papers, the wet fraction m is introduced instead of the air fraction S r . The wet fraction is the relative portion of the tube section occupied by the liquid. These two parameters are linked by the formula: m = 1 -Sr .

3.2.6. Quasi-static and dynamic wet fraction An air finger pushing a liquid plug in a cylindrical tube adopts a cylindrical shape closely fitting the shape of the cylindrical tube. Such close fitting is not possible in rectangular channels due to the presence of sharp corners. In the low Bond number (Bo 1) and low capillary number limit (Ca 1), the section of the bubble far behind the rear meniscus can adopt two configurations [START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF] schematically represented on Fig. 3 depending on the value of the capillary number:

• (i) When Ca → 0, in the quasi-static limit, the liquid only covers the corners of the tube and the limit between the liquid and the air are four quarter circles. In this case the wet fraction tends to an asymptotic value called the quasi-static wet fraction m s , which depends only on the aspect ratio of the channel α: m s = m s (α). Indeed, the liquid deposition in the corners of the tube in this regime mostly relies on the static shape of the rear meniscus, and thus the evolution of the wet fraction with the capillary number is weak.

• (ii) when Ca Ca c , the limit between the liquid and the air becomes two half circle on the side and a liquid film covers the walls in the center of the channel [START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF]. In this case, the wet fraction more strongly depends on the dynamics of the plug and hence on the capillary number.

We will call it the dynamic wet fraction m d = m d (Ca, α).

The transition between these two configurations is progressive (see [START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF]) and occurs at a critical capillary number Ca c whose value was found in our experiments to lie around Ca c = (2 ± 0.1) × 10 -3 , a value coherent with the measurements of [START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF] at similar aspect ratio (α = 15): Ca c ∼ 10 -3 .

The quasi-static wet fraction m s in the absence of gravity (Bo = 0) was theoretically predicted in Wong et al. (1995a):

m s = (4 -π) rs 2 (20)
with rs = r s /l c the dimensionless radius of curvature of the four quarter circles delimiting the liquid and the air in the corners of the tube and rs =

√ α α + 1 + ((α -1) 2 + πα) 1/2 (21)
For an aspect ratio α = 17.5 such as in our experiments, this formula gives the value m s = 0.012. The wet fraction can be estimated experimentally by monitoring the evolution of the plug size as a function of the rear meniscus velocity when Ca Ca c . We found the average estimate of this parameter (over all the experiments performed and described in the next section) to be: m s = 0.03. This value is larger than the theoretical value predicted by Wong et al. (1995a). Nevertheless the theoretical value derived by this author was obtained with zero influence of gravity (Bo = 0). For finite values of the Bond number, it was shown experimentally by [START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF] (see their figure 4) that gravity tends to significantly increase the static wet fraction (these authors found m s ≈ 0.065 for Bo ∼ 1 and α = 15). In our case, though small the Bond number (table 1) is not null, which might explain the larger value of the static wet fraction measured experimentally than expected theoretically. Another tentative explanation would be that, owing to the finite length of the tube, we did not reached the final equilibrium value. For the simulations, we will therefore adopt the constant value:

m s = 1 -Sr = 0.03 for Ca < Ca c ( 22
)
under the critical capillary number. The evolution of the dynamic wet fraction in square and rectangular microfluidic channels was investigated experimentally by [START_REF] Kolb | Coating the inside of a capillary of square cross section[END_REF] 2011) and numerically by [START_REF] Hazel | The steady propagation of a semi-infinite bubble into a tube of elliptical or rectangular cross-section[END_REF]; [START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF]; Kreutzer et al. (2005a). In particular, [START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF][START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF] found that the measured and simulated evolution of the dynamic wet fraction m d as a function of the capillary number collapse for all aspect ratio α providing the introduction of an effective capillary number Ĉa = 1 + α 2 /α 2 t Ca with α t = 6.4. Thus, the dynamic wet fraction in a rectangular channel with any aspect ratio can be inferred from the behaviour in a square tube. This scaling subsist even for finite Bond numbers as demonstrated both theoretically and numerically by [START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF][START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF]. Of course this scaling is only valid for Ca > Ca c since otherwise, the fluid deposition does not depend on Ca but strongly depends on α. Thus, combining (i) the scaling law proposed by [START_REF] Lozàr | The steady propagation of an air finger into a rectangular tube[END_REF] for the effective capillary number, (ii) Aussilous & Quéré law [START_REF] Aussillous | Quick deposition of a fluid on the wall of a tube[END_REF] for the evolution of the wet fraction as a function of the capillary number and (iii) the matching condition at the critical capillary number between the quasi-static and dynamic behaviour gives:

m d = 1 -Sr = m s + G ( Ĉa/ Ĉa c ) 2/3 -1 1 + H(( Ĉa/ Ĉa c ) 2/3 -1) for Ca > Ca c ( 23
)
where the coefficients G = 1.05 and H = 1.75 are obtained from best fit with [START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF] experimental data.

Finally, Sr is taken as Sr = m s -1 when the plug moves at a capillary number lower than Ca c and Sr = m d -1 above. This is of course an approximation since the experimentally observed transition between the dynamic and quasistatic regime at Ca ∼ Ca c is more progressive [START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF].

Numerical resolution of the equations

The closed set of equations ( 16), ( 19), ( 22) and ( 23) are solved using a first order Euler explicit scheme to predict the speed and the evolution of the length of the liquid plug. To cope with the strong acceleration of the liquid plug, an adaptive time step refinement is used: the spatial displacement ∆x = x n+1 r -x n r is kept constant and thus the time step ∆ tn at iteration n is calculated from the formula: ∆ tn = ∆xCa n-1 . Convergence on ∆x has been verified for all the simulations provided in this paper.

4. Effect of the transition between quasi-static and dynamic film deposition on the dynamics of a liquid plug driven by a unidirectional forcing

Direct experimental evidence of the transition

The transition between the quasi-static and dynamic liquid film deposition and the associated changes in the plug dynamics are evidenced on Fig. 4 a-d. These figures illustrate the evolution of a liquid plug of initial size L o = 3.5 mm driven by a unidirectional pressure head ∆P t = 1000e -3e -10t Pa. Blue curves correspond to experimental measurements while red curves correspond to simulations with the model developed in the previous section. Figure 4a is obtained by stacking snapshots of the plug evolution every 8 ms when the capillary number lies below its critical value Ca c and then later on when Ca exceeds Ca c . When Ca < Ca c no liquid film is visible on the tube lateral sides since liquid deposition only occurs in the corners of the channel, while this film is clearly visible when Ca exceeds Ca c . Figure 4b shows the position of the front and rear menisci as a function of time. Figure 4c shows the evolution of the capillary number. The black dashed line (also reported on figure 4d) marks the transition (at time t c ≈ 2.8 s) between the quasi-static and dynamic film deposition regimes. It corresponds to the time when Ca reaches the critical value Ca c = 2 ± 0.1 × 10 -3 . Before t c and after the end of the transient regime (t > t t = 0.3 s) (corresponding to the time required for the pressure controller to achieve a constant value), the increase in the capillary number is very slow. This leads to a quasi-linear variation of the plug size as a function of time as can be seen on figure 4d since the wet fraction m is quasi-constant in the quasi-static regime. Then, when the value of the capillary number overcomes the critical value Ca c , the plug undergoes a strong acceleration leading to more and more fluid deposition and eventually to the plug rupture.

Excellent agreement between the simulations (red) and experiments (blue) is achieved for the evolutions of (i) the position of the menisci (Fig. 4 b), (ii) the plug dimensionless speed (Fig. 4 c) and (iii) the plug length (Fig. 4 d). Our reduced dimension model thus properly captures the main physical ingredients. This model can be used to rationalise the observed tendencies: In the quasi-static film deposition regime, the value of the pressure head prescribes an initial value of the capillary number and the size of the plug diminishes quasi-linearly due to film deposition in the corners of the tube. This regular diminution of the plug size leads to a reduction of the viscous resistance of the plug to motion since the viscous pressure drop depends linearly on L p . This induces a slow increase in the liquid plug speed (since the viscous resistance is weak compared to interfacial resistances). Nevertheless, since the wet fraction m does not depend on Ca, there is no retroaction of the evolution of the plug speed on the liquid film deposition and thus the evolution remains relatively stable. In the dynamic film deposition regime however, the increase in the plug speed leads to more film deposition according to equation ( 23), itself leading to an acceleration of the plug speed. This retroaction is at the origin of the massive acceleration of the plug and rapid decrease in its size when Ca exceeds Ca c . This behaviour is reminiscent of what is observed in cylindrical tubes [START_REF] Magniez | Dynamics of liquid plugs in prewetted capillary tubes: from acceleration and rupture to deceleration and airway obstruction[END_REF][START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF], while the quasi-static film deposition regime only occurs when there is the presence of sharp corners.

Influence of this transition on the plugs rupture time and rupture length

We performed numerous experiments (represented on Fig. 5) for different initial plug lengths L o and two different driving pressures (∆P t = 900e -3e -10t Pa and ∆P t = 1000e -3e -10t Pa) to analyse the evolution of the plug rupture time and rupture length in rectangular microchannels. The rupture time is the time elapsed between the start of the pressure head and the plug rupture (L p = 0), that is to say when the front and rear menisci come into contact at the centerline of the channel and the plug breaks, leaving the air flow freely in the channel. The rupture length is the distance traveled by the liquid plug (D l = max(x f ) -min(x r )) before its rupture. These two parameters quantify the stability of a liquid plug to breaking. The quantitative agreement between experiments (blue stars) and simulations (red lines) enables to validate our model on an an extensive set of experimental data.

Again a transition between two distinct regimes is clearly evidenced on Fig. 5 (for both the rupture time and the rupture length) at a driving pressuredependent critical initial plug length L c o (L c o ≈ 2.6 mm for ∆P t = 900 Pa and L c o ≈ 3.1 mm for ∆P t = 1000 Pa). Under this critical value of the initial plug length L o < L c o , the initial dimensionless plug speed lies above the critical capillary number and thus the dynamics of the liquid plug is only in the dynamic film deposition regime. Thus the plug accelerates rapidly leading to rapid rupture of the plug on a short propagation length scale. Above, this critical initial length

L o > L c
o , the initial capillary number lies under the critical number Ca c and thus the plug dynamics is initially in the quasi-static film deposition regime. This regime leads to larger plug rupture time and thus propagation distance. Moreover, since in this regime the acceleration is weak, the rupture time and rupture length remain relatively linear function of the plug initial length (see Fig. 5). From this analysis, we can infer a theoretical evaluation of the critical initial length L c o , which delimits the transition between these two regimes. The rupture lengths as a function of the rupture times are given in figures (c-c1). The blue stars correspond to experiments and the red lines to simulations with the model developed in section 3. Error bars quantify the error in the determination of the plug initial length owing to the large field of view and the limited resolution of the camera (1024 × 64 pixels). A circle surrounds three data which are out of the global tendency. The most likely reason of this dispersion is that the microfluidic channel was not dried properly and there were some liquid remaining at the front meniscus that lubricated the channel and thus accelerated the plug motion leading to reduced rupture time and length.

Indeed, L c o corresponds to the plug initial length when the initial capillary number is equal to the critical capillary number Ca c . From equation ( 15) and by approximating S r by wh at first order, we obtain:

L c o = h 12Ca c ∆P t h σ -(E 2 + 2Df (α))Ca 2/3 c
This formula gives L c o = 3.1 ± 0.2 mm and L c o = 3.5 ± 0.2 mm for ∆P t = 900 Pa and ∆P t = 1000 Pa respectively. It overestimates by 19% and 16% respectively the critical length but nevertheless remains in good agreement with the experimentally measured values. This formula is also consistent with the increase in L c o as a function of ∆P t observed experimentally. This theoretical prediction of the critical initial length L c o is of the upmost practical interest since it enables to predict in which regime will mainly evolve a liquid plug depending on its initial length. An interesting point is also that despite the regime change, the rupture length and rupture time remain relatively proportional to each other (Fig. 5 c andc1).

Comparison with the dynamics in cylindrical tubes

Figure 6 compares the rupture times and lengths obtained in rectangular and cylindrical tubes in the same range of capillary numbers (5.5×10 -5 Ca 1.2× 10 -2 ). This comparison is not meant to be quantitative since the dimensions of the channels (rectangular: h = 45µm, w = 785µm, cylindrical: diameter D = 470µm) and the driving pressure magnitude (1000 Pa for rectangular channels 80 Pa for cylindrical channels ) are different. The purpose of this figure is only to show that in the same range of capillary numbers, no transition is observed in the case of a cylindrical channel, while a transition is clearly evidenced in a rectangular channel.

5. Response of liquid plugs to periodic pressure forcings in rectangular microfluidic channels 5.1. Detailed analysis of single plug ruptures.

We further investigated the response of liquid plugs to cyclic forcing. For this purpose, liquid plugs are inserted at the center of a rectangular microfluidic channel and a cyclic pressure forcing (represented on Fig. 2) is applied. Fig. 7 illustrates the positions of the rear and front menisci (a-c-e) and the evolution of the plug length (b-d-f) for three different initial plug lengths: L 1 = 3.8 mm (a,b), L 2 = 4.1 mm (c,d) and L 3 = 4.5 mm (e,f). The blue curves correspond to experiments and the red curves to simulations. For these three initial lengths, the plugs undergo oscillations eventually leading to their rupture. The experimental results show that the evolution of the plug length is not monotonic: the plug size first increases and then decreases during each back and forth motion. This is a consequence of the progressive increase in the driving pressure (Fig. 2b): at the beginning the driving pressure is low, the plug moves slowly and leaves less liquid behind it than it recovers from the liquid film lying in front of it. Then, when the driving pressure reaches a critical pressure (derived in [START_REF] Magniez | Dynamics of liquid plugs in prewetted capillary tubes: from acceleration and rupture to deceleration and airway obstruction[END_REF] in the case of cylindrical tubes), the tendency is inverted. The number of oscillations before the plug rupture increases with the initial length of the plug. For the longest plug (L 3 = 4.5 mm), a clear transition can be seen between a first phase where the plug undergoes relatively stable oscillations with weak net evolution of its length from one cycle to another (see Fig. 7f before time t = 15 s) and a second phase with a brutal acceleration of the plug rapidly leading to its rupture (t ≥ 15 s) 5.2. Specificity of the cyclic dynamics of liquid plugs in rectangular channels compared to cylindrical channels.

Such transition is not observed in cylindrical tubes wherein the net variation of the plug size is more regular (see Fig. 8f). [START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF] demonstrated that in cylindrical channels, the two sources of the plug instability leading to its rupture are (i) the cyclic diminution of the plug viscous resistance to motion due to the diminution of its length and (ii) a cyclic reduction of the plug interfacial resistance due to the deposition of a liquid film of increasing thickness at each cycle and lubrication effects. A very interesting point is that these two instability sources rely on the amount of liquid deposited on the walls. If the amount of liquid left on the walls behind the liquid plug would remain constant, there would be no cyclic evolution of the plug size and no instability related to lubrication effects. Thus the plug would undergo stable periodic motion with no remarkable evolution of its size and no rupture. This behaviour is indeed observed for plugs of initial length larger than L c o = 4.7 mm (see Fig. 9a). In this case the plug always moves at a dimensionless speed smaller than the critical capillary number Ca c . The simulations are able (i) to quantitatively reproduce the statistical trends of the evolution of the rupture length and rupture time (80% of the experimental data match with the simulations within the error bar on the determination of the initial plug length, see Fig. 5) and (ii) to qualitatively reproduce individual dynamics of liquid plugs (see Fig. 8). Thus, we use them below to analyse the evolution of .

To understand the transition occurring for liquid plugs smaller than L c o we also plotted the evolution of the wet fraction as a function of time for the initial size L 3 = 4.5 mm (see Fig. 10). In this case, the plug moves initially at a capillary number lying under the critical capillary number Ca c thus leading to the deposition of a liquid film of constant thickness behind the plug. The plug speed increases progressively (see Fig. 8b) due to the increase in the driving pressure (see Fig. 2). At time t ≈ 1.5 s (see Fig. 10) the plug dimensionless speed overcomes Ca c and the wet fractions starts increasing until the direction of motion changes (Fig. 10). During the next cycles, the same behaviour is observed with, at first, the deposition of a film of constant thickness and then the deposition of a film of increasing thickness (Fig. 10). Nevertheless, at each cycle, (i) the plug travels further away, (ii) the plug size decreases, (iii) more and more liquid is left on the walls and (iv) the proportion of the motion above Ca c increases. For time t > 15 s the plug dimensionless speed exceeds Ca c for the most part of the motion and this leads to a rapid evolution of the plug size and speed and eventually its rupture. This second phase is similar to the evolution of liquid plugs in cylindrical tubes. This analysis enables to set a criterion on the stability of a liquid plug driven by a pressure periodic cyclic forcing in a dry rectangular microchannel: If the plug dimensionless speed remains below Ca c during the first cycle, then the plug dynamics will remain stable during the next cycles, while if the plug reaches Ca c during this first cycle, it will accelerate cyclically and eventually rupture. 5.3. Evolution of the rupture time and rupture length and comparison between cyclic and unidirectional forcing.

To get a parametric overview of the liquid plugs dynamics in rectangular microchannels, we performed hundreds experiments with different plug initial lengths and either unidirectional or cyclic pressure drivings (of same maximal amplitude). The measured values of the rupture time and rupture length are represented on and compared to the evolutions in cylindrical tubes (a2-b2). As previously reported in cylindrical channels, we observe a saturation of the rupture length when the plug starts undergoing cycles. Nevertheless a major difference with cylindrical tubes is that the rupture time increases to infinity for a finite value of the critical initial length L c o ≈ 4.7 mm while the increase in the rupture time was shown to follow a more "gradual" exponential trend in [START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF]. This is again a consequence of the existence of the quasi-static deposition regime in rectangular channels which does not exist in cylindrical tubes. Of course, in any case, the rupture time is never really infinite owing to evaporation of the plug occurring in the channel and the weak dependence over the capillary number.

Critical assessment of the experimental dispersion and the model validity.

The experiments presented in this paper were extremely sensitive on the contamination of the channels. To avoid pollution, (i) the channels were fabricated tion, even in the most accelerative phase and the approximation of the viscous pressure drop by a Poiseuille law. 3. The 2D (high aspect ratio) approximation for the estimation of the front interface pressure drop.

4. The approximation of the tube prewetting film by the formula 1 -Sf .

With that said, it is difficult to determine the exact origin of discrepancies between experiments and theory. In particular we are unable at the present state to determine whether observed deviations mostly originate from experimental imperfections or approximations in the theory.

Conclusion

In this work, we studied the dynamics of single liquid plugs in rectangular microfluidic channels under unidirectional and cyclic pressure forcings. First we showed that the transition between quasi-static and dynamic film deposition regimes leads to a dramatic acceleration of the plug rapidly leading to its rupture. A pressure-dependent critical size for the transition between these two regimes is derived analytically. For cyclic periodic pressure forcing, we showed that two regimes can occur depending on the initial size of the plug: the plug can either undergo stable periodic oscillations or cyclically accelerate and eventually rupture. The stable regime is observed when the plug dimensionless speed remains below a critical capillary number during the first cycle, while the second is observed as soon as the plug overcomes this value during the first cycle. We were able to quantitatively reproduce the evolution with a reduced dimension model obtained from the combination of previous elements introduced by [START_REF] Baudoin | Airway reopening through catastrophic events in a hierarchical network[END_REF], [START_REF] Magniez | Dynamics of liquid plugs in prewetted capillary tubes: from acceleration and rupture to deceleration and airway obstruction[END_REF] and [START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF] with additional elements to consider the transition between the quasi-static and dynamic film deposition regimes.

These results are of primary interest since microfluidic channels with rectangular cross sections are widely used in the field of microfluidics owing to their easy fabrication. In particular, for the study of liquid plugs dynamics in complex geometries, such as airway tree, it is extremely difficult to design trees with cylindrical sections. Thus this work also enables to analyse and transpose results obtained in rectangular channels to cylindrical channels and understand the pertinence and limit of such comparison.

Acknowledgment

We acknowledge stimulating discussions with J.C. Magniez. This work was supported by Université de Lille. Anderson, J.R., Chiu, D.T., Wu, H., Schueller, O., Whitesides, G.M., 2000. Fab-rication of microfluidic systems in poly (dimethylsiloxane). Electrophoresis 21, 27-40.

Figure 1 .

 1 Figure1. Sketch of the experimental set-up. A liquid plug is created inside a rectangular PDMS microfluidic channel by pushing alternatively some liquid and some air at a Y-junction. The liquid is pushed by a syringe pump while the air is pushed by a MFCS pressure controller. Then the liquid plug is moved by the pressure controller either in one direction (unidirectional driving) or alternatively in both directions (cyclic periodic driving).

Figure 2 .

 2 Figure 2. (a) Unidirectional pressure forcing measured experimentally (blue solid line) and approximated by the formula (red dotted line): ∆Pt = P 0 e -ae -bt with P 0 = 900 Pa, a = 3 and the growth rate b = 10. (b) Cyclic pressure forcing imposed by the MFCS pressure controller measured experimentally (blue solid line) and approximated by the analytical formula ∆Pt = 1200e -6e -3t Pa for t ∈ [0; T ], and ∆Pt = (-1) n (Pc -P d ) for t ∈ [nT ; (n + 1)T ] with

Figure 3 .

 3 Figure 3. Sketch of the cross-section behind the liquid plug (a) when Ca → 0 (quasi-static film deposition regime) and when (b) CaCac (dynamic film deposition regime). Black: liquid, white: air.

  ; Thulasidas et al. (1995); Lozàr et al. (2007); Fries et al. (2008); Han and Shikazono (2009); Han et al. (

Figure 4 .

 4 Figure 4. Dynamics of a liquid plug of initial length L 0 = 3.55 mm pushed in a dry rectangular microfluidic channel with a unidirectional pressure head ∆Pt = 1000e -3e -10t Pa. (a) Stack showing the evolution of the plug below and above the critical capillary number Cac. Liquid (air) appears light (dark) grey. (b) Position of the rear and front meniscus as a function of time. (c) Evolution of the capillary number as a function of time. (d) Evolution of the plug length as a function of time. Blue solid curves correspond to experimental measurements and red dashed curves to simulations.

Figure 5 .

 5 Figure5. Rupture lengths (a-a1) and rupture times (b-b1) of a set of liquid plugs as a function their initial lengths when they move on a dry rectangular microfluidic channel under two driving pressure ∆Pt = 900e -3e -10t Pa (a-b-c) and ∆Pt = 1000e -3e -10t Pa (a1-b1-c1). The rupture lengths as a function of the rupture times are given in figures (c-c1). The blue stars correspond to experiments and the red lines to simulations with the model developed in section 3. Error bars quantify the error in the determination of the plug initial length owing to the large field of view and the limited resolution of the camera (1024 × 64 pixels). A circle surrounds three data which are out of the global tendency. The most likely reason of this dispersion is that the microfluidic channel was not dried properly and there were some liquid remaining at the front meniscus that lubricated the channel and thus accelerated the plug motion leading to reduced rupture time and length.

Figure 6 .

 6 Figure6. Rupture lengths (a1) and rupture times (b1) of a set of liquid plugs moving in a rectangular microfluidic channel under the pressure driving ∆Pt = 1000e -3e -10t Pa compared to the rupture lengths (a2) and rupture times (b2) of a set of liquid plugs moving in a cylindrical capillary tube of diameter D = 470µm under the pressure driving ∆Pt = 80e -6e -3t Pa (a2-b2) . The blue stars correspond to experiments and red curves to simulations. Error bars quantify the error in the determination of the plug initial length owing to the large field of view and the limited resolution of the camera (1024 × 64 pixels).

Figure 7 .

 7 Figure 7. Spatiotemporal evolution of three liquid plugs of initial lengths (a-b) L 1 = 3.85 mm, (c-d) L 2 = 4.1 mm and (e-f) L 3 = 4.5 mm driven by the cyclic pressure forcing represented on Fig. 2(b). (a-c-e) Positions of the left and right menisci as a function of time. (b-d-f) Evolution of the length of the plug as a function of time.

Figure 8 .Figure 10 .

 810 Figure 8. (a-b-c) Spatiotemporal evolution of a liquid plug of initial length L 1 = 4.5 mm moving in a rectangular microfluidic channel under the driving pressure represented on Fig. 2b. (a) Position of the left and right meniscus. (b) Evolution of the capillary number. (c) Evolution of the plug length. (d-e-f) Spatiotemporal evolution of a liquid plug of initial length L 2 = 3.35 mm moving in a cylindrical capillary tube of diameter D = 470µm driven by a cyclic forcing: ∆Pt = 78e -6e -3t Pa for t ∈ [0, T ], ∆Pt = (-1) n (Pc -P d ) for t ∈ [nT, (n + 1)T ]with Pc = 78e -3e -3(t-nT ) Pa and P d = 78e -1.4(t-nT ) e -0.02e -1.4 * (t-nT ) Pa, T = 2.15 s with 2T = 4 s. (d) Position of the left and right meniscus. (e) Evolution of the capillary number. (f) Evolution of the plug length. For all experiments, the blue curves correspond to experiments and the red curves to simulations with the model presented in this paper for the experiments (a-b-c) in a rectangular tube and the model presented in[START_REF] Signe-Mamba | Dynamics of a liquid plug in a capillary tube under cyclic forcing: memory effects and airway reopening[END_REF] for the experiments (d-e-f) in a cylindrical tube.

) rupture times for a cyclic or unidirectional pressure forcing in a rectangular (a1-b1) or cylindrical (a2-b2) channel. For rectangular channels, the pressure driving corresponds to the one represented on Fig. 2 and for the cyclic pressure driving, it correspond to the one described in Fig. 8. Blue stars correspond to experiments with the cyclic forcing, black squares to experiments with the unidirectional pressure forcing, red solid lines to simulations for the cyclic pressure forcing and black dashed lines to simulations with the unidirectional pressure driving.

in the clean room and stored in a sealed box until their use, (ii) the air injected in the channels by the pressure controller was filtered by several filters, and (iii) we were extremely cautious in the liquid sampling to avoid pollution of Perfluorodecalin. Nevertheless, even with all these precautions, we had to change the channels regularly when contamination was observed to avoid deviation in the results. We believe that this pollution is the main source of discrepancy in the experiments. From the theoretical side, the main shortcomings of the model are:

1. The approximation of the progressive transition between the quasi-static and dynamic deposition regimes (observed in [START_REF] Lozàr | Scaling properties of coating flows in rectangular channels[END_REF]) by a sharp transition between a constant value behind a critical capillary number and equation ( 23) above. 2. The omission of all the unsteady and convective terms from Stokes equa-