domain: sde chim

text:
Asma Tab

Mohamed Dahmane

Belabed Chemseddin

Bachir Bellal

Mohamed Trari

Claire Richard

Mohamad Dahmane

Photocatalytic Degradation of Quinoline Yellow over Ag3PO4

Keywords: Ag 3 PO 4, visible light, dye, Photocatalysis, Photoproducts

come

Introduction

Advanced Oxidation Processes (AOPs), especially heterogeneous photocatalysis, is an environmentally friendly technique [START_REF] Yang | Recent advances in application of graphitic carbon nitride-based catalysts for degrading organic contaminants in water through advanced oxidation processes beyond photocatalysis: A critical review[END_REF] that is used for pollutants elimination in water and air [START_REF] Karpuraranjith | Hexagonal SnSe nanoplate supported SnO 2 -CNTs nanoarchitecture for enhanced photocatalytic degradation under visible light driven[END_REF], bacterial inactivation [START_REF] Yu | Enhanced photocatalytic degradation of tetracycline under visible light by using a ternary photocatalyst of Ag 3 PO 4 /AgBr/g-C 3 N 4 with dual Z-scheme heterojunction[END_REF], reduction of heavy metals to less harmful species and water splitting [START_REF] Wang | Earth-abundant transition-metal-based bifunctional catalysts for overall electrochemical water splitting: A review[END_REF]. Photocatalysis initiates oxidation reactions through the generation of holes (h + ) and electrons (e -) pairs. The resulting radicals act as oxidants of organic pollutants and convert them into CO 2 , H 2 O and inorganic salts [START_REF] Mudhoo | functionalization and performance of doped photocatalysts for dye degradation and mineralization: a review[END_REF]. For this reason, extensive investigation was directed toward semiconducting materials which have positive points like the environmental safety and facile preparation. On the other hand, popular photocatalysts such as ZnO [START_REF] Ong | A review of ZnO nanoparticles as solar photocatalysts: Synthesis, mechanisms and applications[END_REF], TiO 2 [START_REF] Guo | Fundamentals of TiO 2 Photocatalysis: Concepts, Mechanisms, and Challenges[END_REF][START_REF] Vernardou | One pot direct hydrothermal growth of photoactive TiO 2 films on glass[END_REF], and SnO 2 are attractive because of their low cost, environmental friendliness [START_REF] Wang | Earth-abundant transition-metal-based bifunctional catalysts for overall electrochemical water splitting: A review[END_REF], thermal stability and bio safe property [START_REF] Chantana | Impact of Urbach energy on open-circuit voltage deficit of thin-film solar cells[END_REF][START_REF] Fatimah | TiO 2 -pillared saponite and photosensitization using a ruthenium complex for photocatalytic enhancement of the photodegradation of bromophenol blue[END_REF]. However, because they possess a wide gap (Eg > 3 eV), their efficacy and sensitivity only occur under UV-light making solar-light application not enough efficient [START_REF] Chen | Convenient fabrication of Ni-doped SnO 2 quantum dots with improved photodegradation performance for Rhodamine B[END_REF][START_REF] Chakraborty | Averrhoe carrambola fruit extract assisted green synthesis of zno nanoparticles for the photodegradation of congo red dye[END_REF]. In this respect, active research is directed to the elaboration of catalysts active under visible light for the wastewaters treatment [START_REF] Mallakpour | Nanocomposites based on biosafe nano ZnO and different polymeric matrixes for antibacterial, optical, thermal and mechanical applications[END_REF][START_REF] Kadhim | Hydrogen gas sensing based on nanocrystalline SnO 2 thin films operating at low temperatures[END_REF][START_REF] Kumar | Experimental Study on Thermal Properties of Nano-TiO 2 Embedded Paraffin (NEP) for Thermal Energy Storage Applications[END_REF].

Recently, Ag 3 PO 4 has gained great attention as a photocatalyst [START_REF] Liu | P3HT-coated Ag 3 PO 4 core-shell structure for enhanced photocatalysis under visible light irradiation[END_REF]. It exhibits a high photocatalytic activity upon visible-light irradiation for organic pollutants decomposition [START_REF] Borjigin | A solar light-induced photo-thermal catalytic decontamination of gaseous benzene by using Ag/Ag 3 PO 4 /CeO 2 heterojunction[END_REF][START_REF] Sedaghati | Integration of BiOI and Ag 3 PO 4 nanoparticles onto oxygen vacancy rich-TiO 2 for efficient visible-light photocatalytic decontaminations[END_REF][START_REF] Boudiaf | Synthesis and characterization of semiconductor CoAl 2 O 4 for optical and dielectric studies: Application to photodegradation of organic pollutants under visible light[END_REF][START_REF] Tab | Visible light assisted photocatalytic degradation and mineralization of Rhodamine B in aqueous solution by Ag 3 PO 4[END_REF] and selective oxidation of alcohols [START_REF] Sedaghati | Integration of BiOI and Ag 3 PO 4 nanoparticles onto oxygen vacancy rich-TiO 2 for efficient visible-light photocatalytic decontaminations[END_REF]. Its photocatalytic properties depend on morphology, and active surface area [START_REF] Belabed | Optical and dielectric properties for the determination of gap states of the polymer semiconductor: Application to photodegradation of organic pollutants[END_REF][START_REF] Liu | A stable Ag 3 PO 4 @PANI core@shell hybrid: Enrichment photocatalytic degradation with π-π conjugation[END_REF][START_REF] Belabed | Photo-electrochemical characterization of polypyrrol: Application to visible light induced hydrogen production[END_REF]. The photocatalytic degradation of anionic dyes by Ag 3 PO 4 was reported to involve the photogenerated holes that show a high oxidation capability in the valence band [START_REF] Chiu | Mechanistic insights into photodegradation of organic dyes using heterostructure photocatalysts[END_REF]. In addition, the existence of PO [START_REF] Wang | Earth-abundant transition-metal-based bifunctional catalysts for overall electrochemical water splitting: A review[END_REF] 3-makes easy the separation of electrons and holes (e -/h + ) pairs through the interfacial junction electric field which generates a dipolar moment. Ag 3 PO 4 synthesized by a one step process has shown that both the free radicals and the photo holes generated contribute to the dyes mineralization [START_REF] Liu | A stable Ag 3 PO 4 @PANI core@shell hybrid: Enrichment photocatalytic degradation with π-π conjugation[END_REF].

For application of photocatalysis process in the environmental protection, it is important to use cheap, easily synthesized or readily available photocatalyst.

Dyes are widely used in food and cosmetic products as well as in the textile industries for their variety of colors and their chemical stability. The production of synthetic dyes exceeds 700 000 tons per year among which 140000 tons are released in the effluents in a non-controlled way [START_REF] Dias | Residuebased iron catalyst for the degradation of textile dye via heterogeneous photo-Fenton[END_REF]. Given their toxicity [START_REF] Cao | Doping of chlorine from a neoprene adhesive enhances degradation efficiency of dyes by structured TiO 2 -coated photocatalytic fabrics[END_REF], these dyes represent an important source of pollution and a serious threat for the aquatic environment [START_REF] Alahiane | Synchronous role of coupled adsorption-photocatalytic degradation of Direct Red 80 with nanocrystalline TiO 2 -coated non-woven fibres materials in a static batch photoreactor[END_REF]. Many dyes show a high stability against temperature, light and biodegradation [START_REF] Ebrahimi | Photocatalytic degradation of organic dyes using WO 3 -doped ZnO nanoparticles fixed on a glass surface in aqueous solution[END_REF] making their elimination a major challenge in the water treatment processes.

Physico-chemical treatments

(coagulation/filtration, coagulation/flocculation, adsorption, etc.) are currently used for the elimination of industrial effluents [START_REF] Cao | Doping of chlorine from a neoprene adhesive enhances degradation efficiency of dyes by structured TiO 2 -coated photocatalytic fabrics[END_REF]. However, such techniques do not degrade contaminants but only transfer them from a liquid phase to the solid phase, and additional treatments are thus required [START_REF] Molavi | Rapid and tunable selective adsorption of dyes using thermally oxidized nanodiamond[END_REF].

Quinoline Yellow (QY) is a synthetic dye and one of the most employed food additives, even though it has also applications in cosmetic and pharmaceuticals preparations [START_REF] Macioszek | The evaluation of the genotoxicity of two commonly used food colors: Quinoline Yellow (E 104) and Brilliant Black BN (E 151)[END_REF]. This dye may cause some diseases like dermatitis and allergic reactions, and has been banned by the standards of the food industry in countries like USA, Australia, Norway and Iran [START_REF] Lau | Synergistic interactions between commonly used food additives in a developmental neurotoxicity test[END_REF][START_REF] Belabed | Photoelectrochemical properties of doped polyaniline: Application to hydrogen photoproduction[END_REF]. Many N-containing dyes like QY undergo natural anaerobic reductive degradation, which can generate aromatic amines suspected to be carcinogenic [START_REF] Cao | Doping of chlorine from a neoprene adhesive enhances degradation efficiency of dyes by structured TiO 2 -coated photocatalytic fabrics[END_REF]. For these reasons, it is necessary to undertake investigations aiming to eliminate it from water resources. Several studies concerning QY heterogeneous photocatalysis were investigated using ZnO [START_REF] Regulska | Photocatalytic degradation of hazardous Food Yellow 13 in TiO 2 and ZnO aqueous and river water suspensions[END_REF], TiO 2 [START_REF] Regulska | Photocatalytic degradation of hazardous Food Yellow 13 in TiO 2 and ZnO aqueous and river water suspensions[END_REF][START_REF] Gupta | Photodegradation of hazardous dye quinoline yellow catalyzed by TiO 2[END_REF] and TiO 2 /polyaniline composite [START_REF] Salem | Photocatalytic degradation of Allura red and Quinoline yellow with Polyaniline/TiO 2 nanocomposite[END_REF]. However, to our knowledge, product studies are scarce.

Here, we investigated the photocatalytic degradation of QY using Ag 3 PO 4 as a photocatalyst. Ag 3 PO 4 was elaborated by a facile method at room temperature, namely the precipitation. The first part of this paper concerns the physical characterization of the material using XRD, SEM, DRS, BET method, and FTIR spectroscopy. It is complementary to the electrochemical properties previously reported [START_REF] Tab | Visible light assisted photocatalytic degradation and mineralization of Rhodamine B in aqueous solution by Ag 3 PO 4[END_REF]. The second part of this paper is focused on the optimization of operating parameters such as pH, catalyst dose and QY initial concentration. Then, photoproducts were characterized by means of ultrafast high-pressure liquid chromatography (UHPLC) coupled to electrospray ionization high-resolution mass spectrometry (ESI-HRMS), while mineralization and the photocatalyst stability were also carried out.

Experimental

Chemicals

