Menglin Si

Zhao-Liang Li
email: lizhaoliang@caas.cn

Franç Oise Nerry

Xiangyang Liu

Bo-Hui Tang

Spatial heterogeneity of driving factors-induced impacts for global long-term surface urban heat island

Keywords: global surface urban heat island, driving factors, spatial heterogeneity, seasonal variation, day-night contrast

Spatial heterogeneity of driving factors-induced impacts

for global long-term surface urban heat island

Introduction

With the rapid development of global economic integration and population growth, urbanisation has become a worldwide historical process (United-Nations 2019). Urban growth gradually transforms the original natural surfaces of vegetation, water, soil, and so forth into impervious urban areas that are mainly composed of built-up areas and infrastructure, accompanied by a large number of human activities [START_REF] Vitousek | Human domination of Earth's ecosystems[END_REF].

This changes the local meteorology, regional climate conditions, and energy exchange mode at the interface of the surface and the atmosphere between urban and surrounding suburban and rural areas [START_REF] Grimm | Global change and the ecology of cities[END_REF][START_REF] Oke | The energetic basis of the urban heat island[END_REF]. Consequently, urban or metropolitan regions become hotter than surrounding rural areas, leading to the surface urban heat island (SUHI) phenomenon (Oke 1988;[START_REF] Weng | Estimation of land surface temperature-vegetation abundance relationship for urban heat island studies[END_REF]. The advantages and disadvantages of SUHI has been proved hard to discern, though, it can reduce the consumption of energy significantly for mid-high latitude cities in winter while exhibits adverse effect to increase the frequency of heat disasters during heat waves [START_REF] Bründl | Advantages and disadvantages of the urban heat island -an evaluation according to the hygro-thennic Effects[END_REF][START_REF] Mohajerani | The urban heat island effect, its causes, and mitigation, with reference to the thermal properties of asphalt concrete[END_REF]. It is widely known that SUHI can negatively impacts energy use, urban population health, and air quality thus the considerable attention from scientific research teams, urban planners, and government departments [START_REF] Zhao | Prevalent vegetation growth enhancement in urban environment[END_REF][START_REF] Santamouris | On the impact of urban heat island and global warming on the power demand and electricity consumption of buildings-A review[END_REF]. It is necessary to enhance the theoretical research on drivers for global SUHI to support effective mitigation strategies in urban management [START_REF] Zhou | Satellite remote sensing of surface urban heat islands: Progress, challenges, and perspectives[END_REF].

As a vital phenomenon reflecting the urban thermal environment, SUHI is the consequence of the combined effects of surface properties, human activities, and meteorological and climate conditions [START_REF] Li | A review of remotely sensed surface urban heat islands from the fresh perspective of comparisons among different regions (Invited Review)[END_REF]. Urbanisation is the primary cause of alterations in surface attributes. Modifications in the LULC and surface structure constantly alter the surface physical properties and lead to differences in the surface energy redistribution modes between urban and rural areas, which further determines the spatiotemporal pattern of the SUHI [START_REF] Mika | Impact of 1.5 K global warming on urban air pollution and heat island with outlook on human health effects[END_REF][START_REF] Feranec | A review of studies involving the effect of land cover and land use on the urban heat island phenomenon, assessed by means of the MUKLIMO model[END_REF][START_REF] Lee | Observed increase in local cooling effect of deforestation at higher latitudes[END_REF][START_REF] Zhao | Strong contributions of local background climate to urban heat islands[END_REF]. Human activities have dramatically increased energy consumption and generated a large amount of waste heat, leading to additional surface heat flux, which affects the urban climate and near-surface temperature [START_REF] Jin | A new global gridded anthropogenic heat flux dataset with high spatial resolution and long-term time series[END_REF][START_REF] Meng | Quantifying the Anthropogenic Footprint in Eastern China[END_REF][START_REF] Yan | Review of recent studies of the climatic effects of urbanization in China[END_REF]. Background meteorology and climate conditions, such as differences in surface turbulence, precipitation, air temperature, wind, and cloud conditions, over urban and suburban areas often vary across different geographical regions or background climate zones [START_REF] Zhao | Strong contributions of local background climate to urban heat islands[END_REF][START_REF] Zhao | Urban growth and climate adaptation[END_REF][START_REF] Manoli | Magnitude of urban heat islands largely explained by climate and population[END_REF][START_REF] Sun | Understanding the variability of urban heat islands from local background climate and urbanization[END_REF][START_REF] Li | Relationship of surface urban heat island with air temperature and precipitation in global large cities[END_REF], which also has a strong impact on the SUHI.

Although the driving factors for the characteristics of SUHIs have been thoroughly explored for several decades [START_REF] He | Potentials of meteorological characteristics and synoptic conditions to mitigate urban heat island effects[END_REF][START_REF] Deilami | Urban heat island effect: A systematic review of spatio-temporal factors, data, methods, and mitigation measures[END_REF]Zhou et al. 2019), many inconsistencies exist in the current findings owing to different time ranges, regional scopes, and analysis methods [START_REF] Clinton | MODIS detected surface urban heat islands and sinks: Global locations and controls[END_REF][START_REF] Zhou | The role of city size and urban form in the surface urban heat island[END_REF][START_REF] Zhou | Surface urban heat island in China's 32 major cities: Spatial patterns and drivers[END_REF]). Globally, there has been minimal investigations on the disparity in response mechanisms for SUHI over different regions and at different temporal scales out of annual-summer-winter and day-night aspects. Therefore, revealing the spatial heterogeneity of the driving factors of global SUHI from a local perspective is necessary and beneficial in the proposal of targeted mitigation guidance.

To investigate the causality of SUHI, model-based schemes and statistical methods are used universally in the analytical aspect of the driving factors. Based on the surface energy balance (SEB), the underlying mechanisms of SUHI can be scientifically and physically determined [START_REF] Zhao | Strong contributions of local background climate to urban heat islands[END_REF][START_REF] Li | Urban heat island: Aerodynamics or imperviousness?[END_REF]. However, either the model simulation or the numerical deduction process is time-consuming and difficult to conduct on a global scale. Therefore, statistical models are more commonly used in SUHI studies [START_REF] Li | A review of remotely sensed surface urban heat islands from the fresh perspective of comparisons among different regions (Invited Review)[END_REF], wherein a series of analysis tools have been adopted, such as correlation analysis [START_REF] Zhou | Spatiotemporal trends of urban heat island effect along the urban development intensity gradient in China[END_REF][START_REF] Weng | A remote sensing-GIS evaluation of urban expansion and its impact on surface temperature in the Zhujiang Delta, China[END_REF][START_REF] Tran | Assessment with satellite data of the urban heat island effects in Asian mega cities[END_REF], linear regression models [START_REF] Dissanayake | Land-Use/Land-Cover Changes and Their Impact on Surface Urban Heat Islands: Case Study of Kandy City, Sri Lanka[END_REF], spatial association analysis [START_REF] Chen | Surface heat island in Shanghai and its relationship with urban development from 1989 to 2013[END_REF], and geographically weighted regressions [START_REF] Buyantuyev | Urban heat islands and landscape heterogeneity: linking spatiotemporal variations in surface temperatures to land-cover and socioeconomic patterns[END_REF][START_REF] Deilami | Modelling the urban heat island effect of smart growth policy scenarios in Brisbane[END_REF]. However, spatial heterogeneities of relationships generally exist between global SUHI and associated influencing factors, which are hardly detected by traditional statistical methods [START_REF] Chakraborty | A simplified urban-extent algorithm to characterize surface urban heat islands on a global scale and examine vegetation control on their spatiotemporal variability[END_REF][START_REF] Chakraborty | A spatially explicit surface urban heat island database for the United States: Characterization, uncertainties, and possible applications[END_REF][START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF]. Few studies have focused on the global heterogeneity of impacts on SUHI, the driving mechanisms of global SUHI from global to local perspectives still require further investigation. Advanced spatial statistical models such as the commonly used geographically weighted regression (GWR) are recommended to better reveal the spatial heterogeneity of the driving factors for global SUHI (Li, Zha, and Zhang 2020b), however, due to the limitation of fixed bandwidth of the GWR model, improved analytical tools should be introduced to better quantify the global variation of multiple drivers for SUHI.

