Léa Pousse-Beltran

Anne Socquet
email: anne.socquet@univ-grenoble-alpes.fr

Lucilla Benedetti
email: benedetti@cerege.fr

Marie-Pierre Doin
email: marie-pierre.doin@univ-grenoble-alpes.fr

Magali Rizza
email: rizza@cerege.fr

Nicola D ' Agostino
email: nicola.dagostino@ingv.it

Localized afterslip at geometrical complexities revealed by InSAR after the 2016 Central Italy seismic sequence

Keywords:

Key points:

-We monitor pre and post-seismic deformation of the 2016 seismic sequence using twoyear InSAR time-series -Centimetre scale post-seismic surface displacements are detected after October 30, 2016 Mw 6.5 mainshock (Norcia earthquake) -Localized shallow afterslip occurred at structural complexity that may have hindered the propagation of seismic ruptures

Introduction

Monitoring the spatial and temporal variations of the slip on a fault enables researchers to better assess stress build-up on seismic asperities and slip released during seismic cycle [START_REF] Avouac | From Geodetic Imaging of Seismic and Aseismic Fault Slip to Dynamic Modeling of the Seismic Cycle[END_REF][START_REF] Bürgmann | The geophysics, geology and mechanics of slow fault slip[END_REF][START_REF] Chen | Creeping faults: Good news, bad news?[END_REF][START_REF] Harris | Large earthquakes and creeping faults[END_REF]. In the rate and state formulation, rupture propagation can be hindered by rate-strengthening sections of a fault, which tend to slip via creep rather than in seismic rupture (e.g., H. [START_REF] Perfettini | Seismic and aseismic slip on the Central Peru megathrust[END_REF][START_REF] Hirose | Slow Earthquakes Linked Along Dip in the Nankai Subduction Zone[END_REF]. Such barriers are also often associated with structural complexities -such as a change of strike, secondary faulting or interaction with inherited faults (e.g., [START_REF] King | Role of Fault Bends in the Initiation and Termination of Earthquake Rupture[END_REF][START_REF] King | Speculations on the geometry of the initiation and termination processes of earthquake rupture and its relation to morphology and geological structure[END_REF][START_REF] Wesnousky | Seismological and structural evolution of strike-slip faults[END_REF]. These structural complexities can act as a geometrical barrier, and are characterized by an increased equivalent strength [START_REF] Nielsen | The equivalent strength of geometrical barriers to earthquakes[END_REF].

Locating aseismic slip on the fault and comparing these locations with those of seismic slip and fault segmentation is of pivotal importance to better characterize the frictional behavior of a fault system and its relation with structural complexities.

A notable seismic sequence occurred in the Central Apennines (Italy) in 2016-2017 [START_REF] Chiarabba | Faults geometry and the role of fluids in the 2016-2017 Central Italy seismic sequence[END_REF][START_REF] Perouse | Coseismic Slip Vectors of 24 August and 30 October 2016 Earthquakes in Central Italy: Oblique Slip and Regional Kinematic Implications[END_REF][START_REF] Cheloni | Geodetic model of the 2016 Central Italy earthquake sequence inferred from InSAR and GPS data[END_REF][START_REF] Huang | Coseismic deformation and triggered landslides of the 2016 Mw 6.2[END_REF]Cirella et al., 2018;Civico et al., 2018;Scognamiglio et al., 2018;Ragon et al., 2019), with four main events: the 24 th August 2016 MW 6.0 Amatrice event, the 26 th October 2016 MW 5.9 Visso event, the 30 th October 2016 MW 6.5 Norcia event, and the 18 th January 2017 MW 5.5

Campotosto event (Figure 1 and Table S1). This seismic sequence ruptured the complex Mt Vettore fault system (in red in Figure 1) [START_REF] Pizzi | Fault segmentation as constraint to the occurrence of the main shocks of the 2016 Central Italy seismic sequence[END_REF][START_REF] Porreca | Seismic reflection profiles and subsurface geology of the area interested by the 2016-2017 earthquake sequence[END_REF]Villani, Pucci, et al., 2018), and the adjacent Amatrice-Campotosto fault (in orange in Figure 1). During the Norcia earthquake, the rupturing of an antithetic fault on the opposite side of the Castelluccio basin seems necessary to fit geodetic data and is supported by alignments of relocated aftershocks (Chiaraluce et al., 2017;[START_REF] Walters | Dual control of fault intersections on stop-start rupture in the 2016 Central Italy seismic sequence[END_REF]Cheloni et al., 2019). In addition, the role of an inherited west-dipping thrust called OAS (Olevano-Antrodoco-Sibillini thrust) in the Norcia earthquake coseismic rupture geometry has been widely discussed. While some studies suggested that only the Mt Vettore fault system was activated (Chiaraluce et al., 2017;[START_REF] Liu | Rupture features of the 2016 Mw 6.2 Norcia earthquake and its possible relationship with strong seismic hazards[END_REF][START_REF] Papadopoulos | Earthquake Triggering Inferred from Rupture Histories, DInSAR Ground Deformation and Stress-Transfer Modelling: The Case of Central Italy During August 2016-January 2017[END_REF][START_REF] Pavlides | Ground deformation and fault modeling of the 2016 sequence (24 Aug. -30 Oct.) in central Apennines (Central Italy)[END_REF][START_REF] Pizzi | Fault segmentation as constraint to the occurrence of the main shocks of the 2016 Central Italy seismic sequence[END_REF][START_REF] Wang | Source parameters and triggering links of the earthquake sequence in central Italy from 2009 to 2016 analyzed with GPS and InSAR data[END_REF][START_REF] Xu | Source Parameters of the 2016-2017 Central Italy Earthquake Sequence from the Sentinel-1, ALOS-2 and GPS Data[END_REF], others suggested that also the OAS thrust ruptured as a reactivated high-angle normal fault during the event, as suggested by geodetic and seismological observations [START_REF] Cheloni | Geodetic model of the 2016 Central Italy earthquake sequence inferred from InSAR and GPS data[END_REF]Scognamiglio et al., 2018;[START_REF] Walters | Dual control of fault intersections on stop-start rupture in the 2016 Central Italy seismic sequence[END_REF]. Although the reactivation of the thrust is not clearly demonstrated (Cheloni et al., 2019), the OAS appears to have played a role in the aftershock distribution [START_REF] Chiarabba | Faults geometry and the role of fluids in the 2016-2017 Central Italy seismic sequence[END_REF]Chiaraluce et al., 2017;[START_REF] Pizzi | Fault segmentation as constraint to the occurrence of the main shocks of the 2016 Central Italy seismic sequence[END_REF]. The earthquakes of the 2016 sequence appear to have nucleated near crosscutting structures that seem to have been loaded by previous ruptures in the sequence (Chiaraluce et al., 2017;[START_REF] Pino | Clock advance and magnitude limitation through fault interaction: the case of the 2016 central Italy earthquake sequence[END_REF]. This seismic sequence is thus an excellent case study to better understand the link between structural segmentation, aseismic slip and frictional properties that might control this rupture propagation.

Post-seismic processes during this sequence have been observed in the seismicity (e.g., [START_REF] Albano | Aftershock rate and pore fluid diffusion: Insights from the Amatrice-Visso-Norcia (Italy) 2016 seismic sequence[END_REF][START_REF] Tung | Delayed poroelastic triggering of the 2016 October Visso earthquake by the August Amatrice earthquake, Italy[END_REF], but no aseismic slip has been detected with geodetic data so far. Using Interferometric Synthetic Aperture Radar (InSAR) time-series from

Sentinel-1 data, we document a small but significant post-seismic deformation transient that we find is best explained by aseismic slip on faults associated with this sequence. Surface displacements are presented and analysed. We also explore simple modelling schemes that provide a framework for our interpretation and discussion.

Geological setting

The Central Apennines were affected by an extensional phase during the Jurassic, followed by a compressive phase during the Neogene (e.g., [START_REF] Calamita | Contrasting styles of fault reactivation in curved orogenic belts: Examples from the Central Apennines (Italy)[END_REF]. The OAS (Olevano-Antrodoco-Sibillini thrust) is one the main thrusts resulting from the Neogene compressive phase [START_REF] Calamita | Structural styles, chronology rates of deformation, and time-space relationships in the Umbria-Marche thrust system (central Apennines, Italy)[END_REF] and has been interpreted as a transpressive ramp [START_REF] Di Domenica | Relationships between thrusts and normal faults in curved belts: New insight in the inversion tectonics of the Central-Northern Apennines (Italy)[END_REF]. In the area affected by the 2016-2017 seismic sequence, the OAS thrust is expressed by parallel splays associated with fault-bend folds characteristic of structural ramps [START_REF] Calamita | Analysis of thrust shear zones in curve-shaped belts: Deformation mode and timing of the Olevano-Antrodoco-Sibillini thrust (Central/Northern Apennines of Italy)[END_REF] (Figure 1). The ongoing ENE oriented extension of 2 to 4 mm/yr (D' Agostino, 2014;[START_REF] Carafa | Improving deformation models by discounting transient signals in geodetic data: 2. Geodetic data, stress directions, and long-term strain rates in Italy[END_REF][START_REF] Devoti | A Combined Velocity Field of the Mediterranean Region[END_REF], which probably began in the Early Pleistocene (e.g., [START_REF] Galadini | Active Tectonics in the Central Apennines (Italy) -Input Data for Seismic Hazard Assessment[END_REF], is currently accommodated through normal fault systems such as the Monte Vettore fault system that hosted the 2016-2017 seismic sequence.

3

Surface displacements during the seismic sequence

3.1

InSAR processing Synthetic radar interferometry (InSAR) is now systematically used to constrain deformation fields [START_REF] Elliott | The role of space-based observation in understanding and responding to active tectonics and earthquakes[END_REF] and can document centimetre to millimetre scale slow aseismic ground deformation using an adapted processing chain (Hussain et al., 2018;[START_REF] Aslan | Shallow Creep Along the 1999 Izmit Earthquake Rupture (Turkey) From GPS and High Temporal Resolution Interferometric Synthetic Aperture Radar Data (2011-2017)[END_REF]. We used C band (5.5 cm wavelength) images from Sentinel-1A/B images (Figure 1-A) spanning almost two years (July 28, 2015, to June 11, 2017) for the ascending track (A117, subswath IW3). To confirm the main observations made on the ascending track, we processed descending track (D22, subswaths IW2 and IW3) images from October 26, 2016 to February 11, 2017. SAR images were processed in VV polarization. We used the NSBAS processing chain [START_REF] Doin | Presentation of the small baseline NSBAS processing chain on a case example: the Etna deformation monitoring from 2003 to 2010 using ENVISAT data[END_REF][START_REF] Doin | InSAR measurement of the deformation around Siling Co Lake: Inferences on the lower crust viscosity in central Tibet[END_REF] modified for Sentinel data by [START_REF] Grandin | Interferometric Processing of SLC Sentinel-1 TOPS Data[END_REF] to generate differential interferograms. et al., 2003), resampled at 45m resolution, has been used to accurately coregister the focused SAR images and to correct interferograms from the topographic contribution to the interferometric phase. We removed a ramp in range and in azimuth for each interferogram using the methodology of [START_REF] Cavalié | Ground motion measurement in the Lake Mead area, Nevada, by differential synthetic aperture radar interferometry time series analysis: Probing the lithosphere rheological structure[END_REF] included in NSBAS. We de-noised the interferograms before unwrapping using collinearity, a criterion to characterize at each pixel the local spatial variability of the phase [START_REF] Pinel-Puysségur | Multi-link InSAR time series: Enhancement of a wrapped interferometric database[END_REF]. The collinearity was used to adapt the strength of the filter. We filtered for a window of 12 pixels (400m in range and 700m in azimuth); the filter is described in [START_REF] Doin | Presentation of the small baseline NSBAS processing chain on a case example: the Etna deformation monitoring from 2003 to 2010 using ENVISAT data[END_REF] and is based on the collinearity value which weights the complex phase in a sliding window. For the filter we had the option of adapting the weighting of the phase within windows of different sizes. This size adaptation depends on the collinearity within the windows. Unwrapping was performed in 2D with the NSBAS chain [START_REF] Grandin | Long-term growth of the Himalaya inferred from interseismic InSAR measurement[END_REF][START_REF] Doin | InSAR measurement of the deformation around Siling Co Lake: Inferences on the lower crust viscosity in central Tibet[END_REF]. After unwrapping, to account for errors associated with stratified troposphere, we removed a quadratic cross-function of elevation (z)

and azimuth to ramps in azimuth (y) and in range (x) estimation following the function ax+by+c+ez+fz × az+g × (z × az) 2 using a least-square approach [START_REF] Daout | Interseismic and Postseismic Shallow Creep of the North Qaidam Thrust Faults Detected with a Multitemporal InSAR Analysis[END_REF]. Time-series were then calculated following the NSBAS method [START_REF] Doin | Presentation of the small baseline NSBAS processing chain on a case example: the Etna deformation monitoring from 2003 to 2010 using ENVISAT data[END_REF][START_REF] Daout | Alongstrike variations of the partitioning of convergence across the Haiyuan fault system detected by InSAR[END_REF] using an approach based on the Small Baseline Subset time-series Analysis (SBAS) of [START_REF] López-Quiroz | Time series analysis of Mexico City subsidence constrained by radar interferometry[END_REF]'s algorithm. The smoothing of the pixel time-series is performed by minimizing the Laplacian of the temporal evolution of the deformation [START_REF] Cavalié | Ground motion measurement in the Lake Mead area, Nevada, by differential synthetic aperture radar interferometry time series analysis: Probing the lithosphere rheological structure[END_REF]. The final pixel size is 62 m in azimuth and 37 m in range. We removed pixels with an RMS value greater than 0.7. For the ascending track, we build a time series spanning the 2 years (from 28/07/2015 to 11/06/2017 ) (Figure S1 and see missing links in the time series in Figure S3). For this complete 2-year time-series, we encounter problems in unwrapping the co-seismic interferograms due to the large deformation with respect to the Sentinel wavelength in near field (aliasing).

Fringes are too close in space and cannot be unwrapped. At pixels that are incoherent in the co-seismic interferogram (i.e. in near field), this causes gaps in the complete (2 yrs) time-series spanning the main earthquakes (Norcia, Amatrice and Campotosto earthquakes) (Figure S3), or lead to an underestimation of the coseismic displacement (brown circles in Figure 2-D). In addition, we also build three time-series in between earthquakes in order to avoid possible bias: before the 24 th August (Amatrice earthquake), between 27 th August and 26 th October (between the Amatrice and Visso earthquake), and after the 30 th October (Norcia earthquake).

We excluded the SAR data from January 18, 2017, which produced noised interferograms. For the descending track, we build two time-series: between 27 th August and 26 th October (Visso earthquake), and after the 30 th October (Norcia earthquake). To rule out possible bias due to the Campotosto earthquakes affecting the area near Castelluccio and near Arquata del Tronto, we confirmed these previous observations with a shorter time-series calculated between the Norcia and Campotosto earthquakes (November 1 st -January 12 th ) (Figure 4). We prefer to use, for the rest of the manuscript, the longer post-Norcia time-series (November 1 st -June 11 th ) that shows a better signal-noise ratio.

In the post Amatrice earthquake time series (August 27 to October 26) we did not observe any localized surface displacements similar to the pattern observed after October 30 (Norcia earthquake) (Figure S4-B). Yet near the town of Amatrice, we observe diffuse surface displacements (< ~2.5 cm) moving away from the satellite. However, this time-series is constrained by only 9 scenes and 20 interferograms, which prevents from properly (i) constraining a low amplitude signal and (ii) correcting for atmosphere and topography. The surface displacements here have a low signal to noise ratio. The variance of the uncorrelated noise is ~ 40 mm 2 (see Supplementary Text S1 and Figure S5-A), and the standard deviation sigma of the noise is thus 6.3 mm. We take 3*sigma = 20 mm to set our limit of detection. The characteristic length scale of correlated noise is 5.0 km and there is autocovariance for distances smaller than ~15 km (see Figure S5-A). Observed patterns cannot be differentiated from noise, analyses on SAR images from other satellites should be carried out to confirm or not the post-Amatrice (August 24 earthquake) surface displacements.

Descending Track

To confirm the post-October 30 (Norcia earthquake) observations we processed descending interferograms (Figure 2-C). We used a 3-month dataset for the descending track (November 1st, 2016to February 11, 2017). Time-series calculated for the descending track also indicate slow deformation after Norcia earthquake (October, 30), that reaches on average 20.5 ± 2.7 mm on January 24 in the LOS direction near Arquata del Tronto (point 1). The deformation reached on average 10.8 ± 1.5 mm in the Castelluccio basin (in point 2). Assuming negligible north-south displacements, by combining the results from ascending and descending tracks, the displacement is dominated by subsidence in this area (Figure S6).

In the first order, deformation observed in the descending track is compatible with results inferred from the ascending track. Noise and unwrapping issues led us to mask noisy areas and resulted in more blank pixels in the descending track picture which should be used with caution. The relief projected in the descending LOS geometry masked a ridge (Mt Bove -Mt Vettore -Mt Gorzano high massifs) (Figure S7). In addition to noise, this resulted in an incoherent area and made unwrapping difficult in that area. The temporal evolution of the features slightly differs from the one in the ascending track. The removal of the 7 th November SAR images (affected by strong atmosphere conditions) and the poor constraint on the 19 th November SAR images (see the network in Figure S1) could affect the time-series in November. [START_REF] Falcucci | The Campotosto seismic gap in between the 2009 and 2016-2017 seismic sequences of central Italy and the role of inherited lithospheric faults in regional seismotectonic settings[END_REF] had similar difficulties in using descending SAR images to survey the January 2017 Campotosto seismic event. This might be due to the early morning acquisition time, which amplifies decorrelation due to change of moisture level or thawing.

In addition, we produced time-series calculated between the Amatrice and Norcia earthquakes. This time-series is also affected by a short-wavelength atmospheric turbulence that makes it difficult to interpret (Figure S5-B-C). The variance of noise is ~21 mm 2 . The characteristic length scale of correlated noise is 5.0 km and there is autocovariance for distances smaller than ~15 km (see Figure S5). Again, displacement less than 19 mm for the ascending track and under 14 mm for the descending track, cannot be considered as detectable signal, therefore constraining smaller surface displacements between the Amatrice and Norcia events demands complementary time-series with SAR images from other satellite.

Comparison with GNSS

Few GNSS stations exist near the studied area (ARQT, LNSS see station locations in Figure 2-A). Those GNSS time-series agree with InSAR time-series sampled at the same locations (Figure S8).

Pattern temporal evolution

By sampling the time-series, the temporal evolution is shown by averaging pixel displacement in area. To enhance the precision of the displacement temporal evolution without using time-series inversion, we temporally track a stable pattern in our unwrapped interferograms (Grandin, 2009). The method is detailed in Supplementary Text S2 and is illustrated by Figures S9-S10.

This method can only be used on areas where the displacement pattern is well defined, and works well in our case to track the evolution of the area near Arquata del Tronto (pattern in Figure 4 which encompasses the point 1). After the Norcia earthquake, the amplitude ratio is higher, and a logarithmic decay is clearly identified there (Figure 3-B). The pattern tracking shows a temporal evolution similar to the averaging surface displacement method. We are therefore confident that pixel averaging adequately captures the evolution of displacement that we want to describe. For the descending interferograms, the pattern is not precisely trackable and its time evolution is noisy (grey curve in Figure S10-C). This could be due to a poor correlation for the pattern in some interferograms. The evolution of the pattern around Castelluccio Basin (pattern in Figure 4) is poorly resolved too.

4 Modelling 3D deformation field

Comparison with topography and geology

To interpret the displacement pattern near Arquata del Tronto we compared it with the topography and the geology. There is no correlation with either the slope or the relief (Figure S11), thus we can exclude a gravitational process (i.e. landslide) as an explanation for the observed displacements. Concerning geology, the OAS thrusts fractured Meso-Cenozoic carbonate rocks over Neogene pelagic sediments; the carbonates host some of the largest aquifers in the Apennines while the Neogene sediments are thought to be an aquiclude (Figure S12) (Boni et al., 2010). The OAS delimiting these two units therefore acts as an impermeable boundary (Boni et al., 2010). Petitta et al., (2018) and Valigi et al., (2019) show that the 2016-2017 Italian seismic sequence had effects on spring discharge, water-table levels, and streamflow of those aquifers. In our results, the extent of the Maiolica unit, which hosts a shallow aquifer (Boni et al., 2010), corresponds well with the deformation area near Arquata del Tronto (Figure S11-A-B). The long pre-seismic InSAR time-series did not show any seasonal deformation in this aquifer, which means that groundwater seasonal processes can be excluded (point 1 in Figure 2-D). In addition, it seems difficult to get more than 1 cm of subsidence from winter rain in this area (i.e., [START_REF] Silverii | Transient deformation of karst aquifers due to seasonal and multiyear groundwater variations observed by GPS in southern Apennines (Italy)[END_REF].

Poro-elastic modelling

Poro-elastic effects have been studied during the 2016-2017 seismic sequence to explore the aftershocks and earthquakes triggering. [START_REF] Tung | Delayed poroelastic triggering of the 2016 October Visso earthquake by the August Amatrice earthquake, Italy[END_REF] suggest that fluid migration caused by the Amatrice earthquake (August 24, 2016) could have triggered the Visso earthquake (October 26, 2016), as well as some associated aftershocks. They inferred that this poro-elastic triggering should have occurred in an intermediately fractured crust. They also estimated that afterslip and viscoelastic-relaxation were quite negligible with respect to poroelastic effects. [START_REF] Albano | Aftershock rate and pore fluid diffusion: Insights from the Amatrice-Visso-Norcia (Italy) 2016 seismic sequence[END_REF] also suggest that post-seismic fluid diffusion after the Amatrice earthquake is related to aftershocks. According to their pore fluid diffusion model, there could be some associated afterslip (~ 10 cm) after the Amatrice earthquake.

The Maiolica karst aquifer is shallow and unconfined but [START_REF] Roeloffs | Poroelastic Techniques in the Study of Earthquake-Related Hydrologic Phenomena[END_REF] postulates that every aquifer reacts as a confined aquifer to a disturbance at short timescale. Here, to compare with the displacement pattern near Arquata del Tronto, we thus test a forward model of poroelastic rebound using Relax software [START_REF] Barbot | A unified continuum representation of post-seismic relaxation mechanisms: semi-analytic models of afterslip, poroelastic rebound and viscoelastic flow: Semi-analytic models of postseismic transient[END_REF]. We perform a simple poro-elastic model with no lateral variation of diffusivity, that does not account for the spatially heterogeneous distribution of local aquifers. We use the slip distribution models of the Norcia and Visso mainshocks (Maubant et al., 2017), (Figure S13 and Figure 5-A). We do not take into account the Amatrice earthquake since poroelastic rebound from the Amatrice earthquake would likely have stabilized by the time of the Visso and Norcia earthquakes (Figure 4-f in [START_REF] Albano | Aftershock rate and pore fluid diffusion: Insights from the Amatrice-Visso-Norcia (Italy) 2016 seismic sequence[END_REF]). We tested several diffusivities for a shallow layer (0-5 km depth), from 1.5 m².s proposed by the Tung and Masterlark (2018)'s aftershocks analysis to 10 4 m²/s (maximum diffusivity value for karst) [START_REF] Roeloffs | Poroelastic Techniques in the Study of Earthquake-Related Hydrologic Phenomena[END_REF]. Other parameters are described in detail in Figure 5. This model predicts uplift in the Castelluccio Basin and near Arquata del Tronto (Figure 5-D) where we observe subsidence (Figure 5-C). This discrepancy allows us to rule out poro-elastic rebound as the main driver of observed deformation. These poroelastic models are not exhaustive, we only explore simple geometrical configurations using parameter values informed by the local geology, however, other karstic configurations would likely predict uplift as well.

Viscoelastic modelling

Viscoelastic relaxation in the crust and in the mantle is also an important aseismic process, and has been inferred to be the driver of post-seismic deformation following many earthquakes (e.g., [START_REF] Pollitz | Mantle Flow Beneath a Continental Strike-Slip Fault: Postseismic Deformation After the 1999 Hector Mine Earthquake[END_REF][START_REF] Zhao | Dominant Controls of Downdip Afterslip and Viscous Relaxation on the Postseismic Displacements Following the Mw7.9 Gorkha, Nepal, Earthquake[END_REF]. To test whether viscoelastic relaxation is a plausible driving mechanism for post-Norcia deformation, we used Relax to perform forward modelling in a simple layered framework. Our model uses the stress perturbation from the three mainshocks: the Amatrice earthquake using the slip distribution of Ragon et al., (2019) added to the Visso and Norcia earthquakes using the slip distribution model of Maubant et al., (2017). We set the upper, middle and lower crust thickness to the values proposed in [START_REF] Laske | Update on CRUST1.0 -A 1-degree Global Model of Earth's Crust[END_REF] and [START_REF] Verdecchia | The Role of Viscoelastic Stress Transfer in Long-Term Earthquake Cascades: Insights After the Central Italy 2016-2017 Seismic Sequence[END_REF] (see Table 1). We tested a viscoelastic relaxation governed by a Newtonian rheology with 𝛾̇= 𝜏 𝜂 where 𝛾̇ is the viscous strain rate, τ is the deviatoric stress and η is the Newtonian viscosity. We used viscosity values η inferred in studies of viscoelastic relaxation modelling using GPS measurements following 1997 Umbria-Marche earthquakes, or using levelling line measurements following the 1915 Fucino earthquake [START_REF] Amoruso | Clues of postseismic relaxation for the 1915 Fucino earthquake (central Italy) from modeling of leveling data[END_REF][START_REF] Aoudia | Postseismic deformation following the 1997 Umbria-Marche (Italy) moderate normal faulting earthquakes[END_REF][START_REF] Riva | Viscoelastic relaxation and long-lasting after-slip following the 1997 Umbria-Marche (Central Italy) earthquakes[END_REF]. These studies inferred that the ductile structure of the Central Apennines consists of an elastic upper crust overlaying a middle crust (10 18 -10 19 Pa s), a viscous lower crust (10 17 -10 18 Pa s) and the upper mantle (10 21 Pa s) (Table 1).

Although those models are geometrically simple and do not take into account earthquakes prior 2016, they predict at a long-wavelength uplift (>50 km) in the area where we observe subsidence (Figure 6). This discrepancy in both sign and wavelength allows us to rule out a simple visco-elastic relaxation as the main driver of observed deformation. Viscoelastic models are not exhaustive here: nonlinear viscoelastic rheologies (e.g., power law creep) in lower crust [START_REF] Freed | Evidence of power-law flow in the Mojave desert mantle[END_REF]) could also be investigated; however finite element models with creeping lower crust seem also to predict uplift in similar settings of normal fault systems (e.g., [START_REF] Thompson | Vertical deformation associated with normal fault systems evolved over coseismic, postseismic, and multiseismic periods[END_REF]. It is worth also mentioning that the observed logarithmic decay could also potentially be associated with Burgers rheology or shear zone (e.g., [START_REF] Hetland | Effect of shear zones on post-seismic deformation with application to the 1997 Mw 7.6 Manyi earthquake[END_REF]) and could also be investigated in future work. Although, as with the poroelastic modeling, it seems unlikely that a different geometry or rheology could completely reverse the sign of uplift and subsidence (for example, a different rheology would only affect the temporal evolution/spatial distribution).

Modelling the temporal evolution of afterslip

The logarithmic-like temporal evolution of displacement (Figure 3-A), and the disagreement between the observations and the predictions of poroelastic and viscoelastic models, suggest that afterslip may have been the main driver of postseismic deformation [START_REF] Marone | On the mechanics of earthquake afterslip[END_REF]. Thus to characterize the temporal evolution of the deformation, we fitted the deformation decay (of the raw time-series and of the pattern tracking evolution)

with the logarithmic function (Figure 3 In which U: afterslip in time t, α: characteristic length scale, β× 𝑉 𝑖 : initial rate at the beginning of the post-seismic period, with 𝛽 a scaling vector by which the sliding rate evolves in response to the stress and 𝑉 𝑖 the pre-seismic slip rate.

At point 1 (near Arquata del Tronto), the good fit of the afterslip model suggests that the main post-seismic process is likely an afterslip phenomenon. Since we found similar results for c value (describing the temporal decay) using the raw pixel time series and using the pattern tracking (independent for the time-series inversion), we are confident in our fitting (Table 2).

The first post-Norcia Sentinel image has been acquired 2.6 days after the Norcia earthquake; therefore those 2.6 days of early afterslip are missing in our data. We compare the temporal decay of the deformation with the cumulative moment and number of aftershocks (blue and green curves in Figure S14). Although the completeness magnitude is high (ML 3 in Figure S5 in Chiaraluce et al., (2017)), those two curves obtained with less than 50 events are not following the same trends, aftershocks do not seem to be induced by afterslip. At point 2 (near Castelluccio di Norcia), the temporal decay does not fit with the afterslip law as the rate at the beginning of the post-seismic period is close to zero (c value for the purple curve fit in Figure 3-A). Other post-seismic processes could be a complementary driver of the deformation here and should be taken into account in further modelling (i.e. finite element forward model accounting for afterslip superimposed with poro-elastic rebound and viscoelastic relaxation).

Afterslip modelling on faults 4.5.1 Inversion strategy and fit to the data

To obtain the afterslip distribution we use the Classic Slip Inversion (CSI) Python tools [START_REF] Elliott | Himalayan megathrust geometry and relation to topography revealed by the Gorkha earthquake[END_REF]Jolivet et al., 2015) to invert for the slip. We use the constrained least-squares formula of [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF] to solve the inverse problem (see details in Supplementary Text S3 and Figures S15 to S18). We chose a dip of 40° for the Mt Vettore fault (as [START_REF] Cheloni | Geodetic model of the 2016 Central Italy earthquake sequence inferred from InSAR and GPS data[END_REF] see Figure S15) and project the fault down dip from the mapped fault at the surface (modeled fault geometry in yellow Figure 2-A). We discretize the fault into 88 rectangular patches. The smoothed cumulative surface displacement on February 11, 2017, measured in both ascending and descending tracks are inverted to obtain the afterslip distribution. The resolution is good for short-wavelength features at shallow depths (<5 km depth) (Figure S16). As several geometries have been used to model the mainshock rupture (OAS reactivation, antithetic fault), we explored three cases: (1) slip only on the Mt Vettore Fault, (2) slip on the Mt Vettore Fault and on the OAS and (3) slip on the Mt Vettore Fault and an antithetic fault with a geometry similar to the one proposed by Cheloni et al., (2019) and Maubant et al., (2017) (dip 65°) (Figures S17 to S20). We also explored the rake: (1) dip-slip only and (2) variable rake. A rigorous statistical comparison between the cases is difficult due to variation in model parameters and number of degrees of freedom, as explained by Cheloni et al., (2019). We thus used the RMS values to compare the different cases. Here the case that reduces the residuals the most is the inversion with an antithetic fault. This case is in addition consistent with geological and seismological observations (i.e. Cheloni et al., 2019). The different cases were compared using the RMS (root mean square), and allows us to choose the inversion with an antithetic fault (Table S2).

Both ascending and descending tracks are well modelled (Figure 7). It is worth noting that the resolution is low at depth (Figure S16). The total geodetic moment released by the afterslip model is equivalent to Mw ~5.78 (without taking into account the Campotosto earthquake) which corresponds to ~8.7% of the geodetic moment released during the Norcia earthquake (Mw 6.5).

Afterslip distribution 4.5.2.1 Near Arquata del Tronto

We obtain a maximum afterslip of ~100 mm below Arquata del Tronto at shallow depth (0-2 km depth) (Figure 7 Figure 8). The shallow part depth < 3-4 km toward the south is associated with very low coseismic slip for both the Norcia and Amatrice earthquakes and is located at the edge of the coseismic asperities. This is also consistent with an afterslip process:

the afterslip is located where the slip gradient increased the shear stress on the unruptured portions. This feature has also been observed after the L'Aquila earthquake [START_REF] D'agostino | Space-time distribution of afterslip following the 2009 L'Aquila earthquake[END_REF]. We calculate a geodetic moment released equivalent to Mw ~5.1 in this area (patches above 2km depth and south of Arquata del Tronto). The cumulative moment released by the seismicity (Figure S14-D) at this date corresponds to ~2.0% of the moment released by the afterslip model. The deformation is thus mainly aseismic. [START_REF] Brozzetti | High-Resolution Field Mapping and Analysis of the August-October 2016 Coseismic Surface Faulting (Central Italy Earthquakes): Slip Distribution, Parameterization, and Comparison With Global Earthquakes[END_REF] found some surface ruptures south of the 30 th October Norcia rupture (42°47'N in their Fig. 1-C) mapped by Villani et al., (2018). The ruptures mapped by [START_REF] Brozzetti | High-Resolution Field Mapping and Analysis of the August-October 2016 Coseismic Surface Faulting (Central Italy Earthquakes): Slip Distribution, Parameterization, and Comparison With Global Earthquakes[END_REF] are subtle and according to us may correspond to post-seismic deformation, as was observed after the L'Aquila earthquake [START_REF] D'agostino | Space-time distribution of afterslip following the 2009 L'Aquila earthquake[END_REF].

Castelluccio Basin

We observe a maximum slip of ~ 170 mm below Castelluccio at ~5 km depth. Some afterslip overlaps with the coseismic rupture area, but not in the area of maximum coseismic slip (Figure 7). As also the shape of the measured time-series did not fit with the afterslip's law (Figure 3-A), a sole afterslip process is unlikely. It is possible that poroelastic and fluid flow processes could be at work here since there is a large basal aquifer (Boni et al., 2010).

Modelling accounting for both processes could be performed, for example fully coupled poroelastic finite element numerical modelling with spatially variable material properties (e.g., [START_REF] Albano | Aftershocks, groundwater changes and postseismic ground displacements related to pore pressure gradients: Insights from the 2012 Emilia-Romagna earthquake[END_REF].

Discussion

Afterslip mechanism near Arquata del Tronto

Based on Dietrich (1979), [START_REF] Ruina | Slip instability and state variable friction laws[END_REF] and [START_REF] Marone | On the mechanics of earthquake afterslip[END_REF] equations, we estimated 𝛼 (the characteristic length scale over which the elastic stress changes by order of the frictional stress) in Figure 3 and Table 2. We can now estimate the friction parameter (𝑎 -𝑏) since according to [START_REF] Perfettini | Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan[END_REF] and [START_REF] Zhou | Time-dependent postseismic slip following the 1978 Mw 7.3 Tabas-e-Golshan, Iran earthquake revealed by over 20 years of ESA InSAR observations[END_REF]:

𝛼 = (𝑎 -𝑏) × 𝜎 𝑘
with k: the effective stiffness, and 𝜎 the effective normal stress. This relation can be applied to our case since the duration of our analysis (~200 days) is much shorter than the characteristic time td [START_REF] Gualandi | Space-time evolution of crustal deformation related to the Mw 6.3, 2009 L'Aquila earthquake (central Italy) from principal component analysis inversion of GPS position time-series[END_REF]. Where td= 𝛼 𝑉𝑝𝑙 (>12 years) where Vpl is the local plate loading rate (<2.1 mm/yr (Puliti et al., 2020)). The fit to the displacement of the timeseries sampled in point 1 by an afterslip law (Figure 3) leads to α= 30.9±4.0 mm. To estimate σ, we assume that the mean normal stress may vary from the hydrostatic to lithostatic pressure, with a rock density of 2.5 kg.m -3 [START_REF] Albano | Aftershock rate and pore fluid diffusion: Insights from the Amatrice-Visso-Norcia (Italy) 2016 seismic sequence[END_REF] and a depth of the slipping area of 3-5 km. This leads to σ = 44-124 MPa. For k =G/h with G=30 GPa (the shear modulus near the surface) and h=8 km (rate-strengthening depth based on the coseismic slip models), we obtain (a-b)=3 x 10 -3 -7.8 x 10 -4 . More complex models (e.g. finite element model, heterogeneous (a-b) values) could be performed by future studies in order to reproduce the surface displacements and to propose a more precise value of (a-b). However, our estimation is in agreement with experiments on carbonates. [START_REF] Scuderi | The role of fluid pressure in induced vs. triggered seismicity: insights from rock deformation experiments on carbonates[END_REF] found values of (a - b) evolve from velocity strengthening behaviour (a - b ≈ 0.005) at fluid pressure condition of sub-hydrostatic to a velocity neutral behaviour (a - b approaching 0), when the fault is at near lithostatic fluid pressure. [START_REF] Pluymakers | Frictional properties of simulated anhydrite-dolomite fault gouge and implications for seismogenic potential[END_REF] found that wet anhydrite and dolomite gouges at depths <6 km, exhibit (a-b) values ranging from 10 -2 to 10 -4 .

In Arquata del Tronto, the obtained (a-b) positive value (0.0026) is in agreement with a velocity-strengthening area. The obtained (a-b) value is also very low and corresponds to a quasi-neutral rate-dependency of friction implying a high sensitivity to stress perturbations.

Potential stress perturbations needed to reach strengthening may involve either a decrease in effective normal stress and / or an increase in shear stress [START_REF] Scholz | Earthquakes and friction laws[END_REF]. [START_REF] Walters | Dual control of fault intersections on stop-start rupture in the 2016 Central Italy seismic sequence[END_REF] calculated the Coulomb stress change (CFF) after the Norcia earthquake. At shallow depths, they found a positive CFF change below Arquata del Tronto, which might have triggered the observed afterslip in this area. An increase in pore fluid pressure would have a similar effect.

The afterslip is located near to the Olevano-Antrodoco-Sibillini (OAS) thrust. [START_REF] Chiarabba | Faults geometry and the role of fluids in the 2016-2017 Central Italy seismic sequence[END_REF] detected a high Vp/Vs near the OAS thrust (around 3 km depth), indicating high pore pressure that could have been further increased by nearby major shocks (Amatrice and Norcia earthquake). This high pore pressure could have induced the observed afterslip allowing for the release of the residual accumulated stress. This scenario is also supported by the fact that pore fluid pressure variation is common in this area, where the inherited structures seem to control the rupture pattern and could also concentrate high shear stress. Future work could examine detailed Coulomb Stress modelling with high spatial resolution taking all fault complexity and pore fluid effects into account. These works could be useful, despite some limitations related to the resolution of the coseismic slip and by the ambiguity about which faults were involved.

If this overlapping of coseismic slip and afterslip is real below 4-5 km depth (below Arquata del Tronto) could be an artefact caused by the smoothing used to regularize the inversion of the slip distribution, combined with the fact that the resolution at depth is limited.

If this overlapping of co-seismic and slow slip was confirmed, a simple interpretation of rate and state friction cannot be proposed here. Overlaps between coseismic slip and afterslip have been observed after several earthquakes. For example after 2004 MW 6 Parkfield earthquake, afterslip was inferred to overlap the coseismic slip [START_REF] Freed | Afterslip (and only afterslip) following the 2004 Parkfield, California, earthquake[END_REF]. [START_REF] Johnson | Frictional Properties on the San Andreas Fault near Parkfield, California, Inferred from Models of Afterslip following the 2004 Earthquake[END_REF] proposed frictional spatial heterogeneity to explain the afterslip distribution, and proposed (ab) values on the order of 10 -4 -10 -3 . After the 2015 Ilapel MW 8 megathrust earthquake, [START_REF] Barnhart | Coseismic slip and early afterslip of the 2015 Illapel, Chile, earthquake: Implications for frictional heterogeneity and coastal uplift[END_REF] also observed afterslip and coseismic slip overlapping, and advocated that stress heterogeneities likely provide the primary control on the afterslip distribution. After the 1978 MW 7.3 Tabas-e-Golshan earthquake, [START_REF] Zhou | Time-dependent postseismic slip following the 1978 Mw 7.3 Tabas-e-Golshan, Iran earthquake revealed by over 20 years of ESA InSAR observations[END_REF] also observe such overlapping and found (a-b)~3.10 -3 . They proposed that the fault is creeping during the whole interseismic period and that the earthquake propagates through the rate-strengthening region although they also stated that a change in frictional properties is also possible due to shear heating. [START_REF] Thomas | Rate-and-state friction properties of the Longitudinal Valley Fault from kinematic and dynamic modeling of seismic and aseismic slip[END_REF], after the 2003 MW 6.8 Chengkung earthquake estimated that (a-b) varies with depth from 0.018 near the surface to less than 0.001 at depth larger than 19 km, to explain areas showing both seismic and aseismic behavior they advocated shear heating processes. Using numerical modelling, [START_REF] Noda | Stable creeping fault segments can become destructive as a result of dynamic weakening[END_REF] argued that coseismic slip could propagate into velocity-strengthening regions of a fault, and that these regions can appear locked or creeping during the interseismic period. They observed this behavior due to rapid shear heating. For our case it is difficult to pinpoint the particular fault zone mechanism here, but we can propose a switch to (a-b) positive value below 4-5 km depth, due to coseismic shear heating processes [START_REF] Thomas | Rate-and-state friction properties of the Longitudinal Valley Fault from kinematic and dynamic modeling of seismic and aseismic slip[END_REF] requiring temperature changes of several hundred degrees (e.g., [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]. This could be a potential explanation since [START_REF] Smeraglia | Field-to nano-scale evidence for weakening mechanisms along the fault of the 2016 Amatrice and Norcia earthquakes[END_REF] suggested that nanostructures from the Mt Vettoreto Fault rocks were generated by coseismic shear heating and grain comminution (reduction of particle sizes).

Barrier of rupture propagation

Near Arquata del Tronto, at shallow depth (< 3-4 km) we observe that the aseismic slip is localized on an area associated with limited coseismic slip, for both the Amatrice and Norcia earthquakes. The observed aseismic slip seems located in an area where coseismic ruptures have not propagated during the Norcia or Amatrice earthquakes. This slowly slipping zone is located near a structural complexity (inherited thrust OAS), that Chiaraluce et al., (2017[START_REF] Pizzi | Fault segmentation as constraint to the occurrence of the main shocks of the 2016 Central Italy seismic sequence[END_REF] and Puliti et al. (2020) interpreted as a barrier that concentrated stress after the Amatrice earthquake. We propose that maybe as a result of this geometric complexity (that may concentrate high stress and high porosity resulting from fracturation), this zone might be also characterized by frictional properties at shallow depth that could have favoured the arrest of the ruptures (e.g., [START_REF] Kaneko | Towards inferring earthquake patterns from geodetic observations of interseismic coupling[END_REF]. [START_REF] Scholz | Earthquakes and friction laws[END_REF] shows that when an earthquake propagates into a velocity strengthening field it will produce a negative stress drop, that rapidly terminates propagation. Our simple calculations yield here a slightly positive (a-b) value, consistent with this model.

Conclusion

The Sentinel-1 InSAR time-series show cm-scale post-seismic displacements after the Norcia earthquake (October 30). The displacements are more obvious in the ascending track than in the descending track, and correspond to ~1 to 5 cm of subsidence in the Castelluccio Such small and localized slow deformation due to afterslip may be best detected with InSAR time-series in regions of sparse GNSS coverage. It might actually be quite common and could be an underestimated phenomenon. This likely points toward a bias in the literature in favour of high and wider afterslip, that is easier to detect. However detection of small afterslip transient is crucial to further understand the physics of earthquakes, and the link between slow slip and seismic rupture (e.g., [START_REF] Kaneko | Towards inferring earthquake patterns from geodetic observations of interseismic coupling[END_REF][START_REF] Rousset | A geodetic matched filter search for slow slip with application to the Mexico subduction zone: GEODETIC MATCHED FILTER FOR SLOW SLIP[END_REF]. 2: Parameters obtained from the fit of the afterslip model in Figure 3 using a non-linear least squares method [START_REF] Moré | The Levenberg-Marquardt algorithm: Implementation and theory[END_REF]. R2 is the coefficient of determination of the regression.

The pattern encompasses the point 1 (see Figure 4) 2).

(B) Temporal variation of pattern amplitude. Pattern is indicated in Figure 4. In red, the fit with an afterslip model (see results in Table 2). Figure 4: Mean velocity map between Mw 6.5 Norcia and Mw 5.7 Campotosto earthquakes (November 1 st -January 12 th ). This map is computed using 35 ascending interferograms. The swath A-B is used in Figure S11 to compare InSAR displacements and topography. Circles represent the earthquakes (see color code in Figure 1).

Figure 5

Figure 5: Poro-elastic rebound modelled with Relax software [START_REF] Barbot | A unified continuum representation of post-seismic relaxation mechanisms: semi-analytic models of afterslip, poroelastic rebound and viscoelastic flow: Semi-analytic models of postseismic transient[END_REF]. We used the coseismic slip models of the Norcia and Visso earthquakes from Maubant et al., (2017).

3. 2

 2 Time-series results and description of the main features 3.2.1 Ascending Track The ascending time-series built before the seismic sequence (from July 28 th 2015 to August 21 th 2016) does not show significant nor localized surface displacements along the main faults (Figure S4-A,B). After October 30 th , the (post-Norcia) time-series shows centimetre-scale displacements going away from the satellite in the LOS direction in three areas (Figure 2-B):  South of Amatrice: In the area affected by the Campotosto earthquake (January 18 th , 2017 Mw 5.0 -5.5 EQ) the ground surface moved away from the satellite by more than 60 mm in LOS (Figure 2-B). This coseismic displacement results in a step function in the time series, and can therefore be easily separated from any gradual post-seismic deformation.  Near Arquata del Tronto: At the southern extremity of surface rupture of the Mw 6.5 October 30 th Norcia earthquake (red faults in Figure 2), surface displacements are detected over an area of ~12 km 2 , and follow logarithmic evolution (in Figure 2-D and Figure 3-A see time-series at point 1 where the cumulative post-Norcia displacements in LOS on January 06, 2017 is in average ~35.5 ± 1.7 mm and reach 50.5 ± 2.1 mm on 30 April 2017)  Castelluccio Basin: On the hanging wall of the Mt Vettore fault, slow deformation affects an area of ~50 km 2 (Figure 2), and is associated with displacements in the LOS direction that are 13.2 ± 1.4 mm in average at Point 2 on 6 January and reach 37.9 ± 1.3 mm on 30 April 2017 (Figure 2-D and Figure 3-A).

  -A-B respectively) from the Marone et al. (1991) model and Zhou et al. (2018) reformulation for rate-strengthening afterslip:

  basin and at the southern tip of the Norcia surface rupture. In the area of the Castelluccio basin, whether the deformation can be explained by afterslip only is less clear and would require additional modelling to include the effects of poro-elastic release or shear zone deformation. At the southern edge of the Norcia coseismic rupture, deformation evolves with time following a logarithmic decay consistent with afterslip triggered by the Norcia earthquake. Although good observations of the 2.6 first days of afterslip are not available, after 3 months the slow deformation evidenced in this study released a geodetic moment of Mw ~5.78, which corresponds to ~8.7 % of the geodetic moment released during the Norcia earthquake. Cumulated moment released by aftershocks in the southern tip of the Norcia rupture accounts for ~2.0 % of the moment released by the modelled deformation, the deformation is thus here mainly aseismic. This afterslip takes place at the southern tip of the Norcia rupturing patch, and seems to have acted as a barrier to the propagation of Norcia and Amatrice ruptures at shallow depth (<3-4 km). Our results suggest that the observed post-Norcia earthquake afterslip might have been triggered in response to heterogeneities of pore fluid pressure, potentially facilitated by structural complexity and intense faulting in the area.

Figure 1 Figure 1 :

 11 Figure 1: Neotectonic framework of the Central Apennines. Black lines are main normal active faults compiled from Benedetti (1999), Tesson et al., (2016) and references therein (VBFS: Mt Vettore -Mt Bove Fault System, A-C F: Amatrice-Campotosto Fault). In red, coseismic surface ruptures observed in the field after the mainshocks of the Central Italy 2016-2017 seismic sequence (Villani, Civico, et al., 2018). (A) Blue and red rectangles are, respectively, the descending (D22 subswaths IW2 and IW3) and ascending (A117, subswath IW3) Sentinel tracks processed in this study. (B) Black dashed line represents the Neogene OAS thrust (Olevano-Antrodoco-Sibillini thrust) (Di Domenica et al., 2012). In orange the Amatrice-Campotosto fault which was activated during the 2016-2017 seismic sequence. ARQT and LNSS (in B) are GNSS stations available in this region fromhttp://geodesy.unr.edu/NGLStationPages/gpsnetmap/GPSNetMap_MAG.html(Blewitt et al., 2018).

( A )

 A Figure6: Surface displacements predicted by visco-elastic models, we used the Amatrice, Norcia and Visso earthquakes coseismic slip distribution modelled byRagon et al., (2019) andMaubant et al., (2017).(A) Coseismic slip distribution of the Amatrice earthquake from Ragon et al., (2019) (B-C-E-F) Surface displacements predicted by four visco-elastic models on 24 th February 2017 using Relax software (Barbot & Fialko, 2010). Mantle viscosity is fixed for all models at 10 21 Pas. (D) Post-Norcia earthquake (30 th October) cumulative displacement map (11 th February 2017) observed by InSAR (uplift component). (VBFS: Mt Vettore -Mt Bove Fault System, A-C F: Amatrice-Campotosto Fault, OAS: Olevano-Antrodoco-Sibillini thrust)

  The interferogram network, and examples of unfiltered and uncorrected interferograms, are provided in Figures S1 and S2. The Shuttle Radar Topography Mission digital elevation model (DEM) at 3 arc sec resolution (Rabus

Table 1 :

 1 Crust and mantle viscosities (η) tested in Figure6.

	7 Tables
	7.1 Table 1
		Depths (km)	Model 1 η	Model 2 η	Model 3 η	Model 4 η
	Upper crust	0-10.7	Elastic	Elastic	Elastic	Elastic
	Middle crust	10.7-21.7	Elastic	10 19 Pa s	10 18 Pa s	10 18 Pa s
	Lower crust	21.7-33.1	10 18 Pa s	10 18 Pa s	10 18 Pa s	10 17 Pa s
	Upper mantle	33.1-below	10 21 Pa s	10 21 Pa s	10 21 Pa s	10 21 Pa s
	7.2 Table 2
		Fit	α	C (days -1 )	γ	R2
	Point 1 near Arquata del Tronto	30.9 ± 4.0	0.26 ± 0.26	-6.4 ± 6.2	0.87
	Point 2 near Castelluccio di Norcia	156.4 ± 230.7	0.003 ± 0.007	2.54 ± 3.5	0.81
	Pattern tracking	1.6 ± 0.1	0.24 ± 0.1	-0.37 ± 0.21	0.95
	Table

between drained and undrained conditions. Poisson drained is 0.26 (cf. layer 1 in Albano et al. 2018), and for undrained conditions is 0.34 (from beta=0.4 in the table 1 and equation 15 in Barbot et al., (2010)). The layer between (0 and 5 km) is characterized by a diffusivity equals to 1.5 m²/s (Tung & Masterlark, 2018). (E) Modelled poro-elastic rebound time-series at point A using several diffusivities, from 1.5 m²/s (Tung & Masterlark, 2018) to 10 4 m²/s (Roeloffs, 1996). Because the delay time is inversely proportional to the assumed diffusivity, the response for high diffusivity is shorter. (VBFS: Mt Vettore -Mt Bove Fault System, A-C F: Amatrice-Campotosto Fault, OAS: Olevano-Antrodoco-Sibillini thrust). Cirella, A., Pezzo, G., & Piatanesi, A. (2018). Rupture Kinematics and Structural -Rheological Control of the 2016 Mw6.1 Amatrice (Central Italy) Earthquake from Joint Inversion of Seismic and Geodetic Data. Geophysical Research Letters, 0(ja). https://doi.org/10.1029/2018GL080894 Civico, R., Pucci, S., Villani, F., Pizzimenti, L., Martini, P. M. D., Nappi, R., & Group, the O. E. W. (2018). Surface ruptures following the 30 October 2016 Mw 6.5 Norcia earthquake, central Italy. Journal of Maps, 14(2), Application to the June 2000 Kleifarvatn Earthquake, Iceland." Geophysical Journal International 176 (2): 389-404. https://doi.org/10.1111/j.1365-246X.2008.03989.x. Tarantola, Albert. 2005. Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics.

Acknowledgments

We would like to thank the Editor, Associate Editor, Mong-Han Huang and two anonymous reviewers for their constructive suggestions, which helped to improve the manuscript substantially. We (Chiaraluce et al., 2017) and in http://cnt.rm.ingv.it/. The GNSS stations are downloaded from http://geodesy.unr.edu/NGLStationPages/stations/


Sébastien Gauthier
email: sebastien.gauthier@univ-rennes1.fr

Françoise Robin-Le Guen

Laurianne Wojcik

Nicolas Le Poul

Aurélien Planchat

Yann Pellegrin
email: yann.pellegrin@univ-nantes.fr

P G Level

Nadine Szuwarski

Mohammed Boujtita

Denis Jacquemin
email: denis.jacquemin@univ-nantes.fr

Patricia Guevara Level

Fabrice Odobel
email: fabrice.odobel@univ-nantes.fr

Synthesis and properties of novel pyranylidene-based organic

Keywords:

come

Introduction

The increasing demand for energy along with the numerous environmental pollution issues resulting from burning fossil fuel have triggered the expansion of research into clean, affordable and renewable energy sources to accommodate and maintain the needs for a modern life-style of a growing worldwide population [1]. In this respect, sunlight is one of the most widely available natural renewable energy sources that can globally contribute to resolving the environmental and energy issues of the present century [2]. Over the past thirty years, several types of solar cells capable of efficiently harvesting the incident sunlight have been developed [3]. Currently, traditional silicon-based solar cells involving long energy payback time still dominate the photovoltaic market. However, new solar cell technologies based on thin films have emerged as promising efficient substitutes [4].

As alternatives to silicon-based solar cells, Dye-Sensitized Solar Cells (DSSCs) embody a new generation of photovoltaic devices exhibiting high potential in solar power conversion efficiency (PCE) [5]. DSSCs are easy to fabricate, colorful, transparent and more importantly they work well in low light conditions [5a, 6]. A typical DSSC is made of a photoanode consisting of a crystalline mesoporous semiconductor (n-type TiO 2 ) on which sensitizer dye molecules are adsorbed, an electrolyte acting as redox mediator (I -/I 3 -), and a counter electrode (Pt). Among these components, the sensitizer dye is the key element defining the efficiency and the lifetime of the solar cell since it captures the sun light and initiates the electron transfer. Conventionally, transition metal complexes typically based on ruthenium and osmium metals have been used as efficient sensitizers with remarkable Power Conversion Efficiency (PCE) exceeding 10% [7].

However, these transition metals are on the verge of disappearing due to their limited natural availability combined with extensive recent uses. Consequently, they have become quite expensive, making them inappropriate for the development of cost-effective and environmentally friendly photovoltaic devices. Therefore, low-cost copper and zinc complexes have been investigated as DSSC dyes [8]. While such complexes exhibit good conversion efficiency comparable to those of ruthenium complexes, their practical applications in actual devices are hampered by their complicated and low-yield syntheses.

To circumvent the issues associated with transition metal complexes, metal-free organic dye sensitizers, usually based on a donor-(π-linker)-acceptor type push-pull molecules have attracted considerable attention due to their simple and affordable synthetic process, their tunable structure, their high molar extinction coefficient, and their excellent intramolecular charge transfer (ICT) properties [9]. Numerous metal-free organic dyes have been investigated, and some have been found to exhibit conversion efficiency performances similar to those displayed by ruthenium-based complexes [10]. Various chromophores containing conjugated electron donor groups, such as triphenylamine [START_REF] Wang | Triarylamine: versatile platform for organic, dye-sensitized, and perovskite solar cells[END_REF], dialkylaminophenyl [START_REF] Liang | Arylamine organic dyes for dye-sensitized solar cells[END_REF],

indoline In this paper, we describe the syntheses, characterizations, electronic properties and theoretical investigations of this series of new dyes. The measurements of the photovoltaic performances in TiO 2 -based DSSCs, coupled with Intensity-Modulated Photovoltage Spectroscopy (IMVS) experiments enable us to establish structure-property relationships for these dyes by correlating the changes of their efficiencies to their structural modifications, and to show that the photovoltaic performances are predominantly controlled by the -conjugated spacer as well as by the substituents placed on the exocyclic carbon on the pyranylidene electron donor.

Chart. Structures of the dyes investigated in this work. The dye 20a has been drawn as reference to show the moieties corresponding to the investigated parameters within each series (highlighted in yellow, blue and green).

Results and Discussion

Synthesis of the Dyes

The key intermediate for obtaining the push-pull dyes is the alkynes 5b-d, 7a-b, 13a, and 14a

(see Schemes 1-3). First, the two alkynes 13a and 14a are prepared following known synthetic routes [19c, 21a]. All the target sensitizer dyes have been characterized by 1 H and 13 C NMR, and highresolution mass spectroscopy, see Supporting Information (SI). The prepared dyes are also perfectly stable in their solid state form and did not require any specific storage conditions.

Electronic UV-visible absorption

The electronic spectra of compounds 10b-d, 12a-b, 19a, 20a, 21a, and 22a recorded in dichloromethane solution (1.1 × 10 -5 M) in the 300-800 nm wavelength range are displayed in Figure 1 and the positions of the bands are listed in Table 1. For all compounds, the least energetic absorption band in the visible region can be assigned to an ICT process between the donor and the acceptor moieties (see theoretical calculations below). The molar extinction coefficients (ε) of compounds 10b, 12a, 19a, 20a, 21a, and 22a exceed 20,000 M -1 cm -1 with a maximum at 53,000 M -1 cm -1 for compound 22a. These coefficients exceed those of the standard ruthenium N3 [START_REF] Horiuchi | Highly-efficient metal-free organic dyes for dye sensitized solar cells[END_REF] and N719 [START_REF] Wang | A high molar extinction coefcient sensitizer for stable dye-sensitized solar cells[END_REF] dyes (13,900 and 14,000 M -1 cm -1 , respectively), indicating that most of the dyes of this series are good candidates for light harvesting. The UV-visible spectra also show that the modification of the aromatic fragment on the exocyclic carbon on 20a (phenyl) and 12a (thienyl) resulted in red-shifted ICT absorption bands. The groups placed in the 2-and the 6-positions of the pyranylidene fragment, i.e., phenyl, thienyl, naphtyl and biphenyl have a limited influence (10 nm) on the absorption spectra. The increase of the intensity (molar extinction coefficients) of the maximum absorption of 21a and 22a, as compared to that of 20a, can be explained by the enhanced electron delocalization over the oligothiophene spacers of increasing length. Overall, among all the dyes, 12a and 12b demonstrate the broader absorption spectra, while 10d and 10c appear as being the poorest light collectors. Dye  max /nm (10 4 M -1 .cm -1 )

 onset /nm a 10b 389 (5.3),472sh (3.0) 568 10c 391 (5.0),488sh (2.4) 581 10d 396 (3.0),490sh (1.0) 568 12a 399 (3.5),530 (1.9) 653 12b 406 (3.4),520 (2.3 a Estimated from the wavelength with the intersections tangent method (onset).

Electrochemical study and electron transfer driving forces

The synthesized dyes were studied by cyclic voltammetry (CV) at a Pt working electrode in CH 2 Cl 2 /NBu 4 PF 6 0.1 M. The electrochemical data are gathered in Table 2. All dyes, except 19a, display two successive oxidation processes, as previously observed with pyranylidene analogues [20a, 21a]. As shown in Figure 2 for compound 20a, the first process at E 0 (1) is fully reversible at v = 0.1 V s -1 (Figure 2A, black curve) and is ascribed to the mono-oxidation of the pyranylidene moiety, leading to the formation of a pyrylium compound. The second process at E 0 (2) is quasi-reversible for all compounds excepted for 19a and leads to the formation of an unidentified species which is reducible at ca. -0.1 V vs. NHE on the backscan (Figure 2A, red curve). Scanning towards negative potentials (Figure 2B) shows the presence of broad and irreversible peaks in the -0.4 V to -0.8 V potential range, corresponding to the reduction of the cyanoacrylic acid anchoring group [21a]. The specific redox behavior of compound 19a is consistent with previous studies on H-pyranylidenes, which showed that the oxidation of the pyranylidene is followed by fast dimerization yielding bis-pyrylium species [25]. Such dimerization process does not occur when a substituent different than H (19a) is linked to the exocyclic carbon. Analysis of the data for the oxidation processes indicates that the nature of the phenyl or thienyl spacer has a significant effect on the redox properties. For example, the simple switch from a phenyl to a thienyl spacer in close position to the pyranylidene (20a to 12a) decreases both E 0 (1) and E 0 (2) values by 60 mV and 200 mV, respectively (see Table 2). This effect probably results from an increase of the electron density delocalization, (see computational part below). Hence, compound 12b displays the lowest E 0 (1) value amongst the investigated compounds. Moreover, as previously concluded for Pt complexes [21a], the increase of the chain length between the cyanoacrylic acid group and the phenyl-substituted pyranylidene moiety from one to three thienyl groups (compounds 20a, 21a and 22a) does not induce any major changes on the voltammetric response, neither in oxidation nor in reduction (Table 2). Finally, CV studies emphasize that groups in the 2,6 positions of the pyranylidene ring do not influence the redox behavior (compounds 20a, 10bd). e Irreversible anodic peak value. f Values estimated by assuming that E 0 (1) = (E pa (1)

The measured oxidation potentials were used to estimate the injection (G inj ) and dye regeneration (G reg ) free energies (Table 2). As these dyes are not luminescent, the 0-0 energy ( -0 ) was estimated from  onset of the absorption spectra. The results of Table 2 clearly indicate that all the dyes display very significant Gibbs free energies, which is more than sufficient to obtain efficient electron injection in the conduction band of TiO 2 (G inj ≈ -0.5 to -0.8 eV) and to oxidize iodide into triiodide (G reg ≈ -0.6 eV) (Table 2). In fact, these values are quite constant within this series of dyes, suggesting that the photovoltaic performances will not likely be controlled by these energetic parameters. Accordingly, ΔG reg and ΔG inj values are slightly over-estimated for 19a, because E 0 (1) corresponds to an anodic peak potential, not a half-wave potential, due to the irreversibility of the oxidation process for this dye.

Quantum chemical calculations

To gain more insights into the nature of the electronic excited states, we performed Time-Dependent Density Functional Theory Calculations on all dyes, with an approach taking into account solvation effects (CH 2 Cl 2 ) and detailed in the SI together with all HOMO and LUMO plots and energies. One representative example is given on Figure 3. As illustrated in Figure 3, the HOMOs are systematically centered on the pyranylidene segment, whereas the LUMO is located on the cyanoacrylic moiety and vicinal thienyl ring.

HOMO LUMO

This represents a very favorable configuration for the efficient injection of electrons in the TiO 2 conduction band, since, upon light absorption, the electron density undergoes a significant shift from one extremity of the molecule to the other one at the site of the anchoring group. More quantitatively, the computed lowest dipole-allowed transitions are listed in Table 3 together with the ICT parameters as obtained by the so-called Le Bahers ' d CT approach [START_REF] Bahers | A qualitative index of spatial extent in charge-transfer excitations[END_REF][START_REF] Jacquemin | What is the "best" atomic charge model to describe through-space charge-transfer excitations ?[END_REF].

Although TD-DFT vertical transition energies do not strictly correspond to the absolute experimental  max , we note a good qualitative match between the computed values (Table 3) and the measured ones (Table 1), with, in particular: i) almost no change of the absorption band positions within a series of homologous dyes (10b, 10c, 10d and 12a, 12b); ii) a significant bathochromic shift when going from 10 to 12 (computed value: ca. +60 nm, measured value: ca. +40 nm), which can be explained by the more delocalized nature of the HOMO in the latter compounds (see SI); and iii) 19a-22a transition energies quite similar to those of 10b-10d. In addition, the energies of the frontier MOs are not strongly affected by the substitutions on the pyranylidene moiety (see SI), consistently with the electrochemical study described above. The data listed in table 3 show that the ICT parameters are large for all dyes at the exception of the dyes 12a, 12b, and 22a, due to the absorption of light inducing a charge transfer of ca. 0.6-0.7 electron over more than 5 Å.

This ICT is illustrated in Figure 4 with the electron density difference plots clearly showing strong variations of the densities. The cyanoacrylic moiety, mostly in red, exhibits the characteristic electron density of an acceptor and the pyranylidene group, mostly in blue, that of a donor, whereas the side rings of the donor do not play any significant role in the charge transfer. In both 12a and 12b, the addition of thienyl groups on the exocyclic carbon on the pyranylidene group results in an increased electron density delocalization, inducing absorption redshifts. However, the theory shows that, as a consequence of this delocalization, a smaller ICT is obtained, probably due to the diffuse nature of the HOMO (vide supra see SI)

For 22a, the data in table 3 also show a less effective ICT, which is clearly explainable from the plots of Figure 4: in this dye, the donating and accepting moieties are so far apart that the charge transfer actually occurs from the -linker to the accepting moiety, the donor playing almost no role in the lowest electronic transition. This illustrates the well-known fact that longer linkers do not necessarily yield strong CT [START_REF] Ciofini | Through-space charge transfer in rod-like molecules: lessons from theory[END_REF].

Photovoltaic measurements and Intensity-Modulated Photovoltage Spectroscopy (IMVS) experiments

The synthesized push-pull dyes were used to prepare dye sensitized solar cells made of one layer of nanocrystalline TiO 2 particles (12 m thick) and one layer of scattering TiO 2 particles (4 m thick) and containing an iodide/triiodide-based electrolyte (see experimental part for details). The naked electrodes were dyed overnight at room temperature in a terbutanol/acetonitrile solvent mixture. All the cells were also studied by IMVS in order to determine the electron lifetime (τ e-) for each dye. The electron lifetime (τ e-) corresponds to the average time that the injected electrons in TiO 2 last before charge recombination either with the oxidized dye or with triiodide redox mediator in the electrolyte takes place.

It is well known that organic dyes tend to aggregate after their chemisorption on TiO 2 surface.

This process induces an excited state quenching and is therefore detrimental to the electron injection efficiency [START_REF] Zhang | Dye aggregation in dye-sensitized solar cells[END_REF]. Accordingly, we first investigated the impact of the increasing concentration of chenodeoxycholic acid (CDCA) in the dye bath for chemisorption of dye 20b and found that the optimal concentration was 1 mM (Table S1). Indeed, in absence of CDCA the aggregation adversely impacts the performances of the solar cell, but, since CDCA is a competitor to the dye adsorption a large concentration of CDCA decreases the dye loading and consequently the short-circuit current (Jsc). These optimal chemisorption conditions (dye at 0.25 mM in CH 2 Cl 2 /EtOH:1/1 with 1 mM of CDCA) and electrolyte composition ([tBuPy] = 0.1 M) were used to determine the photovoltaic performances of the entire series of dyes.

The UV-Vis spectra of the photoanodes are given in the SI (Figure S60), and show that all electrodes strongly absorb light in the visible and thus display more than 90% light harvesting efficiency until 550 nm. Moreover, 12a and 12b dyed photoanodes exhibit particularly redshifted absorption, in agreement with their absorption spectra recorded in solution. The metrics of the solar cells compiled in Table 4 are the average of at least 6 cells made with the same dye. The incident photon-to-current efficiency (IPCE) spectra of the solar cells are shown in Figure 5. Overall, the best performing dyes are 20a and 10c, while the least efficient ones are 12a, 12b and 19a. The dyes 12a and 12b were found to show a less effective CT than most dyes investigated, whereas 20a and 10c showed amongst the strongest CT character (Table 3). In this series of dyes, there is almost no correlation between the Light Harvesting Efficiency (LHE) (broad and intense absorption spectra) and the Jsc delivered by the sensitizers.

However, it is worth noting that the small d CT estimated by TD-DFT in the case of 12a and 12b correlates well with the quite low Jsc measured for those two dyes.

Conversely, the most efficient dyes are characterized by the longest electron lifetime in TiO 2 , meaning that this parameter is the major key player for the efficiency. Below, we discuss the influence of the dye structural modifications (parameters), as described on the chart given in the introduction, on their photovoltaic performances.

The first parameter is the length of the spacer (i.e. the number of the thienyl units in the conjugated spacer), which connects the pyranylidene moiety to the anchoring cyanoacrylic acid group (20a versus 21a versus 22a). Although the Light Harvesting Efficiency (LHE) increases, as expected, with the number of thienyl units (Figure 2 and table 1), since the enhanced-conjugation decreases the energy gap between the HOMO and LUMO, and consequently, red-shifts the absorption bands (Figures 1 andS58) but without improving the CT character (vide supra), the Power Conversion Efficiency (PCE) steadily decreases when the spacer length increases, due to a decrease of both Jsc and Voc.

These results point to an increase of the charge recombination process, which is known to reduce both the charge collection efficiency and the concentration of electrons in the conduction band and thereby the Voc. Upon inspection, the electron lifetime of these three dyes (20a, 21a and 22a) in TiO 2 can be ranked as follows: 20a ≈ 21a > 22a, which is consistent with the photovoltaic results of these three molecules. Indeed, when the electron lifetime increases, the Voc improves (Figure S59). Logically and assuming a dye standing upright, when the length of the spacer increases, the distance between the hole on the dye and the electron in TiO 2 should increase and consequently the charge recombination should diminish. However, the opposite result is obtained with the dyes 20a, 21a and 22a for which an increase of the spacer length leads to an increase of the charge recombination process. As shown in table 3, the calculated d CT for these dyes decreases when the length of the linker increases, which is in agreement with a lower injection efficiency. However, it is also possible that the 20a, 21a and 22a dyes are bent on the TiO 2 surface. Since the oligothiophene is not a linear molecule, these dyes may not be able to stand perpendicular to the TiO 2 surface, thus resulting is a distance between the hole on the dye and the TiO 2 surface being similar for all three dyes. Conversely, the dye loading for 20a is significantly smaller than those of 21a and 22a (table S3) indicating that those dyes could be more aggregated at the surface of TiO 2 than the 20a dye, thus resulting in smaller Voc and Jsc values.

The second parameter investigated was the influence of the nature of the substituents in the 2 and 6 positions of the pyranylidene moiety on the photovoltaic properties of the 20a, 10b, 10c and 10d dyes. Within this series of dyes, the nature of these substituents does not seem to significantly alter the electronic properties of the dyes, since the redox potentials (Table 2), the absorption spectra (Figure 1 and Table 1) and the CT character (Table 3) of these dyes (20a, 10b, 10c and 10d) are all quite similar. This is consistent with the topology of the HOMO and LUMO that almost do not involve these substituents (Figure S59).

Unsurprisingly, the photovoltaic performances of all these dyes are also quite similar, at the exception of the dye 10b containing thienyl groups as 2 and 6 substituents and displaying a significantly lower PCE (see 20a versus 10b) due to a lower Jsc and Voc. These results can be explained by a lower electron lifetime in TiO 2 (Figure S61) observed for 10b compared to those observed for 20a. Since both dyes exhibit similar electronic properties and similar loading capacities (table S3), we infer that the thienyl groups (whose affinity for iodine has been reported) [START_REF] Clifford | Dye mediated charge recombination dynamics in nanocrystalline TiO 2 dye sensitized solar cells[END_REF] tend to increase the local I 3 -concentration and thus favor interfacial charge recombination. This is in line with the lower Voc observed for 10b-based devices.

Finally, the influence of the nature of the substituents placed on the lower part of the pyranylidene moiety (substituents on the exocyclic carbon) on the photovoltaic performances of 20a, 12a, 19a, 12b and 10b dyes was examined. For those substituents, the replacement of a phenyl by a thienyl (20a versus 12a or 10b versus 12b) naturally leads to a red shift of the absorption band, but the Jsc of the DSSC made with 12a is lower than that made with 20a

(the same observation can be made when comparing 10b to 12b), indicating that either electron injection quantum yield or charge recombination come into play. The Gibbs free energy of the electron injection of this series of dyes is ranked according to the same order than their Jsc values (20a > 19a > 12a ≈ 12b) suggesting that the differences of Jsc could be partly related to the driving forces (Table 2). In addition, the Voc values follow the same trend in relation with the electron lifetime in TiO 2 , implying that charge recombination processes also significantly control the Jsc and the Voc values (Figure S63). Interestingly, the dye loading for 19a and 12a is drastically higher than that of 20a, and it could be once again inferred that 19a and 12a are more prone to aggregate, leading to overall poor photoconversion efficiencies as well as very poor Voc and lower Jsc values. As mentioned above, however, the thienyl groups can contribute to increase the I 3 -concentration close to TiO 2 and this will also lead to low Voc. Overall, the replacement of the phenyl substituents on the lower part of the pyranylidene or their substitution on the upper part by a thienyl seems to negatively impact the organization of the dyes on the TiO 2 surface, and consequently increases the charge recombination, as demonstrated by the decreasing electron lifetime measurements.

Conclusion

In this study, we present the synthesis of nine new organic push-pull dye sensitizers using pyranylidene fragments, which have been scarcely investigated as potential electron donors in push-pull chromophores so far. The dye electrochemical and light absorption properties were studied and the excited state properties modeled by TD-DFT calculations. These analyses show that all the dyes fulfill the criteria to act as efficient sensitizers for TiO 2 in DSSCs (large Gibbs free energy of the electron injection and dye regeneration, presence of strong charge transfer absorption bands and significant electron density of the LUMO on the anchoring groups). We also found that small structural changes significantly impact the photovoltaic performances, since the PCEs of this series of dyes vary from 1.5% to 5.5% depending on the nature of the substituents on the electron donor group. Together, the measurements of the photovoltaic performances and the IMVS experiment results show that the linker between the pyranylidene and the cyanoacrylic acid plays an important role to optimize these dyes for DSSC applications. The oligothiophene spacer brings -conjugation, which is favorable to Light Harvesting Efficiency (LHE), but is detrimental to the charge recombination processes.

Accordingly, new rod like molecules, such as oligophenylene of ethynylene or fused dithienylene, could be attractive -conjugated spacer candidates to engineer new push-pull dye sensitizers. All the substituents on the upper part (2 and 6 positions) of the pyranylidene electron donor moiety, except the thienyl group, poorly affect the electronic properties of the resulting dyes and consequently have a limited influence on the photovoltaic performances.

The major role of these upper part substituents is most likely to control the packing arrangement of the dyes on the TiO 2 surface. However, a variety of these substituents can be used without affecting the photovoltaic properties of the dyes, since we found that the replacement of a phenyl substituent by a biphenyl or a naphthyl substituent on the 2 and 6 positions of the electron donor group does not significantly change the PCE. On the other hand, the nature of the substituents on the lower part (substituents on the exocyclic carbon) of the pyranylidene moiety significantly affects the photovoltaic performances of the dyes, since the replacement of a phenyl substituent by a thienyl substituent leads to a noticeable PCE decrease.

This study shows that obtaining an optimal sensitizer dye exhibiting the highest photovoltaic performance requires phenyl substituents being placed in the lower part (substituents on the exocyclic carbon) of the pyranylidene electron donor group, since these groups provide the most suitable dye packing arrangement and allow to control charge recombination reactions.

Overall, this study provides a better understanding of the factors controlling the efficiency of these types of dyes and offers some molecular engineering rules to design dyes with improved photovoltaic performances.

  [13], carbazole[14], dithiafulvene[15] and tetrathiafulvalene (TTF)[16] groups have been investigated. While a variety of -linkers have been incorporated into these organic dyes to foster the ICT and improve the light-harvesting ability, cyanoacrylic acid remains the most popular due to its dual role as electron acceptor group (cyano) and an anchoring (carboxylic acid) moiety, resulting in excellent electron injection properties on titanium oxide[17].γ-Pyranylidene fragments are pro-aromatic electron-donating groups widely developed byGarín et al. [18] and our group[19]. Their electron-donating ability is based on the formation of an aromatic pyrylium fragment resulting from an intramolecular charge transfer process. In previous studies, the γ-pyranylidene moiety has been incorporated as an electron-donating group in various platinum-and ruthenium-based organometallic complexes and organic molecules for the preparation of push-pull structures with large second-order Non-Linear Optical (NLO) properties[20] or valuable behavior in DSSCs[21]. Pyranylidenes are easily synthesized molecules and they can be directly modified by the addition of different substituents in the 2 and 6 or on the exocyclic carbon, which allows to incorporate them as electron donor groups in push-pull chromophores. Very recently, the modification and optimization of pyranylidene derivatives, investigated as sensitizers for DSSCs, have led to encouraging Power Conversion Efficiency results reaching up to 6%[START_REF] Marco | Pyranylidene/thienothiophene-based organic sensitizers for dye-sensitized solar cells[END_REF]. Herein, we present nine new push-pull organic dyes 10b-d, 12a-b, 19a, 20a, 21a and 22a (Chart) based on the pyranylidene group as push (electron donor) unit. Our goal was to investigate the structure-property relationships within this class of dyes. The compound 20a was taken as reference. First, we studied the impact of the length of the thienyl -conjugated spacer (Chart Parameter I) on the photovoltaic properties (comparing 20a to both 21a and 22a). Then we evaluated how the nature of the substituents on the upper part (positions 2 and 6, Chart Parameter II) or on the lower part (substituents on the exocyclic carbon, Chart Parameter III) of the pyranylidene moiety affects these properties by comparing 20a to 10b to 10c to 10d, and 20a to 12a to 12b to 19a, respectively.

(

  Scheme 1. Synthetic route of alkyne compounds 5b-d and 7a-b. Reagents and conditions : a) THF, n-BuLi, -78°C to r.t., 2h b) K 2 CO 3 , CH 2 Cl 2 /MeOH, r.t., overnight.

Figure 1 .

 1 Figure 1. UV-Vis absorption spectra of dyes 10b-d, 12a-b, 19a, 20a, 21a, and 22a in CH 2 Cl 2 solution (1.1 × 10 -5 M) at room temperature.

Figure 2 .

 2 Figure 2. Cyclic voltammetry at a Pt electrode of compound 20a in CH 2 Cl 2 /NBu 4 PF 6 0.1 M (v = 0.1 V s -1 , C = 0.5 mM); A) in positive initial scanning; B) in negative initial scanning. Experimental curves were measured against Fc and recalibrated by taking E 0 (Fc + /Fc) = 0.71 V vs. NHE in those conditions.

  ) and E 0 (SCE) = 0.24 V vs. NHE. b E(S*/S + ): Fermi level (in V vs NHE) of excited molecules, calculated from : E(S*/S + ) = E 0 (1) -E 0-0 (S/S), where  -0 is the stored excitation energy, obtained from  onset . c G in eV, G reg = E 0 (I -/I 3 -) -E 0 (1) with E 0 (I -/I 3 -) = 0.35 V vs NHE; G inj = E(S*/S + ) -E CB (TiO 2 ) with E CB (TiO 2 ) = -0.5 V vs NHE.d Peak not detected.

Figure 3 .

 3 Figure 3. Contour plots of the HOMO and -LUMO orbitals of dye 20a.

Figure 4 .

 4 Figure 4. EDD plots obtained for all dyes. The blue (red) lobes represent regions of density decrease (increase) upon excitation, that is donor (acceptor) regions. Contour threshold: 0.0008 au.

Figure 5 .

 5 Figure 5. IPCE spectra of the dyes.

Table 1 .

 1 Absorption

properties of the dyes 10b-d, 12a-b, 19a, 20a, 21a, and 22a recorded in CH 2 Cl 2 (1.1 × 10 -5 M) at room temperature.

  Values calibrated vs. NHE from the experimental curves obtained vs. Fc + /Fc, by taking E 0 (Fc + /Fc) = 0.47 V vs. SCE in CH 2 Cl 2 /NBu 4 PF 6 (measured experimentally

				d	-0.65	-1.20	-0.63	-0.70
	10c	0.95	1.28	-0.51	d	-1.19	-0.60	-0.69
	10d 0.97	1.29	-0.44	-0.83	-1.21	-0.62	-0.71
	12a	0.91	1.10	-0.53	-0.78	-0.99	-0.56	-0.49
	12b 0.90	1.08	-0.55	-0.79	-0.99	-0.55	-0.49
	19a	1.02 e	d	-0.39	-0.61	-1.14 f	-0.67 f	-0.64 f
	20a	0.97	1.30	-0.47	-0.83	-1.21	-0.62	-0.71
	21a	0.96	1.27	-0.47	-0.75	-1.29	-0.61	-0.69
	22a	0.96	1.24	d	-0.84	-1.26	-0.61	-0.76

a

Table 3 .

 3 Theoretical vertical excitation wavelength [PCM-CAM-B3LYP/6-311+G(2d,p)] in nm, oscillator strength (between brackets), together with the CT parameters of all dyes: d CT and q CT given respectively in Å and in e.

	Dye	 vert /nm	d CT /Å	q CT /e
	10b	485 (1.6)	5.44	0.73
	10c	485 (1.7)	5.42	0.75
	10d	485 (1.7)	5.40	0.73
	12a	547 (1.4)	4.49	0.67
	12b	554 (1.4)	4.55	0.68
	19a	489 (2.0)	5.33	0.68
	20a	480 (1.6)	5.50	0.75
	21a	485 (2.3)	5.29	0.63
	22a	495 (2.4)	4.56	0.58

Table 4 .

 4 Metrics of the solar cells fabricated illuminated with a calibrated sunlight AM1.5 (100 mW/cm 2 ). The electron lifetime was measured by IMVS at a power of 150W/m 2 .

	Dye	Jsc (mA/cm²)	Voc (mV)	FF (%)	PCE (%)	Electron Lifetime (ms)
	10b	10.8 (±0.1)	588 (±4)	70 (±1)	4.44 (±0.03)	24.19
	10c	12.9 (±0.2)	605 (±1)	70 (±1)	5.41 (±0.09)	23.96
	10d	11.9 (±0.2)	594(±1)	71(±1)	5.05 (±0.10)	29.47
	12a	4.6 (±0.1)	481 (±5)	66 (±1)	1.47 (±0.11)	14.23
	12b	6.7 (±0.1)	505 (±1)	68 (±1)	2.31 (±0.05)	13.70
	19a	7.3 (±0.3)	480 (±2)	63 (±1)	2.19 (±0.10)	25.50
	20a	13.1 (±0.1)	603 (±1)	70 (±1)	5.52 (±0.01)	43.67
	21a	10.5 (±0.2)	590 (±1)	72 (±1)	4.44 (±0.08)	44.29
	22a	11.1 (±0.3)	566 (±2)	70 (±1)	4.41 (±0.09)	17.62

Acknowledgments

This research used resources of the GENCI-CINES/IDRIS, of the CCIPL and the Troy cluster, the two latter being installed in Nantes. D.J. is indebted to the Région des Pays de la Loire, for continuous support.