Quinoline Yellow was a mixture of the mono-and disulfonic acids (acid yellow 3, Sigma-Aldrich, Fig. 1). The other chemicals were: titanium dioxide Degussa P25 (TiO acetonitrile gradient grade for HPLC ACS (≥ 99.9 %, Carlo Erba). All reagents were used without any purification process. Water was purified using a RIOS 5 reverse osmosis and synergy (Millipore) device (resistivity 18 MΩ.cm, DOC < 0.1 mg.L -1 ).

Synthesis of Ag 3 PO 4

Ag 3 PO 4 was prepared by precipitation as previously described [START_REF] Tab | Visible light assisted photocatalytic degradation and mineralization of Rhodamine B in aqueous solution by Ag 3 PO 4[END_REF]. Briefly, 15 mL of a The synthesized phase was confirmed by X-ray diffraction using Co K line (λ= 1.790300 Å) at a 2θ scan rate of 2 per min. The micrograph of the sample was examined by SEM analysis using a JEOL microscope (Model JSM 6360-LV). The BET surface area and pore volume were calculated from the N 2 adsorption/desorption isotherm at 77K using an ASAP 2010

Micrometrics equipment. The external surface area, micropore area and micropore volume were computed by the t-plot technique. The total pore volume was determined from the liquid N 2 volume at high relative pressure (0.99). The subtraction of the micropore volume from the total volume gave the volume of mesopore. The DRS was measured using a spectrophotometer (Specord 200 Plus). The concentration of Ag before and after recycling in hydrolyzed samples was determined using Atomic Absorption Spectrophotometry (AAS) (Perkin-Elmer A700). A calibration curve was created by analyzing Ag (1, 2, 3, 4 μg.mL -1 ) standard solutions, and the instrument was adjusted to 328.1 nm.

Photocatalytic experiments

Irradiations were conducted in two devices; one was equipped with fluorescent tubes (Philips HPK 15 W) emitting in the near UV region (maximum at 365 nm) while the other with LED (Delleled 5W) emitting within the wavelength range 400-800 nm. Pyrex reactors were connected to a thermostatically controlled water bath. For irradiation under visible 50 mL of suspensions containing Ag 3 PO 4 and QY were added to the reactor of 100 mL of capacity while for UVA irradiation 15 mL of suspensions were placed in the reactor of 30 mL capacity. Solution pH was adjusted using NaOH or HClO 4 . In both systems, the suspensions were kept at 25 °C under constant agitation (600 rpm). Aliquots of 0.5 mL were regularly withdrawn, vigorously centrifuged and filtered through a 0.45 μm filter (Waters, Millipore) before analyses. The degradation yield (%) was evaluated from the relationship:

(

C o being the initial QY concentration and C t the concentration at time (t).

Analytical Analyses

A UV-vis spectrophotometer (Model Specord 200 plus) was used to monitor the degradation of QY (λ max = 434 nm). QY and its photoproducts were separated using a Shimadzu 8040 HPLC system, equipped with a SPD-M30A DAD detector, and a Nucleodur 100-3 C 18 endcapped column (125 mm×4.0 mm, 3 µm particle size). The mobile phase was a mixture of (A) 30 mmol.L -1 sodium acetate and (B) acetonitrile in the proportions of 80% A and 20% B, at a constant flow rate of 0.75 mL.min -1 . The detection wavelength was set at 290 nm and the temperature oven at 30 °C.

UHPLC-ESI-HRMS was performed on an Orbitrap Q-Exactive (Thermo Scientific) mass spectrometer coupled to a Ultimate 3000 RSLC (Thermo Scientific) apparatus. The mass spectrometer with heated electrospray ionization was operated in both positive and negative ESI modes with a spray voltage of 3 kV for both modes. The column was a Waters analytical column (2.1 mm×100 mm, 1.7 μm particle size) and the column oven temperature was regulated at 30 °C. The aqueous solvent (A) was a mixture of formic acid (0.1%) and acetonitrile as organic phase (B). The separation was achieved using a gradient program consisting of 0-7.5 min 5%, 7.5-8.5 min 99% of the mobile phase B. The injection volume was taken at 10 μL while the flow rate was fixed at 0.450 mL.min -1 .

A Shimadzu 5050 TOC analyzer was used to monitor the mineralization of QY. The percentage of mineralization was calculated using the following equation:

(

where OC 0 is the organic carbon present in the solution at time 0 and OC t at time t. OC is obtained by subtracting the inorganic carbon from total carbon.

Results and discussion

Characterization of Ag 3 PO 4

The XRD pattern of Ag 3 PO 4 had narrow peaks, characteristic of a high crystallization (Fig. 2 a). The simulated and the measured diffractograms showed a good fit of 1.6 with R wp and R exp (%) less than 10% (Fig. 2b). The microstructural parameters obtained after the Rietveld refinement were used to determine the dislocation density after elaboration. The coherent domains of diffraction (Dv) had an average value of 178 nm and the corresponding microstrain (ε) a percentage of 2.5×10 -3 . It is well known that the presence of crystallographic defects, such as dislocations, has a strong influence on the properties of materials. The dislocation density (ρ) calculated from equations (3, 4, 5) was ρ= 2.6×10 12 in this work. This value is acceptable and confirm the crystal quality. During the elaboration, the dislocation assisted slip activities occured in the studied material. In addition, the processes used in this study had no effect on the lattice parameters (a= 6.0316 nm in this work), variation of which is only sensitive to the effects of electronic linkages of the chemical species.

ρ=(ρ D ρ ε ) 1/2 (3)
Where ρ D and ρ ε are calculated using eq 4 and 5:

ρ D (4) ρ ε (5)
b the Burger vector, and K is a constant equal to 10 for a Gaussian distribution of the deformations. nature. In addition, the desorption isotherm showed a H3-hysteresis, characteristic of a capillary condensation observed in mesoporous structures [START_REF] Salem | Photocatalytic degradation of Allura red and Quinoline yellow with Polyaniline/TiO 2 nanocomposite[END_REF]. The detailed textual parameters of Ag 3 PO 4 are gathered in Table 1. The specific surface area and total pore volume were equal to 6 m 2 .g -1 and 0.0192 cm 3 .g -1 , respectively. The formation of Ag 3 PO 4 was further supported by the FTIR spectroscopy, using the routine KBr technique in the 400 -4500 cm -1 region. The peaks centered at 545 and 1082 cm - 1 (Fig. 5) were assigned to molecular vibrations of P-O of the phosphate (PO 4 3-) unit [START_REF] Sedaghati | Integration of BiOI and Ag 3 PO 4 nanoparticles onto oxygen vacancy rich-TiO 2 for efficient visible-light photocatalytic decontaminations[END_REF], and the bands at 3401 and 1641 cm -1 were attributed to the stretching vibration of O-H groups of adsorbed water. The H 2 O molecules and OH groups adsorbed on the surface of Ag 3 PO 4 are oxidized into free hydroxyl radicals that achieve the photodegradation of organic pollutants [START_REF] Shen | Synthesis of Z-scheme g-C 3 N 4 /Ag/Ag 3 PO 4 composite for enhanced photocatalytic degradation of phenol and selective oxidation of gaseous isopropanol[END_REF]. The intense peak at 1383 cm -1 corresponded to the nitrate impurity, AgNO 3 having been used for the preparation of Ag 3 PO 4 [START_REF] Du | Corn-silk-templated synthesis of TiO 2 nanotube arrays with Ag 3 PO 4 nanoparticles for efficient oxidation of organic pollutants and pathogenic bacteria under solar light[END_REF].

The knowledge of the gap band (E g ) is of great utility in the solar energy conversion; it is evaluated from DRS by using the known relationship:

(  h) k = Constant × (h -E g ) (6)
α  (cm -1 ) is the optical adsorption coefficient and hν (eV) the energy of the incident photon.

The nature of the optical allowed transition is deduced from the exponent k (= 2, direct) or (0.5, indirect). The gap E g is deduced from the intersection of the plot (αh) 2 with the abscissa-axis. The direct transition at 2.52 eV (Fig. 6) was in conformity with the yellow color of Ag 3 PO 4 that is a promising photocatalyst under visible light [START_REF] Tab | Visible light assisted photocatalytic degradation and mineralization of Rhodamine B in aqueous solution by Ag 3 PO 4[END_REF]. In contrast, when the suspension containing Ag 3 PO 4 and QY was exposed to UVA or visible light, an almost complete degradation was achieved within 40 min and 70 min under visible light and UVA, respectively (Fig. 7). Experiments were then undertaken to optimize parameters.

Parameters optimization

In this set of experiments, Ag ) was studied by varying pH in the range of 4-11. As shown in Fig. 8, the degradation rate of QY was pH-dependent. The highest rate constant (0.094 min -1 ) was observed at pH 7. The rate constant drastically decreased below pH 7 by 65-90% while by 40-50% above (Fig. 8). This highest efficiency at neutral pH, close to the natural environment, is a very positive point in terms of application.

The reaction was also studied at QY concentrations in the range of 5 -40 ppm, at a fixed dose of Ag 3 PO 4 (0.5 g.L -1 ) and at pH 7. Between 5 and 20 ppm, the rate constant laid between 0.065 and 0.072 min -1 , while at 30 and 40 ppm the value dropped to 0.02 and 0.011 min -1 , respectively. At 40 ppm, the absorbance of the QY solution was equal to 3.2 at 434 nm, QY had therefore a very strong filter effect and prevents Ag 3 PO 4 to absorb light.

The Ag 3 PO 4 photocatalyst dose was varied in the range of 0.25 -1 g.L -1 to avoid the photocatalyst excess and ensure a total absorption of efficient incident photons. Fig. 10 shows that the degradation rate increased with the catalyst dose and reached a maximum for 0.5 g.L -

1 . The rate constant k was equal to 8.6 min -1 . An increase of the catalyst dose up to 1.5 g. L -1 led to a reduced photoactivity due to a lower availability of photoactive sites on the Ag 3 PO 4 surface and a weak light utilization [START_REF] Gaya | Heterogeneous photocatalytic degradation of organic contaminants over titanium dioxide: A review of fundamentals, progress and problems[END_REF]. The optimal conditions (Ag 3 PO 4 dose = 0.5 g.L -1 ,

[QY]=20 ppm and pH~7) were kept for experiments under UVA irradiation.

Compared photocatalytic activity Ag 3 PO 4 /TiO 2 Degussa P25

The photocatalytic efficiency Ag 3 PO 4 and TiO 2 Degussa P25 were compared at the dose of 0.5 g.L -1 . Under visible light Ag 3 PO 4 was active while TiO 2 non active and under UVA, Ag 3 PO 4 was 1.6-times more efficient that TiO 2 . This demonstrates the interest of Ag 3 PO 4 as a photocatalyst.

Table 2. QY degradation rate constants (k) using Ag 3 PO 4 (0.5 g.L -1 ) or TiO 2 Degussa P25 (0.5 g.L -1 ) as a photocatalyst Conditions k( min -1 ) R 2 Reference

Product studies

HPLC and UHPLC-ESI-HRMS analyses were conducted for products identification.

QY is sold by Sigma-Aldrich as a mixture of mono and disulfonic acids C 18 H 12 NO 5 SNa - TP7 TP6 ₃ =0.5 g.L -1 and pH 7) (Fig. 13). After 70 min of irradiation under UVA corresponding to the full degradation of QY, only 5% of OC were transformed into CO 2 . Total OC removal was reached after 48 h of irradiation time. This long time to achieve complete mineralization indicates that intermediates were difficult to oxidize.

Catalyst stability and reusability

To assess the photostability and reusability of Ag 3 PO 4 under visible irradiation for the photodegradation of QY, three photodegradation experiments were successively performed.

After each cycle, the catalyst was centrifuged, filtered, thoroughly rinsed with water and ethanol and finally dried at 70°C. After three recycles, 87, 78 and 64% of QY were photodegraded (Fig. 14 a) within 40 min. These results revealed some deactivation effects.

This reduction can be explained by the reduction of Ag + to metallic state Ag in the presence of photogenerated electrons. To quantify the percentage of photoreduced Ag + both XRD and AAS analyses were carried out (Table 4). XRD results after the three cycles are shown in Fig. 14 b. Small peaks of Ag appeared (< 5%) [START_REF] Tab | Visible light assisted photocatalytic degradation and mineralization of Rhodamine B in aqueous solution by Ag 3 PO 4[END_REF][START_REF] Zhang | Carbon quantum dots/Ag 3 PO 4 complex photocatalysts with enhanced photocatalytic activity and stability under visible light[END_REF]. On the other hand, results of the literature showed that when dispersed in water and irradiated under visible light for a long time, Ag 3 PO 4 remained stable (dissolution <1 ppm). Our analyses in atomic absorption spectroscopy confirmed this result (Table 4).

Na 2

 2 HPO 4 solution (3 mM) was mixed dropwise to 30 mL of AgNO 3 solution (3 mM) under magnetic stirring during 3 h at room temperature. The precipitate with a yellow color characteristic of Ag 3 PO 4 , was collected and washed with deionized water and ethanol. The yellow powder was recovered after drying at 70 °C under static vacuum overnight.

Fig. 1 2 . 3 .

 123 Fig. 1 Quinoline Yellow structure using ChemSketch software

Fig. 2 .Fig. 3 .

 23 Fig. 2. XRD spectrum: Identified atomic plane of the Ag 3 PO 4 phase (a). Rietveld refinement of the Ag 3 PO 4 phase (b)

Fig. 4 . 1 )Fig. 5 .Fig. 6 . 4 3. 2 .

 415642 Fig. 4. BET of Ag 3 PO 4 at 77 K