To this end, this study aims to explore the potential driving factors of the dynamic surface properties (surface vegetation condition and albedo), climatic conditions (precipitation and wind), and urbanisation processes (population and urban size) for the spatial variation of global long-term SUHI with an improved spatial statistical model. The two analytical perspectives of global and local aspects were considered simultaneously as a primitive attempt to reveal the driving mechanisms of SUHI. Meanwhile, the spatial heterogeneity of the driving factors for global daytime/nighttime SUHIs on the annual, summer, and winter average scales were also disclosed by an advanced spatial statistical model.

Data

The land surface temperature was extracted from the LST products (MOD11A1 and MOD11A1, Version 6) onboard the Terra and Aqua satellites, either of which is a daily surface temperature/emissivity product dataset with a spatial resolution of 1 km. Annual land cover types were extracted from the MODIS land cover product (MCD12Q1, Version 6). The Enhanced Vegetation Index (EVI) was extracted from the vegetation index product (MOD13A2 and MYD13A2, Version 6) onboard the Terra and Aqua satellites, with a spatial resolution of 1 km. Surface albedo was extracted from the MCD43A3 (version 6) dataset. Daily images were generated using 16 days of data (centred on a given date) with a spatial resolution of 500 m. The white-sky albedo (WSA) was selected in this study as it is linearly related to black-sky albedo and demonstrates a similar impact on SUHI [START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF]). The surface elevation was extracted from the GTOPO30 dataset, with a spatial resolution of 1 km, from the US Geological Survey. The meteorological elements were extracted from the monthly climate and climatic water balance datasets for global terrestrial surfaces -TerraClimate with a spatial resolution of 2.5 arc minutes (~ 5 km). Monthly averaged precipitation and wind speed datasets were selected for this study. The population density in the years 2000, 2005, 2010, 2015, and 2020 was obtained from the Gridded Population of the World (GPW, v4.11) population grid density estimation data with a spatial resolution of 30 arc seconds (~ 1 km). For all available datasets, the time span was from January 2003 to December 2019, and the spatial resolution was resampled to 1 km with a sinusoidal projection. For years with missing data between 2003 and 2019, the annual global population density was constructed by borrowing the dataset of its neighbouring years.

Methodology

Definition of SUHII

The instantaneous SUHI intensity (SUHII) was defined as the LST difference between the urban and rural regions as in Equation (1).

__ t t urban t rural SUHII T T  (1)
where the subscript t represents the four time points when Terra transitions at 10:30 and 22:30 and Aqua transitions at 13:30 and 01:30, respectively. _ t urban T and _ t rural T were the mean LSTs for urban and rural regions, respectively. Subsequently, the SUHIIs from the MOD and MYD were averaged to obtain the daily mean daytime and nighttime SUHIIs. Surface elevation data were used to filter the rural pixels which had an elevation difference of more than 50 m from the mean urban elevation. Urban and rural boundaries were defined by a dynamic urban-extent method, where the urban regions were identified by city cluster algorithms with yearly land cover data, and the rural boundaries are defined as the equal-areal buffers after excluding water and impervious surfaces [START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF][START_REF] Si | Spatiotemporal pattern and long-term trend of global surface urban heat islands characterized by dynamic urban-extent method and MODIS data[END_REF]. Finally, the SUHII at the annual, summer, and winter scales was averaged from daily quantification. The summer months of summer (June, July, and August) and winter (December, January, and February) for the Northern Hemisphere are opposite to that for the Southern Hemisphere.

Driving factor indices

Compared with the natural surface in rural areas, the surface energy sources in urban areas have extra radiation from anthropogenic heat, in addition to solar radiation and atmospheric downward radiation. The energy received by the earth's surface is used for heat exchange between the earth and atmosphere in the form of turbulent flow, as well as for soil heat exchange. Taking rural LST as a reference, changes in the net surface shortwave radiation will cause changes in the emission of surface long-wave radiation, which in turn will cause the redistribution of surface energy. The LST change in urban areas can be deduced by a combination of the urban-rural difference of each term related to Equation (1) [START_REF] Lee | Observed increase in local cooling effect of deforestation at higher latitudes[END_REF]. Correspondingly, the difference in surface radiation input and energy distribution directly leads to LST differences between urban and rural areas, that is, the SUHI phenomenon. However, as surface energy flux data with high spatial and temporal resolution and accuracy are difficult to retrieve, several proxy indicators are commonly used to reflect energy variance. Generally, changes in urban atmospheric conditions, surface properties, and climate factors cause differences in net surface shortwave net radiation and surface downward longwave radiation between urban and rural regions.

Based on the findings of previous studies [START_REF] Li | A review of remotely sensed surface urban heat islands from the fresh perspective of comparisons among different regions (Invited Review)[END_REF][START_REF] Sun | Contribution of urbanization to warming in China[END_REF][START_REF] Yang | The footprint of urban heat island effect in 302 Chinese cities: Temporal trends and associated factors[END_REF][START_REF] Sun | Understanding the variability of urban heat islands from local background climate and urbanization[END_REF][START_REF] Manoli | Magnitude of urban heat islands largely explained by climate and population[END_REF], several potential factors associated with surface properties, background climate conditions, and urbanisation were selected to analyse the driving mechanisms of global SUHII at multiple temporal scales. The surface EVI and WSA differences between the urban and rural regions (ΔEVI, ΔWSA) were calculated using the same method as that in Equation (1).

The monthly accumulated precipitation (MAP) and wind speed (WS) in the urban region were averaged to represent the climatic conditions. As for the factors related to urbanisation, the mean urban population density (PD) was assumed to be invariant within a year, and the PD for the years lacking the GPWv411 data was derived from the nearest year, while the urban area (UA) was calculated according to the previously recognised urban clusters. Thereafter, ΔEVI, ΔWSA, MAP, WS, PD, and UA were averaged at annual, summer, and winter scales.

Analytical method

From a global perspective, multiple regression analysis was used to analyse the driving factors of SUHI across cities on annual, summer, and winter averaged scales. In a traditional global regression model, ordinary least squares (OLS) methods are typically used, thus the OLS regression was used by a step-wise process to obtain the influencing factors of multiple drivers at global scale. The long-term average SUHII and each indicator were averaged by time and the regression were conducted by annual, summer, and winter scales, respectively.

Spatial regression models

Considering the spatial variability of the driving mechanism for SUHI from a local perspective, the relationship between SUHI and the related driving factors in different regions also varies with geographical location, which is called spatial non-stationarity. An advanced spatial statistical model can fully consider the first law of geography and detect the non-stationarity of the spatial relationship between variables. Therefore, GWR and multi-scale geographically weighted regression (MGWR) models were simultaneously introduced to further explore the driving factors of global SUHI at multiple spatial and temporal scales.

GWR is a modelling method for local relationships that helps explain spatial relationships under spatial heterogeneity or spatial non-stationary conditions [START_REF] Fotheringham | Geographically weighted regression: the analysis of spatially varying relationships[END_REF]. In contrast to OLS, which can only estimate coefficients in the global sense, GWR uses the local weighted least-squares method to estimate point-by-point parameters by incorporating the spatial position of the observation into the equation. The weight is a function of the distance between the regression point and other observation points, resulting in a continuous and smooth parameter estimation surface [START_REF] Brunsdon | Geographically weighted regression: A method for exploring spatial nonstationarity[END_REF]. The GWR model extends the traditional global regression model and adds geographic location parameters, as follows:

    0 1 ,, p j j j i j j ij j i y u v u v x         (2)
where,   Previous studies have proven that GWR is not sensitive to the choice of weight function, but is very sensitive to the bandwidth of the weight function. In view of the heterogeneous distribution of global SUHI and its drivers, the adaptive variable bandwidth was employed in this study, and the corrected Akaike Information Criterion ( c AIC ) was introduced to determine the optimal model [START_REF] Fotheringham | Geographically weighted regression: the analysis of spatially varying relationships[END_REF]. To evaluate the credibility of the model, it was necessary to further calculate the spatial autocorrelation coefficient of the residuals. The global Moran's I index with a value range of [-1,1] was selected. A negative value indicates a negative spatial correlation, that is, high and low values will accumulate in a cluster; a positive value indicates a positive spatial correlation, that is, a cluster of high and high values, or low and low values, will accumulate in a cluster. A higher absolute value of the spatial autocorrelation coefficient indicates stronger spatial autocorrelation [START_REF] Moran | A test for the serial independence of residuals[END_REF].

With respect to the fact that the GWR model has a scale effect, that is, the effects of certain variables may vary on different scales, while others at the same location may be spatially stable [START_REF] Goodchild | Models of scale and scales of modelling[END_REF]. Accordingly, [START_REF] Fotheringham | Multiscale geographically weighted regression (MGWR)[END_REF] proposed the MGWR model whose essence was designed by relaxing the assumption that all the spatially varying processes in a model operate at the same spatial scale. An optimal bandwidth vector takes place for each particular geographical process considering in reality the discrepancy occurs between micro-and macroprocesses or between local and global processes [START_REF] Mcmaster | Introduction: Scale and Geographic Inquiry[END_REF]. Specifically, the spatial scale at which the conditional relationship between the dependent and each predictor was confirmed to produce a more powerful spatial model MGWR. The mathematical theory behind is to replace the fixed bandwidths in GWR as flexible for each independent variable indicating different processes as follows:

    0 1 ,, p j j j bwi j j ij j i y u v u v x         (3)
where the subscript bwi in bwi  indicates the bandwidth used for calibration of the ith conditional relationship.

Whereas the GWR model constrain the local relationships within each model to vary at the same spatial scale, MGWR allows the conditional relationships between the response variable and the different predictor variables to vary at different spatial scales, which is closer to the real-world geographic processes [START_REF] Yang | An extension of geographically weighted regression with flexible bandwidths[END_REF]. Therefore, the MGWR model was introduced to compare with the GWR model which were exploited within the MGWR software [START_REF] Oshan | mgwr: A Python implementation of multiscale geographically weighted regression for investigating process spatial heterogeneity and scale[END_REF] to analyse the spatial non-stationarity of the driving factors for SUHII, and the results were compared with the findings from traditional multivariate regression schemes.

Results

Model diagnosis

Collinearity analysis

In statistical models, the linear correlation between two or more predictors is called multicollinearity [START_REF] Alin | Multicollinearity[END_REF]. Due to the non-independence of predictor variables, multiple variables may synchronously explain a certain underlying mechanism, resulting in larger deviations in model parameter estimates. In order to avoid the influence of multicollinearity, the Variance Inflation Factor (VIF) is usually used to judge whether multicollinearity exists between the selected predictor variables before the regression analysis. Generally, the closer the VIF is to 1, the smaller the collinearity [START_REF] Dormann | Collinearity: a review of methods to deal with it and a simulation study evaluating their performance[END_REF]. The VIF of the aforementioned six driving factors (ΔEVI, ΔWSA, MAP, WS, PD, and UA) are listed in Table 1. Because there is no day-night variation in the selected variable, the day-night VIF is equal. The VIF value of each variable was close to 1 and far below the threshold value of 10. Therefore, it can be considered that there is no prominent multicollinearity between the selected factors, which can be used for the subsequent analysis of the multiple regression model.

Model comparison for regression coefficients

Table 3 shows a comparison of the regression coefficients of the OLS, GWR, and For the annual daytime SUHII in Figure 1, ΔEVI has a significant negative contribution in most cities around the world, indicating that the increase in vegetation activity means a greater daytime evapotranspiration, which leads to a decrease in SUHII.

MGWR

Spatial non-stationarity

In arid regions, such as North Africa and the Middle East, where cities with cold islands are distributed [START_REF] Si | Spatiotemporal pattern and long-term trend of global surface urban heat islands characterized by dynamic urban-extent method and MODIS data[END_REF], the negative contribution of ΔEVI is more evident. There is a positive correlation between ΔWSA and annual daytime SUHII, and the coefficients vary across different latitude zones. The MAP is positively correlated with the annual daytime SUHII for most cities, and South America displays the largest coefficient magnitude. The WS shows a significant negative contribution in some cities. PD has the largest positive correlation in cities on the east coast of Australia, followed by cities in Southeast Asia and northeastern China. Significant negative correlations appear in some cities in North America, Western Europe, and Africa. UA has significant positive correlations mainly in North America, South America, Africa, and Western Europe in the Western Hemisphere, where PD shows a large proportion of negative contributions.

Figure 1 The spatial variation of coefficients of multi-independents for annual daytime SUHII.

According to Figure 2, for the annual nighttime SUHII, the correlation of ΔEVI in most cities is significantly reduced. ΔWSA has a negative correlation with the annual nighttime SUHII, whose magnitude gradually decreases from the equatorial region to high latitudes, and a few positive contributions arise in the Southern Hemisphere.

According to Equation ( 4), ΔWSA and solar radiation regulate nighttime SUHII by influencing the net surface shortwave radiation during the daytime. Normally, the intensity of solar radiation in the equatorial region is higher than that at high latitudes;

thus, equatorial SUHII is more sensitive to changes in ΔWSA. MAP is negatively correlated with annual nighttime SUHII, while the number of cities with a significant correlation decreased compared with daytime cases. WS shows negative correlations, while the absolute value of the coefficients is lower than that during the daytime. It is inferred that the WS has a smaller impact on the impedance at night, which leads to a weaker regulation effect on the sensible heat flux than during the day. For PD in developed countries such as the United States, Europe, and other regions, the positive correlation at night is greater than in developing countries, while it is relatively minimal in India. The positive effect of UA in North America and South America is weaker at night than during the day, whereas the impact on the SUHII in the Middle East and India has increased significantly, which may be related to the emission of more heat fluxes from human activities in these areas.

Figure 2 The spatial variation of coefficients of multi-independents for annual nighttime SUHII.

Figure 3 and Figure 4 show the spatial distribution of the regression coefficients of the summer daytime and nighttime SUHII, respectively. The negative correlation of ΔEVI with summer daytime SUHII is more evident than that in the annual cases. The latitudinal variations of the correlation with ΔWSA are not that obvious during the daytime, and the negative correlation at nighttime is also lower than the annual average cases. The daytime positive correlation and nighttime negative correlation with MAP in the Northern Hemisphere are significantly higher than those in the Southern Hemisphere. This may be related to the stronger vegetation activities in the Northern Hemisphere during summer.

Precipitation not only regulates the evaporation ratio by affecting the surface soil moisture but also promotes the growth of vegetation, making summer SUHII in the Northern Hemisphere more sensitive to changes in MAP. The negative influence of WS on summer SUHII at night is stronger than that during the daytime. PD had a significantly higher impact on SUHII at night than during the day. UA is significantly positively correlated with summer daytime SUHII in North America and Western Europe in the Western Hemisphere, while in India, the impact of UA is stronger on nighttime SUHII.

Figure 3 The spatial variation of coefficients of multi-independents for summer daytime SUHII.

Figure 4 The spatial variation of coefficients of multi-independents for summer nighttime SUHII.

Figure 5 and Figure 6 show the spatial distribution of the regression coefficients of the winter daytime and nighttime SUHII, respectively. Compared with the results in summer, the negative contribution of ΔEVI in the daytime is relatively weaker, which may be related to the weakening of vegetation activities in winter. The positive correlation of ΔWSA with winter SUHII shows limited significance during the daytime, while the negative correlation at nighttime shows a decreasing magnitude from the equator to higher latitudes. In winter, the positive correlation between MAP and daytime SUHII is more obvious in the equator and arid zone, while the magnitude of the negative correlation between MAP and nighttime SUHII decreases from the equator to the poles.

WS shows a negative correlation, which is more evident during the daytime than at nighttime. PD and UA had significant positive contributions, especially at night. In most cities in the Northern Hemisphere, PD is more positively correlated with nighttime SUHII.

The positive contribution of UA at night is greater in areas such as East Africa and India.

Figure 5 The spatial variation of coefficients of multi-independents for winter daytime SUHII.

Figure 6 The spatial variation of coefficients of multi-independents for winter nighttime SUHII.

Discussion

As one of the most severe urban thermal environmental problems, the SUHI phenomenon based on remote sensing technology acts as a bridge between research on satellite product

Impact of multiple factors