3 PO 4

 34 suspensions were irradiated with visible light. The degradation profile of QY was monitored until 40 min. The QY decrease rate, r, obeyed an apparent pseudo-first order model (r = k × C), therefore the degradation rate constants, k, were obtained from the slope of the plots ln(C/C 0 ) vs irradiation time. The photocatalytic degradation of QY (20 ppm) by Ag 3 PO 4 (0.5 g.L -1

Fig. 7 Fig. 8 .Fig. 9 .Fig. 10 .

 78910 Fig. 7 Loss of QY (20 ppm) in the absence of Ag 3 PO 4 under visible light (▲) or UVA (), in the presence of Ag 3 PO 4 in the dark (▼), upon UVA irradiation () or upon visible

C 18 H

 18 [START_REF] Chen | Convenient fabrication of Ni-doped SnO 2 quantum dots with improved photodegradation performance for Rhodamine B[END_REF] NO 8 S 2 Na 2 . Yet, we only detected two peaks at m/z = 353.0352 in ES -. They correspond to two isomeric forms of the monosulfonic acid (Fig.11a, Fig.11 b). Three samples were analyzed after 44, 70 and 175 min of irradiation under UVA. In the first two samples the QY conversion extent was of 50% and 80% and in the third one QY was fully degraded (Fig.11 c). A typical HPLC chromatogram is shown Fig.11cand Table3gives the accurate and experimental m/z, the chemical formula in the limit of 5 ppm and the possible structure of the 8 detected photoproducts (TPs) by UHPLC-ESI-HRMS. In general, TPs were detected by negative electro spray ionization. TP1 showed m/z = 208.0074 corresponding to C 9 H 7 O 3 NS for the neutral molecule. This TP was likely quinoline sulfonic acid formed by C-C cleavage and its detection allowed us to locate the sulfonic group of QY on the quinoline ring. With m/z = 251.9977, TP2 had the molecular formula of C 10 H 7 O 5 N and corresponded to TP1+CO 2 . It was therefore likely a carboxylated

Fig. 11 2DFig. 12

 1112 Fig. 11 2D HPLC chromatogram (a) and 3D HPLC chromatogram (b) of QY of monosulfonic acid of QY isomers. HPLC chromatogram of photodegraded QY solution (c) ([QY] 0 = 20 ppm, Ag 3 PO 4 dose = 0.5 g.L -1 and pH ∼ 7)

Fig. 13

 13 Fig. 13 Organic carbon evolution during QY photocatalytic degradation in presence of Ag 3 PO 4 under UVA . [QY] 0 = 20 ppm, Ag 3 PO 4 dose = 0.5 g.L -1 and pH ∼ 7)

Fig. 14 Table 4 .

 144 Fig. 14 Reusability of Ag 3 PO 4 photocatalyst. XRD pattern after three cycles (a). Degradation yield after 40 min of irradiation under visible light for three runs (b). ([QY] 0 = 20 ppm, Ag 3 PO 4 dose= 0.5 g.L -1 and pH ∼ 7)

Table 1 .

 1 Textural characteristics of Ag 3 PO 4

		S BET	S ext	S mic	V T	V mes	M mic
	Ag 3 PO 4	6	5.18 0.82	0.0192 0.0188 0.0004

20

Sample

Surface area (m 2 .g -1 ) Pore volume (cm 3 .g -1 )

Table 3 .

 3 UHPLC-HRMS data of QY-Ag 3 PO 4 irradiated under UVA

	TiO 2 Degussa P25 + QY + UVA	0.018	0.998	This work
	Ag 3 PO 4 + QY + UVA	0.029	0.987	This work
	0.099 TiO 2 Degussa P25 + QY + Visible light <1.810 -4 Ag 3 PO 4 + QY + Visible light	0.992	This work This work

-

  and nitrate ions in the samples irradiated for 175 min. The degradation pathways are summarized in Fig.12.

Acknowledgments: A.T. thanks the Franco-Algerien "PROFAS B+" program for having granted her stay in the ICCF. The authors also thank Martin Leremboure (Engineer CNRS) and Frederic Emmenegger (Tech CNRS) for the UHPLC-HRMS analyses.


Te-Sheng Lin
email: tslin@math.nctu.edu.tw

Wei-Fan Hu
email: wfhu@nchu.edu.tw

Chaouqi Misbah
email: chaouqi.misbah@univ-grenoble-alpes.fr

A direct Poisson solver in spherical geometry with an application to diffusiophoretic problems

Keywords: Fast Poisson solver, spherical harmonics expansion, diffusiophoresis, microswimmer 2010 MSC: 65N06, 65N35, 76D07

We propose a simple and efficient class of direct solvers for Poisson equation in finite or infinite domains related to spherical geometry. The solver was developed based on truncated spherical harmonics expansion, where the differential mode equations were solved by second-order finite difference method without handling coordinate singularities. The solver was further extended to study the dynamics of a diffusiophoretic particle suspended in Stokes flow. Numerical experiments suggested that the particle can achieve a self-sustained unidirectional motion at moderate Péclet numbers, whereas the particle motion becomes chaotic in high Péclet number regimes. The statistical analysis illustrates the run-and-tumble-like nature at short times and diffusive nature at long times without any source of noise.

Introduction

The design of artificial microswimmers has become an active research area in the past decade. A prototypical example, Janus particle, named after the two-faced Roman god, uses the asymmetry property of the active surface to generate motions [START_REF] Paxton | Catalytic nanomotors: autonomous movement of striped nanorods[END_REF]. It was shown by Michelin et al. [START_REF] Michelin | Spontaneous autophoretic motion of isotropic particles[END_REF] that, under the framework of axisymmetric assumption, the anisotropic design is not necessary for locomotion, whereas a spontaneous autophoretic directed motion can be activated for an isotropic particle when the Péclet number is above a critical number. In this Note, we numerically examine the dynamics of the diffusiophoretic system in the full three-dimensional model. Interestingly, we observe a transition to chaotic dynamics by increasing the Péclet number further. Such results show evidence that the complicated motions may already be hidden in purely isotropic media due to the intrinsic nonlinearities of the problem.

To investigate the dynamics of diffusiophoretic particles, we first develop a direct numerical solver for Poisson equation in a unit sphere based on spherical harmonics finite difference method. It is known that, in spherical coordinates, appropriate conditions need to be imposed at the coordinate singularities to have the desired accuracy (see, for example, the discussion by Shen [START_REF] Shen | Efficient spectral-galerkin methods iv. spherical geometries[END_REF]). We avoid placing grid points directly on the coordinate singularities by arranging properly the grid points and thus the numerical scheme does not require any pole conditions. A similar trick based on the combination of Fourier series expansions and finite difference approaches was proposed by Lai et al. [START_REF] Lai | A fast spectral/difference method without pole conditions for poisson-type equations in cylindrical and spherical geometries[END_REF]. We then solve Poisson equation in an infinite domain by using Kelvin's inversion which maps the domain into a finite one with the perfect preservation of the Laplace operator. As a result, the proposed fast solver can be straightforwardly applied.

This Note is organized as follows: The description of the solver is presented in Section 2 and in Section 3 we show our findings on the dynamics of a diffusiophoretic particle. Some concluding remarks and future works are given in Section 4.

Fast Poisson solver in spherical geometry

Under the spherical coordinate system (radial distance r, polar angle θ, and azimuthal angle φ), the Poisson equation is written as

∆c(r, θ, φ) = 1 r 2 ∂ ∂r r 2 ∂c ∂r + 1 r 2 sin θ ∂ ∂θ sin θ ∂c ∂θ + 1 r 2 sin 2 θ ∂ 2 c ∂φ 2 = f (r, θ, φ), (1)
where c is an unknown function and f is a given smooth function in which θ ∈ [0, π], φ ∈ [0, 2π). Two types of domains are considered in this Note, namely, the finite spherical domain (r ∈ [0, 1]) and the infinite domain (r ∈ [1, ∞)). For both cases, the boundary condition should be imposed on the spherical surface (r = 1); three different boundary conditions are considered, that is,

Dirichlet (c(1, θ, φ) = c D (θ, φ)), Neumann ( ∂c ∂r (1, θ, φ) = c N (θ, φ)), and Robin type ( ∂c ∂r (1, θ, φ) + αc(1, θ, φ) = c R (θ, φ))
, where c D , c N and c R are given functions and α is a constant. Notice that for the infinite domain case, we impose a vanishing condition c(r → ∞) = 0 to ensure the uniqueness of the solution. It is important to mention that our numerical scheme can treat all types of boundary conditions equally well as we can see from the implementation later.

Finite domain r ∈ [0, 1]

Let us first consider the finite domain case. Due to the spherical coordinate representation (1), the singularities occur not only at the origin (r = 0), but also at the north (θ = 0) and the south poles (θ = π). Therefore, how to handle those singularities properly becomes an important issue in numerical discretizations. Since the considered domain is related to spherical geometry, it is convenient to represent the solution by spherical harmonics expansion

c(r, θ, φ) = p n=0 n m=0 P m n (cos θ) (c mn (r) cos mφ + cmn (r) sin mφ) , (2)
where P m n denotes the associated Legendre polynomial of order m and degree n [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF]. Notice that the solution is truncated up to the spherical harmonics of order p for numerical purpose. The coefficients in the spherical harmonics expansion (2) are given by