Based on the regression from the three schemes, the results prove that the average contributions of the different statistical models show consistent evidence. Generally, ΔEVI is negatively correlated with daytime SUHII, which is mainly caused by the effect of vegetation activity on the latent heat flux. ΔWSA mainly regulates nighttime SUHII by affecting the net surface shortwave radiation during the day and surface heat storage at night. The MAP was positively correlated with daytime SUHII and negatively correlated with nighttime SUHII, which mainly affects the specific heat capacity, surface impedance, and other parameters related to the surface heat flux in rural natural land surface, resulting in different rates of LST increase during the day and LST decrease at night to regulate the daytime and nighttime SUHII. WS mainly contributes negatively to SUHII by disturbing the surface impedance and regulating the sensible heat flux. The positive correlation between PD and UA with SUHII was more significant at nighttime, mainly due to anthropogenic heat emissions amplifying SUHII.

Although we have compared three different analytical methods, the execution of each model was simultaneously used one by one in order to keep the consistency for impact analysis on global cities. Previously, most SUHI studies have focused on a single city or representative cities and urban clusters on a national or regional scale [START_REF] Sun | Estimating local-scale urban heat island intensity using nighttime light satellite imageries[END_REF]Li, Zha, and Zhang 2020a;[START_REF] Hu | Comparison of surface and canopy urban heat islands within megacities of eastern China[END_REF][START_REF] Yue | Impacts of urban configuration on urban heat island: An empirical study in China mega-cities[END_REF][START_REF] Sun | Understanding the variability of urban heat islands from local background climate and urbanization[END_REF][START_REF] Lai | Statistical estimation of next-day nighttime surface urban heat islands[END_REF], whereas the impact analysis on long-term SUHI at a global scale has been less conducted [START_REF] Wu | A bibliometric review of past trends and future prospects in urban heat island research from 1990 to 2017[END_REF]. Specifically, the relationship between spatially distributed SUHII and associated drivers may not necessarily be consistent due to the different definition of urban and rural regions, various analytical approach, and even different data source [START_REF] Li | Satellite remote sensing of global land surface temperature: definition, methods, products, and applications[END_REF]. For example, studies at the global scale show that city size has a positive impact on SUHI [START_REF] Zhou | The role of city size and urban form in the surface urban heat island[END_REF], while studies at the regional scale showed no significant correlation between them [START_REF] Zhou | Surface urban heat island in China's 32 major cities: Spatial patterns and drivers[END_REF]. To avoid this inconsistency as much as possible, this work presented a comprehensive study across global cities with long-term time span and based on identical analysis method.

Spatial non-stationarity

The limited global studies [START_REF] Zhang | Characterizing urban heat islands of global settlements using MODIS and nighttime lights products[END_REF][START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF][START_REF] Clinton | MODIS detected surface urban heat islands and sinks: Global locations and controls[END_REF] did not completely analyze the spatial variation of impacts on SUHI. In this study, we revealed the spatial heterogeneity of multiple factor-induced impact on global longterm SUHI. Under the combined effect of multiple factors, the driving mechanism of the global SUHII is spatially nonstationary for different seasons. During the day, larger differences in vegetation activity between urban and rural areas generally produce a stronger intensity of cold islands in arid regions, where the evaporative heat flux dominates more significantly and shows a larger negative coefficient on ΔEVI. At night, the influence of ΔEVI was weak across the globe. The contribution of ΔEVI is mainly from the summer scale.

As for ΔWSA, the insignificant positive correlations with daytime SUHII were not apparent globally. The negative magnitude of its coefficient at night has a latitudinal effect, indicating global solar radiation distribution. The magnitude of the contribution is larger in equatorial regions, where higher solar radiation exists and decreases at higher latitudes.

Generally, the ΔWSA in winter shows a more obvious negative contribution than in summer, and the nighttime SUHII is dominated by the negative effect of ΔWSA in the higher latitudinal zone where snow cover appears more frequently.

The positive correlations between the daytime SUHII and MAP showed latitudinal and longitudinal variations. At the annual scale, there is a trend of decreasing correlations from low latitudes to high latitudes and from the Western to Eastern hemisphere. The highest contributions were distributed in the equatorial regions of South America. At night, the negative contribution of MAP also shows a larger magnitude in the Northern Hemisphere, with a limited number of significant correlations compared to daytime conditions. Specifically, in summer, the disparity between the Northern and Southern

Hemispheres is obvious both during the day and night. In the Northern Hemisphere with higher vegetation activity, the LST is more sensitive to surface soil moisture and vegetation growth, thus has a larger contribution to SUHII from MAP. In winter, the influence of MAP is weak compared with in summer.

As for WS, its negative contribution in some cities also displays spatial variability, and the magnitude of the effect is relatively larger in coastal cities. Regarding the summer-winter contrast, the WS in winter daytime shows significant and larger cooling effects on SUHI in most cities across the globe. At night, the WS shows negative correlations, while the absolute value of the coefficients is lower than that of during the day. It is inferred that the WS has a smaller impact on the impedance at night, which leads to a weaker regulation effect on the sensible heat flux at night than during the day.

PD showed significant positive correlations with daytime SUHII in cities in the Eastern hemisphere, while negative correlations appeared in regions with larger positive UA contributions, mainly distributed in some Western cities. At night, the sign and magnitude of PD are even larger in Western developed countries, where the positive contribution of UA becomes weaker, which shows contrary regulation as that at daytime.

There was no significant summer-winter difference in the influencing effect of PD and UA on its spatial heterogeneity.

Implications and prospects

The driving factors for SUHI have been fully investigated for several decades; however, SUHI is typically quantified using land surface temperature, which is a spatial

Conclusions

Revealing the driving mechanisms for the global SUHI effect, especially at multiple This study mainly draws the following findings:

(1) The results from the MGWR model provide an optimal estimation of the contribution coefficients compared to the GWR and OLS models. For MGWR diagnosis, the R 2 values were the highest (0.87, 0.73, 0.90, 0.74, 0.85, and 0.76), and the absolute values of Moran's I were closest to 0 for different temporal scales. The contribution of multiple factors to global SUHI is spatially non-stationary.

(2) ΔEVI was negatively correlated with daytime SUHII, which is more sensitive in arid regions. The summer-winter contrast of its influence is evident.

(3) The day-night difference for ΔWSA mainly comes from the insignificant positive and significant negative coefficients of the daytime and nighttime SUHIIs, respectively.

The latitudinal variation in the negative contribution of ΔWSA is more obvious in winter and more sensitive to larger solar radiation in the lower latitudinal zone compared with in summer.

(4) Evident latitudinal and longitudinal variation appears in the positive correlations between MAP and annual and summer daytime SUHII, which decreases from lower to higher latitude zones and from the Western to Eastern Hemisphere. At night, the negative contribution of MAP is superior in the Northern Hemisphere, but with limited significant cases.

(5) The WS is mainly a negative contribution to daytime and nighttime SUHII, which shows a more evident effect along coastal cities and a significant summer-winter contrast with a larger contribution during winter daytime.

(6) The positive correlations between PD and UA with SUHII are more significant at night, with less seasonal contrast. The day-night disparity mainly originates from the daytime negative correlations from PD, accompanied by a larger positive UA contribution in some Western cities, and vice versa at nighttime.

  estimation parameter of the independent variable ij x , p is the number of variables, j  is the residual at position j . Considering distance attenuation, the GWR model executes the modification by weighting the observations around all sample points. Commonly used weight matrix models include the distance threshold, k-neighbour, quadratic kernel function, and Gaussian kernel function matrices.The most frequently used Gaussian kernel function is chosen as the weight function which embodies the bandwidth for controlling the local regression.

  methods. The regression coefficients of the OLS model are unique on a global scale. For the GWR and MGWR models, the global average and variance of the local coefficients are provided. The OLS model has an insignificant parameter fitting at a 95% confidence interval on several scales. Given the results of the three models, the regression coefficients of ΔEVI are all negative, the regression coefficients of ΔWSA and MAP are positive during the day and negative at night, and the regression coefficients of WS are all negative, while the regression coefficients of PD and UA are positive in most cases.The sign of the influence of each factor on the SUHII is consistent with the results of previous studies, while the magnitudes vary between different models due to different fitting process, which indicates the reasonability and reliability of the regression models for characterising the drivers for SUHII[START_REF] Si | Spatiotemporal pattern and long-term trend of global surface urban heat islands characterized by dynamic urban-extent method and MODIS data[END_REF][START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF][START_REF] Chakraborty | A simplified urban-extent algorithm to characterize surface urban heat islands on a global scale and examine vegetation control on their spatiotemporal variability[END_REF]) .

Figure 1 -

 1 Figure 1-Figure 6 show the spatial distribution of the regression coefficients of each

  algorithms and urban remote sensing applications. In view of the uncertainty of the driving mechanisms and the deficiency of heterogeneity of the influencing effect for global long-term SUHI, this study proposes the necessity of exploring the spatial and temporal variation of influencing magnitude and signs by introducing advanced spatial statistical models. First, the global daily instantaneous daytime and nighttime SUHII dataset from 2003 to 2019 were generated and averaged over different seasons. Second, the driving factors for the global SUHII were analysed from a global perspective, which was further compared with the GWR and MGWR models at the local scale. The results and findings are expected to provide target indications for mitigation policies in urban management in different geographic regions.

  schemes exhibited a better model fitting effect, and the latter performed best and

Table 1

 1 Variance inflation factor (VIF) of explanatory variablesThe model fitting and diagnosis information for the multivariate regression, GWR, and MGWRs models at the average annual, summer, and winter scales are shown in Table2. The diagnosis results of the OLS, GWR, and MGWR models change with day-night and seasonal variations. For each temporal scale, the AICs values of the OLS, GWR, and MGWR models gradually decreased, and the lower AICs of GWR and MGWR models indicate a reduction of model complexity which can better avoid the model overfitting than that of OLS. The coefficients of determination (R 2 ) of the OLS regression models are 0.57, 0.24, 0.67, 0.24, 0.34, and 0.31, respectively, which are worse than those of the GWR and MGWR models for the different temporal scales. The R 2 values of the MGWR

		ΔEVI	ΔWSA	MAP	WS	PD	UA
	Annual daytime	1.38	1.25	1.47	1.24	1.10	1.04
	Annual nighttime	1.38	1.25	1.47	1.24	1.10	1.04
	Summer daytime	1.11	1.22	1.30	1.29	1.06	1.03
	Summer nighttime	1.11	1.22	1.30	1.29	1.06	1.03
	Winter daytime	1.39	1.21	1.22	1.16	1.12	1.05
	Winter nighttime	1.39	1.21	1.22	1.16	1.12	1.05

model are 0.87, 0.73, 0.90, 0.74, 0.85, and 0.76, respectively. Moreover, compared with the OLS model, the absolute values of Moran's I for the GWR and MGWR models are closer to 0, indicating that they can effectively capture the spatial non-stationarity of the driving factors, and that MGWR performs better than GWR. In general, the explanatory factors of the SUHII have spatial non-stationarity, and the MGWR model has the optimal performance for the interpretation of global SUHI from a local perspective.

Table 2

 2 Diagnostic Information of OLS, GWR and MGWR models

		Annual	Annual	Summer	Summer	Winter	Winter
		daytime	nighttime	daytime	nighttime	daytime	nighttime
	OLS diagnostics
	AICs	3407.272	4398.981	2975.548	4401.896	4170.915	4235.052
	R 2	0.575	0.240	0.670	0.237	0.335	0.310
	R 2 adjusted	0.573	0.238	0.669	0.234	0.333	0.307
	Moran's Index	0.521	0.440	0.668	0.451	0.572	0.489
	GWR diagnostics
	AICs	2029.597	3415.593	1697.603	3264.442	2404.918	3147.025
	R 2	0.867	0.731	0.889	0.760	0.856	0.757

Table 3

 3 Regression coefficients of OLS, GWR and MGWR models

		Annual	Annual	Summer	Summer	Winter	Winter
		daytime	nighttime	daytime	nighttime	daytime	nighttime
	OLS global regression
	Intercept	-0.000	-0.000	0.000	0.000	0.000	0.000
	ΔEVI	-0.544	-0.182	-0.719	-0.216	-0.469	-0.111
	ΔWSA	0.152	-0.184	0.241	-0.068	-0.024*	-0.392

  variable at the local scale. Assessing the spatial heterogeneity of the driving mechanisms for global long-term SUHI at multiple temporal scales is usually not a point of focus. In this study, we presented a comprehensive investigation on the influencing effects of multiple factors on global SUHI from global to local perspective and distinguished the day-night and summer-winter contrast with an improved MGWR model. To the best of our knowledge, this work has been paid little intention in previous global SUHI studies. This proves to be a significant consideration in further detection of global SUHI and its underlying mechanisms. Meanwhile, our findings exhibit magnitudes and signs of drivers consistent with the results from a global perspective in previous studies, indicating that different quantification schemes of driving mechanism assessment may not affect the average pattern of influencing factors for global SUHI. The findings on global SUHI variation and its associated factors are assumed to have implications for urban planning. factors by distinguished temporal scales, that is, day-night and summer-winter contrast. With the development of more advanced simultaneous spatiotemporal GWR models in the future, the findings of this study are expected to be further validated using other schemes.

	Furthermore, it is necessary to assess additional potential driving factors of anthropogenic
	activities (e.g. air pollution and human heat flux) and urban landscapes (Cao et al. 2016;
	Li et al. 2018; Liu and Weng 2009; Lu and Weng 2006; Weng and Lu 2008; Huang and
	Wang 2019).
	Generally, spatial non-stationarity is assessed using advanced spatial statistical
	models, which can also be applied to SUHI research. To further reveal the temporal
	heterogeneity of global long-term SUHI, we analysed the disparity of spatial non-
	stationarities of driving

Acknowledgements

We would like to acknowledge the Google Earth Engine platform for providing MODIS, Terra Climate, and topographic data, and the cloud computing capability free of charge. This work was supported by the National Natural Science Foundation of China (Grant No. 41921001 and No. 42201404), and the China Postdoctoral Science Foundation (No. 2022M723439). The author Menglin Si would like to thank the China Scholarship Council for the financially support of her stay in ICube, France.


E Centeno

J Cellier

C Bernard

R Rapegno

F Réveret

E Tomasella

A Farhaoui

A Bousquet

R Smaali

Reactive Gas Pulsing Sputtering Process, a promising technique to elaborate silicon oxynitride multilayer nanometric antireflective coatings

Keywords: Oxynitride, Reactive sputtering, antireflective coating

come

Reactive gas pulsing sputtering process, a promising technique to elaborate silicon oxynitride multilayer nanometric antireflective coatings

Introduction

Over last few decades, oxynitride materials have been extensively studied as potential candidates to tune mechanical, optical or electrical properties between oxide and nitride ones [START_REF] Bedjaoui | Physico-chemical, structural and physical properties of hydrogenated silicon oxinitride films elaborated by pulsed radiofrequency discharge[END_REF][START_REF] Martin | Enhanced sputtering of titanium oxide, nitride and oxynitride thin films by the reactive gas pulsing technique[END_REF][START_REF] Cristea | Properties of tantalum oxynitride thin films produced by magnetron sputtering: The influence of processing parameters[END_REF]. Among them, silicon oxynitrides present a great interest for numerous applications, such as an intermediate dielectric constant for microelectronic [4] or a diffusion barrier for boron, phosphor [5] and also for gas [6,7].

More recently, they are increasingly used as antireflective coatings and passivation layers to improve efficiency of solar cell [8,9,10]. In this domain, a recent study has for instance shown that multilayered antireflective coatings using SiO x N y give better transmission coefficients than single layered silicon nitride, and therefore they increasing the penetration of photons into the solar cell [11]. These silicon oxynitrides are commonly prepared by low-pressure plasma enhanced chemical vapor deposition (PECVD) [8][9][10]. However, recently, the sputtering process has started to be used in photovoltaic industries thanks to its several advantages the process avoids using toxic gases as silane is used as the silicon source, and is of its low cost and high versatility in the deposition of film composition.

The silicon oxynitride deposition from the silicon target sputtering in the Ar/O 2 /N 2 atmosphere was shown to be possible [12,13]. However, the reactive sputtering process, especially with two reactive gases, is known to be very difficult to control. Indeed, target poisoning by reactive gas induces a memory effect, which leads to a nonlinear process behavior. This phenomenon is even more complex with two reactive gases because of the competitive reactions between both gases [14,15,16]. Different techniques can be used to overcome this problem, for instance, a process mapping [16] or a feedback control of reactive gas partial pressure [17]. Another technique was proposed by Aronson et al. which and consists of pulses to the reactive gas flow rate to switch between the elemental sputtering mode and the compound sputtering one [18]. This technique is called Reactive Gas Pulsing Process (RGPP). It was successfully used to deposit titanium, iron or silicon oxynitride films [START_REF] Martin | Enhanced sputtering of titanium oxide, nitride and oxynitride thin films by the reactive gas pulsing technique[END_REF]19,20] by firstly using a rectangular shape for pulsed flow rate. Later, this technique reached an even better control by using an exponential shape of pulsed flow rate [21].