c mn (r) = (2n + 1)(n -m)! 2(n + m)! N (m) 2π 0 π 0 c(r, θ, φ)P m n (cos θ) cos mφ sin θ dθdφ, cmn (r) = (2n + 1)(n -m)! 2(n + m)! N (m) 2π 0 π 0 c(r, θ, φ)P m n (cos θ) sin mφ sin θ dθdφ, (3)
where N (0) = 2π, and N (m) = π for m = 0. One can verify that the spherical harmonic functions Y m n (θ, φ) ≡ P m n (cos θ)e imφ are the eigenfunctions of the Laplace-Beltrami operator with

∆ s Y m n = -n(n + 1)Y m n , where ∆ s = 1 sin θ ∂ ∂θ sin θ ∂ ∂θ + 1 sin 2 θ ∂ 2 ∂φ 2 . ( 4
)
Substituting the spherical harmonics expansion (2) into the Poisson equation (1) and using Eq. ( 4), consequently, one can derive a set of one-dimensional ordinary differential equations

where f mn (r) and fmn (r) are obtained in a similar fashion to Eq. ( 3). Accordingly, the associated boundary conditions take the form: Dirichlet c mn (1) = c D mn and cmn (1) = cD mn ; Neumann c mn (1) = c N mn and c mn (1) = cN mn ; or Robin type c mn (1) + αc mn (1) = c R mn and c mn (1) + αc mn (1) = cR mn . The above boundary values can also be found by using Eq. (3).

The discretized variables are defined on (p + 1) nodes along the θ-direction and 2(p + 1) uniformly distributed nodes φ k = k 2π 2(p+1) in the φ-direction, so that θ j = cos -1 (s j ), where s j 's are (p + 1)-point Gauss-Legendre quadrature nodes. Under such a θ-φ grid layout, the numerical quadrature for the integrations in Eq. (3) converges superalgebraically via the combination of Gaussian quadrature (in θ-direction) and trapezoidal rule (in φ-direction). This step can be performed efficiently using the public software SPHEREPACK [START_REF] Adams | Spherepack[END_REF]. There are several advantages of using spherical harmonics expansion: (i) The number of Gauss-Legendre nodes can be taken small since the truncated series ( 2) is spectrally convergent [START_REF] Orszag | Fourier series on spheres[END_REF]. (ii) The usage of Gauss-Legendre nodes naturally avoids placing the points directly at the poles (θ = 0, π) where the singularities of the Poisson equation ( 1) taken place. (iii) The three-dimensional Poisson equation is converted to one-dimensional differential mode equations ( 5) that significantly reduces the numerical complexity of the problem. Furthermore, parallel computations on solving those mode equations can be implemented naturally. Now, the remaining issue is how to handle the coordinate singularity at r = 0 for the mode equations [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF].

To proceed, we choose a grid in the radial direction by setting the uniform grid half mesh away from the origin

r i = (i -1/2)∆r, i = 1, 2, • • • , M, (6)
with the mesh width ∆r = 1/M . We deliberately avoid putting grid point at the origin directly. For succinctness, we will denote the discretized variables by C i ≈ c mn (r i ) and F i ≈ f mn (r i ). By multiplying the term r 2 on both sides of the differential mode equations in Eq. ( 5), and then apply second-order finite difference discretization, we obtain

r 2 i-1/2 ∆r 2 C i-1 + - r 2 i-1/2 + r 2 i+1/2 ∆r 2 -n(n + 1) C i + r 2 i+1/2 ∆r 2 C i+1 = r 2 i F i , i = 1, 2, • • • M. (7)
The resulting matrix formed by the above discretizations has two advantages. One is being symmetric, which is consistent with the self-adjointness of the original differential operator. The other one is being tridiagonal so that the linear system can be solved easily via the Thomas algorithm.

When solving the linear system [START_REF] Orszag | Fourier series on spheres[END_REF], the ghost values C 0 and C M +1 remain to be specified. For i = 1, one can immediately find that the coefficient for C 0 vanishes thanks to the setup of the radial grid ( 6), hence no extra condition at r = 0 is required. The other ghost value C M +1 can be found by the secondorder approximation for Dirichlet

C M +1 +C M 2 = c D mn , Neumann C M +1 -C M ∆r = c N mn , or Robin type boundary condition C M +1 -C M ∆r + α C M +1 +C M 2 = c R mn .
A similar procedure can be straightforwardly applied for solving the coefficients cmn (r) and thus the repeated scheme is omitted here. We emphasize that since the finite difference discretization is adopted here, in the case where the solution has sharp transitions such as boundary layers or interior layers, to achieve the desired accuracy one can easily extend the present scheme to adaptive mesh with the same computational cost.

Infinite domain r ∈ [1, ∞)

Now we turn our attention to tackle the infinite domain problem. Here, we employ Kelvin's inversion [START_REF] Wermer | Potential Theory[END_REF] by introducing the radial coordinate transformation r = 1/r and further define the associated variables by c(r, θ, φ) = 1 r c( 1 r , θ, φ) and f (r, θ, φ) = 1 r f ( 1 r , θ, φ) so that the considered domain will be finite r ∈ (0, 1]. After some careful calculations, the Poisson equation (1) becomes ∆c(r, θ, φ) = 1 r2

∂ ∂ r r2 ∂c ∂ r + 1 r2 sin θ ∂ ∂θ sin θ ∂c ∂θ + 1 r2 sin 2 θ ∂ 2 c ∂φ 2 = 1 r4 f (r, θ, φ). (8)
Remarkably, we see that the Laplace operator is perfectly preserved through this coordinate transformation. Moreover, we should point out that the transformed solution c exists only when the solution admits the behavior c(r → ∞) ≈ r -β with β ≥ 1. As a consequence c is bounded at r = 0 and it is not necessary to impose any additional condition here. Hence, to solve [START_REF] Wermer | Potential Theory[END_REF] we can just follow exactly the same procedure as presented in the previous subsection.

Numerical verification

In this subsection we perform the convergence study of the proposed scheme. For the finite domain case (r ∈ [0, 1]) we solve Poisson equation with exact solution c e (r, θ, φ) = sin(10r sin θ cos φ+π/4) sin(5r sin θ sin φ+ π/3) sin(5r cos θ + π/2) and set the parameter α = 1 in the cases of Robin boundary condition (notice that the chosen solution is significantly oscillatory). For the infinite domain case (r ∈ [1, ∞)) we set c e (r, θ, φ) = r 3 cos θ sin 2 θ cos φ sin φ/ exp(r 2 ) and α = -1. We denote the numerical solutions by c h and show the results in Table 1, where the upper panel gives the mesh refinement results for the finite domain and the lower panel for the infinite domain. In each run we fix 32 Gauss-Legendre nodes for polar angle and 64 nodes for azimuthal angle. Notice that the discretization [START_REF] Orszag | Fourier series on spheres[END_REF], with truncation error O(∆r 2 ), introduces an error c h -c e ∞ with order O(c rr (r = 0)∆r 2 log ∆r)+O(∆r 2 ). As expected, we obtain a roughly secondorder convergence for the finite domain case while a clean second-order of accuracy for the infinite domain case. In addition, these tests show the robustness of our proposed Poisson solver for dealing with oscillatory solutions. On the other hand, the computational cost for the present method is only O(M p 2 (log p) 2 ); the CPU time for the finest mesh case M = 1024 is just within one second on a MacBook Pro (2016), showing the efficiency of the present scheme.

M

Dirichlet

Neumann Robin

c h -c e ∞ Rate c h -c e ∞ Rate c h -c e

Dynamics of a diffusiophoretic particle

We consider the dynamics of a spherical shaped colloidal particle immersed in a Newtonian fluid. The fluid is assumed to be incompressible and obeys the Stokes equations. Surrounding solute is isotropically emitted from the particle surface. The solute interacts with the particle through a short-range potential with a characteristic length that is much less than the particle radius, leading to the so-called sharp interface limit amounting to a slip tangential velocity on the particle surface [START_REF] Anderson | Colloid transport by interfacial forces[END_REF] (see below). The solute concentration is assumed to diffuse and be advected by the surrounding fluid flow.

We use the same scalings to non-dimensionalize the governing equations as in [START_REF] Michelin | Spontaneous autophoretic motion of isotropic particles[END_REF]. The dimensionless governing equations are described in the following: In the co-moving frame attached to the particle center, the fluid flow and solute concentration coupled system is modeled by

-∇p + ∆u = 0, ∇ • u = 0, ∂c ∂t + u • ∇c = 1 Pe ∆c, (9)
where c stands for the solute concentration, the velocity field u = u r e r + u θ e θ + u φ e φ (e r , e θ and e φ stand for the basis vectors in spherical coordinates, u r , u θ and u φ are scalar functions), p is the pressure and Pe is the Péclet number. The particle is suspended in an unbounded fluid domain

Ω = {(r, θ, φ)| r ∈ [1, ∞), θ ∈ [0, π], φ ∈ [0, 2π)}.
The tangential slip velocity on the particle surface (r = 1) is given by u(1, θ, φ) = ∇ s c(1, θ, φ), where ∇ s denotes the surface gradient operator. The fixed-flux boundary condition of surface activity reads ∂c ∂r (1, θ, φ) = -1, where the negative sign denotes the emission mechanism. In the far-field limit, the velocity field converges to the phoretic velocity u(r → ∞) = -(U + Ω × r) and the concentration is assumed to attenuate by c(r → ∞) = 0. The phoretic kinematics of the particle can be arrived by using the reciprocal theorem [START_REF] Stone | Propulsion of microorganisms by surface distortions[END_REF], which allows the relation for the translational velocity U and rotational velocity Ω to the surface slip velocity via the surface integral on the unit sphere S (in Cartesian components):

U = - 1 4π u(r = 1) dS, Ω = - 3 8π e r × u(r = 1) dS. ( 10
)
It can be deduced that the particle motion is always irrotational (Ω = 0) due to the imposed surface slip flow.

It is interesting to see that, a stationary solution where there is no net flow and zero phoretic velocity U = 0 exists at all Péclet numbers with the solute concentration c 0 (r) = 1/r. Whereas at high Péclet numbers, the coupled nonlinear advective term u•∇c would lead to a polarization of the concentration around the particle surface, thus any infinitesimal perturbation to the immobile state will trigger a spontaneous autophoretic swimming motion. Although this symmetry-breaking behavior was reported in [START_REF] Michelin | Spontaneous autophoretic motion of isotropic particles[END_REF] where an axisymmetric constraint is imposed, we need to emphasize that, we have released such a constraint and attempt to explore all possible equilibriums of the entire system in the full three-dimensional simulations.

Since the considered problem is relevant to spherical geometry, it is ideal to use the Lamb's general solution in Stokes flow [START_REF] Pak | Generalized squirming motion of a sphere[END_REF]. By applying the present setup, the velocity field can be analytically expressed through spherical harmonics expansion:

where

ψ 1 (r) = 2 3 (1 -1 r 3 ), ψ n (r) = 1 r n (1 -1 r 2 ) for n ≥ 2, χ 1 (r) = 2 3 (-2 -1 r 3 ), and χ n (r) = 1 r n ( n-2 n -1 r 2 )
for n ≥ 2, and the prime in P m n represents the differentiation. Representing the solute concentration on the surface in Fourier modes along the azimuthal direction by