The present paper aims to use the RGPP for silicon oxynitrides film deposition with a fine tune of elemental composition in order to control their optical properties, especially their refractive indexes.

Having materials with a controlled refractive index in a large range, allow us ti synthesize antireflective multilayer coatings for silicon wafer. Our goal is to grow this multilayer system by only changing the pulsed gas flow parameters during silicon target sputtering, contrary to a system where different targets are sputtered to obtain different materials with different refractive index, for instance in bilayer of MgF 2 -Al 2 O 3 or trilayer of MgF 2 -ZrO 2 -CeF 3 [22]. Furthermore, we aim to obtain this antireflective system not only for normal incident light but also for a large range of incident angles.

In this paper, we will first check, as already demonstrated in other paper [Erreur ! Signet non défini.0,1], that the RGPP technique allows us to deposit silicon oxynitride films with a controlled elemental composition. Deposition rates and film structure will also be studied. After that, the optical properties of these films will be more deeply investigated. Finally, an antireflective multilayer system will be simulated, optimized and realized thanks to the RGPP technique.

Experimental details

SiO x N y thin film deposition

Silicon oxynitride films were deposited by radiofrequency reactive magnetron sputtering. A pure Silicon target (99.99% and 100 mm in diameter) is placed at 95 mm from the substrate and sputtered in Ar/O 2 /N 2 atmosphere. Before the deposition, the sputtering chamber was evacuated to 10 -5 Pa and then the silicon target was sputtered for 30 min at 250 W in pure argon atmosphere at 1 Pa to clean up the surface. Then, the gas mixture is introduced for 2 min in order to prepare the target to the reactive conditions. The deposition is then started by removing the shutter in front of the substrate. The Ar, N 2 and O 2 gas flow rates,  Ar ,  N2 and  O2 respectively, are controlled by mass flowmeters. For all the study, the Ar flow rate is fixed at 7.8 sccm. The radiofrequency power was fixed at 250 W. The total gas pressure was 1 Pa. The substrate neither biased nor heated. The deposition time was modified to obtain 200 nm thin films for monolayer material characterization or a lower thickness for multilayer elaboration.

Two sets of samples were deposited with Conventional (CP) and Pulsed processes (PP). The CP referred to a continuous flow rates injection during the film deposition. The oxynitride films were obtained by changing both the O 2 and the N 2 flow rates. In order to keep constant the total pressure, we fixed the total reactive flow rate,  O2+N2 , and we varied the O 2 one in the total reactive flow rate ratio, R F , from 0 (Ar/N 2 atmosphere) to 1 (Ar/O 2 atmosphere), where R F is defined by:

(

In this conventional process mode, we also deposited 3 films illustrating the deposition in an Elemental Sputtering Mode, when the target is a pure silicon, and in the Compound Sputtering Mode, where the target was covered by an oxide and by a nitride. For the first sample, only Ar ( Ar = 7.8 sccm) was injected; whereas both others samples were deposited with Ar/N 2 ( Ar = 7.8 sccm and  N2 = 3.0 sccm) and Ar/O 2 ( Ar = 7.8 sccm and  O2 = 2.2 sccm) mixtures, respectively.

In Pulsed Process, Ar and N 2 gases were introduced continuously at constant mass flow rates (7.8 and 1.4 sccm for Ar and N 2 respectively), whereas O 2 was periodically supplied using a rectangular pulsed flow rate, from 0 sccm during the off-time, T off , to a maximum flow rate,  max (O 2 ), of 0.4 or 0.8 sccm during the on-time, T on . The period, T, was fixed at 20 or 40 s and the duty cyle (DC = T on /T) was varied from 0.1 to 0.9. From these pulse conditions ( max (O 2 ), T and DC), the average injected O 2 flow rate, < O2 > , was defined following this equation:

(2)

SiO x N y monolayer characterization

The film elemental composition was determined by Rutherford Backscattering Spectrometry (RBS) using 2 MeV apha particles and 50 nA current intensity. A charge of 40 C was collected at a detection angle of 172°. Since carbon appears at lower energy than silicon, oxygen and nitrogen, the films were deposited on a vitreous carbon substrate in order to avoid overlapping between the substrate and film elements. The experimental spectra were simulated using the SIMNRA program [START_REF] Layer | SMINRA user's guide[END_REF].

A qualitative bond analysis was performed employing Fourier Transformed Infrared spectroscopy (FTIR) for absorption measurements in a transmission mode in the range of 400 -4000 cm -1 , with a resolution of 4 cm -1 and 250 scans per spectrum. To limit the infrared light intensity loss, films were deposited on a double side polished intrinsic silicon wafer.

The optical properties of monolayers deposited on silicon substrates were analyzed by spectroscopic ellipsometry in a range of 0.59 to 4.79 eV at a 70° fixed incident angle. The measurements were then simulated with a "New Amorphous" model. In this model, the dispersion formula, which was derived from the Forouhi-Bloomer formulation, was established to give a Lorentzian shape to the expression of the extinction coefficient and the refractive index, and then this was better for the inter-band absorption [START_REF] Forouhi | [END_REF]25,26]. In this model, the extinction coefficient and the refractive index are described by the following equations:

and Where

In this model, f j is related to the strength (and so, to the amplitude) of the extinction coefficient peak,  i is the broadening term of the peak absorption,  i is the light pulsation, at which the extinction coefficient is maximum,  g is the pulsation corresponding to the band gap energy.

This model is particularly consistent for amorphous materials, which present an absorption in the visible or UV range and was already used to simulate ellipsometric measurements of silicon oxynitride thin films [8].

The optical bandgap was determined from UV-Visible spectroscopy performed with a double beam Perkin Elmer spectrometer in 200-1100 nm range for films deposited on quartz substrates.

Multilayer antireflective optical properties simulation and characterization

Among the several gradual index functions (linear, quintic, Gaussian, etc) [27], we chose a linear graded refractive index profile. Ideally the spatial index variation is continuous and varies adiabatically from the lower to the higher optical indices of the targeted materials. This strategy demands a thick antireflective coating (ARC) and imposes important process constraints to realize the refractive index close to 1 at the air interface. Here, we focus on unmatched refractive index gradients which ranges from n min = 1.53 to n max = 1.80, which are boundary values imposed by the available fabrication process conditions. This approach allows us to produce thin ARCs presenting nanogradients of sub-wavelength spatial variations. Here, the ARCs are backed by c-Si substrate. The linear gradient is spatially discretized by a step function consisting of a stack of p-layers of equal thicknesses h and refractive indices n q :

Linear profile:

Where q is an integer that varies between 1 and p-1. The ARCs are optically optimized according to the number of the layers, p, and to the thickness h to reduce the system's reflectance averaged over wavelengths between 300 nm to 750 nm and for incident angle ranging from 0° to 80°. The refractive indices of the p-layers were chosen for those that obtained the wavelength 600 nm since the optical dispersion of SiO x N y materials can be neglected for visible frequencies. However, the optical properties of the c-Si substrate were taken into account in the calculations by integrating spectroscopic ellipsometry measurements retrieved in the spectral domain. These electromagnetic simulations are obtained with a homemade code based on a rigorous S-matrix method [28].

To compare with simulation results, the realized multilayer ARCs are characterized by reflectivity measurements. The structures are deposited on a double side polished silicon wafer. Then, measurements are performed at room temperature under normal incidence. A xenon lamp is used as a source light to illuminate the sample. The reflected signal is directed on the slits entrance of a 64 cm focal length monochromator. The optical signal is detected by a CCD camera for the 230-750 nm range, and a germanium detector cooled at liquid nitrogen temperature is used with a lock in amplification for larger wavelengths. Then the recorded signal is normalized with a measurement performed on a standard aluminum mirror.

Results and comments

Composition and deposition rate

The elemental composition is determined by Rutherford Backscattering Spectrometry for both film sets deposited with CP and PP. These elemental contents are determined at ± 1 at.% due to RBS measurements incertitude [15]. On one hand, the Figure 1a presents the Si, O and N content (in atomic percent, at.%) for films deposited in the Continuous Process as a function of R F . We remark that the silicon content is almost constant around 35 at. %, whereas O and N contents vary in a larger range: from 28 to 63 and from 0.1 to 34 at.% respectively. The film composition can then be characterized by O/(O+N) content ratio, which is plotted in Figure 1b for films deposited in CP and PP modes. This ratio changes from 0.34 to 0.94 with < O2 > in PP, and from 0.45 to 1.00 with R F in CP. This Figure then highlights < O2 > is the main parameter that controls the film composition in Pulsed Process. For < O2 > ≤ 0.4 sccm, the deposited films are oxynitrides with a linearly increasing O/(O+N) ratio, whereas for < O2 > ≥ 0.55 sccm, the film elemental composition is close to silicon oxide one (O/(O+N) ratio > 0.8). In the first region (<(O 2 )> ≤ 0.4 sccm or R F ≤ 0.45), the films deposited with T = 20 s in PP mode present higher O/(O+N) ratios than those obtained in the CP mode, without being majorly affected by  max (O 2 ). Indeed, reducing  max (O 2 ) from 0.8 to 0.4 sccm allows the injection of fine amounts of < O2 > and then reaching smaller modification in oxynitride composition. However, for  max (O 2 ) = 0.4 sccm and DC = 0.9, the O/(O+N) ratio decreases, which could be due to different target surface for longer on-time in this condition.

In the second region (<(O 2 )> > 0.5 sccm or R F > 0.6), the higher O/(O+N) ratios are reached in the CP mode than in PP mode. Moreover, the O/(O+N) ratio for T = 20 s (full triangle) is higher than for T = 40 s (full square). This second minor effect on film composition could again be due to the influence of target oxidizing and nitriding speeds depending on T on and T off pulse times. This difference in speed was already studied in the reactive gas pulsed sputtering process [Erreur ! Signet non défini.]. However, further study in plasma kinetics is needed for a better understanding. Nevertheless, we can emphasize the interest of the pulsed process to adjust finely the injected O 2 flow rate to control the O and N contents in SiO x N y films. From thicknesses determined by ellipsometry measurements, we calculate the film deposition rates.

The Figure 2a presents values obtained with CP. The growth rates increase slightly from 25 to 35 nm.min -1 with R F . These values are close to the ones we classically observed in silicon reactive radiofrequency sputtering [12]. Because of the deposition of additional samples in CP, which correspond to different sputtering modes, we may also notice that these rates are lower than those for the silicon-rich film (at 40 nm/min) but higher than values when sputtering an oxidized or a nitrided target, respectively at 8 and 12 nm/min. This indicates that the deposition does not occur in a compound sputtering mode but in a transition one. For the Pulsed Process, the deposition rates are plotted in Figure 2b as function of < O2 >. For < O2 > ≤ 0.4 sccm, the deposition rates are slowly increasing from 20 to 25 nm.min -1 . For < O2 > ≥ 0.55 sccm, the deposition rates depend on the pulse period: they stay around 20 nm.min -1 for T = 40 s and present a maximum for T = 20 s. Overall, all these values are slightly lower than those obtained in CP, but they are clearly higher than those observed in the Compound Sputtering Mode. We may consider that the target is not completely poisoned in these conditions, but in the transition mode. Again, the pulse period seems to be a second order influencing parameter. Indeed, deposition rates are slightly higher for shorter pulse period. Since the films deposited with T = 20 s are a little bit richer in oxygen, this is consistent with the small increase of deposition rate with higher < O2 >.

To better understand the growth rate variation, we determine the silicon flux incorporated into the film per time and area unit. It is well known that the poisoned target has a reduced sputtering yield which leads to lower deposition rate [29]. Because of the RBS measurements, the total number of atoms number deposited by area unit, N total , can be known and converted to he silicon flux incorporated into the film per area and time unit, F(Si), using the following equation :

Where N total : is the total number of atoms incorporated into the film by area unit, %Si : is the silicon content, t deposition : is the deposition time.

This silicon flux F(Si) is added to Figure 2a and presents the same evolution than the deposition rate in CP. Indeed, a slight increase of F(Si) for higher R F explains the slight increase of the deposition rate. This rise could be due to a difference of the target's partial coverage by an oxide and/or a nitride, which influences the total sputtered silicon flux from the target, and/or due to an ion bombarding flux modification with atmosphere change [14]. In PP, F(Si) is ranging from around 4-5.10 16 at.cm -2 .min -1 . These values of incoming silicon fluxes are relatively close to the one observed in CP (6-8.10 16 at.cm - 2 .min -1 ), which explains why the deposition rates in PP remain a little bit lower than those in CP. This is also coherent with the sputtering of the target in the elemental or the transition mode, but not with the compound sputtering mode.

In However, the elemental composition determined from RBS measurements is a global composition. So for a better understanding of the films structure, FTIR measurements have been performed.

Film structure

The Figures 3 a andb show the IR spectra for  max (O 2 ) = 0.8 and 0.4 sccm respectively, focused on the molecular bonding vibrations zone of the Si bonds in the 600-1600 cm -1 region. These spectra have been normalized by film thickness in order to plot the absorption. In this region, we may observe the absorption bonds of oxide with Si-O and Si-O-Si bonds (around 1050-1070 and 1100-1250 cm -1 respectively), nitride with Si-N (around 790-890 and 960-970 cm -1 ) and oxynitride with O-Si-N bond (around 850-1000 and at 1100 cm -1 ) [30,31]. For  max (O 2 ) = 0.8 sccm, the peak maximum  0 shifts towards higher wavenumbers, when the duty cycle increases. For  max (O 2 ) = 0.4 sccm, this shift to a higher wavenumber is again observed when DC rise to 0.7. However for DC = 0.9, a more important nitride contribution is observed. This unexpected behavior could be due to a specific interaction between both the reactive gases for a long O 2 injection time, but it is also in agreement with the low O/(O+N) ratio observed in Figure 1b. Further studies on the process of the kinetics are needed to explain this point.

The peak maximum  0 shifts is seen more clearly in Figure 4 as a function of the film composition ratio, demonstrating that it is linked to the film nature. Hence, for O/(O+N) ≤ 0.6, the  0 increases from 891 to 943 cm -1 , corresponds to the major bonding variation from nitride to oxynitride. For O/(O+N) ≥ 0.7,  0 rises to values higher than 1050 cm -1 , indicating there is mainly oxide bonding present. To go and check the influence of both < O2 > and the pulse period, all the absorption spectra were decomposed into Gaussian contributions corresponding to the different bonding types, previously reported. The relative contributions of different bonds are plotted in Figure 4b as a function of < O2 > for T = 20 and 40 s. These relative contributions correspond to the contribution of a specific bond peak to the total absorption band. Since the IR absorption peak area depends not only of the bond density, but also on the oscillator strength, we need to pay attention with this plot as the higher relative area does not necessary mean higher bonding content in the film. However, we might observe their qualitative evolution. For < O2 > ≤ 0.4 sccm, the O-Si-N contribution decreases, whereas the oxide bond contribution increases and the Si-N bond contribution stays low. This is in agreement with the  0 evolution, indicating oxynitride films. For < O2 > ≥ 0.55 sccm, the oxide bond contribution increases, whereas O-Si-N and Si-N are lower. Again, this is consistent with the  0 position around 1050 cm -1 and O/(O+N) content > 0.8, corresponding to a film composition close to silicon oxide. Finally, tuning < O2 >, thanks to the chosen pulse parameters, controls not only the global atomic concentration in O and N but also the element bonded to silicon.

Optical properties