c(1, θ, φ, t) = ∞ m=0 (D m (θ, t) cos mφ + D m (θ, t) sin mφ) (12)
and using the prescribed slip boundary condition, we can obtain the coefficients in Eq. ( 11) through Eqs. (82, 90-93) in Ref. [START_REF] Pak | Generalized squirming motion of a sphere[END_REF] as

B mn = 1 4 (2n + 1)(n -m)! (n + 1)(n + m)! 1 -1 n(n + 1)D m (µ, t)P m n (µ) dµ, Bmn = 1 4 (2n + 1)(n -m)! (n + 1)(n + m)! 1 -1 n(n + 1) Dm (µ, t)P m n (µ) dµ, (13)
where we have used the notation µ = cos θ. From the above equations ( 10)-(13), we can see that the velocities u and U are given in terms of the concentration field. Substituting Eq. ( 11) into the concentration equation ( 9) yields a closed nonlinear equation for c.

By plugging-in the analytic velocity solution (11) at r = 1 to Eq. ( 10) gives the phoretic translational speed U = 4 3 (B 11 e x + B11 e y -B 01 e z ). One can see that only the first harmonic modes B 11 , B11 , and B 01 are responsible for the swimming velocity in the x-,y-, and z-directions, respectively, and this allows the phoretic particle to swim freely in the three-dimensional domain.

Numerical algorithm

As aforementioned, by using the Lamb's solution for the velocity field, we only need to solve the nonlinear transportation equation for the concentration field. Since the considered domain is unbounded (r ∈ [1, ∞)), we can map the domain into a finite one (r ∈ (0, 1]) by applying Kelvin's inversion with r = 1/r. The new concentration c(r, θ, φ)

= 1 r c( 1 r , θ, φ) fulfills 1 r5 r ∂c ∂t + u(c) • ∇c = 1 Pe ∆c, (14)
where

∇ = -r 2 -r3 ∂ ∂ r , r2 ∂ ∂θ , r2 sin θ ∂ ∂φ
and the Laplace operator ∆ is preserved (see Eq. ( 8)). We employ a first-order implicit-explicit time advancing scheme to discretize Eq. ( 14), that is, we treat the nonlinear advective term explicitly while the diffusive component is handled implicitly:

1 r5 r cn+1 -cn ∆t + u(c n ) • ∇c n = 1 Pe ∆c n+1 , (15)
where the superscript on a variable represents the time level index and ∆t is the time step size. The spatial discretization for all numerical variables is defined as stated in Section 2. For the computation of the velocity u(c n ), firstly we need to calculate D m and Dm in Eq. ( 12) by applying Fast Fourier Transform (in a truncated sense with the highest frequency mode (p + 1)), and then use them to find B mn and Bmn in Eq. ( 13) in which the integrations can be accurately approximated by Gaussian quadrature. Finally, those coefficients B mn and Bmn are substituted into Eq. ( 11) to obtain the velocity (truncated at the highest harmonic mode of order p) with a change of variable from r to r. The spatial derivatives in the gradient term ∇c n is approximated by central difference in the radial direction while a spectral method is used in the polar and azimuthal directions (see [START_REF] Adams | Spherepack[END_REF] for detail). Equation (15) forms a Helmholtz equation for cn+1 which can be efficiently solved by the method proposed in Section. 2.

Numerical simulations

To investigate the spontaneous locomotion of the phoretic particle, we perform the simulations at various Péclet numbers where the concentrations being initially set as c 0 (r) = 1/r (immobile solution) with a small perturbation on the particle surface. For all the present simulations, we set 32 grid points in the radial direction, 32 and 64 points in the polar and azimuthal directions, respectively; the time step size ∆t is set to be of the same order as ∆r to maintain the numerical stability.

Our numerical results reveal that, for small Pe numbers, the system recovers to the immobile solution regardless of the imposed perturbation, whereas when Pe is sufficiently large, the interplay between the advection and diffusion terms triggers a spontaneous symmetry-breaking swimming motion along a given direction (fixed by the initial perturbation) with a constant speed. As an example, the snapshots of the concentration distribution at different times for P e = 10 are demonstrated in Fig. 1. After a long time, a comet-like pattern is formed (see T = 2000) and the concentration field is axisymmetric about the direction of the phoretic velocity U. In addition, it is found that the symmetry-breaking solution occurs in the form of a supercritical bifurcation when Pe > Pe c = 4; the particle sets into a directed motion with the long-time swimming velocity U . We show the velocity as a function of the Péclet number in Fig. 2. The nonlinear variation of the phoretic velocity agrees quantitatively with the ones obtained under axisymmetric assumption in [START_REF] Michelin | Spontaneous autophoretic motion of isotropic particles[END_REF].

Interestingly, by increasing further the Péclet number reveals a new instability: the particle quits the unidirectional swimming motion and enters into an apparently chaotic regime. As an example, see Fig. 3 for the numerical results at Pe = 24.5. It can be seen that the particle exhibits an apparently random motion (see the trajectory in Fig. 3(a)) and we quantify this random walk by measuring the mean square displacement MSD(τ ) =< r(t + τ ) -r(t) 2 >, where r(t) means the location of the particle at time t and the notation < • > computes the average along the entire trajectory. The MSD result is given in Fig. 3(b) as the solid blue line. For short times the particle attains a persistent run-and-tumble-like motion where the MSD is quadratic in time, while at longer times, due to the random turns of swimming directions a de-correlation process takes place, the MSD shows a linear behavior. This demonstration gives evidence that the complicated swimming motion can arise under the minimal version of dynamics such as in this system, a strong nonlinearity induced by the advective component. There are still interesting questions which require further investigation and we leave them as future works.

Conclusions

In this note, a fast direct solver for Poisson equation was developed in spherical geometry with finite and infinite domains. Representing a solution by spherical harmonics expansion, we only need to solve a set of one-dimensional mode differential equations. From the numerical point of view, this step reduces significantly the numerical complexity of the problem. We use finite difference method to discretize those mode differential equations with a proper choice of the radial grid layout that is shifted half mesh away from the origin. In this way, the pole singularity can be naturally avoided. We deal with Poisson equation in an infinite domain by using Kelvin's inversion which maps the domain into a finite one with the perfect preservation of the Laplace operator. As a result, fast solvers can be easily applied. A numerical verification is conducted to illustrate the accuracy and efficiency of the method. The proposed fast Poisson solver is then applied to study the diffusiophoretic system. The emergence of complex dynamics going from straight swimming to chaotic motion for a rigid isotropic particle in a purely phoretic fluid medium has been identified. The present results may serve as a guide for future systematic analyses for diffusiophoretic problems.

  r)(n + 1)P m n (cos θ)(B mn cos mφ + Bmn sin mφ), r) sin θP m n (cos θ)(B mn cos mφ + Bmn sin mφ), mφ -B mn sin mφ),

Figure 1 :Figure 2 :

 12 Figure 1: The snapshots of the concentration distribution on the surface of the particle and on the planes of symmetry at different times for Pe = 10. At T = 2000, the particle attains a symmetry-breaking self-propelled motion along a straight direction (indicated by the yellow arrow) with the constant phoretic speed U = 8.634 × 10 -2 .

Figure 3 :

 3 Figure 3: (a) The erratic particle trajectory at Pe = 24.5 (the red dot marks the origin). (b) The corresponding MSD (solid line) shows a classical random walk as normal diffusion.

Table 1 :

 1 Mesh refinement results with different types of boundary conditions for the finite domain (upper panel) and infinite domain cases (lower panel).

	∞	Rate

Acknowledgement

We would like to thank S. Rafai for the inspiring discussions. T.-S. Lin and W.-F. Hu acknowledge support by Ministry of Science and Technology, Taiwan, under research grant MOST-107-2115-M-009-008-MY2 and MOST-107-2115-M-005-004-MY2, respectively. C. Misbah thanks CNES (Centre National d'Etudes Spatiales) and the French-German university program "Living Fluids" (grant CFDA-Q1-14).


Thanh T Nguyen

Charles Soussen

Jérôme Idier

El-Hadi Djermoune

Exact Recovery Analysis of Non-Negative Orthogonal Greedy Algorithms

It is well-known that Orthogonal Matching Pursuit (OMP) recovers the exact support of K-sparse signals under the condition µ < 1/(2K -1) where µ denotes the mutual coherence of the dictionary. In this communication, we show that under the same condition and if the unknown K-sparse signal is non-negative, the weights of the atoms selected by OMP are nonnegative at any of the first K iterations. Therefore, the generalized version of OMP to the non-negative setting (NNOMP) identifies with OMP, which allows us to establish an exact recovery analysis of NNOMP under the mutual coherence condition. We further establish a similar analysis of the non-negative extension of Orthogonal Least Squares (OLS), and discuss open issues related to the elaboration of weakened guarantees as compared to mutual coherence guarantees.

Introduction

Greedy algorithms are popular iterative schemes for sparse signal reconstruction. Their principle is to sequentially select atoms in a given dictionary and to update the sparse approximation coefficients by solving a least-square problem whenever a new atom is selected. Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] is a well-known greedy algorithm in which at each iteration, the atom having the largest inner product with the current residual is selected. Then, an Unconstrained Least Squares (ULS) problem is solved to update the sparse approximation coefficients.