After succeeding to finely vary the composition of the SiO x N y films, we performed spectroscopic ellipsometry analysis to study their optical properties. The Figure 5a presents the refractive index in the studied energy range for films deposited with the PP for  max (O 2 ) = 0.8 sccm, T = 20 s and for various DC. All of the films refractive indices have the same evolution with their wavelength; they only shift to lower n values for higher DC. The Figure 5b summarizes the refractive index values at 1.95 eV for films deposited with the PP as a function of < O2 >. The < O2 > was shown previously to control the films elemental composition, and here we can obviously notice its effect on the refractive indices control. Indeed, for the PP, when < O2 > increases up to 0.4 sccm, n decreases linearly from 1.83 to 1.55, which in agreement with the film's composition changing from N-rich to O-rich SiO x N y [32]. Then, for < O2 > ≥ 0.55 sccm, the film's composition is stabilized close to an oxide one, and the refractive index slowly decreases from 1.5 to 1.45. Finally, the n range is reached with the pulsed mode is the same than with the conventional mode (plotted as stars) between the nitride (2.03 at 1.95 eV [START_REF] Palik | Handbook of optical Constants of Solid[END_REF]) and the oxide ones (1.46 at 1.95 eV [START_REF] Palik | Handbook of optical Constants of Solid[END_REF]). In conclusion, by only changing the pulse parameters of one reactive gas, it is possible to tune the refractive index of the deposited material, which could be particularly interesting for optical devices such as multilayer antireflective coatings. However, before studying the antireflective multilayer optimized system, it is necessary to check the monolayer transparency. A first indication is given by the extinction coefficient k as function of the wavelength, plotted in Figure 6a for several pulse conditions. When the nitrogen content increases while < O2 > decreases, the extinction coefficient rises especially for the photon energy that is higher than 2.5 eV. However, even for the lowest DC, for the highest N-content film, k stays lower than 0.005 in the visible range (1.6 to 3.1 eV) and the films IR and UV light absorption seems to be limited.

To confirm this, optical bandgaps are determined from UV-Visible spectroscopic measurements, considering an indirect gap in Tauc relationship [START_REF] Tauc | [END_REF], and they are plotted in Figure 6b. The optical bandgap increases from 4.55 to 4.83 eV when O/(O+N) content ratio increases from 0.34 to 0.64. For the higher O/(O+N) ratio, the optical bandgap can no more be determined from our measurements performed in the 1.1 to 6.0 eV range, meaning that values are higher than 5.0 eV. The referenced bandgaps for silicon nitride and oxide are respectively about 4.6 and 8.0 eV [START_REF] Cristea | Properties of tantalum oxynitride thin films produced by magnetron sputtering: The influence of processing parameters[END_REF]. In our case, the values lower than 5.0 eV are consistent with oxynitride amorphous materials, as reported in [35]. In general, the optical absorption edge shape of amorphous films is directed by two phenomena: at high energy, absorption varies according to Tauc-type expression, whereas at low energy absorption, it follows an exponential function. This latter is due to a broadening of states because of a disordered structure besides the presence of exponentially distributed states below the absorption edge. Indeed, the films defect level can be characterized by the exponential absorption edge's slope, called the Urbach tail energy, E u , which is plotted in Figure 6b. For all films, E u stays roughly between 0.7 and 0.5 eV, which is in agreement with values observed for amorphous oxynitride films [START_REF] Layer | SMINRA user's guide[END_REF]36], taking into account that the values for O-rich films cannot be determined. Finally, because of the large variation of refractive indices, in addition to the good transparency properties, we can use the SiO x N y films obtained with PP for antireflective multilayer systems.

Simulation and realization of antireflective system

The ARC was optimized by computing the reflectivity for a different number of layers. We found that the optimal structure is achieved with 3 layers presenting a linear spatial variation. In that case, the step refractive index function takes the following values [1.53, 1.66, 1.80] from the air-ARC interface to the Si-ARC interface. Furthermore, for the tri-layer ARC, the thickness of each layer was varied. In Figure 7, the influence of light incident angle on ARC is studied. The averaged reflectivity over the wavelengths between 300 nm to 750 nm is then represented as a function of the incident angle for several individual thicknesses h. The lowest reflectivity (= 18 % for incident angle up to 45) is reached for an ARC of 75 nm total thickness, which corresponds to 25 nm per layer, and stays constant to incident angle of 45°. This optimal graded refractive index ARC was then deposited thanks to RGPP technique and characterized. The measured reflectivity for a normal incidence is represented in figure 8a for a bare polished Si substrate taken as a reference, and this is then covered with an AR system: either a monolayer, which is expected to be the first layer (25 nm, n = 1.8) of the tri-layer ARC system, or the optimized tri-layer. We might notice that the monolayer allows a decrease of the silicon reflectivity by a third; whereas as predicted with the model calculations, low reflectivity values are reached with the tri-layer. This trilayer presents a minimum reflectivity of 3.0 % around 510 and 520 nm and an averaged reflectivity in the visible range (400 -750 nm range) of 8.1 %. This result is comparable to the one obtained in the literature for the same refractive index range and without substrate texturation [37,38]. In addition, in order to check the refractive index and thickness of each real layer in tri-layer system, we performed a comparison between the experimental and theoretical reflectivity spectra. To do so, we studied the tri-layer system by ellipsometry at 70° (as the only available incident angle with this apparatus) and fit the experimental data. An almost perfect agreement is achieved for thicknesses [25, 34.5, 30] nm, which corresponds to the respective refraction index variation [1.53, 1.66, 1.80] and is presented on Figure 8b. The small discrepancy is due to a difference either in thicknesses or refractive indices between the model results and the realized structure. This could be explained by a slight variation in the estimated deposition rate for very thin films. But, we think that this difference is more likely linked to a small difference in refractive indices. Indeed, the optimized refractive indices, as mentioned previously, are [1.53, 1.66, 1.80] from the air-ARC interface to the Si-ARC interface, whereas the refractive indices of realized films are [1.55, 1.66, 1.80]. This small difference of the upper layer refractive index could be the reason of such a small observed difference.

Finally, this optimized ARC system is realized easily and quickly thanks to RGPP technique. Hence, with the deposition rates higher than 20 nm/min obtained with this process, the tri-layer system is deposited in less than 4 min. Furthermore,

Conclusion

In this paper, the RGPP is used as a powerful technique to control finely the composition of the SiO x N y films and then the refractive indexes, only by tuning the pulsing conditions. We demonstrate that a large range of refractive indices can be reached, in comparison to the one reached in CP, without drastically reducing the deposition rate. Furthermore, it is possible to go forward and use it to obtain the antireflective multilayer system. For this goal, an antireflective coating structure is optimized using a genetic algorithm and then realized with the RGPP technique. A reflectivity reduction by a factor of 6 is observed compared to bared polished silicon wafer, reaching the averaged reflectivity lower than 8.1 % (in 400 -750 nm visible light range). Nevertheless, this reflectivity result could further be improved. If we go forward, a more efficient graded refractive index system could be reached by enlarging the refractive index window which matches the silicon and the air indices. Furthermore, we could imagine that, in addition to the interferential effect, the contribution of the graded structure on the reflectivity.

Figure 1 :

 1 Figure 1: (a) Elemental composition of films deposited in Continuous Process and (b) O/(O+N) ratio for films deposited in continuous (star) and in pulsed process with different periods and maximum O 2 flow rates.

Figure 2 :

 2 Figure 2: (a) Deposition rates (open) and incorporated silicon flux (full) calculated in CP and, (b) deposition rates in PP for both periods T = 20 (open) and 40 s (full).

Figure 3 :Figure 4 :

 34 Figure 3: IR spectra for films deposited in PP with various duty cycles, for T = 20 s (full line) and T = 40 s (dash line), for (a)  max (O 2 ) = 0.8 sccm and (b)  max (O 2 ) = 0.4 sccm.

Figure 5 :

 5 Photon energy (eV)

Figure 6 :

 6 Photon energy (eV)

Figure 7 :

 7 Figure 7: Simulated average reflectance as function of the incident angle for various h.

Figure 8 :

 8 Figure 8: (a) Measured reflectivity for a normal incidence on a bare polished silicon substrate and covered by optimized mono-and tri-layer ARC, (b) measured and simulated reflectivity for a 70° incident angle.

  the literature, similar or larger composition ranges can be reached: O/(O+N) from 0.31-0.95 and 0.15-0.99 for FeO x N y and SiO x N y respectively in[19,Erreur ! Signet non défini.0], but this is often detrimental to the growth rate: which decreases from 20 to 5 nm.min -1 or stays around 5-6 nm.min -1 for FeO x N y and SiO x N y respectively in [21,Erreur ! Signet non défini.2]. Finally, in the present paper, we confirm that by using the reactive gas pulsing technique, it was possible to deposit films with a relatively large composition variation, depending on pulse parameters, without losing in the deposition speed.

Acknowledgments: This work was founded by the French National Research Agency, the European Commission (Feder funds) and the Région Auvergne in the Framework of the LabEx IMobS 3 .