In the theoretical viewpoint, OMP is guaranteed to exactly recover the support of K-sparse representations in K steps (irrespective of the magnitude of the non-zero coefficients in the K-sparse representation) when the mutual coherence of the dictionary is lower than 1/(2K -1) [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. Moreover, this condition was proved to be sharp [START_REF] Cai | Stable recovery of sparse signals using an oracle inequality[END_REF].

Sign preservation with OMP

Our first contribution states that under the same condition and for K-sparse representations with non-negative weights, the weights of the atoms selected by OMP are non-negative at any iteration.

Theorem 1. [4, Corollary 3.1] Consider a dictionary H whose mutual coherence satisfies µ < 1/(2K -1). Let y = Hx be a K-sparse representation with non-negative weights x. Then, OMP recovers the support of x in K iterations, and at each iteration k = 1, . . . , K, the weights of selected atoms are nonnegative.

K-step analysis of NNOMP

Non-Negative OMP (NNOMP) was first introduced by Bruckstein et al. as a generalization of OMP to address sparse reconstruction under non-negativity constraints [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equation[END_REF]. Similar to OMP, the empty support is used for initialization of NNOMP. At each iteration, NNOMP selects the atom having the most positive inner product with the current residual, the atoms yielding negative inner products being ignored. Then, the non-negative weights are updated by solving a Non-Negative Least Squares (NNLS) problem. It is noticeable that unlike ULS problems, NNLS problems do not have a closed-form solution. Therefore, a subroutine must be called at each iteration of NNOMP. Note also that since NNLS are constrained least-square problems, some non-negativity constraints may be activated, corresponding to the cancellation of some weights in the sparse approximation. Hence, more than K iterations might be necessary to reach a solution that is exactly K-sparse. The reader is referred to [START_REF] Nguyen | Nonnegative orthogonal greedy algorithms[END_REF] for further developments about fully recursive implementations of NNOMP.

On the theoretical side, Bruckstein et al. conjectured that NNOMP is guaranteed to exactly recover the support of nonnegative K-sparse representations in K steps when µ < 1/(2K -1) [5, Theorem 3], but the rigorous proof was not provided. In [START_REF] Nguyen | Sign preservation analysis of orthogonal greedy algorithms[END_REF], we have proved this result as a corollary of Theorem 1. In a nutshell, Theorem 1 states that the ULS solution related to the support found at each iteration of OMP is non-negative, which implies that the NNLS solution related to the same support coincides with this ULS solution. Since the NNOMP selection rule is close to that of OMP, we could thus deduce that the NNOMP and OMP iterates coincide, which yields the following K-step recovery result. Theorem 2. [4, Corollary 3.2] Consider a dictionary H whose mutual coherence satisfies µ < 1/(2K -1). Let y = Hx be a K-sparse representation with non-negative weights x. Then, the NNOMP iterates identify with those of OMP. Thus, NNOMP recovers the support of x in K iterations.

Details of these results can be found in our technical report [START_REF] Nguyen | Sign preservation analysis of orthogonal greedy algorithms[END_REF]. Note that in [START_REF] Nguyen | Sign preservation analysis of orthogonal greedy algorithms[END_REF], we also carry out unified analysis of the non-negative extensions of Orthogonal Least Squares (OLS) [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF] proposed in [START_REF] Yaghoobi | Fast non-negative orthogonal least squares[END_REF]. Open issues related to the elaboration of weakened theoretical guarantees as compared to mutual coherence guarantees will be addressed at the iTWIST workshop.

Acknowledgment

This work was supported by the project ANR BECOSE (ANR-15-CE23-0021).


Alexandre Canitano
email: canitano@earth.sinica.edu.tw

Maxime Godano

Marion Y Thomas

Inherited state of stress as a key factor controlling slip and slip mode: inference from the study of a slow slip event in the Longitudinal Valley, Taiwan

Using borehole strainmeters, we detected a 13-day long slow slip event in 2011 with M w = 5.45 on the Longitudinal Valley Fault, Taiwan. It has been likely promoted by the significant Coulomb stress changes (∼ 0.5-1 MPa) imparted by a combination of coseismic and postseismic slip of the 2003 Chengkung earthquake. Using kinematic slip models, we infer that the slow event has accomodated about 15 to 35% of the slip deficit accumulated in 2011 in its source region. Further, we find that inherited state of stress may have controlled the slow event southward propagation. We also highlight a spatiotemporal correlation between the slow event and a cluster of repeating microearthquakes, suggesting a complex interplay between seismic and aseismic processes on the fault, where transient aseismic regions are adjacent to highly coupled zones.

D R A F T

November 6, 2020, 11:55am D R A F T

Introduction

Over the last two decades, the growing development of geodetic and seismological monitoring arrays in active regions has revealed episodic aseismic slip in the crust, spanning timescales from seconds to years [e.g., [START_REF] Peng | An integrated perspective of the continuum between earthquakes and slow-slip phenomena[END_REF]. These slow slip events (SSEs) play an important role in redistributing stress in the Earth's crust [e.g., Linde et al., 1996] and are now observed in various tectonic regions [e.g., [START_REF] Bürgmann | The geophysics, geology and mechanics of slow fault slip[END_REF].

SSEs are often accompanied by earthquake swarms [Vallée et al., 2013;[START_REF] Gualandi | Aseismic deformation associated with an earthquake swarm in the northern Apennines (Italy)[END_REF]Fasola et al., 2019] or nonvolcanic tremors [e.g., [START_REF] Beroza | Slow earthquakes and nonvolcanic tremors[END_REF], and together, represent an important mechanism of strain release in active regions. Therefore, investigating the stress conditions, the faulting mechanisms of slow slip and what role they play in the earthquake cycle is fundamental to determine time-dependent earthquake hazard.

In eastern Taiwan, the Longitudinal Valley (LV) is an active collision boundary between the Eurasian and Philippine Sea plates [START_REF] Barrier | Active collision in eastern Taiwan: the Coastal Range[END_REF], and accounts for more than half of the 9 cm.yr -1 of oblique plate convergence [START_REF] Yu | Velocity field of GPS stations in the Taiwan area[END_REF]. The Longitudinal Valley Fault (LVF), which runs along the eastern side of the LV, represents the major active structure in the region and accounts for about 4.5 cm.yr -1 of total plate convergence [Thomas et al., 2014a] (Figure 1). The fault is creeping at the surface at the rate of 1-6 cm/yr between latitudes 23 • 00 ′ and 23 • 30 ′ [Thomas et al., 2014a] and also experiences seasonal and transient creep episodes [START_REF] Lee | Active fault creep variations at Chihshang, Taiwan, revealed by creep meter monitoring, 1998-2001[END_REF][START_REF] Murase | An episodic creepslip event detected by precise levelling surveys in the central part of the Longitudinal Valley Fault, eastern Taiwan, in 2011-2012[END_REF].

Despite significant historical earthquakes, there is a paucity of large shocks along the fault relatively to the high convergence rate, suggesting that a significant fraction of the D R A F T November 6, 2020, 11:55am D R A F T long-term slip rate, in the seismogenic depth range, is released aseismically [START_REF] Liu | Slow earthquakes triggered by typhoons[END_REF]. Indeed, based on the analysis of geodetic data for the 1992-2010 period, Thomas et al. [2014a] demonstrate that 80 to 90% of the long-term slip budget on the southern section of the LVF is the result of aseismic slip. Following the 2003 M w 6.8 Chengkung earthquake, a 7-year long afterslip has been detected by Global Positioning System (GPS) stations along the Chihshang Fault (CF) [Thomas et al., 2014a], a 30-km long section of the southern LVF. Borehole strainmeters have captured a very shallow SSE (2 to 4 km depth) with geodetic moment magnitude M w ∼ 4.5 in central LV [Canitano et al., 2019].

A 1-month long afterslip following a M w 4.6 earthquake on the CF was also shown to control the rate of aftershocks near the earthquake source region [Canitano et al., 2018a].

However, the dearth of aseismic transient observations in the region strongly limits our ability to investigate the mechanisms of slow deformation along the fault and to further interpret the interplay between seismic and aseismic processes.

In this study, we document a 2-week long SSE with an equivalent moment magnitude of Sacks-Evertson [START_REF] Sacks | Sacks-Evertson strainmeter, its installation in Japan and some preliminary results concerning strain steps[END_REF] borehole strainmeters along the LVF (Figure 1). They monitor rock volume change (dilatation ǫ v ) and complement the GPS measurements for detecting crustal transients at short to intermediate periods (minutes to weeks). We calibrate the dilatometers using solid-Earth and ocean tides [Canitano et al., 2018b] (Figure S1), and correct the strain data for borehole relaxation, for solid-Earth and ocean tidal strain and air pressure induced strain. 2. We process the GPS data (Figure 1b) with the GAMIT10.42/GLOBK5.16 software packages [START_REF] Herring | Introduction to GAMIT/GLOBK release 10[END_REF] using the 2005 International Terrestrial Reference Frame (ITRF2005) [START_REF] Altamini | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF] coordinates in GLOBK processing (Figure S2).

D R A F T November 6, 2020, 11:55am D R A F T

Detection and characterization of the SSE

Volumetric strain response to rainfall represents the largest signal observed in strainmeter time-series. Hours to days following the rainfall, a contraction signal proportional to the amount of cumulative rainfall is recorded (Figure S3a). A rainfall-strain admittance of -0.51 nǫ/mm is found for the central LV [START_REF] Hsu | Revisiting borehole strain, typhoons, and slow earthquakes using quantitative estimates of precipitation-induced strain changes[END_REF]. At longer periods (weeks to months), an alternance of phases of contraction and expansion is observed, which suggests that volumetric crustal strain is likely modulated by the cumulative rainfall loading (Fig-

ure S3b). Such deformation is complex, highly site-dependent (e.g., local topography, rock elastic/poroelastic properties) [START_REF] Mouyen | Typhoon-induced ground deformation[END_REF], and its mechanisms remain poorly constrained in the region.

Starting on 29 January 2011 (day 29), we observe an expansion of +50 nanostrain (nǫ) at SSTB station, synchronous with a contraction of -39 nǫ at ZANB station and a moderate contraction (∼ -10 nǫ) at HGSB station (Figure 2a). FBRB do not record any relevant signal (SSNB and CHMB experience a power outage). On day 28.2, a light rainfall episode (12 mm of rainfall within 12 hours) (Figures 2b andS4) occurs in the CR (Figure 1b, station C0T9M0 near ZANB) inducing a contraction of -7 nǫ at ZANB, which is about the value expected based on the rainfall-strain admittance in the region [START_REF] Hsu | Revisiting borehole strain, typhoons, and slow earthquakes using quantitative estimates of precipitation-induced strain changes[END_REF]. Rainfall has likely stopped when SSTB signal initiates its large expansion (Figure S4a), that is unlikely to be the strain response to rainfall. Besides, to produce the contraction detected by ZANB and HGSB (-39 nǫ and -10 nǫ, respectively), it would require continuous rainfall from days 29 to 42 with an amount of ∼ 80 mm and of 20 mm

D R A F T November 6, 2020, 11:55am D R A F T X -7
of water, respectively, which is not observed (Figure S4b). Therefore, strain variations are not associated with precipitation and neither are they induced by hydrology as no transient change in groundwater level is detected during January-February 2011 (Figure 2b).

There are also no appreciable sea level changes, the largest variations occurring during day 34 (∼ 0.1 m) and remain undetected by near-coastal stations (Figure 2c). Finally, no signal is detected by the GPS stations (Figure S2).

Strain signals with nearly similar temporal evolution are detected from days 29 to 41.8 in two networks distant by 35 km (Figure 2d). They are unrelated to environmental changes and represent unambiguous evidence for a slow slip episode. Signals follow a logarithmic time-dependence, a pattern observed on borehole strainmeters during a previous SSE [Canitano et al., 2019]. Additionally, following day 41.8, ZANB records a moderate expansion (+20 nǫ) during 3 days (Figure 2a), while a contraction of ∼ -5 nǫ is expected due to 10 mm of cumulative rainfall (Figure S4, station C0T9M0). HGSB records a contraction of -10 nǫ that is about 5 times larger than the expected changes due to 3-4 mm of rainfall (Figure S4, station C0Z070), while FBRB do not record any variation. Since these variations are unrelated to environmental perturbations or to sensor drift or relaxation, they may be the signature of a propagating source [e.g., [START_REF] D R A F T Linde | A slow earthquake sequence on the San Andreas fault[END_REF]. However, it remains difficult to characterize this possible second slow phase because of the lack of data (no signal for SSTB), which limits our ability to fully decipher the complexity of this event. Consequently, we only characterize the 13-day long main phase and thus estimate D R A F T November 6, 2020, 11:55am D R A F T a minimum moment magnitude for the sequence.

We search for a source compatible with the LVF (strike = 23 • NE) and a geologic rake of 70 • , which corresponds to the mean slip vector direction in the Yuli-Fuli region [START_REF] Peyret | Present-day interseismic surface deformation along the Longitudinal Valley, eastern Taiwan, from a PS-InSAR analysis of the ERS satellite archives[END_REF]. We estimate the optimal source location and magnitude for a rectangular fault with uniform slip embedded in an elastic homogeneous half-space [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF] with a rigidity G = 30 GPa using a grid search approach derived from GPS source modeling [START_REF] Lin | Quick determination of earthquake source parameters from GPS measurements: a study of suitability for Taiwan[END_REF] (Section S1 and Figure S5). A minimum misfit is achieved (rms < 1 nǫ) for a source with length of L = 16-24 km, width of W = 3.8-12 km and aseismic slip of U = 2-10 cm at the depth of 10.5-11.5 km. The fault surface is constrained in the range 65-288 km 2 , the source aspect ratio (L/W ) is ranging from 1.8 to 4.6 and the geodetic moment magnitude M w = 5.41-5.45 (Figure S6a). Our preferred source model (rms misfit = 0.36 nǫ) is located at 10.75 km depth, it has length and width of 19 km x 6.4 km, respectively, and a total displacement of 5 cm (M w = 5.44) (Figure S6b). The SSE geodetic moment magnitude is strongly constrained, it represents a typical value for a 2-week long SSE [START_REF] Bürgmann | The geophysics, geology and mechanics of slow fault slip[END_REF][START_REF] Michel | Similar scaling laws for earthquakes and Cascadia slow-slip events[END_REF] and agrees with an earthquakelike cubic moment-duration scaling, as reported in Cascadia [START_REF] Michel | Similar scaling laws for earthquakes and Cascadia slow-slip events[END_REF]Dal Zilio et al., 2020] and Mexico [Frank and Brodsky, 2019]. We compute the static stress drop ∆σ following [START_REF] Kanamori | Theoretical basis of some empirical relation in seismology[END_REF]:

∆σ = 8 3π GU √ LW (1) D R A F T November 6, 2020, 11:55am D R A F T X -9
The average stress drop ranges from 0.03 MPa to 0.32 MPa with ∆σ = 0.115 MPa for our best source model, which is consistent with values from [START_REF] Gao | Scaling relationships of source parameters for slow slip events[END_REF].

We search for recurrent events for the 2004-2019 period using a geodetic template matching [e.g., [START_REF] Rousset | A geodetic matched filter search for slow slip with application to the Mexico subduction zone[END_REF] but found no additional events (Section S2 and Figure S7).

At the minimum, the recurrence time of such event is about 7-8 years. Our search is however limited by the strain templates allowing to detect only events with nearly similar duration and location as the 2011 SSE. We cannot exclude either that our method failed to detect recurrent events with lower magnitude. Indeed, such signals are below the GPS ambient noise level and frequent rainfall strongly impact strainmeter records, potentially concealing transient signals.

Analyze of the seismicity during the SSE episode

We analyze the seismicity from days 28 to 43 in 2011 but find no evidence for the occurrence of nonvolcanic tremors or for a temporary increase in seismicity (Figure S8).

We also search for repeating earthquakes (REs), as a possible seismological signature of aseismic slip [Uchida and Bürgmann, 2019]. They represent nearly identical asperities that rupture repeatedly under the influence of aseismic slip [e.g., [START_REF] Beeler | A simple stick-slip and creep slip model for repeating earthquakes and its implication for micro-earthquakes at Parkfield[END_REF].

We calculate the cross-correlation coefficients (ccc) of the vertical velocity signals for all earthquake pairs (13 events) recorded by 10 broadband seismic stations distant from 5 to 65 km from the SSE source (Figure 1a). We filter the 100-Hz sampling signals between includes P -and S -wave arrivals. Window lengths are ranging from 4 s to 13 s between nearby and remote stations. We find high waveform coherence for 3 events (ccc > 0.92) (Figure 3a). They occur during a 40-min period on 5 February 2011 (day 36) and have local magnitude M L < 2 (events 1-3, Table S2).

To further characterize the earthquake multiplet, we search for additional events for the 2003-2019 period using waveform template matching. We use event 1 as a template and perform sliding-window cross correlations [START_REF] Yang | Fault-plane determination of the 18 April 2008 Mount Carmel, Illinois, earthquake by detecting and relocating aftershocks[END_REF] for the same previous stations and same P -S phase time-windows. A new event is added to the multiplet if its crosscorrelation coefficient with event 1 is ≥ 0.80, at minimum, for three stations. We find a total of 11 new events (Table S2 and Figure 3b) with M L = 1.43-2.90. The multiplet is relocated using the double-difference algorithm HypoDD [Waldhauser and [START_REF] D R A F T Waldhauser | A double-difference earthquake location algorithm: method and application to the northern Hayward fault[END_REF] with manually picked P -and S -wave absolute arrival times and relative cross-correlation delay times. The relative times are determined by cross-correlating all earthquake pairs in the multiplet for both P and S waves using 2-s length windows beginning 0.1 s before manual wave picks. All pairs satisfying ccc ≥ 0.80 are retained for the inversions (85 pairs). To ensure a good constraint on the absolute location of the multiplet, we jointly relocate it with 100 earthquakes selected near and within the SSE source. For these earthquakes, only manually picked P -and S -arrival times are inverted.

As a result, 10 events of the multiplet (M L = 1.43-2.76) are relocated at the depth of about 14.2 km along a structure dipping about 55 • southeastward (Figure S9), which is S10). The source radius R of each relocated event is estimated based on the relationship between stress drop ∆σ and seismic moment M 0 [Eshelby, 1957]. Local magnitude M L is converted to moment magnitude M w using an empirical formula derived for Taiwan [START_REF] Huang | The determination of M L -M w in Taiwan[END_REF][START_REF] Wu | Revised M L determination for crustal earthquakes in Taiwan[END_REF] and [START_REF] Kanamori | Theoretical basis of some empirical relation in seismology[END_REF] (Section S3). Source radii range between 20 and 90 m (30-40 m for the 3 events in 2011) with uncertainties varying between -40 and +250%. Then, using a stochastic approach, we find a maximum of 4 cumulative ruptures inside the cluster, suggesting limited source overlapping and an asperity that likely ruptured by small parts (Figure 3c).

D R A F
M 0 is estimated from M w following

Influence of static stress changes from the 2003 Chengkung earthquake

The SSE source is adjacent to the rupture area of the 2003 M w 6.8 Chengkung earthquake. The mainshock was followed by an afterslip with an aseismic moment equivalent is defined as:

δCF F = δτ + µ ′ δσ n (2)
where µ ′ is the effective friction coefficient, here taken as 0. We propose that the SSE occurrence has been promoted by the significant Coulomb stress changes (δCF F > 0.5-1 MPa, about 5-10 times greater than the SSE average stress drop)

imparted by a combination of coseismic and postseismic slip from the earthquake, and which persists several years after the mainshock. δCF F distribution onto the SSE fault plane reveals a complex pattern (Figures 4 andS11). First, the preseismic stress (1997 to 2003) is essentially negative, with a minimum level similar to the positive coseismic δCF F . Therefore, it is possible that coseismic Couloumb stress changes could not overwhelm Coulomb pre-stress [e.g., [START_REF] Mildon | Coulomb pre-stress and fault bends are ignored yet vital factors for earthquake triggering and hazard[END_REF], which may explain the lack of detected recurrent SSEs right after the mainshock, despite the favourable stress conditions.

Conversely, we simply missed an event since no strain data are available within 6 months following the mainshock. Then, the positive coseismic δCF F on the southern part of the SSE area may have induced the postseismic slip (0.4-0.6 m) which in turn generated significant Coulomb stress changes (δCF F > 0.5 MPa) in the northern section. The latter could be the trigger of the SSE. However, δCF F only gives the relative changes of stress for a defined time interval, here the 7-yr preseismic, coseismic and 7-yr postseismic periods.

We now analyze the location of the SSE with respect to the total slip budget using

Thomas et al.

[2014a] kinematic study. We plot the along-dip variations of slip through time, at the SSE latitude (Figure 4e, area A4) and compare the accumulated slip over ∼ 14 years with the total motion of the fault if it had crept at the long term slip rate. We observe a slip deficit between 0 and 16 km depth, which encompasses the SSE source area.

Analyzing the along-strike variation of slip within area A4, we see that for the southern-D R A F T November 6, 2020, 11:55am D R A F T most section (column c12, Figure 4f), there is almost no slip deficit unlike for the columns c13 to c18 (Figures S13). seismic ruptures can partially or completely penetrate it [e.g., [START_REF] Lin | Overlapping regions of coseismic and transient slow slip on the Hawaiian décollement[END_REF]. Finally, the SSE region also lies in an area of moderate afterslip, suggesting that sections of the LVF experiencing afterslip can also host SSEs [START_REF] Yarai | Quasi-periodic slow slip events in the afterslip area of the 1996 Hyuga-nada earthquakes, Japan[END_REF][START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF].

In addition to triggering aseismic slip, the Chengkung earthquake had also impacted the dynamics of repeater sequences along the LVF, notably halving their recurrence interval [Chen et al., 2020]. During the 2011 SSE, we observe burst-type of nearly identical REs which are events that have extremely short recurrence intervals [START_REF] Igarashi | Repeating earthquakes and interplate aseismic slip in the northeastern Japan subduction zone[END_REF]Templeton et al., 2008]. However, 8 events (out of 10) of the cluster have M L < 2 for which distinguishing REs requires high-frequency signal contents and may remain poorly resolved in the seismic waveforms [Uchida and Bürgmann, 2019]. Magnitude limitations as well as uncertainties on the source stress drops and relocations strongly limit our ability to resolve the source overlapping and to assess the source mechanisms of the REs. We therefore propose two possible scenarii based on the previous analysis:

i : Since burst-type REs tend to occur off the major fault planes [START_REF] Igarashi | Repeating earthquakes and interplate aseismic slip in the northeastern Japan subduction zone[END_REF]Templeton et al., 2008] and Coulomb stress changes caused by the SSE in the RE region is negative (δCF F ∼ -0.1 MPa) (Figure S14), the REs may have occurred spontaneously with no mechanical link with the SSE. Besides, no SSE is detected during other RE occurrences, and particularly during other burst-type REs (Figure S15), suggesting that aseismic slip is unlikely to be the dominant mechanism of REs in the region [Chen et al., 2009,2020] (Figure 3d).

D R A F T November 6, 2020, 11:55am D R A F T
ii : Conversely, since the fault plane solution of the REs is unknown, our Coulomb stress model, which consider the LVF as receiver fault, may not be adequate. Moreover, given the uncertainties on the RE relocations (Section 4) and that our best SSE source locations are within a range of a kilometer (Section S1), aseismic and seismic sources can be spatially colocated. This clustering may suggest a dominant control of the SSE on the repeating microearthquakes. However, given that source overlapping is not observed, the 40-min RE burst is unlikely to be the result of recurrent seismic ruptures to accommodate aseismic slip in the surrounding area [e.g., [START_REF] Beeler | A simple stick-slip and creep slip model for repeating earthquakes and its implication for micro-earthquakes at Parkfield[END_REF]. It rather suggests a triggering cascade of ruptures [e.g., [START_REF] Lengliné | Inferring the coseismic and postseismic stress changes caused by the 2004 M w =6 Parkfield earthquake from variations of recurrence times of microearthquakes[END_REF], in which event 1, likely triggered by aseismic slip, has then ruptured close-by, but distinct, asperities (as events 2 and 3), together partially rupturing the main asperity (event 54). Overall, burst-type REs are likely caused by local, transient stress concentration from nearby earthquakes and/or aseismic slip events [Templeton et al., 2008] and they are also observed during postseismic slip on the LVF [Canitano et al., 2018a].

To conclude, while a SSE was detected at the transition between the creeping section of the LVF and the locked zone [Canitano et al., 2019], we now show that the fault can also

5. 45 .

 45 It occurred in January-February 2011 on the central section of the LVF, northeast of the source rupture of the 2003 Chengkung earthquake. This event occurred between about 8 to 14 km depth and was detected by borehole strainmeters deployed in two networks distant by 35 km. Using static Coulomb stress modeling, we investigate a possible contribution of coseismic and postseismic slip of the Chengkung event to the SSE occurrence. We also analyze the spatiotemporal pattern of a cluster of earthquake multiplets occurring during the SSE episode and investigate its relationship with aseismic slip. Fi-D R A F T November 6, 2020, 11:55am D R A F T X -5 nally, we investigate the behaviour of the SSE region during the Chengkung earthquake and discuss the SSE contribution to the release of slip deficit along the LVF. 2. Near-fault instrumentation and data processing Beginning in 2003, to shed additional light on the nature of the deformation in the LV, the Institute of Earth Sciences (IES) Academia Sinica, in cooperation with the Department of Terrestrial Magnetism, Carnegie Institution of Washington, has deployed 11

3

  and 20 Hz to suppress microseismic noise and select a cross-correlation window which D R A F T November 6, 2020, 11:55am D R A F T

  to 0.8 times the seismic moment at the end of the study period (end of November 2010)[Thomas et al., 2014a[Thomas et al., , 2017]]. The SSE occurred two months after the end of the study D R A F T November 6, 2020, 11:55am D R A F T period while some parts of the fault were still creeping at higher rate than during the preseismic period [Thomas et al., 2014a, 2017]. Given that SSEs are highly sensitive to small stress perturbations (a few kPa) [e.g., Hawthorne and Rubin, 2010], we estimate the contribution of coseismic and postseismic slip to the SSE occurrence using the Coulomb failure criterion. If δσ n and δτ represent respectively the changes in normal stress (tensile stress is positive) and shear stress on the fault plane (positive in the direction of the longterm fault slip) induced by the cumulative slip, the static Coulomb stress change δCF F

  4. The fault moves closer to failure if δCF F > 0 and away from failure if δCF F < 0. δCF F is resolved onto the LVF plane using Coulomb 3.3[Toda et al., 2005;[START_REF] Lin | Stress triggering in thrust and subduction earthquakes and stress interaction between the southern San Andreas and nearby thrust and strike-slip faults[END_REF] for an assigned rake, we further argue for.Unlike the simple geometry we use to compute the SSE source (justified by the strain data resolution), the fault geometry in the kinematic study followed the local variation of strike of the surface fault trace (to be consistent with the high InSAR spatial resolution) as well as the variation of the dip angle, highlighted by the microseismicity[Thomas et al., 2014a]. Thus, even the projection onto the fault of a smooth displacement field leads to a corresponding rake that varies from one patch to another. It can therefore be misleading to compute the δCF F on each individual patch for an imposed rake corresponding to that D R A F T November 6, 2020, 11:55am D R A F T of the SSE. We consequently plot the cumulative displacement as well as the δCF F along the direction of long-term displacement of the LVF. The latter is given by the projection onto the fault of the long-term displacement of the Central Range relative to the Coastal Range[Thomas et al., 2014a]. This further allows to discuss if the SSE accommodates some slip deficit. For the sake of the discussion we also compute the δCF F for a patch corresponding to the inferred SSE source (70 • ).The time-evolution of postseismic slip is modeled using 7 years of geodetic data following the mainshock[Thomas et al., 2014a]. Thus, modeled δCF F represents a very close estimate of the static stress changes on the fault at the SSE onset timing but predicts the cumulative stress as if it had occurred instantaneously. It does not account for the complex interactions between fault plane patches that occurred over the postseismic period.The SSE source is located in an area of very low coseismic slip, ≤ 0.1 m (Figure4a), which results in positive Coulomb stress changes of 0.5-0.7 MPa, essentially in the southern SSE region (Figure4c). Moderate cumulative postseismic slip (about 0.2-0.6 m during the 7-yr period) (Figure4b) is observed in the source region with significant positive (δCF F ∼ 0.5-0.7 MPa) and negative (δCF F ∼ -0.5 MPa) Coulomb stress changes in the northern and southern SSE area, respectively (Figure4d). Cumulative slip estimated for the preseismic period (01/1997-12/2003) is about 20-30 cm in the SSE source region (Figures4e and S11c). It is however more subject to caution given that data are more sparse and as a direct consequence the applied smoothing is stronger (see Thomas et al. [2014a] for details), which inevitably impacts the δCF F computation. However, δCF F is essentially D R A F T November 6, 2020, 11:55am D R A F T null and negative in the southern and central sections (δCF F ∼ -0.5 MPa) and positive (δCF F ∼ 0.5 MPa) in a limited northermost section. Note that since Coulomb pre-stress changes are only estimated over a 7-year period, they may represent a minimum amount in the SSE area at the onset of the Chengkung rupture. Finally, for a rake of 70 • on the patch corresponding to the SSE source, δCF F ∼ 0.1-0.2 MPa for all models (Figures 4c & dand S11d).6. Discussion and concluding remarksCoseismic static stress changes are permanent perturbations, with effects lasting from days to years [e.g.,[START_REF] Segou | The role of seismic and slow slip events in triggering the 2018 M7.1 Anchorage earthquake in the Southcentral Alaska subduction zone[END_REF], and can therefore explain an extended period for triggering of aftershocks[START_REF] King | Static stress changes and the triggering of earthquakes[END_REF] and SSEs[START_REF] Hayes | Triggered aseismic slip adjacent to the 6 February 2013 M w 8.0 Santa Cruz Islands megathrust earthquake[END_REF][START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF]. Besides, long-lasting stress effects from postseismic deformation can also play a significant role in promoting delayed rupture[START_REF] Segou | Testing earthquake links in Mexico from 1978 to the 2017 M = 8.1 Chiapas and M = 7.1 Puebla shocks[END_REF]. The 2017 M w 7.1 Puebla earthquake, Mexico, was likely triggered by postseismic stress changes (up to 0.1 MPa) following the 2012 M w 7.5 Oaxaca event [Segou and Parsons, 2018], while a combination of postseismic stress changes from the 1964 M w 9.2 Alaska earthquake and a 12 year-long SSE possibly promoted the 2018 M w 7.1 Anchorage earthquake [Segou and Parsons, 2020]. Here, the 2003 Chengkung earthquake represents the largest event impacting the LVF during the past decades. It generated maximum static Coulomb stress perturbations in the SSE source region two orders of magnitude larger than any other event (Figure S12). D R A F T November 6, 2020, 11:55am D R A F T

  Combining δCF F and the kinematic analysis we propose the following scenario: 7 years of afterslip leads to a δCF F > 0 in the north, which triggers the SSE. It then propagates southward and was further stopped by a stress barrier, i.e. in an area where all the stress has been released by a combination of coseismic and postseimic slip. We also observe a slip deficit near the surface (Figure4e) which is not covered by the inferred SSE source model. We still know little about the physical mechanisms and conditions that lead to SSE generation near the surface, where deformation is more distributed[Thomas et al., 2014b]. Thus, while most of the deformation occurs on the LVF, we could have failed to detect slip on secondary faults, underestimating the deformation near the surface, therefore leading to an apparent slip deficit. This is supported by the evidence that the shallow creeping portion of the fault (at the depths < 7 km) had acted as a rather effective barrier during the 2003 seismic rupture[Thomas et al., 2017]. This means that no remaining stress (therefore no slip deficit) was available to help the rupture propagates to the surface. These observations raise the question of the role of the SSE transient zone during the 2003 Chengkung rupture. Since insignificant coseismic slip is infer in the transient region(Figures 4 and S13), the latter may have acted as a barrier, limiting the rupture propagation northeastward. Rupture penetration in transient regions being mainly controlled by effective stress differences between "fast" and "slow" slip regions[START_REF] Lin | Overlapping regions of coseismic and transient slow slip on the Hawaiian décollement[END_REF], further analysis are needed to decipher whether a dominant aseismic slip mode prevails throughout the earthquake cycle for the transient slow slip zone [e.g.,[START_REF] Perfettini | Seismic and aseismic slip on the Central Peru megathrust[END_REF][START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF] (i.e., a permanent barrier), or if D R A F T November 6, 2020, 11:55am D R A F T

Figure 1 .Figure 2 .Figure 4 .

 124 Figure 1. (a) Map of southern Taiwan. Inverted triangles denote broadband seismome-

  compatible with the geometry of the LVF. We assess the relative location uncertainties by a bootstrap resampling method following Waldhauser and Ellsworth[2000] for two types of uncertainties. First we assess the relative uncertainty inside the multiplet by bootstrapping residual times from the relocation of the multiplet without the additional earthquakes. Second is absolute uncertainty of the multiplet location estimated from the bootstrapping performed on the joint relocation. We obtain relative errors of about 24, 19 and 42 m in EW, NS, and vertical directions, and absolute errors of 178, 124 and 243 m in EW, NS, and vertical directions, respectively (Figure

	T	November 6, 2020, 11:55am	D R A F T

Acknowledgments

We are grateful to Associate Editor German Prieto and two anonymous reviewers for their constructive suggestions and comments allowing to improve the manuscript. We thank Alan Linde, Selwyn Sacks and the support staff of the Carnegie Institution of Washington for the construction, installation, and maintenance of the dilatometers, and Hsin-Ming Lee for collecting the strainmeter data. The authors acknowledge Yi-Chuen Tsai and Mei-Jun Yu for processing the GPS data and Ya-Ju Hsu for providing codes for visualizing seismicity. We are grateful to Kate Chen and Yaochieh Chen for providing the repeater catalog. We thank Wen-Tzong Liang and Wen-Hui Lee for the access to the seismological data. Some figures were drawn using the Generic Mapping Tools [Wessel and Smith, 1998]. Strainmeter, hydrological, tide-gauge and broadband seismological data (CWB and BATS network [Institute of Earth Sciences, Academia Sinica, Taiwan, 1996]) are archived in PANGAEA (https://issues.pangaea.de/browse/PDI-25970). This research is supported by the Ministry of Science and Technology grant MOST 108-2116-M-001-027-MY2. This is the contribution of the Institute of Earth Sciences, Academia D R A F T November 6, 2020, 11:55am D R A F T

Templeton, D. C., R. M. Nadeau and R. Bürgmann (2008), Behavior of repeating earthquake sequences in Central California and the implications for subsurface fault creep, Bull. Seismol. Soc. Am., 98 (1), 52-65.

Thomas, M. Y., J. P. Avouac, J. Champenois, J. C. Lee and L. C. Kuo (2014a), Spatiotemporal evolution of seismic and aseismic slip on the Longitudinal Valley Fault, Taiwan, J. Geophys. Res., 119 (6), 5114-5139.

Thomas, M. Y., J. P. Avouac, J. P. Gratier and J. C. Lee (2014b)